Pyomo Model of the Double Pipe Heat Exchanger
Contents
1.6. Pyomo Model of the Double Pipe Heat Exchanger¶
This is additional material regarding the modeling and analysis of the double pipe heat exchanger. If you are using this notebook in Google Colab, run the following cell to import needed libraries to run the notebook code.
import sys
if "google.colab" in sys.modules:
!pip install -q pyomo
!wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"
!unzip -o -q ipopt-linux64
1.6.1. Temperature Dependence Heat Transfer Coefficient¶
1.6.1.1. Dittus-Boelter Equation¶
The Prandtl number for water is a reasonably strong function of temperature over the operating range of this heat exchanger.
where
1.6.1.2. Whitaker Correlation¶
From a meta-analysis of data for well developed turbulent flow in a pipe, Stephen Whitaker proposed a correlation
where
where
For water in the temperature range of interest, both
where
1.6.2. Colburn’s Analogy¶
Colburn’s analogy between momentum and heat transfer in turbulent flow yields a correlation
where properties are evaluated at the bulk flow temperature.
1.6.2.1. Physical property equations from¶
Pátek, J., Hrubý, J., Klomfar, J., Součková, M., & Harvey, A. H. (2009). Reference correlations for thermophysical properties of liquid water at 0.1 MPa. Journal of Physical and Chemical Reference Data, 38(1), 21-29. https://aip.scitation.org/doi/10.1063/1.3043575
import numpy as np
import math
import matplotlib.pyplot as plt
cp = 4183.0 # J/kg/K
def viscosity(T):
ac = [280.68, 511.45, 61.131, 0.45903]
bc = [-1.9, -7.7, -19.6, -40.0]
Tr = (T + 273.15)/300
return 1e-6*sum(a*Tr**b for a, b in zip(ac, bc))
def thermal_conductivity(T):
Tr = (T + 273.15)/300
cc = [0.80201, -0.25992, 0.10024, -0.032005]
dc = [-0.32, -5.7, -12.0, -15.0]
return sum(c*Tr**d for c, d in zip(cc, dc))
def temperature_dependence(T_centigrade):
return thermal_conductivity(T_centigrade)**(2/3) * viscosity(T_centigrade)**(-7/15)
def phi(T_centigrade):
return temperature_dependence(T_centigrade)/temperature_dependence(20)
T = np.linspace(18, 55)
fig, ax = plt.subplots(4, 1, figsize=(10, 8), sharex=True)
ax[0].plot(T, list(map(viscosity, T)))
ax[0].set_title("Viscosity")
ax[0].set_ylabel("Pascal-seconds")
ax[0].grid(True)
ax[1].plot(T, list(map(thermal_conductivity, T)))
ax[1].set_title("Thermal Conductivity")
ax[1].set_ylabel("J/kg/K")
ax[1].grid(True)
ax[2].plot(T, list(map(temperature_dependence, T)))
ax[2].set_title("Temperature Dependence")
ax[2].set_xlabel("Temperature deg C")
ax[2].grid(True)
ax[3].plot(T, list(map(phi, T)))
ax[3].set_title("phi(T)")
ax[3].set_xlabel("Temperature deg C")
ax[3].grid(True)
plt.tight_layout()

1.6.3. Pyomo DAE Model for Double Pipe Heat Exchanger¶
Let
A model for the double-pipe heat exchanger is given by
where
From the Colburn analogy, we assume
where
Two parameters,
Parameter
is used to normalize
import pyomo.environ as pyo
import pyomo.dae as dae
import pandas as pd
import matplotlib.pyplot as plt
def build_hx(y_flow_config=0, qh=600, qc=600, Th_feed=55.0, Tc_feed=18.0):
"""
Parameters:
y_flow_config: 0 for counter-current, 1 for co-current
qh (float): hot water flowrate in liters/hour
qc (float): cold water flowrate in liters/hour
Th_feed (float): hot water feed temperature in C
Tc_feed (float): cold water feed temperature in C
Returns a Pyomo model:
m.z: set of positions, 0 to 1.
m.Th[z]: hot water temperature profile
m.Tc[z]: cold water temperature profile
"""
assert (y_flow_config >= 0) and (y_flow_config <= 1), "Unrealizalbe flow configuration"
# known parameter values
A = 0.5 # square meters
Cp = 4184 # Joules/kg/deg C
rho = 1.0 # 1 kg/liter
m = pyo.ConcreteModel()
# estimated model parameters
m.Ch = pyo.Param(mutable=True, default=5000.0)
m.Cc = pyo.Param(mutable=True, default=5000.0)
# operating inputs
m.y = pyo.Param(mutable=True, initialize=y_flow_config)
m.qh = pyo.Param(mutable=True, initialize=qh)
m.qc = pyo.Param(mutable=True, initialize=qc)
m.z = dae.ContinuousSet(bounds=(0, 1))
m.Th = pyo.Var(m.z, initialize=Th_feed)
m.Tc = pyo.Var(m.z, initialize=Tc_feed)
m.Tw = pyo.Var(m.z, initialize=(Th_feed + Tc_feed)/2)
m.Uh = pyo.Var(m.z, initialize=1500)
m.Uc = pyo.Var(m.z, initialize=1500)
m.dTh = dae.DerivativeVar(m.Th, wrt=m.z)
m.dTc = dae.DerivativeVar(m.Tc, wrt=m.z)
# feed conditions
m.Tc[0].fix(Tc_feed)
@m.Constraint()
def flow_configuration(m):
return m.y*m.Th[0] + (1 - m.y)*m.Th[1] == Th_feed
# local heat transfer coefficients
@m.Constraint(m.z)
def heat_transfer_ho(m, z):
return m.Uh[z] == m.Ch * (m.qh/3600)**0.8 * phi(m.Th[z])
@m.Constraint(m.z)
def heat_transfer_cold(m, z):
return m.Uc[z] == m.Cc * (m.qc/3600)**0.8 * phi(m.Tc[z])
# stream energy balances
@m.Constraint(m.z)
def dThdz(m, z):
return rho*Cp*(m.qh/3600)*m.dTh[z] == (1 - 2*m.y) * m.Uh[z]*A*(m.Th[z] - m.Tw[z])
@m.Constraint(m.z)
def dTcdz(m, z):
return rho*Cp*(m.qc/3600)*m.dTc[z] == m.Uc[z]*A*(m.Th[z] - m.Tw[z])
# wall energy balance
@m.Constraint(m.z)
def wall(m, z):
return m.Uh[z]*(m.Th[z] - m.Tw[z]) == m.Uc[z]*(m.Tw[z] - m.Tc[z])
pyo.TransformationFactory('dae.finite_difference').apply_to(m, nfe=20, wrt=m.z)
# function to solve the model
m.solve = lambda: pyo.SolverFactory("ipopt").solve(m)
# function to plot the results
def _visualize(ax):
df = pd.DataFrame({
"Th": [m.Th[z]() for z in m.z],
"Tw": [m.Tw[z]() for z in m.z],
"Tc": [m.Tc[z]() for z in m.z],
"Uh": [m.Uh[z]() for z in m.z],
"Uc": [m.Uc[z]() for z in m.z],
}, index=m.z)
if ax is None:
fig, ax = plt.subplots(1, 1)
df.plot(y=["Th", "Tc", "Tw"], ax=ax, grid=True, lw=3)
df.plot(y=["Uh", "Uc"], ax=ax, secondary_y=True, grid=True, lw=3, style=['--', '--'])
ax.set_xlabel("Position")
ax.set_ylabel("deg C")
ax.set_title(f"Qc = {qc}, Qh = {qh}")
plt.tight_layout()
return ax
m.visualize = lambda ax=None: _visualize(ax)
return m
hx_model = build_hx(y_flow_config=0, qc=500, qh=500)
hx_model.solve()
hx_model.visualize()
<AxesSubplot:title={'center':'Qc = 500, Qh = 500'}, xlabel='Position', ylabel='deg C'>

1.6.4. Comparing simulations to experimental data¶
# guess parameter values for heat transfer model
Cc = 4000
Ch = 5000
# set grid of flowrates
Qh = [300, 450, 600]
Qc = [250, 400, 550]
# set grid of plots
fig, ax = plt.subplots(len(Qc), len(Qh), figsize=(15, 10))
# do simulations
for i, qc in enumerate(Qc):
for j, qh in enumerate(Qh):
hx_model = build_hx(y_flow_config=0, qc=qc, qh=qh)
hx_model.Cc = Cc
hx_model.Ch = Ch
hx_model.solve()
hx_model.visualize(ax[i, j])
# overlay temperature measurements
for i, qc in enumerate(Qc):
for j, qh in enumerate(Qh):
z_pos = []
Th_data = []
Tc_data = []
ax[i, j].plot(z_pos, Th_data, 'r.', z_pos, Tc_data, 'b.', ms=10)

1.6.5. Parameter Estimation¶
The model is parameterized by two coefficients
Operationally, the model wilk be used to determine the hot water flow necessary to temper the cold stream to a desired temperature at a desired flowrate. That is, given feed temperatures, flow configuration, desired cold stream flow and cold exit temperature, the model will be used to compute the required hot stream flowrate. The quality of the model will be measured by comparing the model’s prediction for cold exit temperature to actual measurements.
hx_model = build_hx()
hx_model.y = 0
hx_model.solve()
hx_model.visualize()
<AxesSubplot:title={'center':'Heat Transfer Coefficient'}, xlabel='Position', ylabel='kW/m**2/K'>
