diff --git a/docs/reference_guides/model_libraries/power_generation/unit_models/cross_flow_heat_exchanger_1D.rst b/docs/reference_guides/model_libraries/power_generation/unit_models/cross_flow_heat_exchanger_1D.rst new file mode 100644 index 0000000000..ed995a84c8 --- /dev/null +++ b/docs/reference_guides/model_libraries/power_generation/unit_models/cross_flow_heat_exchanger_1D.rst @@ -0,0 +1,323 @@ +CrossFlowHeatExchanger1D +======================== + +.. index:: + pair: idaes.models_extra.power_generation.unit_models.cross_flow_heat_exchanger_1D;CrossFlowHeatExchanger1D + +.. module:: idaes.models_extra.power_generation.unit_models.cross_flow_heat_exchanger_1D + +This model is for a cross flow heat exchanger between two gases. The gas in the shell has a straight path, +while the gas in the tubes snakes back and forth across the shell's path. Note that the ``finite_elements`` +option in the shell and tube control volume configs should be set to an integer factor of ``number_passes`` +in order for the discretization equations to make sense as a cross-flow heat exchanger. + +Example +------- + +.. code-block:: python + + import pyomo.environ as pyo + + from idaes.core import FlowsheetBlock + import idaes.core.util.scaling as iscale + from idaes.models.unit_models import HeatExchangerFlowPattern + from idaes.models.properties.modular_properties import GenericParameterBlock + from idaes.models_extra.power_generation.properties.natural_gas_PR import ( + get_prop, + EosType, + ) + from idaes.models_extra.power_generation.unit_models import CrossFlowHeatExchanger1D + from idaes.core.util.model_statistics import degrees_of_freedom + + optarg = { + "constr_viol_tol": 1e-8, + "nlp_scaling_method": "user-scaling", + "linear_solver": "ma57", + "OF_ma57_automatic_scaling": "yes", + "max_iter": 350, + "tol": 1e-8, + "halt_on_ampl_error": "no", + } + + m = pyo.ConcreteModel() + m.fs = pyo.FlowsheetBlock(dynamic=False) + m.fs.properties = GenericParameterBlock( + **get_prop(["H2", "H2O", "Ar", "N2"], {"Vap"}, eos=EosType.IDEAL), + doc="H2O + H2 gas property parameters", + ) + m.fs.heat_exchanger = CrossFlowHeatExchanger1D( + has_holdup=True, + dynamic=False, + cold_side={ + "property_package": m.fs.h2_side_prop_params, + "has_holdup": False, + "dynamic": False, + "has_pressure_change": pressure_drop, + "transformation_method": "dae.finite_difference", + "transformation_scheme": "FORWARD", + }, + hot_side={ + "property_package": m.fs.h2_side_prop_params, + "has_holdup": False, + "dynamic": False, + "has_pressure_change": pressure_drop, + "transformation_method": "dae.finite_difference", + "transformation_scheme": "BACKWARD", + }, + shell_is_hot=True, + flow_type=HeatExchangerFlowPattern.countercurrent, + finite_elements=12, + tube_arrangement="staggered", + ) + + hx = m.fs.heat_exchanger + + hx.hot_side_inlet.flow_mol.fix(2619.7) + hx.hot_side_inlet.temperature.fix(971.6) + hx.hot_side_inlet.pressure.fix(1.2e5) + hx.hot_side_inlet.mole_frac_comp[0, "H2"].fix(0.79715) + hx.hot_side_inlet.mole_frac_comp[0, "H2O"].fix(0.20177) + hx.hot_side_inlet.mole_frac_comp[0, "Ar"].fix(0.00086358) + hx.hot_side_inlet.mole_frac_comp[0, "N2"].fix(0.00021589) + + hx.cold_side_inlet.flow_mol.fix(2619.7) + hx.cold_side_inlet.temperature.fix(446.21) + hx.cold_side_inlet.pressure.fix(1.2e5) + hx.cold_side_inlet.mole_frac_comp[0, "H2"].fix(0.36203) + hx.cold_side_inlet.mole_frac_comp[0, "H2O"].fix(0.63689) + hx.cold_side_inlet.mole_frac_comp[0, "Ar"].fix(0.00086358) + hx.cold_side_inlet.mole_frac_comp[0, "N2"].fix(0.00021589) + + hx.di_tube.fix(0.0525018) + hx.thickness_tube.fix(0.0039116) + hx.length_tube_seg.fix(4.3) + hx.number_passes.fix(12) + hx.number_columns_per_pass.fix(50) + hx.number_rows_per_pass.fix(25) + + hx.pitch_x.fix(0.1) + hx.pitch_y.fix(0.1) + hx.therm_cond_wall = 43.0 + hx.rfouling_tube = 0.0001 + hx.rfouling_shell = 0.0001 + hx.fcorrection_htc_tube.fix(1) + hx.fcorrection_htc_shell.fix(1) + + hx.cp_wall.value = 502.4 + + pp = m.fs.properties + pp.set_default_scaling("enth_mol_phase", 1e-3) + pp.set_default_scaling("pressure", 1e-5) + pp.set_default_scaling("temperature", 1e-2) + pp.set_default_scaling("flow_mol", 1e-3) + + _mf_scale = { + "H2": 1, + "H2O": 1, + "N2": 10, + "Ar": 10, + } + for comp, s in _mf_scale.items(): + pp.set_default_scaling("mole_frac_comp", s, index=comp) + pp.set_default_scaling("mole_frac_phase_comp", s, index=("Vap", comp)) + pp.set_default_scaling("flow_mol_phase_comp", s * 1e-3, index=("Vap", comp)) + + shell = hx.hot_side + tube = hx.cold_side + iscale.set_scaling_factor(shell.area, 1e-1) + iscale.set_scaling_factor(shell.heat, 1e-6) + iscale.set_scaling_factor(tube.area, 1) + iscale.set_scaling_factor(tube.heat, 1e-6) + iscale.set_scaling_factor(shell._enthalpy_flow, 1e-8) # pylint: disable=W0212 + iscale.set_scaling_factor(tube._enthalpy_flow, 1e-8) # pylint: disable=W0212 + iscale.set_scaling_factor(shell.enthalpy_flow_dx, 1e-7) + iscale.set_scaling_factor(tube.enthalpy_flow_dx, 1e-7) + iscale.set_scaling_factor(hx.heat_holdup, 1e-8) + + iscale.calculate_scaling_factors(m) + + initializer = m.fs.heat_exchanger.default_initializer( + solver="ipopt", + solver_options=optarg + ) + initializer.initialize(m.fs.heat_exchanger) + + + + +Heat Exchanger Geometry +----------------------- +=========================== =========== =================================================================================================== +Variable Index Sets Doc +=========================== =========== =================================================================================================== +``number_columns_per_pass`` None Number of columns of tube per pass +``number_rows_per_pass`` None Number of rows of tube per pass +``number_passes`` None Number of tube banks of ``nrow_tube * ncol_inlet`` tubes +``pitch_x`` None Distance between tubes parallel to shell flow, measured from center-of-tube to center-of-tube +``pitch_y`` None Distance between tubes perpendicular to shell flow, measured from center-of-tube to center-of-tube +``length_tube_seg`` None Length of tube segment perpendicular to flow in each pass +``area_flow_shell`` None Reference to flow area on shell control volume +``length_flow_shell`` None Reference to flow length on shell control volume +``area_flow_shell_min`` None Minimum flow area on shell side +``di_tube`` None Inner diameter of tubes +``thickness_tube`` None Thickness of tube wall. +``area_flow_tube`` None Reference to flow area on tube control volume +``length_flow_tube`` None Reference to flow length on tube control volume +=========================== =========== =================================================================================================== + +============================ =========== =========================================================================== +Expression Index Sets Doc +============================ =========== =========================================================================== +``nrow_tube`` None Total number of rows of tube +``do_tube`` None Outer diameter of tube (equal to ``di_tube+2*thickness_tube``) +``pitch_x_to_do`` None Ratio of ``pitch_x`` to ``do_tube`` +``pitch_y_to_do`` None Ratio of ``pitch_y`` to ``do_tube`` +``area_wall_seg`` None Total cross-sectional area of tube per pass +``total_heat_transfer_area`` None Total heat transfer area, as measured on outer surface of tubes +============================ =========== =========================================================================== + +=========================== =========== ================================================================================================= +Constraint Index Sets Doc +=========================== =========== ================================================================================================= +``length_flow_shell_eqn`` None Constrains shell flow length from control volume to equal value implied by geometry +``area_flow_shell_eqn`` None Constrains shell flow cross-sectional area from control volume to equal value implied by geometry +``area_flow_shell_min_eqn`` None Constraints ``area_flow_shell_min`` to equal value determined by geometry +``length_flow_tube_eqn`` None Constrains tube flow length from control volume to equal value implied by geometry +``area_flow_tube_eqn`` None Constrains tube flow cross-sectional area from control volume to equal value implied by geometry +=========================== =========== ================================================================================================= + +Performance Equations +----------------------- + +================================== ============ ================================================================================= +Variable Index Sets Doc +================================== ============ ================================================================================= +``fcorrection_htc_shell`` time, length Correction factor for shell convective heat transfer +``conv_heat_transfer_coeff_shell`` time, length Shell-side convective heat transfer coefficient +``temp_wall_shell`` time, length Shell-side wall temperature of tube +``temp_wall_center`` time, length Temperature at center of tube wall +``v_shell`` time, length Flow velocity on shell side through minimum area +``N_Re_shell`` time, length Reynolds number on shell side +``N_Nu_shell`` time, length Nusselt number on shell side +``heat_transfer_coeff_tube`` time, length Tube-side heat transfer coefficient +``temp_wall_tube`` time, length Tube-side wall temperature of tube +``v_tube`` time, length Flow velocity of gas in tube +``N_Re_tube`` time, length Reynolds number in tube +``N_Nu_tube`` time, length Nusselt number on tube side +================================== ============ ================================================================================= + +=========================== =========== ================================================================================= +Parameter Index Sets Doc +=========================== =========== ================================================================================= +``therm_cond_wall`` None Thermal conductivity of tube wall +``density_wall`` None Mass density of tube wall metal +``cp_wall`` None Tube wall heat capacity (mass basis) +``rfouling_shell`` None Fouling resistance on shell side +``f_arrangement`` None Adjustment factor depending on ``tube_arrangement`` in config +=========================== =========== ================================================================================= + +====================================== ============ ================================================================================= +Constraint Index Sets Doc +====================================== ============ ================================================================================= +``v_shell_eqn`` time, length Calculates velocity of flow through shell using ``area_flow_shell_min`` +``N_Re_shell_eqn`` time, length Calculates the shell-side Reynolds number +``conv_heat_transfer_coeff_shell_eqn`` time, length Calculates the shell-side convective heat transfer coefficient +``v_tube_eqn`` time, length Calculates gas velocity in tube +``N_Re_tube_eqn`` time, length Calculates the tube-side Reynolds number +``heat_transfer_coeff_tube_eqn`` time, length Calculates the tube-side heat transfer coefficient +``N_Nu_shell_eqn`` time, length Calculate the shell-side Nusselt number +``N_Nu_tube_eqn`` time, length Calculate the tube-side Nusselt number +``heat_tube_eqn`` time, length Calculates heat transfer per unit length +``heat_shell_eqn`` time, length Calculates heat transfer per unit length +``temp_wall_tube_eqn`` time, length Calculate the wall temperature of the inner tube +``temp_wall_shell_eqn`` time, length Calculate the wall temperature of the outer tube +``temp_wall_center_eqn`` time, length Overall energy balance on tube metal +====================================== ============ ================================================================================= + +====================================== ============ =================================================================================== +Expression Index Sets Doc +====================================== ============ =================================================================================== +``total_heat_transfer_coeff_shell`` time Returns ``conv_heat_transfer_coeff_shell``. Could be extended to include radiation. +``total_heat_duty`` time Created only if not ``dynamic``. Gives total heat transferred over entire exchanger +``log_mean_delta_temperature`` time Created only if not ``dynamic``. Gives the log mean temperature difference (LMTD). +``overall_heat_transfer_coefficient`` time Created only if not ``dynamic``. Calculated from total heat transfer, area, and LMTD. +====================================== ============ =================================================================================== + +Pressure Change Equations +------------------------- + +=========================== ============ ================================================================================= +Parameter Index Sets Doc +=========================== ============ ================================================================================= +``fcorrection_dp_shell`` None Correction factor for shell side pressure drop +``kloss_uturn`` None Loss coefficient of a tube u-turn +``fcorrection_dp_tube`` None Correction factor for tube side pressure drop +=========================== ============ ================================================================================= + +=========================== ============ ================================================================================= +Variable Index Sets Doc +=========================== ============ ================================================================================= +``fcorrection_dp_shell`` None Correction factor for shell side pressure drop +``friction_factor_shell`` time, length Friction factor on shell side +``friction_factor_tube`` time, length Friction factor on tube side +``deltaP_tube_friction`` time, length Change of pressure in tube due to friction +``deltaP_tube_uturn`` time, length Change of pressure in tube due to U turns +=========================== ============ ================================================================================= + +================================== ============ ================================================================================== +Constraint Index Sets Doc +================================== ============ ================================================================================== +``friction_factor_shell_eqn`` time, length Calculates the shell-side friction factor +``deltaP_shell_eqn`` time, length Sets ``deltaP_shell`` based on the friction factor and shell properties +``friction_factor_tube_eqn`` time, length Calculates the tube-side friction factor +``deltaP_tube_friction_eqn`` time, length Sets ``deltaP_tube_friction`` based on friction factor +``deltaP_tube_uturn_eqn`` time, length Sets ``deltaP_tube_uturn`` based on ``kloss_uturn`` +``deltaP_tube_eqn`` time, length Sets ``deltaP_tube`` by summing ``deltaP_tube_friction`` and ``deltaP_tube_uturn`` +================================== ============ ================================================================================== + + +Holdup Equations +---------------- + +Created when ``has_holdup=True`` in the config. + +=========================== ============ ================================================================================= +Variable Index Sets Doc +=========================== ============ ================================================================================= +``heat_holdup`` time, length Energy holdup per unit length of shell flow path +=========================== ============ ================================================================================= + +=========================== ============ ================================================================================= +Constraint Index Sets Doc +=========================== ============ ================================================================================= +``heat_holdup_eqn`` time, length Defines heat holdup in terms of geometry and physical properties +=========================== ============ ================================================================================= + +Dynamic Equations +----------------- + +Created when ``dynamic=True`` in the config. + +=========================== ============ ================================================================================= +Derivative Variable Index Sets Doc +=========================== ============ ================================================================================= +``heat_accumulation`` time, length Energy accumulation in tube wall per unit length of shell flow path per unit time +=========================== ============ ================================================================================= + +CrossFlowHeatExchanger1D Class +------------------------------ + +.. autoclass:: CrossFlowHeatExchanger1D + :members: + +CrossFlowHeatExchanger1DData Class +---------------------------------- + +.. autoclass:: CrossFlowHeatExchanger1DData + :members: + +CrossFlowHeatExchanger1DInitializer Class +----------------------------------------- + +.. autoclass:: CrossFlowHeatExchanger1DInitializer + :members: \ No newline at end of file diff --git a/docs/reference_guides/model_libraries/power_generation/unit_models/heater_1D.rst b/docs/reference_guides/model_libraries/power_generation/unit_models/heater_1D.rst new file mode 100644 index 0000000000..89b757ffc7 --- /dev/null +++ b/docs/reference_guides/model_libraries/power_generation/unit_models/heater_1D.rst @@ -0,0 +1,276 @@ +Heater1D +======== + +.. index:: + pair: idaes.models_extra.power_generation.unit_models.heater_1D;Heater1D + +.. module:: idaes.models_extra.power_generation.unit_models.heater_1D + +This model is for a gas trim heater modeled as gas being blown perpendicularly across banks of hollow tubes, +which are heated by resistive heating. Note that the ``finite_elements`` option in the control +volume config should be set to an integer factor of ``number_passes`` in order for the +discretization equations to make sense as a cross-flow heat exchanger. + +Example +------- + +.. code-block:: python + + import pyomo.environ as pyo + + from idaes.core import FlowsheetBlock + import idaes.core.util.scaling as iscale + from idaes.models.properties.modular_properties import GenericParameterBlock + from idaes.models_extra.power_generation.properties.natural_gas_PR import ( + get_prop, + EosType, + ) + from idaes.models_extra.power_generation.unit_models import Heater1D + from idaes.core.util.model_statistics import degrees_of_freedom + + optarg = { + "constr_viol_tol": 1e-8, + "nlp_scaling_method": "user-scaling", + "linear_solver": "ma57", + "OF_ma57_automatic_scaling": "yes", + "max_iter": 350, + "tol": 1e-8, + "halt_on_ampl_error": "no", + } + + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.h2_side_prop_params = GenericParameterBlock( + **get_prop(["H2", "H2O", "Ar", "N2"], {"Vap"}, eos=EosType.IDEAL), + doc="H2O + H2 gas property parameters", + ) + m.fs.heater = Heater1D( + property_package=m.fs.h2_side_prop_params, + has_holdup=True, + dynamic=False, + has_fluid_holdup=False, + has_pressure_change=pressure_drop, + finite_elements=4, + tube_arrangement="in-line", + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + ) + + heater = m.fs.heater + + heater.inlet.flow_mol.fix(5102.5) + heater.inlet.temperature.fix(938.83) + heater.inlet.pressure.fix(1.2e5) + heater.inlet.mole_frac_comp[0, "H2"].fix(0.57375) + heater.inlet.mole_frac_comp[0, "H2O"].fix(0.42517) + heater.inlet.mole_frac_comp[0, "Ar"].fix(0.00086358) + heater.inlet.mole_frac_comp[0, "N2"].fix(0.00021589) + + heater.di_tube.fix(0.0525018) + heater.thickness_tube.fix(0.0039116) + heater.pitch_x.fix(0.1) + heater.pitch_y.fix(0.1) + heater.length_tube_seg.fix(10) + heater.number_passes.fix(1) + heater.rfouling = 0.0001 + heater.fcorrection_htc_shell.fix(1) + heater.cp_wall = 502.4 + if pressure_drop: + heater.fcorrection_dp_shell.fix(1) + + heater.number_columns_per_pass.fix(40) + heater.number_rows_per_pass.fix(40) + heater.electric_heat_duty.fix(3.6504e06) + + pp = m.fs.h2_side_prop_params + pp.set_default_scaling("enth_mol_phase", 1e-3) + pp.set_default_scaling("pressure", 1e-5) + pp.set_default_scaling("temperature", 1e-2) + pp.set_default_scaling("flow_mol", 1e-3) + + _mf_scale = { + "H2": 1, + "H2O": 1, + "N2": 10, + "Ar": 10, + } + for comp, s in _mf_scale.items(): + pp.set_default_scaling("mole_frac_comp", s, index=comp) + pp.set_default_scaling("mole_frac_phase_comp", s, index=("Vap", comp)) + pp.set_default_scaling("flow_mol_phase_comp", s * 1e-3, index=("Vap", comp)) + + shell = heater.control_volume + iscale.set_scaling_factor(shell.area, 1e-1) + iscale.set_scaling_factor(shell.heat, 1e-6) + iscale.set_scaling_factor(shell.enthalpy_flow_dx, 1e-7) + iscale.set_scaling_factor(heater.heat_holdup, 1e-8) + + iscale.calculate_scaling_factors(m) + + initializer = m.fs.heat_exchanger.default_initializer( + solver="ipopt", + solver_options=optarg + ) + initializer.initialize(m.fs.heat_exchanger) + + + + +Heater Geometry +--------------- +=========================== =========== ============================================================================================= +Variable Index Sets Doc +=========================== =========== ============================================================================================= +``number_columns_per_pass`` None Number of columns of tube per pass +``number_rows_per_pass`` None Number of rows of tube per pass +``number_passes`` None Number of tube banks of ``nrow_tube * ncol_inlet`` tubes +``pitch_x`` None Distance between tubes parallel to flow, measured from center-of-tube to center-of-tube +``pitch_y`` None Distance between tubes perpendicular to flow, measured from center-of-tube to center-of-tube +``length_tube_seg`` None Length of tube segment perpendicular to flow in each pass +``area_flow_shell`` None Reference to flow area on control volume +``length_flow_shell`` None Reference to flow length on control volume +``area_flow_shell_min`` None Minimum flow area on shell side +``di_tube`` None Inner diameter of tubes +``thickness_tube`` None Thickness of tube wall. +=========================== =========== ============================================================================================= + +============================ =========== =========================================================================== +Expression Index Sets Doc +============================ =========== =========================================================================== +``nrow_tube`` None Total number of rows of tube +``do_tube`` None Outer diameter of tube (equal to ``di_tube+2*thickness_tube``) +``pitch_x_to_do`` None Ratio of ``pitch_x`` to ``do_tube`` +``pitch_y_to_do`` None Ratio of ``pitch_y`` to ``do_tube`` +``area_wall_seg`` None Total cross-sectional area of tube per pass +``total_heat_transfer_area`` None Total heat transfer area, as measured on outer surface of tubes +============================ =========== =========================================================================== + +=========================== =========== ================================================================================================= +Constraint Index Sets Doc +=========================== =========== ================================================================================================= +``length_flow_shell_eqn`` None Constrains flow length from control volume to equal value implied by geometry +``area_flow_shell_eqn`` None Constrains flow cross-sectional area from control volume to equal value implied by geometry +``area_flow_shell_min_eqn`` None Constraints ``area_flow_shell_min`` to equal value determined by geometry +=========================== =========== ================================================================================================= + +Performance Equations +----------------------- + +================================== ============ ================================================================================= +Variable Index Sets Doc +================================== ============ ================================================================================= +``electric_heat_duty`` time Electric heat duty supplied to entire heater unit +``fcorrection_htc_shell`` time, length Correction factor for convective heat transfer +``conv_heat_transfer_coeff_shell`` time, length Convective heat transfer coefficient +``temp_wall_shell`` time, length Wall temperature of tube +``temp_wall_center`` time, length Temperature at center of tube wall +``v_shell`` time, length Flow velocity through minimum area +``N_Re_shell`` time, length Reynolds number +``N_Nu_shell`` time, length Nusselt number +================================== ============ ================================================================================= + +=========================== =========== ================================================================================= +Parameter Index Sets Doc +=========================== =========== ================================================================================= +``therm_cond_wall`` None Thermal conductivity of tube wall +``density_wall`` None Mass density of tube wall metal +``cp_wall`` None Tube wall heat capacity (mass basis) +``rfouling_shell`` None Fouling resistance on shell side +``f_arrangement`` None Adjustment factor depending on ``tube_arrangement`` in config +=========================== =========== ================================================================================= + +====================================== ============ ================================================================================= +Constraint Index Sets Doc +====================================== ============ ================================================================================= +``v_shell_eqn`` time, length Calculates velocity of flow through shell using ``area_flow_shell_min`` +``N_Re_shell_eqn`` time, length Calculates the Reynolds number +``conv_heat_transfer_coeff_shell_eqn`` time, length Calculates the convective heat transfer coefficient +``N_Nu_shell_eqn`` time, length Calculate the Nusselt number +``heat_shell_eqn`` time, length Calculates heat transfer per unit length +``temp_wall_shell_eqn`` time, length Calculate the wall temperature of the outer tube +``temp_wall_center_eqn`` time, length Overall energy balance on tube metal +====================================== ============ ================================================================================= + +====================================== ============ =================================================================================== +Expression Index Sets Doc +====================================== ============ =================================================================================== +``total_heat_transfer_coeff_shell`` time Returns ``conv_heat_transfer_coeff_shell``. Could be extended to include radiation. +====================================== ============ =================================================================================== + +Pressure Change Equations +------------------------- + +=========================== ============ ================================================================================= +Parameter Index Sets Doc +=========================== ============ ================================================================================= +``fcorrection_dp_shell`` None Correction factor for pressure drop +=========================== ============ ================================================================================= + +=========================== ============ ================================================================================= +Variable Index Sets Doc +=========================== ============ ================================================================================= +``fcorrection_dp_shell`` None Correction factor for pressure drop +``friction_factor_shell`` time, length Friction factor +=========================== ============ ================================================================================= + +================================== ============ ================================================================================= +Constraint Index Sets Doc +================================== ============ ================================================================================= +``friction_factor_shell_eqn`` time, length Calculates the friction factor +``deltaP_shell_eqn`` time, length Sets ``deltaP_shell`` based on the friction factor and physical properties +================================== ============ ================================================================================= + + +Holdup Equations +---------------- + +Created when ``has_holdup=True`` in the config. + +=========================== ============ ================================================================================= +Variable Index Sets Doc +=========================== ============ ================================================================================= +``heat_holdup`` time, length Energy holdup per unit length of flow path +=========================== ============ ================================================================================= + +=========================== ============ ================================================================================= +Constraint Index Sets Doc +=========================== ============ ================================================================================= +``heat_holdup_eqn`` time, length Defines heat holdup in terms of geometry and physical properties +=========================== ============ ================================================================================= + +Dynamic Equations +----------------- + +Created when ``dynamic=True`` in the config. + +=========================== ============ ================================================================================= +Derivative Variable Index Sets Doc +=========================== ============ ================================================================================= +``heat_accumulation`` time, length Energy accumulation in tube wall per unit length of shell flow path per unit time +=========================== ============ ================================================================================= + + +Initialization +-------------- + +A simple initialization method that first initializes the control volume without heat transfer, +then adds heat transfer in and solves it again, then finally solves the entire model. + + +Heater1D Class +-------------- + +.. autoclass:: Heater1D + :members: + +Heater1DData Class +------------------ + +.. autoclass:: Heater1DData + :members: + +Heater1DInitializer Class +------------------------- + +.. autoclass:: Heater1DInitializer + :members: \ No newline at end of file diff --git a/docs/reference_guides/model_libraries/power_generation/unit_models/index.rst b/docs/reference_guides/model_libraries/power_generation/unit_models/index.rst index 1a10eb0952..b226b2817a 100644 --- a/docs/reference_guides/model_libraries/power_generation/unit_models/index.rst +++ b/docs/reference_guides/model_libraries/power_generation/unit_models/index.rst @@ -26,3 +26,5 @@ Unit Models waterpipe boiler_heat_exchanger_3streams feedwater_heater_0D_dynamic + cross_flow_heat_exchanger_1D + heater_1D diff --git a/idaes/models_extra/power_generation/unit_models/__init__.py b/idaes/models_extra/power_generation/unit_models/__init__.py index abd6fbb5d6..92cdde53a4 100644 --- a/idaes/models_extra/power_generation/unit_models/__init__.py +++ b/idaes/models_extra/power_generation/unit_models/__init__.py @@ -15,10 +15,15 @@ from .boiler_fireside import BoilerFireside from .boiler_heat_exchanger import BoilerHeatExchanger from .boiler_heat_exchanger_2D import HeatExchangerCrossFlow2D_Header +from .cross_flow_heat_exchanger_1D import ( + CrossFlowHeatExchanger1D, + CrossFlowHeatExchanger1DInitializer, +) from .downcomer import Downcomer from .drum import Drum from .drum1D import Drum1D from .feedwater_heater_0D_dynamic import FWH0DDynamic +from .heater_1D import Heater1D, Heater1DInitializer from .heat_exchanger_3streams import HeatExchangerWith3Streams from .steamheater import SteamHeater from .waterpipe import WaterPipe diff --git a/idaes/models_extra/power_generation/unit_models/cross_flow_heat_exchanger_1D.py b/idaes/models_extra/power_generation/unit_models/cross_flow_heat_exchanger_1D.py new file mode 100644 index 0000000000..d4615698d5 --- /dev/null +++ b/idaes/models_extra/power_generation/unit_models/cross_flow_heat_exchanger_1D.py @@ -0,0 +1,1058 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, 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.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +1-D Cross Flow Heat Exchanger Model With Wall Temperatures + +Discretization based on tube rows +""" + + +# Import Pyomo libraries +from pyomo.environ import ( + assert_optimal_termination, + Block, + value, + Var, + log, + Param, + Reference, + units as pyunits, +) +from pyomo.common.config import ConfigValue, In, Bool +from pyomo.network import Port + +# Import IDAES cores +from idaes.core import declare_process_block_class +from idaes.core.util.constants import Constants as const +import idaes.core.util.scaling as iscale +from idaes.core.solvers import get_solver +from idaes.core.util.exceptions import InitializationError +from idaes.core.util.misc import add_object_reference +from idaes.core.util.exceptions import ConfigurationError, BurntToast +import idaes.logger as idaeslog +from idaes.core.util.tables import create_stream_table_dataframe +from idaes.models.unit_models.heat_exchanger import ( + HeatExchangerFlowPattern, +) +from idaes.models.unit_models.heat_exchanger_1D import HeatExchanger1DData +from idaes.models_extra.power_generation.unit_models import heat_exchanger_common +from idaes.core.initialization import SingleControlVolumeUnitInitializer + + +__author__ = "Jinliang Ma, Douglas Allan" + + +class CrossFlowHeatExchanger1DInitializer(SingleControlVolumeUnitInitializer): + """ + Initializer for Cross Flow Heat Exchanger 1D units. + + First, the shell and tube control volumes are initialized without heat transfer. Next + the total possible heat transfer between streams is estimated based on heat capacity, + flow rate, and inlet/outlet temperatures. The actual temperature change is set to be + half the theoretical maximum, and the shell and tube are initialized with linear + temperature profiles. Finally, temperatures besides the inlets are unfixed and + the performance equations are activated before a full solve of the system model. + """ + + def initialize_main_model( + self, + model: Block, + copy_inlet_state: bool = False, + ): + """ + Initialization routine for the main Cross Flow Heat Exchanger 1D model + (as opposed to submodels like costing, which presently do not exist). + + 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. + + """ + # 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" + ) + + solver_obj = get_solver(self.config.solver, self.config.solver_options) + + hot_side = model.hot_side + cold_side = model.cold_side + t0 = model.flowsheet().time.first() + if not ( + "temperature" in hot_side.properties[t0, 0].define_state_vars().keys() + and "temperature" in cold_side.properties[t0, 0].define_state_vars().keys() + ): + raise NotImplementedError( + "Presently, initialization of the CrossFlowHeatExchanger1D requires " + "temperature to be a state variable of both hot side and cold side " + "property packages. Extension to enth_mol or enth_mass as state variables " + "is straightforward---feel free to open a pull request implementing it." + ) + + hot_units = model.config.hot_side.property_package.get_metadata().derived_units + cold_units = ( + model.config.cold_side.property_package.get_metadata().derived_units + ) + + if model.config.shell_is_hot: + shell = model.hot_side + tube = model.cold_side + shell_has_pressure_change = model.config.hot_side.has_pressure_change + tube_has_pressure_change = model.config.cold_side.has_pressure_change + shell_units = ( + model.config.hot_side.property_package.get_metadata().derived_units + ) + tube_units = ( + model.config.cold_side.property_package.get_metadata().derived_units + ) + else: + shell = model.cold_side + tube = model.hot_side + shell_has_pressure_change = model.config.cold_side.has_pressure_change + tube_has_pressure_change = model.config.hot_side.has_pressure_change + shell_units = ( + model.config.cold_side.property_package.get_metadata().derived_units + ) + tube_units = ( + model.config.hot_side.property_package.get_metadata().derived_units + ) + + # Trigger creation of cp for use in future initialization + # Important to do before initializing property packages in + # case it is implemented as Var-Constraint pair instead of + # an Expression + hot_side.properties[t0, 0].cp_mol # pylint: disable=pointless-statement + cold_side.properties[t0, 0].cp_mol # pylint: disable=pointless-statement + + # --------------------------------------------------------------------- + # Initialize shell block + self.initialize_control_volume(tube, copy_inlet_state) + self.initialize_control_volume(shell, copy_inlet_state) + + init_log.info_high("Initialization Step 1 Complete.") + + # Set tube thermal conductivity to a small value to avoid IPOPT unable to solve initially + therm_cond_wall_save = model.therm_cond_wall.value + model.therm_cond_wall.set_value(0.05) + # In Step 2, fix tube metal temperatures fix fluid state variables (enthalpy/temperature and pressure) + # calculate maximum heat duty assuming infinite area and use half of the maximum duty as initial guess to calculate outlet temperature + + for t in model.flowsheet().config.time: + mcp_hot_side = value( + pyunits.convert( + hot_side.properties[t, 0].flow_mol + * hot_side.properties[t, 0].cp_mol, + to_units=shell_units["power"] / shell_units["temperature"], + ) + ) + T_in_hot_side = value( + pyunits.convert( + hot_side.properties[t, 0].temperature, + to_units=shell_units["temperature"], + ) + ) + P_in_hot_side = value(hot_side.properties[t, 0].pressure) + if model.config.flow_type == HeatExchangerFlowPattern.cocurrent: + mcp_cold_side = value( + pyunits.convert( + cold_side.properties[t, 0].flow_mol + * cold_side.properties[t, 0].cp_mol, + to_units=shell_units["power"] / shell_units["temperature"], + ) + ) + T_in_cold_side = value( + pyunits.convert( + cold_side.properties[t, 0].temperature, + to_units=shell_units["temperature"], + ) + ) + P_in_cold_side = value(cold_side.properties[t, 0].pressure) + + T_out_max = ( + mcp_cold_side * T_in_cold_side + mcp_hot_side * T_in_hot_side + ) / (mcp_cold_side + mcp_hot_side) + + q_guess = mcp_cold_side * (T_out_max - T_in_cold_side) / 2 + + temp_out_cold_side_guess = T_in_cold_side + q_guess / mcp_cold_side + + cold_side.properties[t, 1].temperature.fix( + pyunits.convert_value( + temp_out_cold_side_guess, + from_units=shell_units["temperature"], + to_units=cold_units["temperature"], + ) + ) + + temp_out_hot_side_guess = T_in_cold_side - q_guess / mcp_hot_side + hot_side.properties[t, 1].temperature.fix( + pyunits.convert_value( + temp_out_hot_side_guess, + from_units=shell_units["temperature"], + to_units=hot_units["temperature"], + ) + ) + + elif model.config.flow_type == HeatExchangerFlowPattern.countercurrent: + mcp_cold_side = value( + pyunits.convert( + cold_side.properties[t, 1].flow_mol + * cold_side.properties[t, 1].cp_mol, + to_units=shell_units["power"] / shell_units["temperature"], + ) + ) + T_in_cold_side = value( + pyunits.convert( + cold_side.properties[t, 1].temperature, + to_units=shell_units["temperature"], + ) + ) + P_in_cold_side = value(cold_side.properties[t, 1].pressure) + + if mcp_cold_side < mcp_hot_side: + q_guess = mcp_cold_side * (T_in_hot_side - T_in_cold_side) / 2 + else: + q_guess = mcp_hot_side * (T_in_hot_side - T_in_cold_side) / 2 + + temp_out_cold_side_guess = T_in_cold_side + q_guess / mcp_cold_side + cold_side.properties[t, 0].temperature.fix( + pyunits.convert_value( + temp_out_cold_side_guess, + from_units=shell_units["temperature"], + to_units=cold_units["temperature"], + ) + ) + + temp_out_hot_side_guess = T_in_hot_side - q_guess / mcp_hot_side + hot_side.properties[t, 1].temperature.fix( + pyunits.convert_value( + temp_out_hot_side_guess, + from_units=shell_units["temperature"], + to_units=hot_units["temperature"], + ) + ) + + else: + raise BurntToast( + "HeatExchangerFlowPattern should be limited to cocurrent " + "or countercurrent flow by parent model. Please open an " + "issue on the IDAES Github so this error can be fixed." + ) + + for z in cold_side.length_domain: + hot_side.properties[t, z].temperature.fix( + value( + (1 - z) * hot_side.properties[t, 0].temperature + + z * hot_side.properties[t, 1].temperature + ) + ) + cold_side.properties[t, z].temperature.fix( + value( + (1 - z) * cold_side.properties[t, 0].temperature + + z * cold_side.properties[t, 1].temperature + ) + ) + model.temp_wall_center[t, z].fix( + value( + pyunits.convert( + hot_side.properties[t, z].temperature, + to_units=shell_units["temperature"], + ) + + pyunits.convert( + cold_side.properties[t, z].temperature, + to_units=shell_units["temperature"], + ) + ) + / 2 + ) + + model.temp_wall_shell[t, z].set_value( + model.temp_wall_center[t, z].value + ) + model.temp_wall_tube[t, z].set_value( + pyunits.convert_value( + model.temp_wall_center[t, z].value, + from_units=shell_units["temperature"], + to_units=tube_units["temperature"], + ) + ) + + if model.config.cold_side.has_pressure_change: + cold_side.properties[t, z].pressure.fix(P_in_cold_side) + if model.config.hot_side.has_pressure_change: + hot_side.properties[t, z].pressure.fix(P_in_hot_side) + + model.temp_wall_center_eqn.deactivate() + if tube_has_pressure_change: + model.deltaP_tube_eqn.deactivate() + if shell_has_pressure_change: + model.deltaP_shell_eqn.deactivate() + model.heat_tube_eqn.deactivate() + model.heat_shell_eqn.deactivate() + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(model, tee=slc.tee) + try: + assert_optimal_termination(res) + except AssertionError: + raise InitializationError("Initialization solve failed.") + init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res))) + + # In Step 3, unfix fluid state variables (enthalpy/temperature and pressure) + # keep the inlet state variables fixed, otherwise, the degree of freedom > 0 + hot_side.properties[:, :].temperature.unfix() + hot_side.properties[:, :].pressure.unfix() + hot_side.properties[:, 0].temperature.fix() + hot_side.properties[:, 0].pressure.fix() + + cold_side.properties[:, :].temperature.unfix() + cold_side.properties[:, :].pressure.unfix() + if model.config.flow_type == HeatExchangerFlowPattern.cocurrent: + cold_side.properties[:, 0].temperature.fix() + cold_side.properties[:, 0].pressure.fix() + elif model.config.flow_type == HeatExchangerFlowPattern.countercurrent: + cold_side.properties[:, 1].temperature.fix() + cold_side.properties[:, 1].pressure.fix() + else: + raise BurntToast( + "HeatExchangerFlowPattern should be limited to cocurrent " + "or countercurrent flow by parent model. Please open an " + "issue on the IDAES Github so this error can be fixed." + ) + + if tube_has_pressure_change: + model.deltaP_tube_eqn.activate() + if shell_has_pressure_change: + model.deltaP_shell_eqn.activate() + model.heat_tube_eqn.activate() + model.heat_shell_eqn.activate() + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(model, tee=slc.tee) + try: + assert_optimal_termination(res) + except AssertionError: + raise InitializationError("Initialization solve failed.") + + init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res))) + + model.temp_wall_center.unfix() + model.temp_wall_center_eqn.activate() + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(model, tee=slc.tee) + try: + assert_optimal_termination(res) + except AssertionError: + raise InitializationError("Initialization solve failed.") + + init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res))) + + # set the wall thermal conductivity back to the user specified value + model.therm_cond_wall.set_value(therm_cond_wall_save) + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(model, tee=slc.tee) + init_log.info_high("Initialization Step 5 {}.".format(idaeslog.condition(res))) + init_log.info("Initialization Complete.") + return res + + +@declare_process_block_class("CrossFlowHeatExchanger1D") +class CrossFlowHeatExchanger1DData(HeatExchanger1DData): + """Standard Cross Flow Heat Exchanger Unit Model Class.""" + + default_initializer = CrossFlowHeatExchanger1DInitializer + + CONFIG = HeatExchanger1DData.CONFIG() + CONFIG.declare( + "shell_is_hot", + ConfigValue( + default=True, + domain=Bool, + description="Shell side contains hot fluid", + doc="""Boolean flag indicating whether shell side contains hot fluid (default=True). + If True, shell side will be the hot_side, if False shell side will be cold_side.""", + ), + ) + CONFIG.declare( + "tube_arrangement", + ConfigValue( + default="in-line", + domain=In(["in-line", "staggered"]), + description="tube configuration", + doc="tube arrangement could be in-line or staggered", + ), + ) + + def _process_config(self): + # Copy and pasted from ShellAndTube1D + super()._process_config() + + # Check for custom names, and if not present assign defaults + if self.config.hot_side_name is None: + if self.config.shell_is_hot: + self.config.hot_side_name = "Shell" + else: + self.config.hot_side_name = "Tube" + + if self.config.cold_side_name is None: + if self.config.shell_is_hot: + self.config.cold_side_name = "Tube" + else: + self.config.cold_side_name = "Shell" + + def build(self): + """ + Begin building model (pre-DAE transformation). + + Args: + None + + Returns: + None + """ + # Call HeatExchanger1DData build to make common components + super().build() + + # The HeatExchanger1DData equates the heat lost by the hot side and heat gained by the cold side. + # That equation is deleted here because heat can accumulate in the wall. + self.del_component(self.heat_conservation) + + # Create aliases for ports + if self.config.shell_is_hot: + self.shell_inlet = Port(extends=self.hot_side_inlet) + self.shell_outlet = Port(extends=self.hot_side_outlet) + self.tube_inlet = Port(extends=self.cold_side_inlet) + self.tube_outlet = Port(extends=self.cold_side_outlet) + else: + self.shell_inlet = Port(extends=self.cold_side_inlet) + self.shell_outlet = Port(extends=self.cold_side_outlet) + self.tube_inlet = Port(extends=self.hot_side_inlet) + self.tube_outlet = Port(extends=self.hot_side_outlet) + + def _make_geometry(self): + """ + Constraints for unit model. + + Args: + None + + Returns: + None + """ + if self.config.shell_is_hot: + shell = self.hot_side + tube = self.cold_side + shell_units = ( + self.config.hot_side.property_package.get_metadata().derived_units + ) + else: + shell = self.cold_side + tube = self.hot_side + shell_units = ( + self.config.cold_side.property_package.get_metadata().derived_units + ) + # Use add_object_reference for scalar vars + add_object_reference(self, "area_flow_shell", shell.area) + add_object_reference(self, "area_flow_tube", tube.area) + add_object_reference(self, "length_flow_shell", shell.length) + add_object_reference(self, "length_flow_tube", tube.length) + + heat_exchanger_common.make_geometry_common(self, shell_units=shell_units) + + # Important that these values about tube geometry are in shell units! + @self.Constraint(doc="Length of tube side flow") + def length_flow_tube_eqn(b): + return ( + pyunits.convert(b.length_flow_tube, to_units=shell_units["length"]) + == b.number_passes * b.length_tube_seg + ) + + @self.Constraint(doc="Total area of tube flow") + def area_flow_tube_eqn(b): + return ( + b.area_flow_tube + == 0.25 + * const.pi + * b.di_tube**2.0 + * b.number_columns_per_pass + * b.number_rows_per_pass + ) + + def _make_performance(self): + """ + Constraints for unit model. + + Args: + None + + Returns: + None + """ + if self.config.shell_is_hot: + shell = self.hot_side + tube = self.cold_side + shell_units = ( + self.config.hot_side.property_package.get_metadata().derived_units + ) + tube_units = ( + self.config.cold_side.property_package.get_metadata().derived_units + ) + else: + shell = self.cold_side + tube = self.hot_side + shell_units = ( + self.config.cold_side.property_package.get_metadata().derived_units + ) + tube_units = ( + self.config.hot_side.property_package.get_metadata().derived_units + ) + if ( + len(self.config.hot_side.property_package.phase_list) != 1 + or len(self.config.cold_side.property_package.phase_list) != 1 + ): + raise ConfigurationError( + "The CrossFlowHeatExchanger1D model is valid only for property packages " + f"with a single phase. Found {len(self.config.hot_side.property_package.phase_list)} " + f"phases on the hot side and {len(self.config.cold_side.property_package.phase_list)} " + "phases on the cold side." + ) + + p_hot = self.config.hot_side.property_package.phase_list.at(1) + pobj_hot = self.config.hot_side.property_package.get_phase(p_hot) + p_cold = self.config.cold_side.property_package.phase_list.at(1) + pobj_cold = self.config.cold_side.property_package.get_phase(p_cold) + if not pobj_hot.is_vapor_phase(): + raise ConfigurationError( + "The CrossFlowHeatExchanger1D model is valid only for property packages " + "whose single phase is a vapor phase. The hot side phase is not a vapor phase." + ) + if not pobj_cold.is_vapor_phase(): + raise ConfigurationError( + "The CrossFlowHeatExchanger1D model is valid only for property packages " + "whose single phase is a vapor phase. The cold side phase is not a vapor phase." + ) + + self.heat_tube = Reference(tube.heat) + self.heat_shell = Reference(shell.heat) + + shell_has_pressure_change = False + tube_has_pressure_change = False + + if self.config.cold_side.has_pressure_change: + if self.config.shell_is_hot: + self.deltaP_tube = Reference(tube.deltaP) + tube_has_pressure_change = True + else: + self.deltaP_shell = Reference(shell.deltaP) + shell_has_pressure_change = True + if self.config.hot_side.has_pressure_change: + if self.config.shell_is_hot: + self.deltaP_shell = Reference(shell.deltaP) + shell_has_pressure_change = True + else: + self.deltaP_tube = Reference(tube.deltaP) + tube_has_pressure_change = True + + heat_exchanger_common.make_performance_common( + self, + shell=shell, + shell_units=shell_units, + shell_has_pressure_change=shell_has_pressure_change, + make_reynolds=True, + make_nusselt=True, + ) + + self.heat_transfer_coeff_tube = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=100.0, + units=tube_units["heat_transfer_coefficient"], + doc="tube side convective heat transfer coefficient", + ) + + # Heat transfer resistance due to the fouling on tube side + self.rfouling_tube = Param( + initialize=0.0, + mutable=True, + units=1 / tube_units["heat_transfer_coefficient"], + doc="fouling resistance on tube side", + ) + # Correction factor for convective heat transfer coefficient on tube side + self.fcorrection_htc_tube = Var( + initialize=1.0, doc="correction factor for convective HTC on tube side" + ) + # Correction factor for tube side pressure drop due to friction + if tube_has_pressure_change: + # Loss coefficient for a 180 degree bend (u-turn), usually related to radius to inside diameter ratio + self.kloss_uturn = Param( + initialize=0.5, mutable=True, doc="loss coefficient of a tube u-turn" + ) + self.fcorrection_dp_tube = Var( + initialize=1.0, doc="correction factor for tube side pressure drop" + ) + + # Boundary wall temperature on tube side + self.temp_wall_tube = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=500, + units=tube_units[ + "temperature" + ], # Want to be in shell units for consistency in equations + doc="boundary wall temperature on tube side", + ) + # Tube side heat transfer coefficient and pressure drop + # ----------------------------------------------------- + # Velocity on tube side + self.v_tube = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=1.0, + units=tube_units["velocity"], + doc="velocity on tube side", + ) + + # Reynalds number on tube side + self.N_Re_tube = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=10000.0, + units=pyunits.dimensionless, + doc="Reynolds number on tube side", + bounds=(1e-7, None), + ) + if tube_has_pressure_change: + # Friction factor on tube side + self.friction_factor_tube = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=1.0, + doc="friction factor on tube side", + ) + + # Pressure drop due to friction on tube side + self.deltaP_tube_friction = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=-10.0, + units=tube_units["pressure"] / tube_units["length"], + doc="pressure drop due to friction on tube side", + ) + + # Pressure drop due to 180 degree turn on tube side + self.deltaP_tube_uturn = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=-10.0, + units=tube_units["pressure"] / tube_units["length"], + doc="pressure drop due to u-turn on tube side", + ) + # Nusselt number on tube side + self.N_Nu_tube = Var( + self.flowsheet().config.time, + tube.length_domain, + initialize=1, + units=pyunits.dimensionless, + doc="Nusselts number on tube side", + bounds=(1e-7, None), + ) + + # Velocity equation + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="tube side velocity equation", + ) + def v_tube_eqn(b, t, x): + return ( + b.v_tube[t, x] + * pyunits.convert(b.area_flow_tube, to_units=tube_units["area"]) + * tube.properties[t, x].dens_mol_phase["Vap"] + == tube.properties[t, x].flow_mol + ) + + # Reynolds number + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="Reynolds number equation on tube side", + ) + def N_Re_tube_eqn(b, t, x): + return ( + b.N_Re_tube[t, x] * tube.properties[t, x].visc_d_phase["Vap"] + == pyunits.convert(b.di_tube, to_units=tube_units["length"]) + * b.v_tube[t, x] + * tube.properties[t, x].dens_mol_phase["Vap"] + * tube.properties[t, x].mw + ) + + if tube_has_pressure_change: + # Friction factor + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="Darcy friction factor on tube side", + ) + def friction_factor_tube_eqn(b, t, x): + return ( + b.friction_factor_tube[t, x] * b.N_Re_tube[t, x] ** 0.25 + == 0.3164 * b.fcorrection_dp_tube + ) + + # Pressure drop due to friction per tube length + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="pressure drop due to friction per tube length", + ) + def deltaP_tube_friction_eqn(b, t, x): + return ( + b.deltaP_tube_friction[t, x] + * pyunits.convert(b.di_tube, to_units=tube_units["length"]) + == -0.5 + * tube.properties[t, x].dens_mass_phase["Vap"] + * b.v_tube[t, x] ** 2 + * b.friction_factor_tube[t, x] + ) + + # Pressure drop due to u-turn + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="pressure drop due to u-turn on tube side", + ) + def deltaP_tube_uturn_eqn(b, t, x): + return ( + b.deltaP_tube_uturn[t, x] + * pyunits.convert(b.length_tube_seg, to_units=tube_units["length"]) + == -0.5 + * tube.properties[t, x].dens_mass_phase["Vap"] + * b.v_tube[t, x] ** 2 + * b.kloss_uturn + ) + + # Total pressure drop on tube side + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="total pressure drop on tube side", + ) + def deltaP_tube_eqn(b, t, x): + return b.deltaP_tube[t, x] == ( + b.deltaP_tube_friction[t, x] + b.deltaP_tube_uturn[t, x] + ) + + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="convective heat transfer coefficient equation on tube side", + ) + def heat_transfer_coeff_tube_eqn(b, t, x): + return ( + b.heat_transfer_coeff_tube[t, x] * b.di_tube + == b.N_Nu_tube[t, x] + * tube.properties[t, x].therm_cond_phase["Vap"] + * b.fcorrection_htc_tube + ) + + # Nusselts number + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="Nusselts number equation on tube side", + ) + def N_Nu_tube_eqn(b, t, x): + return ( + b.N_Nu_tube[t, x] + == 0.023 + * b.N_Re_tube[t, x] ** 0.8 + * tube.properties[t, x].prandtl_number_phase["Vap"] ** 0.4 + ) + + @self.Constraint( + self.flowsheet().config.time, + shell.length_domain, + doc="Nusselts number equation on shell side", + ) + def N_Nu_shell_eqn(b, t, x): + return ( + b.N_Nu_shell[t, x] + == b.f_arrangement + * 0.33 + * b.N_Re_shell[t, x] ** 0.6 + * shell.properties[t, x].prandtl_number_phase["Vap"] ** 0.333333 + ) + + # Energy balance with tube wall + # ------------------------------------ + # Heat to wall per length on tube side + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="heat per length on tube side", + ) + def heat_tube_eqn(b, t, x): + return b.heat_tube[t, x] == ( + b.heat_transfer_coeff_tube[t, x] + * const.pi + * pyunits.convert(b.di_tube, to_units=tube_units["length"]) + * b.number_rows_per_pass + * b.number_columns_per_pass + * (b.temp_wall_tube[t, x] - tube.properties[t, x].temperature) + ) + + # Heat to wall per length on shell side + @self.Constraint( + self.flowsheet().config.time, + shell.length_domain, + doc="heat per length on shell side", + ) + def heat_shell_eqn(b, t, x): + return b.heat_shell[t, x] * b.length_flow_shell == pyunits.convert( + b.length_flow_tube, to_units=shell_units["length"] + ) * b.total_heat_transfer_coeff_shell[ + t, x + ] * const.pi * b.do_tube * b.number_rows_per_pass * b.number_columns_per_pass * ( + b.temp_wall_shell[t, x] - shell.properties[t, x].temperature + ) + + # Tube side wall temperature + # FIXME replace with deviation variables + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="tube side wall temperature", + ) + def temp_wall_tube_eqn(b, t, x): + return b.heat_transfer_coeff_tube[t, x] * ( + tube.properties[t, x].temperature - b.temp_wall_tube[t, x] + ) * ( + pyunits.convert( + b.thickness_tube / 2 / b.therm_cond_wall, + to_units=1 / tube_units["heat_transfer_coefficient"], + ) + + b.rfouling_tube + ) == b.temp_wall_tube[ + t, x + ] - pyunits.convert( + b.temp_wall_center[t, x], to_units=tube_units["temperature"] + ) + + # Shell side wall temperature + @self.Constraint( + self.flowsheet().config.time, + shell.length_domain, + doc="shell side wall temperature", + ) + def temp_wall_shell_eqn(b, t, x): + return ( + b.total_heat_transfer_coeff_shell[t, x] + * (shell.properties[t, x].temperature - b.temp_wall_shell[t, x]) + * (b.thickness_tube / 2 / b.therm_cond_wall + b.rfouling_shell) + == b.temp_wall_shell[t, x] - b.temp_wall_center[t, x] + ) + + def heat_accumulation_term(b, t, x): + return b.heat_accumulation[t, x] if b.config.dynamic else 0 + + # Center point wall temperature based on energy balance for tube wall heat holdup + @self.Constraint( + self.flowsheet().config.time, + tube.length_domain, + doc="center point wall temperature", + ) + def temp_wall_center_eqn(b, t, x): + # heat_shell and heat_tube are positive when heat flows into those + # control volumes (and out of the wall), hence the negative sign + # on heat_accumulation_term + return -heat_accumulation_term(b, t, x) == ( + b.heat_shell[t, x] + * b.length_flow_shell + / pyunits.convert(b.length_flow_tube, to_units=shell_units["length"]) + + pyunits.convert( + b.heat_tube[t, x], + to_units=shell_units["power"] / shell_units["length"], + ) + ) + + if not self.config.dynamic: + z0 = shell.length_domain.first() + z1 = shell.length_domain.last() + + @self.Expression(self.flowsheet().config.time) + def total_heat_duty(b, t): + enth_in = b.hot_side.properties[t, z0].get_enthalpy_flow_terms(p_hot) + enth_out = b.hot_side.properties[t, z1].get_enthalpy_flow_terms(p_hot) + + return pyunits.convert( + enth_in - enth_out, # Hot side loses enthalpy + to_units=shell_units["power"], # Hot side isn't always the shell + ) + + @self.Expression(self.flowsheet().config.time) + def log_mean_delta_temperature(b, t): + dT0 = b.hot_side.properties[t, z0].temperature - pyunits.convert( + b.cold_side.properties[t, z0].temperature, + to_units=shell_units["temperature"], + ) + dT1 = b.hot_side.properties[t, z1].temperature - pyunits.convert( + b.cold_side.properties[t, z1].temperature, + to_units=shell_units["temperature"], + ) + return (dT0 - dT1) / log(dT0 / dT1) + + @self.Expression(self.flowsheet().config.time) + def overall_heat_transfer_coefficient(b, t): + return b.total_heat_duty[t] / ( + b.total_heat_transfer_area * b.log_mean_delta_temperature[t] + ) + + def calculate_scaling_factors(self): + def gsf(obj): + return iscale.get_scaling_factor(obj, default=1, warning=True) + + def ssf(obj, sf): + iscale.set_scaling_factor(obj, sf, overwrite=False) + + def cst(con, sf): + iscale.constraint_scaling_transform(con, sf, overwrite=False) + + sgsf = iscale.set_and_get_scaling_factor + if self.config.shell_is_hot: + shell = self.hot_side + tube = self.cold_side + shell_units = ( + self.config.hot_side.property_package.get_metadata().derived_units + ) + else: + shell = self.cold_side + tube = self.hot_side + shell_units = ( + self.config.cold_side.property_package.get_metadata().derived_units + ) + + shell_has_pressure_change = hasattr(self, "deltaP_shell") + + heat_exchanger_common.scale_common( + self, + shell, + shell_has_pressure_change, + make_reynolds=True, + make_nusselt=True, + ) + + sf_di_tube = iscale.get_scaling_factor( + self.do_tube, default=1 / value(self.di_tube) + ) + + sf_area_per_length_shell = value( + self.length_flow_shell + / ( + pyunits.convert(self.length_flow_tube, to_units=shell_units["length"]) + * const.pi + * self.do_tube + * self.number_rows_per_pass + * self.number_columns_per_pass + ) + ) + + for t in self.flowsheet().time: + for z in shell.length_domain: + # FIXME get better scaling later + ssf(self.v_tube[t, z], 1 / 20) + sf_flow_mol_tube = gsf(tube.properties[t, z].flow_mol) + + cst(self.v_tube_eqn[t, z], sf_flow_mol_tube) + + # FIXME should get scaling of N_Re from defining eqn + sf_N_Re_tube = sgsf(self.N_Re_tube[t, z], 1e-4) + + sf_visc_d_tube = gsf(tube.properties[t, z].visc_d_phase["Vap"]) + cst(self.N_Re_tube_eqn[t, z], sf_N_Re_tube * sf_visc_d_tube) + + sf_k_tube = gsf(tube.properties[t, z].therm_cond_phase["Vap"]) + + sf_N_Nu_tube = sgsf(self.N_Nu_tube[t, z], 1 / 0.023 * sf_N_Re_tube**0.8) + cst(self.N_Nu_tube_eqn[t, z], sf_N_Nu_tube) + + sf_heat_transfer_coeff_tube = sgsf( + self.heat_transfer_coeff_tube[t, z], + sf_N_Nu_tube * sf_k_tube / sf_di_tube, + ) + cst( + self.heat_transfer_coeff_tube_eqn[t, z], + sf_heat_transfer_coeff_tube * sf_di_tube, + ) + + sf_T_tube = gsf(tube.properties[t, z].temperature) + ssf(self.temp_wall_tube[t, z], sf_T_tube) + cst(self.temp_wall_tube_eqn[t, z], sf_T_tube) + + sf_Q_tube = gsf(tube.heat[t, z]) + cst(self.heat_tube_eqn[t, z], sf_Q_tube) + + sf_T_shell = gsf(shell.properties[t, z].temperature) + ssf(self.temp_wall_shell[t, z], sf_T_shell) + cst(self.temp_wall_shell_eqn[t, z], sf_T_shell) + + sf_conv_heat_transfer_coeff_shell = gsf( + self.conv_heat_transfer_coeff_shell[t, z] + ) + s_Q_shell = sgsf( + shell.heat[t, z], + sf_conv_heat_transfer_coeff_shell + * sf_area_per_length_shell + * sf_T_shell, + ) + cst( + self.heat_shell_eqn[t, z], s_Q_shell * value(self.length_flow_shell) + ) + # Geometric mean is overkill for most reasonable cases, but it mitigates + # damage done when one stream has an unset scaling factor + ssf(self.temp_wall_center[t, z], (sf_T_shell * sf_T_tube) ** 0.5) + cst(self.temp_wall_center_eqn[t, z], (sf_Q_tube * s_Q_shell) ** 0.5) + + def _get_performance_contents(self, time_point=0): + var_dict = {} + + expr_dict = {} + expr_dict["HX Area"] = self.total_heat_transfer_area + expr_dict["Delta T Driving"] = self.log_mean_delta_temperature[time_point] + expr_dict["Total Heat Duty"] = self.total_heat_duty[time_point] + expr_dict["Average HX Coefficient"] = self.overall_heat_transfer_coefficient[ + time_point + ] + + return {"vars": var_dict, "exprs": expr_dict} + + def _get_stream_table_contents(self, time_point=0): + return create_stream_table_dataframe( + { + "Hot Inlet": self.shell_inlet, + "Hot Outlet": self.shell_outlet, + "Cold Inlet": self.tube_inlet, + "Cold Outlet": self.tube_outlet, + }, + time_point=time_point, + ) diff --git a/idaes/models_extra/power_generation/unit_models/heat_exchanger_common.py b/idaes/models_extra/power_generation/unit_models/heat_exchanger_common.py new file mode 100644 index 0000000000..fd2932853b --- /dev/null +++ b/idaes/models_extra/power_generation/unit_models/heat_exchanger_common.py @@ -0,0 +1,539 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, 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.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +Methods shared by the CrossFlowHeatExchanger1D and Heater1D models. +""" +# Presently the tube methods are not shared. I'm not sure why I chose to extract +# them here, but I don't want to change things this late into development. --Doug + +# Import Pyomo libraries +from pyomo.environ import ( + Var, + Param, + value, + units as pyunits, +) +from pyomo.util.calc_var_value import calculate_variable_from_constraint + +from pyomo.dae import DerivativeVar + +# Import IDAES cores +from idaes.core.util.exceptions import ConfigurationError +from idaes.core.util.constants import Constants as const +import idaes.core.util.scaling as iscale +from idaes.core.util.misc import add_object_reference + +__author__ = "Jinliang Ma, Douglas Allan" + + +def make_geometry_common(blk, shell_units): + """Function to create variables, constraints, and expressions regarding + unit geometry shared between the CrossFlowHeatExchanger1D and Heater1D models. + + Args: + blk : unit model on which components are being generated + shell_units : derived units for property package of shell control volume + + Returns: + None + """ + # Number of tube columns in the cross section plane perpendicular to shell side fluid flow (y direction) + blk.number_columns_per_pass = Var( + initialize=10.0, doc="Number of tube columns", units=pyunits.dimensionless + ) + + # Number of segments of tube bundles + blk.number_passes = Var( + initialize=10.0, + doc="Number of segments of tube bundles", + units=pyunits.dimensionless, + ) + + # Number of inlet tube rows + blk.number_rows_per_pass = Var( + initialize=1, doc="Number of inlet tube rows", units=pyunits.dimensionless + ) + + # Inner diameter of tubes + blk.di_tube = Var( + initialize=0.05, doc="Inner diameter of tube", units=shell_units["length"] + ) + + blk.thickness_tube = Var( + initialize=0.005, doc="Tube thickness", units=shell_units["length"] + ) + + blk.pitch_y = Var( + initialize=0.1, + doc="Pitch between tubes perpendicular to direction of flow", + units=shell_units["length"], + ) + + blk.pitch_x = Var( + initialize=0.1, + doc="Pitch between tubes in line with direction of flow", + units=shell_units["length"], + ) + + # Length of tube per segment in z direction + blk.length_tube_seg = Var( + initialize=1.0, doc="Length of tube per segment", units=shell_units["length"] + ) + + # Minimum cross-sectional area on shell side + blk.area_flow_shell_min = Var( + initialize=1.0, doc="Minimum flow area on shell side", units=shell_units["area"] + ) + + # total number of tube rows + @blk.Expression(doc="Total number of tube rows") + def nrow_tube(b): + return b.number_passes * b.number_rows_per_pass + + # Tube outside diameter + @blk.Expression(doc="Outside diameter of tube") + def do_tube(b): + return b.di_tube + b.thickness_tube * 2.0 + + # Ratio of pitch_x/do_tube + @blk.Expression(doc="Ratio of pitch in x direction to tube outside diameter") + def pitch_x_to_do(b): + return b.pitch_x / b.do_tube + + # Ratio of pitch_y/do_tube + @blk.Expression(doc="Ratio of pitch in y direction to tube outside diameter") + def pitch_y_to_do(b): + return b.pitch_y / b.do_tube + + # Total cross-sectional area of tube metal per segment + @blk.Expression(doc="Total cross section area of tube metal per segment") + def area_wall_seg(b): + return ( + 0.25 + * const.pi + * (b.do_tube**2 - b.di_tube**2) + * b.number_columns_per_pass + * b.number_rows_per_pass + ) + + @blk.Expression(doc="Total heat transfer area on outer surface of tubes") + def total_heat_transfer_area(b): + return ( + const.pi + * b.do_tube + * b.number_rows_per_pass + * b.number_columns_per_pass + * b.number_passes + * b.length_tube_seg + ) + + # Length of shell side flow + @blk.Constraint(doc="Length of shell side flow") + def length_flow_shell_eqn(b): + return b.length_flow_shell == b.nrow_tube * b.pitch_x + + # Average flow area on shell side + @blk.Constraint(doc="Average cross section area of shell side flow") + def area_flow_shell_eqn(b): + return ( + b.length_flow_shell * b.area_flow_shell + == b.length_tube_seg + * b.length_flow_shell + * b.pitch_y + * b.number_columns_per_pass + - b.number_columns_per_pass + * b.nrow_tube + * 0.25 + * const.pi + * b.do_tube**2 + * b.length_tube_seg + ) + + # Minimum flow area on shell side + @blk.Constraint(doc="Minimum flow area on shell side") + def area_flow_shell_min_eqn(b): + return ( + b.area_flow_shell_min + == b.length_tube_seg * (b.pitch_y - b.do_tube) * b.number_columns_per_pass + ) + + +def make_performance_common( + blk, + shell, + shell_units, + shell_has_pressure_change: bool, + make_reynolds: bool, + make_nusselt: bool, +): + """Function to create variables, constraints, and expressions regarding + performance constraints shared between the CrossFlowHeatExchanger1D and + Heater1D models. + + Args: + blk : unit model on which components are being generated + shell: shell control volume + shell_units : derived units for property package of shell control volume + shell_has_pressure_change: bool about whether to make pressure change components + make_reynolds: bool about whether to create the Reynolds number + make_nusselt: bool about whether to create Nusselt number + + Returns: + None + """ + # We need the Reynolds number for pressure change, even if we don't need it for heat transfer + if shell_has_pressure_change: + make_reynolds = True + + add_object_reference(blk, "heat_shell", shell.heat) + + if shell_has_pressure_change: + add_object_reference(blk, "deltaP_shell", shell.deltaP) + + # Parameters + # Wall thermal conductivity + blk.therm_cond_wall = Param( + initialize=1.0, + mutable=True, + units=shell_units["thermal_conductivity"], + doc="Thermal conductivity of tube wall", + ) + + # Wall heat capacity + blk.cp_wall = Param( + initialize=502.4, + mutable=True, + units=shell_units["heat_capacity_mass"], + doc="Tube wall heat capacity", + ) + + # Wall density + blk.density_wall = Param( + initialize=7800.0, + mutable=True, + units=shell_units["density_mass"], + doc="Tube wall density", + ) + + # Heat transfer resistance due to the fouling on shell side + blk.rfouling_shell = Param( + units=1 / shell_units["heat_transfer_coefficient"], + initialize=0.0001, + mutable=True, + doc="Fouling resistance on tube side", + ) + + # Correction factor for convective heat transfer coefficient on shell side + blk.fcorrection_htc_shell = Var( + initialize=1.0, doc="Correction factor for convective HTC on shell" + ) + + # Correction factor for shell side pressure drop due to friction + if shell_has_pressure_change: + blk.fcorrection_dp_shell = Var( + initialize=1.0, doc="Correction factor for shell side pressure drop" + ) + + # Performance variables + # Shell side convective heat transfer coefficient due to convection only + blk.conv_heat_transfer_coeff_shell = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=100.0, + bounds=(0, None), + units=shell_units["heat_transfer_coefficient"], + doc="Shell side convective heat transfer coefficient due to convection", + ) + + # Boundary wall temperature on shell side + blk.temp_wall_shell = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=500, + units=shell_units["temperature"], + doc="Boundary wall temperature on shell side", + ) + + # Central wall temperature of tube metal, used to calculate energy contained by tube metal + blk.temp_wall_center = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=500, + units=shell_units["temperature"], + doc="Tube wall temperature at center", + ) + + # Tube wall heat holdup per length of shell + if blk.config.has_holdup: + blk.heat_holdup = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=1e6, + units=shell_units["energy"] / shell_units["length"], + doc="Tube wall heat holdup per length of shell", + ) + + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="Heat holdup of tube metal", + ) + def heat_holdup_eqn(b, t, x): + return ( + b.heat_holdup[t, x] + == b.cp_wall + * b.density_wall + * b.area_wall_seg + * pyunits.convert(b.length_flow_tube, to_units=shell_units["length"]) + / b.length_flow_shell + * b.temp_wall_center[t, x] + ) + + # Tube wall heat accumulation term + if blk.config.dynamic: + blk.heat_accumulation = DerivativeVar( + blk.heat_holdup, + initialize=0, + wrt=blk.flowsheet().config.time, + units=shell_units["energy"] / shell_units["length"] / shell_units["time"], + doc="Tube wall heat accumulation per unit length of shell", + ) + + # Pressure drop and heat transfer coefficient on shell side + # ---------------------------------------------------------- + # Tube arrangement factor + if blk.config.tube_arrangement == "in-line": + blk.f_arrangement = Param( + initialize=0.788, doc="in-line tube arrangement factor" + ) + elif blk.config.tube_arrangement == "staggered": + blk.f_arrangement = Param( + initialize=1.0, doc="staggered tube arrangement factor" + ) + else: + raise ConfigurationError() + + if make_reynolds: + # Velocity on shell side + blk.v_shell = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=1.0, + units=shell_units["velocity"], + doc="velocity on shell side", + ) + + # Reynalds number on shell side + blk.N_Re_shell = Var( + blk.flowsheet().config.time, + shell.length_domain, + bounds=(1e-7, None), + initialize=10000.0, + units=pyunits.dimensionless, + doc="Reynolds number on shell side", + ) + + if shell_has_pressure_change: + # Friction factor on shell side + blk.friction_factor_shell = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=1.0, + doc="friction factor on shell side", + ) + + if make_nusselt: + # Nusselt number on shell side + blk.N_Nu_shell = Var( + blk.flowsheet().config.time, + shell.length_domain, + initialize=1, + units=pyunits.dimensionless, + doc="Nusselts number on shell side", + bounds=(1e-7, None), + ) + + if make_reynolds: + # Velocity equation on shell side + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="velocity on shell side", + ) + def v_shell_eqn(b, t, x): + return ( + b.v_shell[t, x] + * shell.properties[t, x].dens_mol_phase["Vap"] + * b.area_flow_shell_min + == shell.properties[t, x].flow_mol + ) + + # Reynolds number + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="Reynolds number equation on shell side", + ) + def N_Re_shell_eqn(b, t, x): + return ( + b.N_Re_shell[t, x] * shell.properties[t, x].visc_d_phase["Vap"] + == b.do_tube + * b.v_shell[t, x] + * shell.properties[t, x].dens_mass_phase["Vap"] + ) + + if shell_has_pressure_change: + # Friction factor on shell side + if blk.config.tube_arrangement == "in-line": + + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="in-line friction factor on shell side", + ) + def friction_factor_shell_eqn(b, t, x): + return ( + b.friction_factor_shell[t, x] * b.N_Re_shell[t, x] ** 0.15 + == ( + 0.044 + + 0.08 + * b.pitch_x_to_do + / (b.pitch_y_to_do - 1.0) ** (0.43 + 1.13 / b.pitch_x_to_do) + ) + * b.fcorrection_dp_shell + ) + + elif blk.config.tube_arrangement == "staggered": + + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="staggered friction factor on shell side", + ) + def friction_factor_shell_eqn(b, t, x): + return ( + b.friction_factor_shell[t, x] * b.N_Re_shell[t, x] ** 0.16 + == (0.25 + 0.118 / (b.pitch_y_to_do - 1.0) ** 1.08) + * b.fcorrection_dp_shell + ) + + else: + raise ConfigurationError() + + # Pressure drop on shell side + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="pressure change on shell side", + ) + def deltaP_shell_eqn(b, t, x): + return ( + b.deltaP_shell[t, x] * b.pitch_x + == -1.4 + * b.friction_factor_shell[t, x] + * shell.properties[t, x].dens_mass_phase["Vap"] + * b.v_shell[t, x] ** 2 + ) + + if make_nusselt: + # The actual Nusselt number correlation needs to be made by the particular heat exchanger + @blk.Constraint( + blk.flowsheet().config.time, + shell.length_domain, + doc="Convective heat transfer coefficient equation on shell side due to convection", + ) + def conv_heat_transfer_coeff_shell_eqn(b, t, x): + return ( + b.conv_heat_transfer_coeff_shell[t, x] * b.do_tube + == b.N_Nu_shell[t, x] + * shell.properties[t, x].therm_cond_phase["Vap"] + * b.fcorrection_htc_shell + ) + + # Total convective heat transfer coefficient on shell side + @blk.Expression( + blk.flowsheet().config.time, + shell.length_domain, + doc="Total convective heat transfer coefficient on shell side", + ) + def total_heat_transfer_coeff_shell(b, t, x): + # Retain in case we add back radiation + return b.conv_heat_transfer_coeff_shell[t, x] + + +def scale_common(blk, shell, shell_has_pressure_change, make_reynolds, make_nusselt): + """Function to scale variables and constraints shared between + the CrossFlowHeatExchanger1D and Heater1D models. + + Args: + blk : unit model on which components are being generated + shell: shell control volume + shell_units : derived units for property package of shell control volume + shell_has_pressure_change: bool about whether to make pressure change components + make_reynolds: bool about whether to create the Reynolds number + make_nusselt: bool about whether to create Nusselt number + + Returns: + None + """ + + def gsf(obj): + return iscale.get_scaling_factor(obj, default=1, warning=True) + + def ssf(obj, sf): + iscale.set_scaling_factor(obj, sf, overwrite=False) + + def cst(con, sf): + iscale.constraint_scaling_transform(con, sf, overwrite=False) + + sgsf = iscale.set_and_get_scaling_factor + + sf_do_tube = iscale.get_scaling_factor(blk.do_tube, default=1 / value(blk.do_tube)) + calculate_variable_from_constraint( + blk.area_flow_shell_min, blk.area_flow_shell_min_eqn + ) + for t in blk.flowsheet().time: + for z in shell.length_domain: + sf_flow_mol_shell = gsf(shell.properties[t, z].flow_mol) + + if make_reynolds: + # FIXME get better scaling later + ssf(blk.v_shell[t, z], 1 / 10) + cst(blk.v_shell_eqn[t, z], sf_flow_mol_shell) + + # FIXME should get scaling of N_Re from defining eqn + sf_N_Re_shell = sgsf(blk.N_Re_shell[t, z], 1e-4) + + sf_visc_d_shell = gsf(shell.properties[t, z].visc_d_phase["Vap"]) + cst(blk.N_Re_shell_eqn[t, z], sf_N_Re_shell * sf_visc_d_shell) + if make_nusselt: + sf_k_shell = gsf(shell.properties[t, z].therm_cond_phase["Vap"]) + + sf_N_Nu_shell = sgsf( + blk.N_Nu_shell[t, z], 1 / 0.33 * sf_N_Re_shell**0.6 + ) + cst(blk.N_Nu_shell_eqn[t, z], sf_N_Nu_shell) + + sf_conv_heat_transfer_coeff_shell = sgsf( + blk.conv_heat_transfer_coeff_shell[t, z], + sf_N_Nu_shell * sf_k_shell / sf_do_tube, + ) + cst( + blk.conv_heat_transfer_coeff_shell_eqn[t, z], + sf_conv_heat_transfer_coeff_shell * sf_do_tube, + ) + + # FIXME estimate from parameters + if blk.config.has_holdup: + s_U_holdup = gsf(blk.heat_holdup[t, z]) + cst(blk.heat_holdup_eqn[t, z], s_U_holdup) diff --git a/idaes/models_extra/power_generation/unit_models/heater_1D.py b/idaes/models_extra/power_generation/unit_models/heater_1D.py new file mode 100644 index 0000000000..e52d5d3d7c --- /dev/null +++ b/idaes/models_extra/power_generation/unit_models/heater_1D.py @@ -0,0 +1,554 @@ +############################################################################## +# Institute for the Design of Advanced Energy Systems Process Systems +# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018-2019, 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.txt and LICENSE.txt for full copyright and +# license information, respectively. Both files are also available online +# at the URL "https://github.com/IDAES/idaes-pse". +############################################################################## +""" +1-D Electric Trim Heater Model With Wall Temperatures + +Discretization based on tube rows +""" +# Import Pyomo libraries +from pyomo.environ import ( + assert_optimal_termination, + Block, + units as pyunits, + Var, + value, +) +from pyomo.common.config import ConfigValue, In +from pyomo.util.calc_var_value import calculate_variable_from_constraint + +# Import IDAES cores +from idaes.core import ( + ControlVolume1DBlock, + declare_process_block_class, + MaterialBalanceType, + EnergyBalanceType, + MomentumBalanceType, + FlowDirection, + UnitModelBlockData, + useDefault, +) +import idaes.core.util.scaling as iscale +from idaes.core.solvers import get_solver +from idaes.core.util.config import is_physical_parameter_block +from idaes.core.util.exceptions import InitializationError +from idaes.core.util.misc import add_object_reference +import idaes.logger as idaeslog +from idaes.core.util.tables import create_stream_table_dataframe +from idaes.models_extra.power_generation.unit_models.heat_exchanger_common import ( + make_geometry_common, + make_performance_common, + scale_common, +) +from idaes.core.initialization import SingleControlVolumeUnitInitializer + +__author__ = "Jinliang Ma, Douglas Allan" + + +class Heater1DInitializer(SingleControlVolumeUnitInitializer): + """ + Initializer for Heater 1D units. + + A simple initialization method that first initializes the control volume + without heat transfer, then adds heat transfer in and solves it again, + then finally solves the entire model. + """ + + def initialize_main_model( + self, + model: Block, + copy_inlet_state: bool = False, + ): + """ + Initialization routine for the main Heater 1D model + (as opposed to submodels like costing, which presently do not exist). + + + 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. + + """ + # 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_obj = get_solver(self.config.solver, self.config.solver_options) + + # --------------------------------------------------------------------- + # Initialize shell block + + self.initialize_control_volume(model.control_volume, copy_inlet_state) + + init_log.info_high("Initialization Step 1 Complete.") + + calc_var = calculate_variable_from_constraint + + calc_var(model.length_flow_shell, model.length_flow_shell_eqn) + calc_var(model.area_flow_shell, model.area_flow_shell_eqn) + calc_var(model.area_flow_shell_min, model.area_flow_shell_min_eqn) + + for t in model.flowsheet().config.time: + for x in model.control_volume.length_domain: + model.control_volume.heat[t, x].fix( + value(model.electric_heat_duty[t] / model.length_flow_shell) + ) + + if model.config.has_pressure_change: + model.control_volume.pressure.fix() + + model.control_volume.length.fix() + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(model.control_volume, tee=slc.tee) + + try: + assert_optimal_termination(res) + except AssertionError: + raise InitializationError("Initialization solve failed.") + + init_log.info_high("Initialization Step 2 Complete.") + model.control_volume.length.unfix() + model.control_volume.heat.unfix() + + for t in model.flowsheet().config.time: + for x in model.control_volume.length_domain: + model.temp_wall_center[t, x].fix( + value(model.control_volume.properties[t, x].temperature) + 10 + ) + calc_var(model.heat_holdup[t, x], model.heat_holdup_eqn[t, x]) + model.temp_wall_center[t, x].unfix() + + if model.config.has_pressure_change: + model.control_volume.pressure.unfix() + model.control_volume.pressure[:, 0].fix() + + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = solver_obj.solve(model, tee=slc.tee) + + try: + assert_optimal_termination(res) + except AssertionError: + raise InitializationError("Initialization solve failed.") + + init_log.info_high("Initialization Step 3 Complete.") + + return res + + +@declare_process_block_class("Heater1D") +class Heater1DData(UnitModelBlockData): + """Standard Trim Heater Model Class Class.""" + + default_initializer = Heater1DInitializer + + CONFIG = UnitModelBlockData.CONFIG() + CONFIG.declare( + "has_fluid_holdup", + ConfigValue( + default=False, + domain=In([False]), + description="Holdup construction flag", + doc="""Indicates whether holdup terms for the fluid should be constructed or not. + **default** - False. + **Valid values:** { + **False** - do not construct holdup terms}""", + ), + ) + CONFIG.declare( + "material_balance_type", + ConfigValue( + default=MaterialBalanceType.componentTotal, + domain=In(MaterialBalanceType), + description="Material balance construction flag", + doc="""Indicates what type of mass balance should be constructed, +**default** - MaterialBalanceType.componentTotal. +**Valid values:** { +**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.}""", + ), + ) + CONFIG.declare( + "energy_balance_type", + ConfigValue( + default=EnergyBalanceType.enthalpyTotal, + domain=In(EnergyBalanceType), + description="Energy balance construction flag", + doc="""Indicates what type of energy balance should be constructed, +**default** - EnergyBalanceType.enthalpyTotal. +**Valid values:** { +**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.}""", + ), + ) + CONFIG.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.}""", + ), + ) + CONFIG.declare( + "has_pressure_change", + ConfigValue( + default=False, + domain=In([True, False]), + 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.}""", + ), + ) + CONFIG.declare( + "property_package", + ConfigValue( + domain=is_physical_parameter_block, + description="Property package to use for control volume", + doc="""Property parameter object used to define property calculations +(default = 'use_parent_value') +- 'use_parent_value' - get package from parent (default = None) +- a ParameterBlock object""", + ), + ) + CONFIG.declare( + "property_package_args", + ConfigValue( + default=None, + description="Arguments for constructing shell property package", + doc="""A dict of arguments to be passed to the PropertyBlockData +and used when constructing these +(default = 'use_parent_value') +- 'use_parent_value' - get package from parent (default = None) +- a dict (see property package for documentation)""", + ), + ) + CONFIG.declare( + "transformation_method", + ConfigValue( + default=useDefault, + description="Discretization method to use for DAE transformation", + doc="""Discretization method to use for DAE transformation. See Pyomo +documentation for supported transformations.""", + ), + ) + CONFIG.declare( + "transformation_scheme", + ConfigValue( + default=useDefault, + description="Discretization scheme to use for DAE transformation", + doc="""Discretization scheme to use when transforming domain. See Pyomo +documentation for supported schemes.""", + ), + ) + + # Common config args for both sides + CONFIG.declare( + "finite_elements", + ConfigValue( + default=5, + domain=int, + description="Number of finite elements length domain", + doc="""Number of finite elements to use when discretizing length +domain (default=5). Should set to the number of tube rows""", + ), + ) + CONFIG.declare( + "collocation_points", + ConfigValue( + default=3, + domain=int, + description="Number of collocation points per finite element", + doc="""If using collocation, number of collocation points to use + per finite element when discretizing length domain (default=3)""", + ), + ) + CONFIG.declare( + "tube_arrangement", + ConfigValue( + default="in-line", + domain=In(["in-line", "staggered"]), + description="tube configuration", + doc="tube arrangement could be in-line or staggered", + ), + ) + + def build(self): + """ + Begin building model (pre-DAE transformation). + + Args: + None + + Returns: + None + """ + # Call UnitModel.build to setup dynamics + super().build() + + # Set flow directions for the control volume blocks and specify + # dicretization if not specified. + set_direction_shell = FlowDirection.forward + if self.config.transformation_method is useDefault: + self.config.transformation_method = "dae.finite_difference" + if self.config.transformation_scheme is useDefault: + self.config.transformation_scheme = "BACKWARD" + + if self.config.property_package_args is None: + self.config.property_package_args = {} + + # Control volume 1D for shell, set to steady-state for fluid + self.control_volume = ControlVolume1DBlock( + dynamic=self.config.dynamic and self.config.has_fluid_holdup, + has_holdup=self.config.has_fluid_holdup, + property_package=self.config.property_package, + property_package_args=self.config.property_package_args, + transformation_method=self.config.transformation_method, + transformation_scheme=self.config.transformation_scheme, + finite_elements=self.config.finite_elements, + collocation_points=self.config.collocation_points, + ) + + self.control_volume.add_geometry(flow_direction=set_direction_shell) + + self.control_volume.add_state_blocks( + information_flow=set_direction_shell, + has_phase_equilibrium=False, + ) + + # Populate shell + self.control_volume.add_material_balances( + balance_type=self.config.material_balance_type, + has_phase_equilibrium=False, + ) + + self.control_volume.add_energy_balances( + balance_type=self.config.energy_balance_type, has_heat_transfer=True + ) + + self.control_volume.add_momentum_balances( + balance_type=self.config.momentum_balance_type, + has_pressure_change=self.config.has_pressure_change, + ) + + self.control_volume.apply_transformation() + + # Populate tube + + # Add Ports for shell side + self.add_inlet_port(name="inlet", block=self.control_volume) + self.add_outlet_port(name="outlet", block=self.control_volume) + + self._make_geometry() + + self._make_performance() + + def _make_geometry(self): + """ + Constraints for unit model. + + Args: + None + + Returns: + None + """ + units = self.config.property_package.get_metadata().derived_units + # Add reference to control volume geometry + add_object_reference(self, "area_flow_shell", self.control_volume.area) + add_object_reference(self, "length_flow_shell", self.control_volume.length) + make_geometry_common(self, shell_units=units) + + @self.Expression( + doc="Common performance equations expect this expression to be here" + ) + def length_flow_tube(b): + return b.number_passes * b.length_tube_seg + + def _make_performance(self): + """ + Constraints for unit model. + + Args: + None + + Returns: + None + """ + self.electric_heat_duty = Var( + self.flowsheet().config.time, + initialize=1e6, + units=pyunits.W, + doc="Heat duty provided to heater through resistive heating", + ) + units = self.config.property_package.get_metadata().derived_units + make_performance_common( + self, + shell=self.control_volume, + shell_units=units, + shell_has_pressure_change=self.config.has_pressure_change, + make_reynolds=True, + make_nusselt=True, + ) + + def heat_accumulation_term(b, t, x): + return b.heat_accumulation[t, x] if b.config.dynamic else 0 + + # Nusselt number, currently assume Re>300 + @self.Constraint( + self.flowsheet().config.time, + self.control_volume.length_domain, + doc="Nusselts number equation", + ) + def N_Nu_shell_eqn(b, t, x): + return ( + b.N_Nu_shell[t, x] + == b.f_arrangement + * 0.33 + * b.N_Re_shell[t, x] ** 0.6 + * b.control_volume.properties[t, x].prandtl_number_phase["Vap"] + ** 0.333333 + ) + + # Energy balance with tube wall + # ------------------------------------ + # Heat to wall per length + @self.Constraint( + self.flowsheet().config.time, + self.control_volume.length_domain, + doc="heat per length", + ) + def heat_shell_eqn(b, t, x): + return b.control_volume.heat[t, x] * b.length_flow_shell == ( + b.total_heat_transfer_coeff_shell[t, x] + * b.total_heat_transfer_area + * ( + b.temp_wall_shell[t, x] + - b.control_volume.properties[t, x].temperature + ) + ) + + # Shell side wall temperature + @self.Constraint( + self.flowsheet().config.time, + self.control_volume.length_domain, + doc="shell side wall temperature", + ) + def temp_wall_shell_eqn(b, t, x): + return ( + b.total_heat_transfer_coeff_shell[t, x] + * ( + b.control_volume.properties[t, x].temperature + - b.temp_wall_shell[t, x] + ) + # Divide thickness by 2 in order to represent center of hollow tube instead of + # interior edge of hollow tube + * (b.thickness_tube / (2 * b.therm_cond_wall) + b.rfouling_shell) + == b.temp_wall_shell[t, x] - b.temp_wall_center[t, x] + ) + + @self.Constraint( + self.flowsheet().config.time, + self.control_volume.length_domain, + doc="wall temperature", + ) + def temp_wall_center_eqn(b, t, x): + return heat_accumulation_term(b, t, x) == ( + -b.control_volume.heat[t, x] + + b.electric_heat_duty[t] / b.length_flow_shell + ) + + def calculate_scaling_factors(self): + def gsf(obj): + return iscale.get_scaling_factor(obj, default=1, warning=True) + + def ssf(obj, sf): + iscale.set_scaling_factor(obj, sf, overwrite=False) + + def cst(con, sf): + iscale.constraint_scaling_transform(con, sf, overwrite=False) + + scale_common( + self, + self.control_volume, + self.config.has_pressure_change, + make_reynolds=True, + make_nusselt=True, + ) + + sf_d_tube = iscale.get_scaling_factor( + self.do_tube, default=1 / value(self.do_tube) + ) + + for t in self.flowsheet().time: + for z in self.control_volume.length_domain: + sf_hconv_conv = gsf(self.conv_heat_transfer_coeff_shell[t, z]) + cst( + self.conv_heat_transfer_coeff_shell_eqn[t, z], + sf_hconv_conv * sf_d_tube, + ) + + sf_T = gsf(self.control_volume.properties[t, z].temperature) + ssf(self.temp_wall_shell[t, z], sf_T) + ssf(self.temp_wall_center[t, z], sf_T) + + s_Q = gsf(self.control_volume.heat[t, z]) + ssf(self.electric_heat_duty[t], s_Q / value(self.length_flow_shell)) + cst(self.heat_shell_eqn[t, z], s_Q * value(self.length_flow_shell)) + ssf(self.temp_wall_center[t, z], sf_T) + cst(self.temp_wall_shell_eqn[t, z], sf_T) + cst(self.temp_wall_center_eqn[t, z], s_Q) + + def _get_performance_contents(self, time_point=0): + var_dict = {} + var_dict["Electric Heat Duty"] = self.electric_heat_duty[time_point] + + expr_dict = {} + expr_dict["HX Area"] = self.total_heat_transfer_area + + return {"vars": var_dict, "exprs": expr_dict} + + def _get_stream_table_contents(self, time_point=0): + return create_stream_table_dataframe( + { + "Inlet": self.inlet, + "Outlet": self.outlet, + }, + time_point=time_point, + ) diff --git a/idaes/models_extra/power_generation/unit_models/tests/test_cross_flow_heat_exchanger_1D.py b/idaes/models_extra/power_generation/unit_models/tests/test_cross_flow_heat_exchanger_1D.py new file mode 100644 index 0000000000..9141f350e6 --- /dev/null +++ b/idaes/models_extra/power_generation/unit_models/tests/test_cross_flow_heat_exchanger_1D.py @@ -0,0 +1,309 @@ +################################################################################# +# 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. +################################################################################# + +import pytest +import pyomo.environ as pyo +from pyomo.util.check_units import assert_units_consistent + +from idaes.core import FlowsheetBlock +import idaes.core.util.scaling as iscale +from idaes.models.unit_models import HeatExchangerFlowPattern +from idaes.models.properties.modular_properties import GenericParameterBlock +from idaes.models_extra.power_generation.properties.natural_gas_PR import ( + get_prop, + EosType, +) +from idaes.models_extra.power_generation.unit_models import ( + CrossFlowHeatExchanger1D, + CrossFlowHeatExchanger1DInitializer, +) +import idaes.core.util.model_statistics as mstat +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.model_diagnostics import DiagnosticsToolbox + +# Set up solver +optarg = { + "constr_viol_tol": 1e-8, + "nlp_scaling_method": "user-scaling", + "linear_solver": "ma57", + "OF_ma57_automatic_scaling": "yes", + "max_iter": 350, + "tol": 1e-8, + "halt_on_ampl_error": "no", +} + + +def _create_model(pressure_drop): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.h2_side_prop_params = GenericParameterBlock( + **get_prop(["H2", "H2O", "Ar", "N2"], {"Vap"}, eos=EosType.IDEAL), + doc="H2O + H2 gas property parameters", + ) + m.fs.heat_exchanger = CrossFlowHeatExchanger1D( + has_holdup=True, + dynamic=False, + cold_side={ + "property_package": m.fs.h2_side_prop_params, + "has_holdup": False, + "dynamic": False, + "has_pressure_change": pressure_drop, + "transformation_method": "dae.finite_difference", + "transformation_scheme": "FORWARD", + }, + hot_side={ + "property_package": m.fs.h2_side_prop_params, + "has_holdup": False, + "dynamic": False, + "has_pressure_change": pressure_drop, + "transformation_method": "dae.finite_difference", + "transformation_scheme": "BACKWARD", + }, + shell_is_hot=True, + flow_type=HeatExchangerFlowPattern.countercurrent, + finite_elements=12, + tube_arrangement="staggered", + ) + + hx = m.fs.heat_exchanger + + hx.hot_side_inlet.flow_mol.fix(2619.7) + hx.hot_side_inlet.temperature.fix(971.6) + hx.hot_side_inlet.pressure.fix(1.2e5) + hx.hot_side_inlet.mole_frac_comp[0, "H2"].fix(0.79715) + hx.hot_side_inlet.mole_frac_comp[0, "H2O"].fix(0.20177) + hx.hot_side_inlet.mole_frac_comp[0, "Ar"].fix(0.00086358) + hx.hot_side_inlet.mole_frac_comp[0, "N2"].fix(0.00021589) + + hx.cold_side_inlet.flow_mol.fix(2619.7) + hx.cold_side_inlet.temperature.fix(446.21) + hx.cold_side_inlet.pressure.fix(1.2e5) + hx.cold_side_inlet.mole_frac_comp[0, "H2"].fix(0.36203) + hx.cold_side_inlet.mole_frac_comp[0, "H2O"].fix(0.63689) + hx.cold_side_inlet.mole_frac_comp[0, "Ar"].fix(0.00086358) + hx.cold_side_inlet.mole_frac_comp[0, "N2"].fix(0.00021589) + + hx.di_tube.fix(0.0525018) + hx.thickness_tube.fix(0.0039116) + hx.length_tube_seg.fix(4.3) + hx.number_passes.fix(12) + hx.number_columns_per_pass.fix(50) + hx.number_rows_per_pass.fix(25) + + hx.pitch_x.fix(0.1) + hx.pitch_y.fix(0.1) + hx.therm_cond_wall = 43.0 + hx.rfouling_tube = 0.0001 + hx.rfouling_shell = 0.0001 + hx.fcorrection_htc_tube.fix(1) + hx.fcorrection_htc_shell.fix(1) + if pressure_drop: + hx.fcorrection_dp_tube.fix(1) + hx.fcorrection_dp_shell.fix(1) + + hx.cp_wall.value = 502.4 + + pp = m.fs.h2_side_prop_params + pp.set_default_scaling("enth_mol_phase", 1e-3) + pp.set_default_scaling("pressure", 1e-5) + pp.set_default_scaling("temperature", 1e-2) + pp.set_default_scaling("flow_mol", 1e-3) + + _mf_scale = { + "H2": 1, + "H2O": 1, + "N2": 10, + "Ar": 10, + } + for comp, s in _mf_scale.items(): + pp.set_default_scaling("mole_frac_comp", s, index=comp) + pp.set_default_scaling("mole_frac_phase_comp", s, index=("Vap", comp)) + pp.set_default_scaling("flow_mol_phase_comp", s * 1e-3, index=("Vap", comp)) + + shell = hx.hot_side + tube = hx.cold_side + iscale.set_scaling_factor(shell.area, 1e-1) + iscale.set_scaling_factor(shell.heat, 1e-6) + iscale.set_scaling_factor(tube.area, 1) + iscale.set_scaling_factor(tube.heat, 1e-6) + iscale.set_scaling_factor(shell.enthalpy_flow_dx, 1e-7) + iscale.set_scaling_factor(tube.enthalpy_flow_dx, 1e-7) + iscale.set_scaling_factor(hx.heat_holdup, 1e-8) + + iscale.calculate_scaling_factors(m) + + return m + + +def _check_model_statistics(m, deltaP): + fixed_unused_var_set = { + "fs.h2_side_prop_params.H2.omega", + "fs.h2_side_prop_params.H2.pressure_crit", + "fs.h2_side_prop_params.H2.temperature_crit", + "fs.h2_side_prop_params.H2.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.H2.cp_mol_ig_comp_coeff_H", + "fs.h2_side_prop_params.H2O.omega", + "fs.h2_side_prop_params.H2O.pressure_crit", + "fs.h2_side_prop_params.H2O.temperature_crit", + "fs.h2_side_prop_params.H2O.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.H2O.cp_mol_ig_comp_coeff_H", + "fs.h2_side_prop_params.H2O.pressure_sat_comp_coeff_A", + "fs.h2_side_prop_params.H2O.pressure_sat_comp_coeff_B", + "fs.h2_side_prop_params.H2O.pressure_sat_comp_coeff_C", + "fs.h2_side_prop_params.Ar.omega", + "fs.h2_side_prop_params.Ar.pressure_crit", + "fs.h2_side_prop_params.Ar.temperature_crit", + "fs.h2_side_prop_params.Ar.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.Ar.cp_mol_ig_comp_coeff_H", + "fs.h2_side_prop_params.N2.omega", + "fs.h2_side_prop_params.N2.pressure_crit", + "fs.h2_side_prop_params.N2.temperature_crit", + "fs.h2_side_prop_params.N2.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.N2.cp_mol_ig_comp_coeff_H", + } + if not deltaP: + fixed_unused_var_set.add("fs.heat_exchanger.delta_elevation") + + for var in mstat.fixed_unused_variables_set(m): + assert var.name in fixed_unused_var_set + + unfixed_unused_var_set = { + "fs.heat_exchanger.hot_side.material_flow_dx[0.0,0.0,Vap,H2]", + "fs.heat_exchanger.hot_side.material_flow_dx[0.0,0.0,Vap,H2O]", + "fs.heat_exchanger.hot_side.material_flow_dx[0.0,0.0,Vap,Ar]", + "fs.heat_exchanger.hot_side.material_flow_dx[0.0,0.0,Vap,N2]", + "fs.heat_exchanger.hot_side.enthalpy_flow_dx[0.0,0.0,Vap]", + "fs.heat_exchanger.hot_side.pressure_dx[0.0,0.0]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,H2]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,H2O]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,Ar]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,N2]", + "fs.heat_exchanger.cold_side.enthalpy_flow_dx[0.0,1.0,Vap]", + "fs.heat_exchanger.cold_side.pressure_dx[0.0,1.0]", + } + + for var in mstat.unused_variables_set(m) - mstat.fixed_unused_variables_set(m): + assert var.name in unfixed_unused_var_set + + assert len(mstat.deactivated_constraints_set(m)) == 0 + + +@pytest.fixture +def model_no_dP(): + m = _create_model(pressure_drop=False) + return m + + +@pytest.mark.component +def test_initialization(model_no_dP): + m = model_no_dP + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=False) + + initializer = m.fs.heat_exchanger.default_initializer( + solver="ipopt", solver_options=optarg + ) + assert isinstance(initializer, CrossFlowHeatExchanger1DInitializer) + initializer.initialize(model=m.fs.heat_exchanger) + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=False) + assert pyo.value( + m.fs.heat_exchanger.hot_side_outlet.temperature[0] + ) == pytest.approx(485.34, abs=1e-1) + assert pyo.value( + m.fs.heat_exchanger.cold_side_outlet.temperature[0] + ) == pytest.approx(911.47, abs=1e-1) + + +@pytest.mark.integration +def test_structural_issues_no_dP(model_no_dP): + dt = DiagnosticsToolbox(model_no_dP) + dt.assert_no_structural_warnings(ignore_evaluation_errors=True) + + +@pytest.mark.integration +def test_numerical_issues_no_dP(model_no_dP): + # Model will already be initialized if the component test is run, + # but reinitialize in case integration tests are run alone + initializer = model_no_dP.fs.heat_exchanger.default_initializer( + solver="ipopt", solver_options=optarg + ) + initializer.initialize(model=model_no_dP.fs.heat_exchanger) + + m_scaled = pyo.TransformationFactory("core.scale_model").create_using( + model_no_dP, rename=False + ) + + dt = DiagnosticsToolbox(m_scaled) + dt.assert_no_numerical_warnings() + + +@pytest.fixture +def model_dP(): + m = _create_model(pressure_drop=True) + return m + + +@pytest.mark.component +def test_initialization_dP(model_dP): + m = model_dP + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=True) + + initializer = m.fs.heat_exchanger.default_initializer( + solver="ipopt", solver_options=optarg + ) + assert isinstance(initializer, CrossFlowHeatExchanger1DInitializer) + initializer.initialize(m.fs.heat_exchanger) + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=True) + + assert pyo.value( + m.fs.heat_exchanger.hot_side_outlet.temperature[0] + ) == pytest.approx(485.34, abs=1e-1) + assert pyo.value( + m.fs.heat_exchanger.cold_side_outlet.temperature[0] + ) == pytest.approx(911.47, abs=1e-1) + assert pyo.value(m.fs.heat_exchanger.hot_side_outlet.pressure[0]) == pytest.approx( + 118870.08569, abs=1 + ) + assert pyo.value(m.fs.heat_exchanger.cold_side_outlet.pressure[0]) == pytest.approx( + 111418.71399, abs=1 + ) + + +@pytest.mark.integration +def test_structural_issues_dP(model_dP): + dt = DiagnosticsToolbox(model_dP) + dt.assert_no_structural_warnings(ignore_evaluation_errors=True) + + +@pytest.mark.integration +def test_numerical_issues_dP(model_dP): + # Model will already be initialized if the component test is run, + # but reinitialize in case integration tests are run alone + initializer = model_dP.fs.heat_exchanger.default_initializer( + solver="ipopt", solver_options=optarg + ) + initializer.initialize(model=model_dP.fs.heat_exchanger) + + m_scaled = pyo.TransformationFactory("core.scale_model").create_using( + model_dP, rename=False + ) + + dt = DiagnosticsToolbox(m_scaled) + dt.assert_no_numerical_warnings() diff --git a/idaes/models_extra/power_generation/unit_models/tests/test_heater_1D.py b/idaes/models_extra/power_generation/unit_models/tests/test_heater_1D.py new file mode 100644 index 0000000000..91c34d2a26 --- /dev/null +++ b/idaes/models_extra/power_generation/unit_models/tests/test_heater_1D.py @@ -0,0 +1,262 @@ +################################################################################# +# 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. +################################################################################# + +import pytest +import pyomo.environ as pyo + +from idaes.core import FlowsheetBlock +import idaes.core.util.scaling as iscale +from idaes.models.properties.modular_properties import GenericParameterBlock +from idaes.models_extra.power_generation.properties.natural_gas_PR import ( + get_prop, + EosType, +) +from idaes.models_extra.power_generation.unit_models import ( + Heater1D, + Heater1DInitializer, +) +import idaes.core.util.model_statistics as mstat +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.model_diagnostics import DiagnosticsToolbox + +# Set up solver +optarg = { + "constr_viol_tol": 1e-8, + "nlp_scaling_method": "user-scaling", + "linear_solver": "ma57", + "OF_ma57_automatic_scaling": "yes", + "max_iter": 350, + "tol": 1e-8, + "halt_on_ampl_error": "no", +} + + +def _create_model(pressure_drop): + m = pyo.ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + m.fs.h2_side_prop_params = GenericParameterBlock( + **get_prop(["H2", "H2O", "Ar", "N2"], {"Vap"}, eos=EosType.IDEAL), + doc="H2O + H2 gas property parameters", + ) + m.fs.heater = Heater1D( + property_package=m.fs.h2_side_prop_params, + has_holdup=True, + dynamic=False, + has_fluid_holdup=False, + has_pressure_change=pressure_drop, + finite_elements=4, + tube_arrangement="in-line", + transformation_method="dae.finite_difference", + transformation_scheme="BACKWARD", + ) + + heater = m.fs.heater + + heater.inlet.flow_mol.fix(5102.5) + heater.inlet.temperature.fix(938.83) + heater.inlet.pressure.fix(1.2e5) + heater.inlet.mole_frac_comp[0, "H2"].fix(0.57375) + heater.inlet.mole_frac_comp[0, "H2O"].fix(0.42517) + heater.inlet.mole_frac_comp[0, "Ar"].fix(0.00086358) + heater.inlet.mole_frac_comp[0, "N2"].fix(0.00021589) + + heater.di_tube.fix(0.0525018) + heater.thickness_tube.fix(0.0039116) + heater.pitch_x.fix(0.1) + heater.pitch_y.fix(0.1) + heater.length_tube_seg.fix(10) + heater.number_passes.fix(1) + heater.rfouling = 0.0001 + heater.fcorrection_htc_shell.fix(1) + heater.cp_wall = 502.4 + if pressure_drop: + heater.fcorrection_dp_shell.fix(1) + + heater.number_columns_per_pass.fix(40) + heater.number_rows_per_pass.fix(40) + heater.electric_heat_duty.fix(3.6504e06) + + pp = m.fs.h2_side_prop_params + pp.set_default_scaling("enth_mol_phase", 1e-3) + pp.set_default_scaling("pressure", 1e-5) + pp.set_default_scaling("temperature", 1e-2) + pp.set_default_scaling("flow_mol", 1e-3) + + _mf_scale = { + "H2": 1, + "H2O": 1, + "N2": 10, + "Ar": 10, + } + for comp, s in _mf_scale.items(): + pp.set_default_scaling("mole_frac_comp", s, index=comp) + pp.set_default_scaling("mole_frac_phase_comp", s, index=("Vap", comp)) + pp.set_default_scaling("flow_mol_phase_comp", s * 1e-3, index=("Vap", comp)) + + shell = heater.control_volume + iscale.set_scaling_factor(shell.area, 1e-1) + iscale.set_scaling_factor(shell.heat, 1e-6) + iscale.set_scaling_factor(shell._enthalpy_flow, 1e-8) # pylint: disable=W0212 + iscale.set_scaling_factor(shell.enthalpy_flow_dx, 1e-7) + iscale.set_scaling_factor(heater.heat_holdup, 1e-8) + + iscale.calculate_scaling_factors(m) + + return m + + +def _check_model_statistics(m, deltaP): + fixed_unused_var_set = { + "fs.h2_side_prop_params.H2.omega", + "fs.h2_side_prop_params.H2.pressure_crit", + "fs.h2_side_prop_params.H2.temperature_crit", + "fs.h2_side_prop_params.H2.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.H2.cp_mol_ig_comp_coeff_H", + "fs.h2_side_prop_params.H2O.omega", + "fs.h2_side_prop_params.H2O.pressure_crit", + "fs.h2_side_prop_params.H2O.temperature_crit", + "fs.h2_side_prop_params.H2O.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.H2O.cp_mol_ig_comp_coeff_H", + "fs.h2_side_prop_params.H2O.pressure_sat_comp_coeff_A", + "fs.h2_side_prop_params.H2O.pressure_sat_comp_coeff_B", + "fs.h2_side_prop_params.H2O.pressure_sat_comp_coeff_C", + "fs.h2_side_prop_params.Ar.omega", + "fs.h2_side_prop_params.Ar.pressure_crit", + "fs.h2_side_prop_params.Ar.temperature_crit", + "fs.h2_side_prop_params.Ar.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.Ar.cp_mol_ig_comp_coeff_H", + "fs.h2_side_prop_params.N2.omega", + "fs.h2_side_prop_params.N2.pressure_crit", + "fs.h2_side_prop_params.N2.temperature_crit", + "fs.h2_side_prop_params.N2.cp_mol_ig_comp_coeff_G", + "fs.h2_side_prop_params.N2.cp_mol_ig_comp_coeff_H", + } + + for var in mstat.fixed_unused_variables_set(m): + assert var.name in fixed_unused_var_set + + unfixed_unused_var_set = { + "fs.heater.control_volume.material_flow_dx[0.0,0.0,Vap,H2]", + "fs.heater.control_volume.material_flow_dx[0.0,0.0,Vap,H2O]", + "fs.heater.control_volume.material_flow_dx[0.0,0.0,Vap,Ar]", + "fs.heater.control_volume.material_flow_dx[0.0,0.0,Vap,N2]", + "fs.heater.control_volume.enthalpy_flow_dx[0.0,0.0,Vap]", + "fs.heater.control_volume.pressure_dx[0.0,0.0]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,H2]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,H2O]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,Ar]", + "fs.heat_exchanger.cold_side.material_flow_dx[0.0,1.0,Vap,N2]", + "fs.heat_exchanger.cold_side.enthalpy_flow_dx[0.0,1.0,Vap]", + "fs.heat_exchanger.cold_side.pressure_dx[0.0,1.0]", + } + + for var in mstat.unused_variables_set(m) - mstat.fixed_unused_variables_set(m): + assert var.name in unfixed_unused_var_set + + assert len(mstat.deactivated_constraints_set(m)) == 0 + + +@pytest.fixture +def model_no_dP(): + m = _create_model(pressure_drop=False) + return m + + +@pytest.mark.component +def test_initialization(model_no_dP): + m = model_no_dP + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=False) + + initializer = m.fs.heater.default_initializer(solver="ipopt", solver_options=optarg) + assert isinstance(initializer, Heater1DInitializer) + initializer.initialize(model=m.fs.heater) + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=False) + assert pyo.value(m.fs.heater.outlet.temperature[0]) == pytest.approx( + 959.55, abs=1e-1 + ) + + +@pytest.mark.integration +def test_structural_issues_no_dP(model_no_dP): + dt = DiagnosticsToolbox(model_no_dP) + dt.assert_no_structural_warnings(ignore_evaluation_errors=True) + + +@pytest.mark.integration +def test_numerical_issues_no_dP(model_no_dP): + # Model will already be initialized if the component test is run, + # but reinitialize in case integration tests are run alone + initializer = model_no_dP.fs.heater.default_initializer( + solver="ipopt", solver_options=optarg + ) + initializer.initialize(model=model_no_dP.fs.heater) + + m_scaled = pyo.TransformationFactory("core.scale_model").create_using( + model_no_dP, rename=False + ) + + dt = DiagnosticsToolbox(m_scaled) + dt.assert_no_numerical_warnings() + + +@pytest.fixture +def model_dP(): + m = _create_model(pressure_drop=True) + return m + + +@pytest.mark.component +def test_initialization_dP(model_dP): + m = model_dP + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=True) + + initializer = m.fs.heater.default_initializer(solver="ipopt", solver_options=optarg) + assert isinstance(initializer, Heater1DInitializer) + initializer.initialize(model=m.fs.heater) + + assert degrees_of_freedom(m) == 0 + _check_model_statistics(m, deltaP=True) + + assert pyo.value(m.fs.heater.outlet.temperature[0]) == pytest.approx( + 959.55, abs=1e-1 + ) + assert pyo.value(m.fs.heater.outlet.pressure[0]) == pytest.approx(119762.3, abs=1) + + +@pytest.mark.integration +def test_structural_issues_dP(model_dP): + dt = DiagnosticsToolbox(model_dP) + dt.assert_no_structural_warnings(ignore_evaluation_errors=True) + + +@pytest.mark.integration +def test_numerical_issues_dP(model_dP): + # Model will already be initialized if the component test is run, + # but reinitialize in case integration tests are run alone + initializer = model_dP.fs.heater.default_initializer( + solver="ipopt", solver_options=optarg + ) + initializer.initialize(model=model_dP.fs.heater) + + m_scaled = pyo.TransformationFactory("core.scale_model").create_using( + model_dP, rename=False + ) + + dt = DiagnosticsToolbox(m_scaled) + dt.assert_no_numerical_warnings()