# -*- coding: utf-8 -*-
"""
Created on Thu Jul 25 16:42:44 2024
heat-exchanger class
@author: welp
"""
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import root
import matplotlib.pyplot as plt
from carbatpy.models.components import heat_transfer
import time
import carbatpy as cb
import os
import pandas as pd
[docs]
class HeatExchanger:
"""
Static heat exchanger model for counterflow configurations.
Supports double-pipe and monoblock geometries. Heat exchanger profiles are
computed by integrating coupled enthalpy and pressure ODEs along the tube
length using ``scipy.integrate.solve_ivp``.
Two inlet-state conventions are available:
- **One-side IVP** (:meth:`calc_IVP_one_side`): Both inlet states are
provided for the *same* side (inlet–outlet pair). Avoids the internal
root-finding step and is the recommended approach.
- **Two-side IVP** (:meth:`calc_IVP`): True counterflow setup — inlet
states are given for *both* fluid inlets (left and right end). Requires
an internal root-find to determine the unknown starting conditions.
Geometry Options
----------------
- **Double-pipe** (:meth:`geo_double_pipe`): Annular secondary-fluid
channel around an inner tube carrying the working fluid.
- **Monoblock** (:meth:`geo_monoblock`): Two parallel tubes embedded in a
solid conductive block. Thermal coupling is via the block material,
using the VDI Heat Atlas shape factor for two eccentric cylinders.
Parameters
----------
fluids : list of fluid objects
Two-element list ``[working_fluid, secondary_fluid]``. Each fluid must
expose a ``set_state`` method compatible with carbatpy.fprop.
inlet_states : list of array_like
Two-element list ``[state_in_wf, state_in_sf]``. Each state is the
standard carbatpy property array (T, p, h, v, s, q, u, …).
mdot : list of float
Two-element list ``[mdot_wf, mdot_sf]``. Mass flow rates in kg/s.
resolution : int
Number of spatial evaluation points along the heat exchanger length.
ht_method : str
Heat transfer correlation for the working fluid passed to
``heat_transfer.alpha_km``. Example: ``"Deng_et_al-2019"``.
pl_method : str
Pressure-loss correlation for the working fluid passed to
``heat_transfer.alpha_km``. Example: ``"Macdonald_et_al-2016"``.
Attributes
----------
state_array_wf : ndarray, shape (n, n_props)
Fluid-property states along the heat exchanger for the working fluid.
state_array_sf : ndarray, shape (n, n_props)
Fluid-property states along the heat exchanger for the secondary fluid.
x : ndarray
Axial positions [m] corresponding to the state arrays.
alpha_wf : list of float
Local heat transfer coefficients of the working fluid [W/(m²·K)].
alpha_sf : list of float
Local heat transfer coefficients of the secondary fluid [W/(m²·K)].
dpdl_wf : list of float
Local pressure-loss gradients of the working fluid [Pa/m].
dpdl_sf : list of float
Local pressure-loss gradients of the secondary fluid [Pa/m].
termination_flag : str or None
Name of the termination-event function that stopped the ODE solver, or
``None`` if the solver reached the end of the integration interval.
termination_position : float or None
Axial position [m] at which termination occurred, or ``None``.
Examples
--------
Double-pipe condenser solved with the one-side IVP method:
>>> import carbatpy as cb
>>> from class_he import HeatExchanger
>>> fluid1 = cb.init_fluid("Butane * CO2", [0.8, 0.2])
>>> fluid2 = cb.fprop.init_fluid("water", [1])
>>> st1_in = fluid1.set_state([3e6, 1], "PQ")
>>> st2_in = fluid2.set_state([st1_in[0] - 10, 3e5], "TP", cb.fprop._TRANS_STRING)
>>> he = HeatExchanger(
... fluids=[fluid1, fluid2],
... inlet_states=[st1_in, st2_in],
... mdot=[0.01, 0.01],
... resolution=40,
... ht_method="Deng_et_al-2019",
... pl_method="Macdonald_et_al-2016",
... )
>>> he.geo_double_pipe(16e-3, 18e-3, 22e-3, 40e-3, 40, 15.4)
>>> he.calc_IVP_one_side(case="wf_condensed")
>>> he.diagram()
See Also
--------
heat_transfer : Heat transfer and pressure-loss correlations used
internally by :meth:`thermal_resistance`.
"""
def __init__(self, fluids, inlet_states, mdot, resolution, ht_method, pl_method):
[docs]
self.resolution = resolution
[docs]
self.ht_method = ht_method
[docs]
self.pl_method = pl_method
[docs]
def geo_monoblock(self, d_wf, d_sf, r, l, lam_block):
"""
Set monoblock heat exchanger geometry.
Two parallel circular tubes are embedded in a solid conductive block.
The thermal coupling between them is computed via the VDI Heat Atlas
(2013) shape factor for two cylinders in an infinite medium
(Table E1-3).
Parameters
----------
d_wf : float
Inner diameter of the working-fluid tube [m].
d_sf : float
Inner diameter of the secondary-fluid tube [m].
r : float
Centre-to-centre distance between the two tubes [m].
l : float
Length of the block (and tubes) [m].
lam_block : float
Thermal conductivity of the block material [W/(m·K)].
Notes
-----
The conduction shape factor for two cylinders in an infinite medium
(VDI Heat Atlas 2013, E1 Table 3) is:
.. math::
S = \\frac{2\\pi}{\\mathrm{arccosh}\\!
\\left(\\frac{r^2 - r_{wf}^2 - r_{sf}^2}
{2\\, r_{wf}\\, r_{sf}}\\right)}
where :math:`r_{wf} = d_{wf}/2` and :math:`r_{sf} = d_{sf}/2`.
"""
self.geo_flag = "monoblock"
self.d_wf = d_wf
self.d_sf = d_sf
self.r = r
self.l = l
self.lam_block = lam_block
self.alpha_wf = []
self.alpha_sf = []
self.dpdl_wf = []
self.dpdl_sf = []
[docs]
def geo_double_pipe(self, di, Di, da, Da, l, lam_tube):
"""
Set double-pipe heat exchanger geometry.
The working fluid flows through the inner tube; the secondary fluid
flows counter-currently through the annular gap between inner and outer
tube.
Parameters
----------
di : float
Inner diameter of the inner tube (working-fluid side) [m].
Di : float
Outer diameter of the inner tube [m].
da : float
Inner diameter of the outer tube (secondary-fluid side) [m].
Da : float
Outer diameter of the outer tube [m].
l : float
Length of the pipe section [m].
lam_tube : float
Thermal conductivity of the inner-tube wall material [W/(m·K)].
Notes
-----
The total heat transfer area referenced to the inner tube surface is
set to :math:`A = \\pi d_i l`.
The thermal resistance network (per unit length) is:
.. math::
R_{tot} = \\underbrace{\\frac{1}{\\alpha_{sf}\\, \\pi D_i}}_{R_{sf}}
+ \\underbrace{\\frac{\\ln(D_i/d_i)}{2\\pi\\lambda_{tube}}}_{R_{tube}}
+ \\underbrace{\\frac{1}{\\alpha_{wf}\\, \\pi d_i}}_{R_{wf}}
Examples
--------
>>> he.geo_double_pipe(16e-3, 18e-3, 22e-3, 40e-3, 40, 15.4)
"""
self.geo_flag = "double_pipe"
self.di = di
self.Di = Di
self.da = da
self.Da = Da
self.l = l
self.A = np.pi * self.di * self.l
self.lam_tube= lam_tube
self.alpha_wf = []
self.alpha_sf = []
self.dpdl_wf = []
self.dpdl_sf = []
[docs]
def calc_IVP_one_side(self, verbose=False, wf_outlet=False, case=None, stop_value=0):
"""
Calculate the heat exchanger by integrating the governing ODEs.
Both inlet states must be provided for the *same* physical side of the
heat exchanger (i.e. an inlet–outlet pair, not two true inlets). The
``wf_outlet`` flag controls which end of the domain is used as the
working-fluid inlet.
This is the **recommended** calculation method because it avoids the
internal root-finding step required by :meth:`calc_IVP`.
Parameters
----------
verbose : bool, default=False
If ``True``, print the solver status message.
wf_outlet : bool, default=False
Flow-direction flag.
- ``False``: working-fluid inlet is at :math:`x = 0`.
- ``True``: working-fluid inlet is at :math:`x = l` (outlet side).
case : str or None, default=None
Optional early-termination criterion. Supported values:
- ``'sf_temperature_min'`` – stop when secondary-fluid temperature
falls below ``stop_value`` [K].
- ``'sf_enthalpy_min'`` – stop when secondary-fluid enthalpy falls
below ``stop_value`` [J/kg].
- ``'wf_enthalpy_min'`` – stop when working-fluid enthalpy falls
below ``stop_value`` [J/kg].
- ``'wf_enthalpy_max'`` – stop when working-fluid enthalpy exceeds
``stop_value`` [J/kg].
- ``'wf_condensed'`` – stop when working fluid is fully condensed
(vapour quality ≤ saturated-liquid enthalpy).
- ``'wf_evaporated'`` – stop when working fluid is fully evaporated
(vapour quality ≥ saturated-vapour enthalpy).
stop_value : float, default=0
Threshold value used by the selected ``case``. Unit depends on the
chosen case (K or J/kg).
Returns
-------
None
Results are stored in ``state_array_wf``, ``state_array_sf``,
``x``, ``alpha_wf``, ``alpha_sf``, ``dpdl_wf``, ``dpdl_sf``,
``termination_flag``, and ``termination_position``.
Examples
--------
>>> he.geo_double_pipe(16e-3, 18e-3, 22e-3, 40e-3, 40, 15.4)
>>> he.calc_IVP_one_side(verbose=True, case="wf_condensed")
>>> print(he.termination_flag)
"""
y0 = [self.input_state_wf[2],
self.input_state_sf[2],
self.input_state_wf[1],
self.input_state_sf[1]]
x = np.linspace(0, self.l, self.resolution)
term_events = self._define_termination_events_per_case(case, stop_value)
res = solve_ivp(self._ODE_ivp, args=[wf_outlet], t_span=(x[0], x[-1]), y0=y0,
t_eval=x, rtol=0.5, max_step=0.5,
events=term_events)
self.termination_flag = None
self.termination_position = None
if res.status == 1: # terminated by event
# collect events
hits = []
for i, ev in enumerate(term_events):
if res.t_events[i].size > 0:
hits.append((res.t_events[i][0], ev.__name__))
# earliest event cased termination
if hits:
t_hit, name = min(hits, key=lambda h: h[0])
self.termination_flag = name
self.termination_position = t_hit
if verbose == True:
print(res.message)
self._calc_final_state_arrays(res)
[docs]
def calc_IVP(self, verbose=False, case=None, stop_value=0):
"""
Calculate the heat exchanger using true counterflow inlet conditions.
.. note::
**Not recommended.** Use :meth:`calc_IVP_one_side` instead,
which avoids the internal root-finding step and is more robust.
Inlet states must be the true inlet states for *both* fluids — one at
each end of the heat exchanger. An internal ``scipy.optimize.root``
call iterates on the unknown secondary-fluid state at :math:`x = 0`
until the prescribed secondary-fluid inlet condition at :math:`x = l`
is matched.
Parameters
----------
verbose : bool, default=False
If ``True``, print the root-finder result and solver status.
case : str or None, default=None
Optional early-termination criterion (see :meth:`calc_IVP_one_side`
for supported values).
stop_value : float, default=0
Threshold for the selected ``case``.
Returns
-------
None
Results are stored in ``state_array_wf``, ``state_array_sf``,
``x``, ``alpha_wf``, ``alpha_sf``, ``dpdl_wf``, and ``dpdl_sf``.
Raises
------
ValueError
If the root-finder does not converge.
"""
T_guess = (self.input_state_wf[0] - self.input_state_sf[0]) * 0.9 + \
self.input_state_sf[0] # guesses start temperature for solver
p_guess = self.input_state_sf[1] # don't change, otherwise calculation without
# pressure loss does not converge
result = root(self._solve_ODE_ivp, [T_guess, p_guess], args=(verbose), method='lm', options={'ftol': 1e-2})
if result.success == 0:
raise ValueError(result.message)
T_correct, p_correct = result.x
state_correct_sf = self.fluids[1].set_state([T_correct, p_correct], "TP")
y0 = [self.input_state_wf[2],
state_correct_sf[2],
self.input_state_wf[1],
state_correct_sf[1]]
x = np.linspace(0, np.sum(self.l), self.resolution)
term_events = self._define_termination_events_per_case(case, stop_value)
res = solve_ivp(self._ODE_ivp, args=[False], t_span=(x[0], x[-1]), y0=y0,
t_eval=x, rtol=0.5, max_step = 1,
events=term_events)
self._calc_final_state_arrays(res)
if verbose == True:
print(f"\nroot result: {result.message}")
print(f"T: {T_correct} K, p: {p_correct} Pa")
def _solve_ODE_ivp(self, y_tp, verbose=False):
"""
Residual function for the root-finder in :meth:`calc_IVP`.
Runs the ODE integrator with a guessed secondary-fluid state at
:math:`x = 0` and returns the difference between the resulting outlet
state and the prescribed secondary-fluid inlet condition.
Parameters
----------
y_tp : array_like, shape (2,)
Guessed ``[T_sf, p_sf]`` at :math:`x = 0` [K, Pa].
verbose : bool, default=False
Print diagnostic output when the guess is rejected or when the
solver fails.
Returns
-------
residuals : list of float, length 2
``[ΔT, 0.0001·Δp]`` between the prescribed inlet and the
calculated outlet of the secondary fluid. Returns ``[1e10, 1e10]``
for invalid guesses or solver failures to guide the root-finder
away from unphysical regions.
"""
T_guess, p_guess = y_tp
if not (min(self.input_state_wf[0], self.input_state_sf[0]) < T_guess <
max(self.input_state_wf[0], self.input_state_sf[0])):
if verbose == True:
print(f"\nguessed T, p: {T_guess} K, {p_guess} Pa")
print("T_guess not valid!")
return [10e9, 10e9]
x = np.linspace(0, np.sum(self.l), self.resolution)
state_guess_sf = self.fluids[1].set_state([T_guess, p_guess], "TP")
y0 = [self.input_state_wf[2],
state_guess_sf[2],
self.input_state_wf[1],
state_guess_sf[1]]
res = solve_ivp(self._ODE_ivp, args=[False], t_span=(x[0], x[-1]), y0=y0,
t_eval=x, rtol=0.5, max_step=0.5,
events=[self._event_freezing,
self._event_neg_pressure,
self._event_fail_secfluid,
self._event_fail_workfluid,
self._event_workfluid_max])
if res.t_events and len(res.t_events[1]) > 0:
raise ValueError("pressure loss in tube exceeds inlet pressure")
if res.t_events and len(res.t_events[3]) > 0:
if verbose == True:
print(f"\nguessed T, p: {T_guess} K, {p_guess} Pa")
print("\nset_state of working fluid failed")
print(res.message)
return [10e9, 10e9]
if res.t_events and len(res.t_events[2]) > 0:
if verbose == True:
print(f"\nguessed T, p: {T_guess} K, {p_guess} Pa")
print("\nset_state of secondary fluid failed")
print(res.message)
return [10e9, 10e9]
if res.t_events and len(res.t_events[0]) > 0:
if verbose == True:
print(f"\nguessed T, p: {T_guess} K, {p_guess} Pa")
print("\nwater freezing")
print(res.message)
return [10e9, 10e9]
if res.message != 'The solver successfully reached the end of the integration interval.':
if verbose == True:
print(f"\nguessed T, p: {T_guess} K, {p_guess} Pa")
print("\nsolver didn't reach end of integration interval")
print(res.message)
return [10e9, 10e9]
state2 = self.fluids[1].set_state([res.y[1][-1], res.y[3][-1]], "HP")
if verbose == True:
print(f"\nguessed T, p: {T_guess} K, {p_guess} Pa")
print(f"delta_T: {self.input_state_sf[0] - state2[0]} K, delta_p: {self.input_state_sf[1] - state2[1]} Pa")
print(res.message)
return [self.input_state_sf[0] - state2[0],
0.0001 * (self.input_state_sf[1] - state2[1])] # factor for pressure deviation!
def _define_termination_events_per_case(self, case, stop_value):
termination_events = [self._event_freezing,
self._event_neg_pressure,
self._event_fail_secfluid,
self._event_fail_workfluid]
if case=='sf_temperature_min':
self.stop_value = stop_value
termination_events.append(self._event_sf_temperature_min)
elif case=='wf_enthalpy_min':
self.stop_value = stop_value
termination_events.append(self._event_wf_enthalpy_min)
elif case =='wf_enthalpy_max':
self.stop_value = stop_value
termination_events.append(self._event_wf_enthalpy_max)
elif case =='sf_enthalpy_min':
self.stop_value = stop_value
termination_events.append(self._event_sf_enthalpy_min)
elif case == 'wf_condensed':
termination_events.append(self._event_wf_condensed)
elif case == 'wf_evaporated':
termination_events.append(self._event_wf_evaporated)
return termination_events
def _calc_final_state_arrays(self, res):
state1 = []
state2 = []
for i in range(len(res.t)):
state1.append(self.fluids[0].set_state([res.y[0][i], res.y[2][i]], "HP"))
state2.append(self.fluids[1].set_state([res.y[1][i], res.y[3][i]], "HP",
cb.fprop._TRANS_STRING))
self.state_array_wf = np.array(state1)
self.state_array_sf = np.array(state2)
self.x = res.t
self.alpha_wf = []
self.alpha_sf = []
self.dpdl_wf = []
self.dpdl_sf = []
for i in range(len(self.state_array_wf)):
self.thermal_resistance(self.state_array_wf[i, :], self.state_array_sf[i, :])
def _ODE_ivp(self, x, y, wf_outlet=False):
"""
Right-hand side of the heat exchanger ODE system for ``solve_ivp``.
Solves the coupled first-order ODEs for specific enthalpy and pressure
of both fluids along the axial coordinate :math:`x`:
.. math::
\\frac{dh_{wf}}{dx} &= \\frac{\\Delta T}{R_{tot}\\,\\dot{m}_{wf}} \\\\
\\frac{dh_{sf}}{dx} &= \\frac{\\Delta T}{R_{tot}\\,\\dot{m}_{sf}} \\\\
\\frac{dp_{wf}}{dx} &= -\\left(\\frac{dp}{dl}\\right)_{wf} \\\\
\\frac{dp_{sf}}{dx} &= +\\left(\\frac{dp}{dl}\\right)_{sf}
where :math:`\\Delta T = T_{sf} - T_{wf}` and :math:`R_{tot}` is the
local total thermal resistance per unit length returned by
:meth:`thermal_resistance`.
Parameters
----------
x : float
Current axial position [m].
y : array_like, shape (4,)
State vector ``[h_wf, h_sf, p_wf, p_sf]``.
wf_outlet : bool, default=False
If ``True``, all derivatives are negated so that integration
proceeds from the working-fluid outlet towards its inlet.
Returns
-------
dy : ndarray, shape (4,)
Derivatives ``[dh_wf/dx, dh_sf/dx, dp_wf/dx, dp_sf/dx]``.
Raises
------
ValueError
If ``p_wf`` becomes negative, or if a fluid-property call fails.
"""
h_wf, h_sf, p_wf, p_sf = y
if p_wf < 0:
raise ValueError("pressure loss exceeds inlet pressure")
try:
state_wf_x = (self.fluids[0].set_state([h_wf, p_wf], "HP"))
except:
raise ValueError(f"At x {x} h_wf {h_wf} and p_wf {p_wf} failed Refprop")
try:
state_sf_x = (self.fluids[1].set_state([h_sf, p_sf], "HP", cb.fprop._TRANS_STRING))
except:
raise ValueError(f"At x {x} h_sf {h_sf} and p_sf {p_sf} failed Refprop")
R_ges, dp_wf, dp_sf= self.thermal_resistance(state_wf_x, state_sf_x)
delta_T = state_sf_x[0] - state_wf_x[0]
dh_wf = 1/R_ges * delta_T / self.mdot_wf
dh_sf = 1/R_ges * delta_T / self.mdot_sf
if np.isnan(dh_wf):
raise ValueError(f"At x {x} h_wf {h_wf} is nan")
if wf_outlet==True:
return np.array([float(-dh_wf), float(-dh_sf), float(dp_wf), float(-dp_sf)])
return np.array([float(dh_wf), float(dh_sf), float(-dp_wf), float(dp_sf)])
# termination events to avoid solver from crashing when running into
# unphysical regime (mainly caused by calling fluid properties)
def _event_freezing(self, x, y, wf_outlet):
# secondary fluid reaches solid state
h_wf, h_sf, p_wf, p_sf = y
state_sf_x = (self.fluids[1].set_state([h_sf, p_sf], "HP"))
if self.fluids[1].fluidmodel.fluid == 'water' or self.fluids[1].fluidmodel.fluid == 'Water':
T_min = 280
elif self.fluids[1].fluidmodel.fluid == "ethylene glycol" or self.fluids[1].fluidmodel.fluid == "Ethylene Glycol":
T_min = 259
elif self.fluids[1].fluidmodel.fluid == "MD4M":
T_min = 250
else:
T_min = 250
if state_sf_x[0] < T_min:
return 0
else:
return 10
def _event_neg_pressure(self, x, y, wf_outlet):
# pressure loss to high
if y[2] < 1e2:
return 0
else:
return 10
def _event_fail_secfluid(self, x, y, wf_outlet):
# set_state for secondary fluid fails
try:
h_wf, h_sf, p_wf, p_sf = y
state_sf_x = (self.fluids[1].set_state([h_sf, p_sf], "HP"))
return 10
except:
return 0
def _event_fail_workfluid(self, x, y, wf_outlet):
# set_state for working fluid fails
try:
h_wf, h_sf, p_wf, p_sf = y
state_sf_x = (self.fluids[0].set_state([h_wf, p_wf], "HP"))
return 10
except:
return 0
def _event_workfluid_max(self, x, y, wf_outlet):
# maximal temperature for working fluid reached
h_wf, h_sf, p_wf, p_sf = y
state_wf_x = (self.fluids[1].set_state([h_wf, p_wf], "HP"))
if state_wf_x[0] > 500:
return 0
else:
return 10
_event_fail_workfluid.terminal = True
_event_fail_secfluid.terminal = True
_event_freezing.terminal = True
_event_neg_pressure.terminal = True
_event_workfluid_max.terminal = True
# termination events for use-cases
def _event_wf_condensed(self, x, y, wf_outlet):
h_wf, h_sf, p_wf, p_sf = y
state_wf_x = (self.fluids[0].set_state([h_wf, p_wf], "HP"))
state_wf_0 = (self.fluids[0].set_state([p_wf, 0], "PQ"))
if state_wf_x[2] <= state_wf_0[2]:
return 0
else:
return 10
def _event_wf_evaporated(self, x, y, wf_outlet):
h_wf, h_sf, p_wf, p_sf = y
state_wf_x = (self.fluids[0].set_state([h_wf, p_wf], "HP"))
state_wf_1 = (self.fluids[0].set_state([p_wf, 1], "PQ"))
if state_wf_x[2] >= state_wf_1[2]:
return 0
else:
return 10
def _event_sf_temperature_min(self, x, y, wf_outlet):
h_wf, h_sf, p_wf, p_sf = y
state_sf_x = (self.fluids[1].set_state([h_sf, p_sf], "HP"))
if state_sf_x[0] <= self.stop_value:
return 0
else:
return 10
def _event_wf_enthalpy_min(self, x, y, wf_outlet):
h_wf, h_sf, p_wf, p_sf = y
state_wf_x = (self.fluids[0].set_state([h_wf, p_wf], "HP"))
if state_wf_x[2] <= self.stop_value:
return 0
else:
return 10
def _event_wf_enthalpy_max(self, x, y, wf_outlet):
h_wf, h_sf, p_wf, p_sf = y
state_wf_x = (self.fluids[0].set_state([h_wf, p_wf], "HP"))
if state_wf_x[2] >= self.stop_value:
return 0
else:
return 10
def _event_sf_enthalpy_min(self, x, y, wf_outlet):
h_wf, h_sf, p_wf, p_sf = y
state_sf_x = (self.fluids[1].set_state([h_sf, p_sf], "HP"))
if state_sf_x[2] <= self.stop_value:
return 0
else:
return 10
# all events defined as termination events
_event_wf_condensed.terminal = True
_event_wf_evaporated.terminal = True
_event_sf_temperature_min.terminal = True
_event_wf_enthalpy_min.terminal = True
_event_wf_enthalpy_max.terminal = True
_event_sf_enthalpy_min.terminal = True
[docs]
def thermal_resistance(self, state_wf_x, state_sf_x):
"""
Calculate the local total thermal resistance and pressure-loss gradients.
Dispatches to the correct resistance model based on ``self.geo_flag``
(``"double_pipe"`` or ``"monoblock"``), appends the local heat transfer
coefficients and pressure-loss gradients to the instance lists, and
updates ``self.dT_wall``.
Parameters
----------
state_wf_x : array_like
Local carbatpy property array of the working fluid (T, p, h, …).
state_sf_x : array_like
Local carbatpy property array of the secondary fluid (T, p, h, …).
Returns
-------
R_ges : float
Total thermal resistance per unit length [(m·K)/W].
dp_wf : float
Pressure-loss gradient of the working fluid [Pa/m].
dp_sf : float
Pressure-loss gradient of the secondary fluid [Pa/m].
Notes
-----
**Double-pipe** resistance network (per unit length):
.. math::
R_{tot} = \\frac{1}{\\alpha_{sf}\\,\\pi D_i}
+ \\frac{\\ln(D_i/d_i)}{2\\pi\\lambda_{tube}}
+ \\frac{1}{\\alpha_{wf}\\,\\pi d_i}
**Monoblock** resistance network (per unit length):
.. math::
R_{tot} = \\frac{1}{\\alpha_{sf}\\,\\pi d_{sf}}
+ \\frac{1}{\\lambda_{block}\\, S}
+ \\frac{1}{\\alpha_{wf}\\,\\pi d_{wf}}
with the VDI shape factor :math:`S` from :meth:`geo_monoblock`.
"""
fluidstate_wf_x = self.fluids[0].set_state([state_wf_x[2],
state_wf_x[1]], "HP",
output="FluidState")
fluidstate_sf_x = self.fluids[1].set_state([state_sf_x[0],
state_sf_x[1]], "TP",
cb.fprop._TRANS_STRING,
output="FluidState")
if self.geo_flag == "double_pipe":
alpha_sf, dp_sf = heat_transfer.calc_alpha_dpdl_tube(self.fluids[1], fluidstate_sf_x, self.mdot_sf,
self.l, self.da, self.Di)
alpha_wf, dp_wf = heat_transfer.calc_alpha_dpdl_tube(self.fluids[0], fluidstate_wf_x,
self.mdot_wf,
self.l, self.di, pl_method=self.pl_method, ht_method= self.ht_method, dT_wall=self.dT_wall)
self.alpha_wf.append(alpha_wf)
self.alpha_sf.append(alpha_sf)
self.dpdl_wf.append(dp_wf)
self.dpdl_sf.append(dp_sf)
# heat transfer resistance
R_sf = 1/(alpha_sf * self.Di * np.pi) # convection secondary fluid
R_tube = np.log(self.Di/self.di) / (2 * self.lam_tube * np.pi) # conduction tube
R_wf = 1/(alpha_wf * self.di * np.pi) # convection working fluid
R_ges = R_sf + R_tube + R_wf
if self.geo_flag == "monoblock":
alpha_sf, dp_sf = heat_transfer.calc_alpha_dpdl_tube(self.fluids[1], fluidstate_sf_x, self.mdot_sf,
self.l, self.d_sf)
alpha_wf, dp_wf = heat_transfer.calc_alpha_dpdl_tube(self.fluids[0], fluidstate_wf_x,
self.mdot_wf,
self.l, self.d_wf, pl_method=self.pl_method, ht_method= self.ht_method, dT_wall=self.dT_wall)
self.alpha_wf.append(alpha_wf)
self.alpha_sf.append(alpha_sf)
self.dpdl_wf.append(dp_wf)
self.dpdl_sf.append(dp_sf)
# heat transfer resistance
R_sf = 1 / (alpha_sf * self.d_sf * np.pi) # convection secondary fluid
S = 2 * np.pi / \
np.arccosh(
(self.r ** 2 - (self.d_wf/2) ** 2 - (self.d_sf/2) ** 2) / \
(2 * self.d_wf / 2 * self.d_sf / 2)) # form factor two tubes in medium, VDI heat atlas 2013 E1-table 3 3)
R_block = 1 / (self.lam_block * S)
R_wf = 1 / (alpha_wf * self.d_wf * np.pi) # convection working fluid
R_ges = R_sf + R_block + R_wf
self.dT_wall = R_wf / R_ges * (fluidstate_wf_x.temperature - fluidstate_sf_x.temperature)
return R_ges, dp_wf, dp_sf
[docs]
def diagram(self, ordinate=0, second_yaxis=False, save_dir=None, show_plot=False):
"""
Plot fluid-property profiles along the heat exchanger length.
The abscissa is the axial coordinate :math:`x` [m]. The ordinate is
selected by index, matching the column order of ``state_array_wf`` /
``state_array_sf``:
- 0: Temperature [K]
- 1: Pressure [Pa]
- 2: Specific enthalpy [J/kg]
- 3: Specific volume [m³/kg]
- 4: Specific entropy [J/(kg·K)]
- 5: Vapour quality [–]
- 6: Specific internal energy [J/kg]
- 7: Dynamic viscosity [Pa·s]
- 8: Thermal conductivity [W/(m·K)]
- 9: Prandtl number [–]
- 10: Kinematic viscosity [m²/s]
- 11: Molar mass [kg/mol]
- 12: Speed of sound [m/s]
Parameters
----------
ordinate : int, default=0
Column index of the property to plot (see list above).
second_yaxis : bool, default=False
If ``True``, working and secondary fluid are plotted on separate
y-axes (useful when the two fluids have very different scales).
save_dir : str or None, default=None
Directory path for saving the figure as a PNG file. If ``None``,
the figure is only displayed.
show_plot : bool, default=False
If "True", plots are shown.
Returns
-------
None
Examples
--------
>>> he.calc_IVP_one_side()
>>> he.diagram() # temperature profile
>>> he.diagram(ordinate=1, second_yaxis=True) # pressure, dual axes
>>> he.diagram(ordinate=2, save_dir=r"C:/results") # save enthalpy plot
"""
fig, ax1 = plt.subplots()
x_vals = self.x
ax1.plot(x_vals, self.state_array_wf[:, ordinate], 'r:o',
label='Working Fluid', markersize=5)
xlabel_string = 'x [m]'
ax1.set_xlabel(xlabel_string)
if second_yaxis == True:
ax1.tick_params(axis='y', labelcolor='r')
ax1.set_ylabel(cb.CB_DEFAULTS["Fluid_Defaults"]['Property_Names'][ordinate],
color='r')
ax2 = ax1.twinx()
# Plot the secondary fluid on the second y-axis
ax2.plot(x_vals, self.state_array_sf[:, ordinate], 'b--s',
label='Secondary Fluid', markersize=5)
ax2.set_ylabel(cb.CB_DEFAULTS["Fluid_Defaults"]['Property_Names'][ordinate],
color='b')
ax2.tick_params(axis='y', labelcolor='b')
else:
ax1.set_ylabel(cb.CB_DEFAULTS["Fluid_Defaults"]['Property_Names'][ordinate])
ax1.plot(x_vals, self.state_array_sf[:, ordinate], 'b--s',
label='Secondary Fluid', markersize=5)
plt.title(cb.CB_DEFAULTS["Fluid_Defaults"]['Property_Names'][ordinate] \
+ ' vs Area')
fig.legend(loc="upper right", bbox_to_anchor=(1,1),
bbox_transform=ax1.transAxes)
ax1.grid(True)
if save_dir is not None:
file_path = os.path.join(save_dir, cb.CB_DEFAULTS["Fluid_Defaults"]['Property_Names'][ordinate] + "_x_"+ ".png")
plt.savefig(file_path)
if show_plot:
plt.show()
[docs]
def calculate_arrays(self, save_dir=None):
"""
Export calculation results to CSV files and a termination-event log.
Writes three files to ``save_dir``:
- ``wf.csv``: Working-fluid state array with columns
``T, p, h, v, s, q, u, alpha_wf, dpdl_wf, tcx, vis``.
``tcx`` and ``vis`` are the saturated-liquid thermal conductivity
and viscosity where the vapour quality is ≥ 0; otherwise 0.
- ``sf.csv``: Secondary-fluid state array with columns
``T, p, h, v, s, q, u, eta, tcx, Pr, cp, ws, M, alpha_sf, dpdl_sf``.
- ``termination_event.txt``: Plain-text log of ``termination_flag``
and ``termination_position``.
Parameters
----------
save_dir : str
Absolute path to an existing directory where the files will be
saved.
Returns
-------
df_wf : DataFrame
df_wf : DataFrame
"""
thermo, trans = heat_transfer.get_correct_trans_string(self.fluids[0])
df_wf = pd.DataFrame(self.state_array_wf)
df_wf.columns = thermo if isinstance(thermo, list) else thermo.split(';')
df_wf['alpha_wf'] = self.alpha_wf
df_wf['dpdl_wf'] = self.dpdl_wf
thermo, trans = heat_transfer.get_correct_trans_string(self.fluids[1])
df_sf = pd.DataFrame(self.state_array_sf)
df_sf.columns = trans if isinstance(trans, list) else trans.split(';')
df_sf['alpha_sf'] = self.alpha_sf
df_sf['dpdl_sf'] = self.dpdl_sf
if save_dir != None:
df_wf.to_csv(save_dir + '\\wf.csv')
df_sf.to_csv(save_dir + '\\sf.csv')
txt_path = save_dir + "\\termination_event.txt"
with open(txt_path, "w", encoding="utf-8") as f:
f.write(f"termination_flag: {self.termination_flag}\n")
f.write(f"termination_position: {self.termination_position}\n")
return df_wf, df_sf
if __name__ == "__main__":
# Create the folder in default results directory
from datetime import datetime
folder_name = now.strftime("%Y-%m-%d_%H-%M-%S")
base_path = cb.CB_DEFAULTS["General"]["RES_DIR"]
new_folder_path = os.path.join(base_path, folder_name)
os.makedirs(new_folder_path)
mdot_1, mdot_2 = 0.02 / 2, 0.01
cells = 40 # resolution
fluid_models = {
"Refprop": {},
#"Gerg-ECS": {"model": "\\gerg-2008", "no": 1, "case": "_ECS"},
#"Gerg-RES": {"model": "\\gerg-2008", "no": 1, "case": "_RES"},
#"PR-ECS": {"model": "", "no": 3, "case": "_ECS"},
#"PR-RES": {"model": "", "no": 3, "case": "_RES"},
#"sim": {"model": "", "no": 1, "case": "_berlin"},
}
ht_method = "Deng_et_al-2019"
pl_method = "Macdonald_et_al-2016"
FLUID1 = "Butane * CO2"
comp1 = [.8, .2]
for key, items in fluid_models.items():
save_dir = new_folder_path + "\\double-pipe\\" + key
os.makedirs(save_dir)
if key == "Refprop":
fluid1 = cb.init_fluid(FLUID1, comp1)
else:
if items["case"] == "_ECS":
cb.CB_DEFAULTS["Fluid_Defaults"]['TRANS_TREND_MIX'] = ['T', 'P', 'H', 'D', 'S', 'QEOS', 'U', 'ETA_ECS',
'TCX_ECS', 'CP', 'WS']
elif items["case"] == "_RES":
cb.CB_DEFAULTS["Fluid_Defaults"]['TRANS_TREND_MIX'] = ['T', 'P', 'H', 'D', 'S', 'QEOS', 'U', 'ETA_RES',
'TCX_RES', 'CP', 'WS']
elif items["case"] == "_berlin":
cb.CB_DEFAULTS["Fluid_Defaults"]['TRANS_TREND_MIX'] = ['T', 'P', 'H', 'D', 'S', 'QEOS', 'U',
'ETA_BERLIN',
'TCX_BERLIN', 'CP', 'WS']
trend_dict = {"trendmodel": items["model"],
"eos_ind": [items["no"], items["no"]],
"mix_ind": items["no"]}
fluid1 = cb.init_fluid(FLUID1, comp1, props="TREND", **trend_dict)
st1_in = fluid1.set_state([3e6, 1], "PQ")
FLUID2 = "water"
comp2 = [1]
fluid2 = cb.fprop.init_fluid(FLUID2, comp2)
st2_in = fluid2.set_state([st1_in[0] - 10, 3e5], "TP", cb.fprop._TRANS_STRING)
print(key)
# for solve IVP
my_doublepipe_he = HeatExchanger([fluid1, fluid2], [st1_in, st2_in],
[mdot_1, mdot_2], cells, pl_method=pl_method, ht_method=ht_method)
my_doublepipe_he.geo_double_pipe(16e-3, 18e-3,
22e-3, 40e-3,
40, 15.4)
st_3 = time.time()
my_doublepipe_he.calc_IVP_one_side(verbose=True, case="wf_condensed")
et_3 = time.time()
my_doublepipe_he.diagram(save_dir=save_dir)
my_doublepipe_he.diagram(1, True, save_dir=save_dir)
my_doublepipe_he.diagram(2, True, save_dir=save_dir)
my_doublepipe_he.calculate_arrays(save_dir)
print(f"solve_ivp-method: {et_3 - st_3} s")
### example monoblock ####
monoblock = True
if monoblock:
save_dir = new_folder_path + "\\monoblock"
os.makedirs(save_dir)
mdot_1, mdot_2 = 0.02/2, 0.04
cells = 40 # resolution
FLUID1 = "Propane * Butane"
comp1 = [.2, .8]
fluid1 = cb.init_fluid(FLUID1, comp1)
st1_in = fluid1.set_state([2.5e5, 1], "PQ")
FLUID2 = "water"
comp2 = [1]
fluid2 = cb.fprop.init_fluid(FLUID2, comp2)
st2_in = fluid2.set_state([st1_in[0]-10, 4e5], "TP", cb.fprop._TRANS_STRING)
my_monoblock_he = HeatExchanger([fluid1, fluid2], [st1_in, st2_in],
[mdot_1, mdot_2], cells, pl_method=pl_method, ht_method=ht_method)
my_monoblock_he.geo_monoblock(16e-3, 18e-3,
22e-3, 60,
15.4)
st_3 = time.time()
my_monoblock_he.calc_IVP_one_side(verbose=True, case="wf_condensed")
et_3 = time.time()
my_monoblock_he.diagram(save_dir=save_dir)
my_monoblock_he.diagram(1, True, save_dir=save_dir)
my_monoblock_he.diagram(2, True, save_dir=save_dir)
my_monoblock_he.calculate_arrays(save_dir=save_dir)
print(f"solve_ivp-method: {et_3-st_3} s")