Source code for idaes.models_extra.column_models.plate_heat_exchanger

#################################################################################
# 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-2024 by the software owners: The Regents of the
# University of California, through Lawrence Berkeley National Laboratory,
# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
# University, West Virginia University Research Corporation, et al.
# All rights reserved.  Please see the files COPYRIGHT.md and LICENSE.md
# for full copyright and license information.
#################################################################################
"""
IDAES plate heat exchanger model (PHE) using effectiveness-NTU method, derived
from the core HeatExchangerNTU model.

Notes:
Let P = number of passes(same for both fluids)
if P is even PHE is in counter-current mode
if P is odd PHE  is in co-current mode
Hot and cold fluids flow alternatively in the channels.
Divider plates may be used to partition the passes in separate sections.
The heat transfer area  of a single plate depends on the total heat transfer
area as specified by the manufacturer and the total number of active plates.

Reference
Detailed model equations can be found in the paper :
[1] Akula, P., Eslick, J., Bhattacharyya, D. and Miller, D.C., 2019.
Modelling and Parameter Estimation of a Plate Heat Exchanger
as Part of a Solvent-Based Post-Combustion CO2 Capture System.
In Computer Aided Chemical Engineering (Vol. 47, pp. 47-52). Elsevier.

"""

# Import Pyomo libraries
from pyomo.environ import (
    Constraint,
    exp,
    Expression,
    NonNegativeIntegers,
    Param,
    PositiveIntegers,
    PositiveReals,
    units as pyunits,
    Var,
)
from pyomo.common.config import ConfigValue, In, Integer
from pyomo.util.calc_var_value import calculate_variable_from_constraint
from pyomo.common.deprecation import deprecated

# Import IDAES cores
from idaes.core import declare_process_block_class
from idaes.models.unit_models.heat_exchanger_ntu import HeatExchangerNTUData
from idaes.core.util.constants import Constants
from idaes.core.solvers import get_solver
import idaes.logger as idaeslog

__author__ = "Paul Akula, Andrew Lee"


# Set up logger
_log = idaeslog.getLogger(__name__)


[docs] @deprecated( "The Plate Heat Exchanger (PHE) model is known to be affected by " "issues causing it to fail to solve on certain platforms starting with Pyomo v6.7.0. " "This might cause the model to be removed in an upcoming IDAES release if these failures are not resolved. " "For more information, see IDAES/idaes-pse#1294", version="2.3.0", ) @declare_process_block_class("PlateHeatExchanger") class PlateHeatExchangerData(HeatExchangerNTUData): """Plate Heat Exchanger(PHE) Unit Model.""" CONFIG = HeatExchangerNTUData.CONFIG() CONFIG.declare( "passes", ConfigValue( default=4, domain=Integer, description="Number of passes", doc="""Number of passes of the fluids through the heat exchanger""", ), ) CONFIG.declare( "channels_per_pass", ConfigValue( default=12, domain=Integer, description="Number of channels for each pass", doc="""Number of channels to be used in each pass where a channel is the space between two plates with a flowing fluid""", ), ) CONFIG.declare( "number_of_divider_plates", ConfigValue( default=0, domain=Integer, description="Number of divider plates in heat exchanger", doc="""Divider plates are used to create separate partitions in the unit. Each pass can be separated by a divider plate""", ), ) # Update config block setting for pressure change to always be true CONFIG.hot_side.has_pressure_change = True CONFIG.hot_side.get("has_pressure_change").set_domain(In([True])) # pylint: disable-next=protected-access CONFIG.hot_side.get("has_pressure_change")._description = ( "Pressure change term construction flag - must be True" ) # pylint: disable-next=protected-access CONFIG.hot_side.get("has_pressure_change")._doc = ( "Plate Heat Exchanger model includes correlations for pressure drop " "thus has_pressure_change must be True" ) CONFIG.cold_side.has_pressure_change = True CONFIG.cold_side.get("has_pressure_change").set_domain(In([True])) # pylint: disable-next=protected-access CONFIG.cold_side.get("has_pressure_change")._description = ( "Pressure change term construction flag - must be True" ) # pylint: disable-next=protected-access CONFIG.cold_side.get("has_pressure_change")._doc = ( "Plate Heat Exchanger model includes correlations for pressure drop " "thus has_pressure_change must be True" )
[docs] def build(self): # Call super.build to setup model # This will create the control volumes, ports and basic equations super().build() # Units will be based on hot side properties units_meta = ( self.hot_side.config.property_package.get_metadata().get_derived_units ) # --------------------------------------------------------------------- # Plate design variables and parameter self.number_of_passes = Param( initialize=self.config.passes, units=pyunits.dimensionless, domain=PositiveIntegers, doc="Number of hot/cold fluid passes", mutable=True, ) # Assuming number of channels is equal in all plates self.channels_per_pass = Param( initialize=self.config.channels_per_pass, units=pyunits.dimensionless, domain=PositiveIntegers, doc="Number of channels in each pass", mutable=True, ) self.number_of_divider_plates = Param( initialize=self.config.number_of_divider_plates, units=pyunits.dimensionless, domain=NonNegativeIntegers, doc="Number of divider plates in heat exchanger", mutable=True, ) self.plate_length = Var( initialize=1.6925, units=units_meta("length"), domain=PositiveReals, doc="Length of heat exchanger plate", ) self.plate_width = Var( initialize=0.6135, units=units_meta("length"), domain=PositiveReals, doc="Width of heat exchanger plate", ) self.plate_thickness = Var( initialize=0.0006, units=units_meta("length"), domain=PositiveReals, doc="Thickness of heat exchanger plate", ) self.plate_pact_length = Var( initialize=0.381, units=units_meta("length"), domain=PositiveReals, doc="Compressed plate pact length", ) self.port_diameter = Var( initialize=0.2045, units=units_meta("length"), domain=PositiveReals, doc="Port diameter", ) self.plate_therm_cond = Var( initialize=16.2, units=units_meta("thermal_conductivity"), domain=PositiveReals, doc="Thermal conductivity heat exchanger plates", ) # Set default value of total heat transfer area self.area.set_value(114.3) # --------------------------------------------------------------------- # Derived geometric quantities total_plates = ( 2 * self.channels_per_pass * self.number_of_passes + 1 + self.number_of_divider_plates ) total_active_plates = 2 * self.channels_per_pass * self.number_of_passes - ( 1 + self.number_of_divider_plates ) self.plate_gap = Expression( expr=self.plate_pact_length / total_plates - self.plate_thickness ) self.plate_area = Expression( expr=self.area / total_active_plates, doc="Heat transfer area of single plate", ) self.surface_enlargement_factor = Expression( expr=self.plate_area / (self.plate_length * self.plate_width) ) # Channel equivalent diameter self.channel_diameter = Expression( expr=2 * self.plate_gap / self.surface_enlargement_factor, doc="Channel equivalent diameter", ) # --------------------------------------------------------------------- # Fluid velocities def rule_port_vel_hot(blk, t): return ( 4 * blk.hot_side.properties_in[t].flow_vol / (Constants.pi * blk.port_diameter**2) ) self.hot_port_velocity = Expression( self.flowsheet().time, rule=rule_port_vel_hot, doc="Hot side port velocity" ) def rule_port_vel_cold(blk, t): return ( 4 * pyunits.convert( blk.cold_side.properties_in[t].flow_vol, to_units=units_meta("flow_vol"), ) / (Constants.pi * blk.port_diameter**2) ) self.cold_port_velocity = Expression( self.flowsheet().time, rule=rule_port_vel_cold, doc="Cold side port velocity", ) def rule_channel_vel_hot(blk, t): return blk.hot_side.properties_in[t].flow_vol / ( blk.channels_per_pass * blk.plate_width * blk.plate_gap ) self.hot_channel_velocity = Expression( self.flowsheet().time, rule=rule_channel_vel_hot, doc="Hot side channel velocity", ) def rule_channel_vel_cold(blk, t): return pyunits.convert( blk.cold_side.properties_in[t].flow_vol, to_units=units_meta("flow_vol") ) / (blk.channels_per_pass * blk.plate_width * blk.plate_gap) self.cold_channel_velocity = Expression( self.flowsheet().time, rule=rule_channel_vel_cold, doc="Cold side channel velocity", ) # --------------------------------------------------------------------- # Reynolds & Prandtl numbers # Density cancels out of Reynolds number if mass flow rate is used def rule_Re_h(blk, t): return ( blk.hot_side.properties_in[t].flow_mass * blk.channel_diameter / ( blk.channels_per_pass * blk.plate_width * blk.plate_gap * blk.hot_side.properties_in[t].visc_d_phase["Liq"] ) ) self.Re_hot = Expression( self.flowsheet().time, rule=rule_Re_h, doc="Hot side Reynolds number" ) def rule_Re_c(blk, t): return ( pyunits.convert( blk.cold_side.properties_in[t].flow_mass / blk.cold_side.properties_in[t].visc_d_phase["Liq"], to_units=units_meta("length"), ) * blk.channel_diameter / (blk.channels_per_pass * blk.plate_width * blk.plate_gap) ) self.Re_cold = Expression( self.flowsheet().time, rule=rule_Re_c, doc="Cold side Reynolds number" ) def rule_Pr_h(blk, t): return ( blk.hot_side.properties_in[t].cp_mol / blk.hot_side.properties_in[t].mw * blk.hot_side.properties_in[t].visc_d_phase["Liq"] / blk.hot_side.properties_in[t].therm_cond_phase["Liq"] ) self.Pr_hot = Expression( self.flowsheet().time, rule=rule_Pr_h, doc="Hot side Prandtl number" ) def rule_Pr_c(blk, t): return ( blk.cold_side.properties_in[t].cp_mol / blk.cold_side.properties_in[t].mw * blk.cold_side.properties_in[t].visc_d_phase["Liq"] / blk.cold_side.properties_in[t].therm_cond_phase["Liq"] ) self.Pr_cold = Expression( self.flowsheet().time, rule=rule_Pr_c, doc="Cold side Prandtl number" ) # --------------------------------------------------------------------- # Heat transfer coefficients # Parameters for Nusselt number correlation self.Nusselt_param_a = Param( initialize=0.4, domain=PositiveReals, units=pyunits.dimensionless, mutable=True, doc="Nusselt parameter A", ) self.Nusselt_param_b = Param( initialize=0.663, domain=PositiveReals, units=pyunits.dimensionless, mutable=True, doc="Nusselt parameter B", ) self.Nusselt_param_c = Param( initialize=0.333, domain=PositiveReals, units=pyunits.dimensionless, mutable=True, doc="Nusselt parameter C", ) # Film heat transfer coefficients def rule_hotside_transfer_coeff(blk, t): return ( blk.hot_side.properties_in[t].therm_cond_phase["Liq"] / blk.channel_diameter * blk.Nusselt_param_a * blk.Re_hot[t] ** blk.Nusselt_param_b * blk.Pr_hot[t] ** blk.Nusselt_param_c ) self.heat_transfer_coefficient_hot_side = Expression( self.flowsheet().time, rule=rule_hotside_transfer_coeff, doc="Hot side heat transfer coefficient", ) def rule_coldside_transfer_coeff(blk, t): return ( pyunits.convert( blk.cold_side.properties_in[t].therm_cond_phase["Liq"], to_units=units_meta("thermal_conductivity"), ) / blk.channel_diameter * blk.Nusselt_param_a * blk.Re_cold[t] ** blk.Nusselt_param_b * blk.Pr_cold[t] ** blk.Nusselt_param_c ) self.heat_transfer_coefficient_cold_side = Expression( self.flowsheet().time, rule=rule_coldside_transfer_coeff, doc="Cold side heat transfer coefficient", ) # Overall heat transfer coefficient def rule_U(blk, t): return blk.heat_transfer_coefficient[t] == ( 1.0 / ( 1.0 / blk.heat_transfer_coefficient_hot_side[t] + blk.plate_gap / blk.plate_therm_cond + 1.0 / blk.heat_transfer_coefficient_cold_side[t] ) ) self.overall_heat_transfer_eq = Constraint( self.flowsheet().time, rule=rule_U, doc="Calculations of overall heat transfer coefficient", ) # Effectiveness based on sub-heat exchangers # Divide NTU by number of channels per pass def rule_Ecf(blk, t): if blk.number_of_passes.value % 2 == 0: return blk.effectiveness[t] == ( 1 - exp(-blk.NTU[t] / blk.channels_per_pass * (1 - blk.Cratio[t])) ) / ( 1 - blk.Cratio[t] * exp(-blk.NTU[t] / blk.channels_per_pass * (1 - blk.Cratio[t])) ) elif blk.pass_num.value % 2 == 1: return blk.effectiveness[t] == ( 1 - exp(-blk.NTU[t] / blk.channels_per_pass * (1 + blk.Cratio[t])) ) / (1 + blk.Cratio[t]) self.effectiveness_correlation = Constraint( self.flowsheet().time, rule=rule_Ecf, doc="Correlation for effectiveness factor", ) # --------------------------------------------------------------------- # Pressure drop correlations # Friction factor calculation self.friction_factor_param_a = Param( initialize=0.0, units=pyunits.dimensionless, doc="Friction factor parameter A", mutable=True, ) self.friction_factor_param_b = Param( initialize=18.29, units=pyunits.dimensionless, doc="Friction factor parameter B", mutable=True, ) self.friction_factor_param_c = Param( initialize=-0.652, units=pyunits.dimensionless, doc="Friction factor parameter C", mutable=True, ) def rule_fric_h(blk, t): return ( blk.friction_factor_param_a + blk.friction_factor_param_b * blk.Re_hot[t] ** (blk.friction_factor_param_c) ) self.friction_factor_hot = Expression( self.flowsheet().time, rule=rule_fric_h, doc="Hot side friction factor" ) def rule_fric_c(blk, t): return ( blk.friction_factor_param_a + blk.friction_factor_param_b * blk.Re_cold[t] ** (blk.friction_factor_param_c) ) self.friction_factor_cold = Expression( self.flowsheet().time, rule=rule_fric_c, doc="Cold side friction factor" ) def rule_hotside_dP(blk, t): return blk.hot_side.deltaP[t] == -( ( 2 * blk.friction_factor_hot[t] * (blk.plate_length + blk.port_diameter) * blk.number_of_passes * blk.hot_channel_velocity[t] ** 2 * blk.hot_side.properties_in[t].dens_mass / blk.channel_diameter ) + ( 0.7 * blk.number_of_passes * blk.hot_port_velocity[t] ** 2 * blk.hot_side.properties_in[t].dens_mass ) + ( blk.hot_side.properties_in[t].dens_mass * pyunits.convert( Constants.acceleration_gravity, to_units=units_meta("acceleration"), ) * (blk.plate_length + blk.port_diameter) ) ) self.hot_side_deltaP_eq = Constraint( self.flowsheet().time, rule=rule_hotside_dP ) def rule_coldside_dP(blk, t): return blk.cold_side.deltaP[t] == -( ( 2 * blk.friction_factor_cold[t] * (blk.plate_length + blk.port_diameter) * blk.number_of_passes * blk.cold_channel_velocity[t] ** 2 * pyunits.convert( blk.cold_side.properties_in[t].dens_mass, to_units=units_meta("density_mass"), ) / blk.channel_diameter ) + ( 0.7 * blk.number_of_passes * blk.cold_port_velocity[t] ** 2 * pyunits.convert( blk.cold_side.properties_in[t].dens_mass, to_units=units_meta("density_mass"), ) ) + ( pyunits.convert( blk.cold_side.properties_in[t].dens_mass, to_units=units_meta("density_mass"), ) * pyunits.convert( Constants.acceleration_gravity, to_units=units_meta("acceleration"), ) * (blk.plate_length + blk.port_diameter) ) ) self.cold_side_deltaP_eq = Constraint( self.flowsheet().time, rule=rule_coldside_dP )
[docs] def initialize( self, hot_side_state_args=None, cold_side_state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None, duty=None, ): """ Heat exchanger initialization method. Args: hot_side_state_args : a dict of arguments to be passed to the property initialization for the hot side (see documentation of the specific property package) (default = None). cold_side_state_args : a dict of arguments to be passed to the property initialization for the cold side (see documentation of the specific property package) (default = None). outlvl : sets output level of initialization routine optarg : solver options dictionary object (default=None, use default solver options) solver : str indicating which solver to use during initialization (default = None, use default solver) duty : an initial guess for the amount of heat transferred. This should be a tuple in the form (value, units), (default = (1000 J/s)) Returns: None """ # Set solver options init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") hot_side = self.hot_side cold_side = self.cold_side # Create solver opt = get_solver(solver, optarg) flags1 = hot_side.initialize( outlvl=outlvl, optarg=optarg, solver=solver, state_args=hot_side_state_args ) init_log.info_high("Initialization Step 1a (hot side) Complete.") flags2 = cold_side.initialize( outlvl=outlvl, optarg=optarg, solver=solver, state_args=cold_side_state_args ) init_log.info_high("Initialization Step 1b (cold side) Complete.") # --------------------------------------------------------------------- # Solve unit without heat transfer equation self.energy_balance_constraint.deactivate() self.effectiveness_correlation.deactivate() self.effectiveness.fix(0.68) # Get side 1 and side 2 heat units, and convert duty as needed s1_units = hot_side.heat.get_units() s2_units = cold_side.heat.get_units() if duty is None: # Assume 1000 J/s and check for unitless properties if s1_units is None and s2_units is None: # Backwards compatibility for unitless properties s1_duty = -1000 s2_duty = 1000 else: s1_duty = pyunits.convert_value( -1000, from_units=pyunits.W, to_units=s1_units ) s2_duty = pyunits.convert_value( 1000, from_units=pyunits.W, to_units=s2_units ) else: # Duty provided with explicit units s1_duty = -pyunits.convert_value( duty[0], from_units=duty[1], to_units=s1_units ) s2_duty = pyunits.convert_value( duty[0], from_units=duty[1], to_units=s2_units ) cold_side.heat.fix(s2_duty) for i in hot_side.heat: hot_side.heat[i].value = s1_duty with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(self, tee=slc.tee) init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res))) cold_side.heat.unfix() self.energy_balance_constraint.activate() for t in self.effectiveness: calculate_variable_from_constraint( self.effectiveness[t], self.effectiveness_correlation[t] ) # --------------------------------------------------------------------- # Solve unit with new effectiveness factor with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(self, tee=slc.tee) init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res))) self.effectiveness_correlation.activate() self.effectiveness.unfix() # --------------------------------------------------------------------- # Final solve of full modelr with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: res = opt.solve(self, tee=slc.tee) init_log.info_high("Initialization Step 4 {}.".format(idaeslog.condition(res))) # --------------------------------------------------------------------- # Release Inlet state hot_side.release_state(flags1, outlvl=outlvl) cold_side.release_state(flags2, outlvl=outlvl) init_log.info("Initialization Completed, {}".format(idaeslog.condition(res)))