Heat Exchanger NTU Unit Model with Aqueous MEA System

../../../_images/heat_exchanger_4.svg

Problem Statement: In this example, we will be transfering heat from hot to cold streams of MEA (monoethanolamine), simulating heat integration of an aqueous MEA system.

Hot Side Inlet

Flow Rate = 60.54879 mol/s

Mole fraction (CO2) = 0.0158

Mole fraction (H2O) = 0.8747

Mole fraction (MEA) = 0.1095

Pressure = 202650 Pa

Temperature = 392.23 K

Cold Side Inlet

Flow Rate = 63.01910 mol/s

Mole fraction (CO2) = 0.0414

Mole fraction (H2O) = 0.8509

Mole fraction (MEA) = 0.1077

Pressure = 202650 Pa

Temperature = 326.36 K

This example will demonstrate the simulation of the NTU heat exchanger by fixing the following degrees of freedom: - heat transfer area - heat transfer coefficient - effectiveness - hot and cold side pressure changes

IDAES documentation reference for heat exchanger NTU model: https://idaes-pse.readthedocs.io/en/latest/reference_guides/model_libraries/generic/unit_models/heat_exchanger_ntu.html

This example utilizes the NTU heat exchanger for a shell and tube system.

Importing necessary tools

First, import the required IDAES and Pyomo modules:

# Import pyomo package
from pyomo.environ import ConcreteModel, Constraint, value, units

# Import idaes logger to set output levels
import idaes.logger as idaeslog

# Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model
from idaes.core import FlowsheetBlock

# Import the MEA property package to create a property block for the flowsheet
from idaes.models_extra.column_models.properties.MEA_solvent \
    import configuration as aqueous_mea

from idaes.models.properties.modular_properties.base.generic_property import (
        GenericParameterBlock)

# Import the degrees_of_freedom function from the idaes.core.util.model_statistics package
from idaes.core.util.model_statistics import degrees_of_freedom

# Import the default IPOPT solver
from idaes.core.solvers import get_solver

# Import a heat exchanger unit
from idaes.models.unit_models.heat_exchanger_ntu import HeatExchangerNTU as HXNTU

Setting up the flowsheet

Then, we will create the ConcreteModel foundation, attach the steady state flowsheet, and declare the property parameter block that will used for the hot and cold sides.

More information on this general workflow can be found here: https://idaes-pse.readthedocs.io/en/stable/how_to_guides/workflow/general.html

m = ConcreteModel()

m.fs = FlowsheetBlock(dynamic=False)

m.fs.hotside_properties = GenericParameterBlock(**aqueous_mea)

m.fs.coldside_properties = GenericParameterBlock(**aqueous_mea)

Next, we specify hot and cold streams; note that heat flow from the hot to the cold stream is calculated using the transfer unit method for the HXNTU unit model.

The unit is created below:

# Create an instance of the heat exchanger unit, attaching it to the flowsheet
# Specify that the property package to be used with the heat exchanger is the one we created earlier.
m.fs.heat_exchanger = HXNTU(
    hot_side_name="shell",
    cold_side_name="tube",
    shell={"property_package": m.fs.hotside_properties,
           "has_pressure_change": True},
    tube={"property_package": m.fs.coldside_properties,
          "has_pressure_change": True})

# Call the degrees_of_freedom function, get initial DOF
DOF_initial = degrees_of_freedom(m)
print("The initial DOF is {0}".format(DOF_initial))
The initial DOF is 15
assert DOF_initial == 15

Fixing input specifications

For this problem, we will fix the inlet conditions, re-evaluate the degrees of freedom to ensure the problem is square (i.e. DOF=0), and run two different options for unit specifications:

m.fs.heat_exchanger.hot_side_inlet.flow_mol[0].fix(60.54879) # mol/s
m.fs.heat_exchanger.hot_side_inlet.temperature[0].fix(392.23) # K
m.fs.heat_exchanger.hot_side_inlet.pressure[0].fix(202650) # Pa
m.fs.heat_exchanger.hot_side_inlet.mole_frac_comp[0, "CO2"].fix(0.0158) # dimensionless
m.fs.heat_exchanger.hot_side_inlet.mole_frac_comp[0, "H2O"].fix(0.8747) # dimensionless
m.fs.heat_exchanger.hot_side_inlet.mole_frac_comp[0, "MEA"].fix(0.1095) # dimensionless

m.fs.heat_exchanger.cold_side_inlet.flow_mol[0].fix(63.01910) # mol/s
m.fs.heat_exchanger.cold_side_inlet.temperature[0].fix(326.36) # K
m.fs.heat_exchanger.cold_side_inlet.pressure[0].fix(202650) # Pa
m.fs.heat_exchanger.cold_side_inlet.mole_frac_comp[0, "CO2"].fix(0.0414) # dimensionless
m.fs.heat_exchanger.cold_side_inlet.mole_frac_comp[0, "H2O"].fix(0.8509) # dimensionless
m.fs.heat_exchanger.cold_side_inlet.mole_frac_comp[0, "MEA"].fix(0.1077) # dimensionless

DOF_initial = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_initial))
The DOF is 3

Option 1: Fix heat transfer coefficient (HTC) and dimensions of each domain

Below, we fix the heat exchanger area, effectiveness factor and heat transfer coefficient, which yields a fully defined problem with zero degrees of freedom that may be initialized and solved:

# Unit design variables
m.fs.heat_exchanger.area.fix(100) # m2
m.fs.heat_exchanger.heat_transfer_coefficient.fix(200) # W/m2/K
m.fs.heat_exchanger.effectiveness.fix(0.7) # dimensionless

m.fs.heat_exchanger.hot_side.deltaP.fix(-2000) # Pa
m.fs.heat_exchanger.cold_side.deltaP.fix(-2000) # Pa

DOF_final = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_final))
The DOF is 0
assert DOF_final == 0

Now that the problem is square (zero degress of freedom), we can initialize and solve the full model:

# Initialize the flowsheet, and set the output at INFO
m.fs.heat_exchanger.initialize(outlvl=idaeslog.INFO)

# Solve the simulation using ipopt
# Note: If the degrees of freedom = 0, we have a square problem
opt = get_solver()
solve_status = opt.solve(m, tee = True)

# Display a readable report
m.fs.heat_exchanger.report()
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side.properties_in: Starting initialization
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side.properties_in: State variable initialization completed.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side.properties_in: Property initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side.properties_out: Starting initialization
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side.properties_out: State variable initialization completed.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side.properties_out: Property initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.hot_side: Initialization Complete
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties_in: Starting initialization
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties_in: State variable initialization completed.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties_in: Property initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties_out: Starting initialization
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties_out: State variable initialization completed.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side.properties_out: Property initialization: optimal - Optimal Solution Found.
2023-03-04 01:48:06 [INFO] idaes.init.fs.heat_exchanger.cold_side: Initialization Complete
2023-03-04 01:48:07 [INFO] idaes.init.fs.heat_exchanger: Initialization Completed, optimal - Optimal Solution Found
Ipopt 3.13.2: nlp_scaling_method=gradient-based
tol=1e-06
max_iter=200


**************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
**************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:      575
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      446

Total number of variables............................:      122
                     variables with only lower bounds:       28
                variables with lower and upper bounds:       76
                     variables with only upper bounds:        0
Total number of equality constraints.................:      122
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 4.33e+02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 8.16e+00 2.63e+01  -1.0 1.01e-02    -  9.90e-01 9.81e-01h  1
   2  0.0000000e+00 7.34e-02 3.03e+01  -1.0 4.29e-04    -  9.90e-01 9.91e-01h  1
   3  0.0000000e+00 1.23e-09 4.64e+02  -1.0 3.96e-06    -  9.91e-01 1.00e+00h  1

Number of Iterations....: 3

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   2.8492324393151369e-12    1.2266809790162370e-09
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.8492324393151369e-12    1.2266809790162370e-09


Number of objective function evaluations             = 4
Number of objective gradient evaluations             = 4
Number of equality constraint evaluations            = 4
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 4
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 3
Total CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

====================================================================================
Unit : fs.heat_exchanger                                                   Time: 0.0
------------------------------------------------------------------------------------
    Stream Table
                                Units       Hot Inlet  Hot Outlet  Cold Inlet  Cold Outlet
    Total Molar Flowrate     mole / second     60.549      60.549      63.019      63.019
    Total Mole Fraction H2O  dimensionless    0.87470     0.87470     0.85090     0.85090
    Total Mole Fraction MEA  dimensionless    0.10950     0.10950     0.10770     0.10770
    Total Mole Fraction CO2  dimensionless   0.015800    0.015800    0.041400    0.041400
    Temperature                     kelvin     392.23      344.00      326.36      374.33
    Pressure                        pascal 2.0265e+05  2.0065e+05  2.0265e+05  2.0065e+05
====================================================================================
from pyomo.environ import assert_optimal_termination
import pytest

# Check if termination condition is optimal
assert_optimal_termination(solve_status)

assert value(m.fs.heat_exchanger.hot_side_outlet.temperature[0]) == pytest.approx(344.00, abs=1e-2)
assert value(m.fs.heat_exchanger.cold_side_outlet.temperature[0]) == pytest.approx(374.33, abs=1e-2)

Option 2: Unfix shell length and fix shell outlet temperatures

In the previous case, we fixed the heat exchanger area and overall heat transfer coefficient. However, given that the models in IDAES are equation oriented, we can fix the outlet variables. For example, we can fix the hot outlet temperature and heat exchanger length, and solve for the heat exchanger area and effectiveness factor that will satisfy that condition.

m.fs.heat_exchanger.effectiveness.unfix()
m.fs.heat_exchanger.area.unfix()
m.fs.heat_exchanger.hot_side_outlet.temperature.fix(344.0)

DOF_final = degrees_of_freedom(m)
print("The DOF is {0}".format(DOF_final))
The DOF is 0
result = opt.solve(m)

print(result)

# Display a readable report
m.fs.heat_exchanger.report()
Problem:
- Lower bound: -inf
  Upper bound: inf
  Number of objectives: 1
  Number of constraints: 122
  Number of variables: 122
  Sense: unknown
Solver:
- Status: ok
  Message: Ipopt 3.13.2x3a Optimal Solution Found
  Termination condition: optimal
  Id: 0
  Error rc: 0
  Time: 0.015180826187133789
Solution:
- number of solutions: 0
  number of solutions displayed: 0


====================================================================================
Unit : fs.heat_exchanger                                                   Time: 0.0
------------------------------------------------------------------------------------
    Stream Table
                                Units       Hot Inlet  Hot Outlet  Cold Inlet  Cold Outlet
    Total Molar Flowrate     mole / second     60.549      60.549      63.019      63.019
    Total Mole Fraction H2O  dimensionless    0.87470     0.87470     0.85090     0.85090
    Total Mole Fraction MEA  dimensionless    0.10950     0.10950     0.10770     0.10770
    Total Mole Fraction CO2  dimensionless   0.015800    0.015800    0.041400    0.041400
    Temperature                     kelvin     392.23      344.00      326.36      374.33
    Pressure                        pascal 2.0265e+05  2.0065e+05  2.0265e+05  2.0065e+05
====================================================================================
# Check if termination condition is optimal
assert_optimal_termination(result)

assert value(m.fs.heat_exchanger.area) == pytest.approx(100, abs=1e-2)
assert value(m.fs.heat_exchanger.effectiveness[0]) == pytest.approx(0.7, abs=1e-2)
assert value(m.fs.heat_exchanger.cold_side_outlet.temperature[0]) == pytest.approx(
    374.33, rel=1e-3)