Finite Element Method for Heat Exchangers¶
Suppose we have a couple of fluid streams exchanging heat (or mass) and want to compute the amount of heat (or mass) exchanged under given conditions of flow configuration and overall heat transfer coefficient U and area through the process. Furthermore, the streams are real -- defined numerically by points (Q,T) or an interpolate, instead of constant stream heat capacity. We want a method that converges with spatial resolution and converges with time for steady-state boundary conditions, and that hopefully gives a reasonable guess for transient conditions. Let's start with some assumptions:
- We want a method that is first order in time, and preferably explicit.
- The method may be spatially implicit.
- The inlets (fluid state and rate of change) are specified as a boundary condition.
- We have all the information we need about U, A, cross-section area, mass flow rate, and density
- It is okay to violate transient conservation of energy insofar as we choose to treat flow as incompressible.
- Assume the fluid state in a given cross-section is uniform, so we are not concerned with spatial resolution of the boundary layer, mixing, etc.
- Assume that only advection acts to carry information downstream, ie. flow is fast enough to neglect heat transfer between points in the same stream.
- Interacting streams are coupled only through local heat transfer, eg. the forcing term.
As a result of some of the assumptions, it is possible for the method to overshoot, so let's keep the time step smallish.
Finite difference schemes¶
Explicit method¶
The obvious starting point is to draw a cell with some $U \Delta A$, write out the transient first law, and extract the relevant equations for the rate of change of state. I'm going to define states at the inlet and outlet to the cell. So we have for one cell
$$ \Delta Q = \dot{m} (h_{out} - h_{in}) + \frac{d}{dt} \int_{cell} u dm $$Now we need to provide an finite element approximation of the mass integral / time rate of change term and estimate $\Delta Q$ resulting from interaction with the other stream. First, let's use a first order central method for the integral term:
\begin{align} \int_{cell} u dm &\approx (\rho A_c \Delta x) \left(\frac{1}{2}\right)\left(u_{in}+u_{out}\right) \\ &= \left(\frac{1}{2}\right) (\alpha \dot{m}) \left(u_{in}+u_{out}\right) \end{align}where for purpose of scaling analysis, we define $\alpha := \rho A_c \Delta x / \dot{m}$, which indicates the thermal inertia of the cell relative to the flow heat capacity, with units of time. We need the whole scheme in terms of one energy variable, so let's substitute in for $u=h-pv$: assuming constant pressure,
\begin{align} \int_{cell} u dm &\approx (\rho A_c \Delta x) \left(\frac{1}{2}\right)\left(h_{in}+h_{out} + p(v_{in} + v_{out}) \right) \\ \end{align}Assume that the cell geometry is constant, although the density may not be. In the time derivative, the specific volume terms cancel, and we have
\begin{align} \frac{d}{dt} \int_{cell} u dm &\approx (\dot{\rho} A_c \Delta x) \left(\frac{1}{2}\right)\left(h_{in}+h_{out} \right) + (\rho A_c \Delta x) \left(\frac{1}{2}\right)\left(\dot{h}_{in}+\dot{h}_{out} \right) \\ \end{align}For this work, let's pretend that the volume change is negligible, ie. we call the fluid incompressible and drop the first term. Let's use an explicit, first order forward difference for the time rate of change:
$$\dot{h} \approx \frac{1}{\Delta t} \left( h(t + \Delta t) - h(t)\right)$$Now we need to write the governing equation with cell and time indices. Let's indicate them like this:
$$ h_{j}^{n} $$where subscript $j$ indexes the node (0 through J) and superscript $n$ indexes the time. Assume the flow and direction of $\Delta x$ is directed from node $j$ to $j+1$. So, associated with each cell from (0) to (J-1) we have
$$ \Delta Q_{(j)}^{n} = \dot{m} (h_{j+1}^{n} - h_{j}^{n}) + \left(\frac{\alpha \dot{m}}{2 \Delta t}\right) (h_{j}^{n+1} - h_{j}^{n} + h_{j+1}^{n+1} - h_{j+1}^{n}) $$The block associated with this cell looks like:
$$ \left[ \Delta Q_{(j)}^n - \dot{m}(h_{j+1}^{n} - h_{j}^{n}) + \left(\frac{\alpha \dot{m}}{2\Delta t}\right) \left(h_{j}^{n} + h_{j+1}^{n}\right) \right] = \left(\frac{\alpha \dot{m}}{2\Delta t}\right) \begin{pmatrix} 1 & 1 \end{pmatrix} \begin{pmatrix} h_{j}^{n+1} \\ h_{j+1}^{n+1} \end{pmatrix} $$Now, we are almost ready to solve. For an inlet boundary condition at $j=0$, we have one more equation that simply specifies $h_0^n$. Now, we start solving with cell j=0, then move to j=1 and solve explicitly, and so on.
Implicit method¶
Following a similar process, we start with the derivative form of the governing equation, and apply the incompressible assumption:
$$ \frac{Q'}{A_c \rho} = \left(\frac{\dot{m}}{A_c \rho}\right) \frac{\partial h}{\partial x} + \frac{\partial h}{\partial t} $$This time, we are write the discrete equations for nodes; thus the cell forcing must be distributed over the nodes in that cell. Also, we use a backwards time approximation for the rate of change term. Anyway,
$$ \frac{Q_{(j)}^{n-1}}{A_c \rho} = \left(\frac{\dot{m}}{A_c \rho}\right) \frac{\left(h_{j}^{n} - h_{j-1}^{n} \right)}{\Delta x} + \frac{\left( h_{j}^{n} - h_{j}^{n-1} \right)}{\Delta t} $$So the block looks like this:
$$ \left[ \frac{\Delta t}{A_c \rho} Q^n_{j} + h_{j}^{n-1} \right] = \begin{pmatrix} -\left(\frac{\dot{m}}{A_c \rho} \frac{\Delta t}{\Delta x} \right) & \left(1 + \frac{\dot{m}}{A_c \rho} \frac{\Delta t}{\Delta x} \right) \end{pmatrix} \begin{pmatrix} h_{j-1}^{n} \\ h_{j}^{n} \end{pmatrix} $$Coupling/Forcing¶
Meanwhile, the local heat transfer rate is
$$ Q = \epsilon Q_{max}, $$where $Q_{max}$ is found for counter-flow by trying to take each stream to equilibrium with the other stream's inlet, and seeing which side is limiting. If we have local estimates of stream heat capacity, so much the better, as we can use the textbook definition. For counter-flow, effectiveness can be computed as
$$ \epsilon(NTU, C_r) = \begin{cases} \frac {1 - \exp[-NTU(1 - C_{r})]}{1 - C_{r}\exp[-NTU(1 - C_{r})]} & C_r < 1 \\[0.5ex] \frac{NTU}{1+NTU} & C_r = 1 \end{cases} $$The pole-zero combination at $C_r = 1$ makes the function less than ideal for numerical work. I will consider this further in another notebook and continue despite the potential issue.
In [1]:
import numpy
from numpy import exp, log, sqrt, power, inf
def counterflow_effectiveness(ntu, c_r):
a = exp(-(ntu * (1.-c_r)))
return numpy.where(c_r == 1,
ntu / (1. + ntu),
(1. - a) / (1. - c_r * a))
def counterflow_effectiveness_bivariate(ntu,c_r):
# As an alternative ...
# This works well for small NTU, and is exact for C_r = 0 or 1.
# However, for large NTU (>= 2), it misses increasingly.
e_c_r_0 = 1 - exp(-ntu)
e_c_r_1 = ntu / (1 + ntu)
alpha = c_r
beta = 1-alpha
return beta * e_c_r_0 + alpha * e_c_r_1
Now let's try some examples. Here are some good cases to test:
- One stream runs against a constant temperature source.
- Two streams start off at uniform (different) temperatures, then come into contact. The inlets are held at the original temperatures.
- Two streams start off at uniform (same) temperatures, then the inlet of one is raised to a different temperature.
- One hot stream runs with fixed inlet, while a cold stream recycles through a mixing tank.
Implementation¶
In [3]:
import CoolProp.CoolProp as CP
import numpy.linalg
import matplotlib
matplotlib.rc_params
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML, Latex
print(matplotlib.__version__)
animation.writers.avail
Out[3]:
In [4]:
matplotlib.rc_params()
?animation.FFMpegWriter
In [5]:
%matplotlib qt5
In [3]:
%matplotlib notebook
In [12]:
for wr in writers.list():
print(wr)
print(writers[wr].supported_formats)
Example 1¶
In [6]:
ac, rho, ua = 1, 1, 10000
dt, dx = 1, 1
m_hot = 1e5
m_cold = 1
cfl_cold = m_cold / (ac * rho) * (dt / dx)
cfl_cold
Out[6]:
In [7]:
p = 1e5
temp_hot_inlet = 350
temp_cold_inlet = 280
h_hot_inlet = CP.PropsSI('H','T',temp_hot_inlet,'P',p,'water')
h_cold_inlet = CP.PropsSI('H','T',temp_cold_inlet,'P',p,'water')
J = 10
ua_cell = ua / J
ua_cell
Out[7]:
In [10]:
# Stiffness matrices
a_cold = numpy.zeros((J+1,J+1))
# Assigning to diagonal doesn't work in numpy v1.9
# a_cold.diagonal(0) = m_cold / dx
# a_cold.diagonal(-1) = ac * rho / dt - m_cold / dx
a_cold[0,0] = 1
for j in range(J):
a_cold[j+1,j+1] += 1 + cfl_cold
a_cold[j+1,j] += -cfl_cold
a_hot = numpy.eye(J+1,J+1)
a_cold_inv = numpy.linalg.inv(a_cold)
a_hot_inv = numpy.linalg.inv(a_hot)
In [11]:
h_hot = numpy.ones(J+1) * h_hot_inlet
h_cold = numpy.ones(J+1) * h_cold_inlet
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
In [13]:
fig=plt.figure()
def init_func():
global fig, line_hot, line_cold
plt.ylim(270,360)
line_hot = plt.plot(temp_hot,'r.-')[0]
line_cold = plt.plot(temp_cold,'b.-')[0]
#for n in range(50):
def update(n):
global h_hot, p, h_cold, temp_hot, temp_cold, ua_cell, ac, rho, dt, line_hot, line_cold
if n < 1:
plt.title("n = 0")
else:
# Compute the forcing terms
c_hot = CP.PropsSI('C','H', h_hot, 'P', p, 'water')
c_cold = CP.PropsSI('C','H', h_cold, 'P', p, 'water')
C_min_cells = m_cold * c_cold[:-1]
C_r_cells = (m_cold * c_cold[:-1]) / (m_hot * c_hot[1:])
Q_max = C_min_cells * (temp_hot[1:] - temp_cold[:-1])
NTU = ua_cell / C_min_cells
effect = counterflow_effectiveness(NTU, C_r_cells)
Q_cells = effect * Q_max
Q_prime_cells = Q_cells / dx
b_cold = numpy.zeros(J+1)
b_cold[0] = h_cold_inlet
b_cold[1:] = Q_prime_cells * dt / (ac * rho) + h_cold[1:]
#b_cold[1:] = Q_prime_cells / (ac * rho) + h_cold[:-1]
b_hot = h_hot
b_hot[-1] = h_hot_inlet
#h_hot = numpy.linalg.solve(a_hot, b_hot)
#h_cold = numpy.linalg.solve(a_cold, b_cold)
h_hot = a_hot_inv.dot(b_hot)
h_cold = a_cold_inv.dot(b_cold)
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
line_cold.set_ydata(temp_cold)
plt.title("n = {}, effectiveness = {:f}"
.format(n, (temp_cold[-1]-temp_cold[0])/(temp_hot_inlet-temp_cold[0])))
#fig.canvas.draw()
plt.draw()
ani = animation.FuncAnimation(fig, update, range(0,21), init_func, repeat=False)
plt.show()
#Writer = animation.FFMpegFileWriter()
#ani.save('fem_hx_example1.mp4',Writer)
#ani.save('fem_hx_example1.gif')
#ani_html5_video = ani.to_html5_video()
ani_js = ani.to_jshtml()
In [ ]:
HTML(ani_js)
In [ ]:
#print('html: ',len(ani_html5_video))
print('js: ',len(ani_js))
import os
print('mp4: ',os.path.getsize('fem_hx_example1.mp4'))
print('gif: ',os.path.getsize('fem_hx_example1.gif'))
In [698]:
total_ntu = ua / c_cold.mean()
total_effect = counterflow_effectiveness(total_ntu, 0)
total_effect
Out[698]:
In [699]:
1-(1-effect).prod()
Out[699]:
In [700]:
a_cold
Out[700]:
Example 2¶
In [53]:
ac, rho, ua = 1, 1, 10000
dt, dx = 5, 10
m_hot = 10
m_cold = 1
cfl_cold = m_cold / (ac * rho) * (dt / dx)
cfl_hot = m_hot / (ac * rho) * (dt / dx)
cfl_cold, cfl_hot
Out[53]:
In [54]:
p = 1e5
temp_hot_inlet = 350
temp_cold_inlet = 280
h_hot_inlet = CP.PropsSI('H','T',temp_hot_inlet,'P',p,'water')
h_cold_inlet = CP.PropsSI('H','T',temp_cold_inlet,'P',p,'water')
J = 5
ua_cell = ua / J
In [55]:
# Stiffness matrices
a_cold = numpy.zeros((J+1, J+1))
# Assigning to diagonal doesn't work in numpy v1.9
# a_cold.diagonal(0) = m_cold / dx
# a_cold.diagonal(-1) = ac * rho / dt - m_cold / dx
a_cold[0,0] = 1
for j in range(J):
a_cold[j+1,j+1] += 1 + cfl_cold
a_cold[j+1,j] += -cfl_cold
a_hot = numpy.zeros((J+1, J+1))
a_hot[-1,-1] = 1
for j in range(J):
a_hot[j,j] += 1 + cfl_hot
a_hot[j,j+1] += -cfl_hot
a_cold_inv = numpy.linalg.inv(a_cold)
a_hot_inv = numpy.linalg.inv(a_hot)
# Do more time steps
#a_cold_inv = numpy.linalg.matrix_power(a_cold_inv, 10)
#a_hot_inv = numpy.linalg.matrix_power(a_hot_inv, 10)
In [62]:
n=0
h_hot = numpy.ones(J+1) * h_hot_inlet
h_cold = numpy.ones(J+1) * h_cold_inlet
#h_hot = numpy.ones(J+1) * h_hot_inlet
#h_cold = numpy.ones(J+1) * h_cold_inlet
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
In [63]:
fig=plt.figure()
plt.ylim(270,360)
line_hot = plt.plot(temp_hot,'r.-')[0]
line_cold = plt.plot(temp_cold,'b.-')[0]
b_cold = numpy.zeros(J+1)
b_hot = numpy.zeros(J+1)
for dn in range(50):
n += 1
# Compute the forcing terms
c_hot = CP.PropsSI('C','H', h_hot, 'P', p, 'water')
c_cold = CP.PropsSI('C','H', h_cold, 'P', p, 'water')
C_cold_cells = m_cold * c_cold[:-1]
C_hot_cells = m_hot * c_hot[1:]
C_min_cells = C_cold_cells
C_r_cells = (m_cold * c_cold[:-1]) / (m_hot * c_hot[1:])
DeltaT_cells = (temp_hot[1:] - temp_cold[:-1])
Q_max = C_min_cells * DeltaT_cells
#NTU = ua_cell / C_min_cells * (1+ 1000*numpy.exp(-1.5*n))
NTU = ua_cell / C_min_cells
effect = counterflow_effectiveness(NTU, C_r_cells)
Q_cells = effect * Q_max
Q_prime_cells = Q_cells / dx
b_cold[0] = h_cold_inlet
b_cold[1:] = Q_prime_cells * dt / (ac * rho) + h_cold[1:]
#b_cold[1:] = Q_prime_cells / (ac * rho) + h_cold[:-1]
b_hot[-1] = h_hot_inlet
b_hot[:-1] = -Q_prime_cells * dt / (ac * rho) + h_hot[:-1]
#b_hot[:-1] = -Q_prime_cells / (ac * rho) + h_hot[1:]
# Update
#h_hot = numpy.linalg.solve(a_hot, b_hot)
#h_cold = numpy.linalg.solve(a_cold, b_cold)
h_hot = a_hot_inv.dot(b_hot)
h_cold = a_cold_inv.dot(b_cold)
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
#plt.plot(temp_cold,'.-')
line_hot.set_ydata(temp_hot)
line_cold.set_ydata(temp_cold)
plt.title("n = {}, effectiveness = {:f}"
.format(n, (temp_cold[-1]-temp_cold[0])/(temp_hot_inlet-temp_cold[0])))
fig.canvas.draw()
fig.canvas.flush_events()
plt.show()
In [712]:
total_ntu = ua / C_cold_cells.mean()
total_cr = C_r_cells.mean()
total_effect = counterflow_effectiveness(total_ntu, total_cr)
total_effect
Out[712]:
In [713]:
DeltaT_cold = numpy.diff(temp_cold)
print(tabulate.tabulate(zip(
range(J),DeltaT_cold,DeltaT_cells,DeltaT_cold/DeltaT_cells,effect)))
#DeltaT_cells
In [255]:
import tabulate
print(tabulate.tabulate(zip(range(J),DeltaT_cells,Q_max, Q_load_cells, Q_prime_cells)))
In [457]:
numpy.linalg.eigvals(a_hot_inv).max()
Out[457]:
Example 3¶
In [64]:
ac, rho, ua = 1, 1, 10000
dt, dx = 0.01, 1
m_hot = 10
m_cold = 1
cfl_cold = m_cold / (ac * rho) * (dt / dx)
cfl_hot = m_hot / (ac * rho) * (dt / dx)
cfl_cold, cfl_hot
Out[64]:
In [65]:
p = 1e5
temp_hot_inlet1 = 280
temp_hot_inlet2 = 350
temp_cold_inlet = 280
h_hot_inlet1 = CP.PropsSI('H','T',temp_hot_inlet1,'P',p,'water')
h_hot_inlet2 = CP.PropsSI('H','T',temp_hot_inlet2,'P',p,'water')
h_cold_inlet = CP.PropsSI('H','T',temp_cold_inlet,'P',p,'water')
J = 20
ua_cell = ua / J
In [66]:
# Stiffness matrices
a_cold = numpy.zeros((J+1, J+1))
# Assigning to diagonal doesn't work in numpy v1.9
# a_cold.diagonal(0) = m_cold / dx
# a_cold.diagonal(-1) = ac * rho / dt - m_cold / dx
a_cold[0,0] = 1
for j in range(J):
a_cold[j+1,j+1] += 1 + cfl_cold
a_cold[j+1,j] += -cfl_cold
a_hot = numpy.zeros((J+1, J+1))
a_hot[-1,-1] = 1
for j in range(J):
a_hot[j,j] += 1 + cfl_hot
a_hot[j,j+1] += -cfl_hot
a_cold_inv = numpy.linalg.inv(a_cold)
a_hot_inv = numpy.linalg.inv(a_hot)
# Do more time steps
a_cold = numpy.zeros((J+1, J+1))
a_cold[0,0] = 1
for j in range(J):
a_cold[j+1,j+1] += 1 + cfl_cold * 10
a_cold[j+1,j] += -cfl_cold * 10
a_hot = numpy.zeros((J+1, J+1))
a_hot[-1,-1] = 1
for j in range(J):
a_hot[j,j] += 1 + cfl_hot * 10
a_hot[j,j+1] += -cfl_hot * 10
a_cold_inv_fast = numpy.linalg.inv(a_cold)
a_hot_inv_fast = numpy.linalg.inv(a_hot)
In [71]:
n=0
h_hot = numpy.ones(J+1) * h_hot_inlet1
h_cold = numpy.ones(J+1) * h_cold_inlet
#h_hot = numpy.ones(J+1) * h_hot_inlet
#h_cold = numpy.ones(J+1) * h_cold_inlet
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
In [73]:
fig=plt.figure()
plt.ylim(270,360)
line_hot = plt.plot(temp_hot,'r.-')[0]
line_cold = plt.plot(temp_cold,'b.-')[0]
b_cold = numpy.zeros(J+1)
b_hot = numpy.zeros(J+1)
for dn in range(50):
n += 1
speed_up = n > 20
h_hot_inlet = h_hot_inlet2 if (n >10) else h_hot_inlet1
# Compute the forcing terms
c_hot = CP.PropsSI('C','H', h_hot, 'P', p, 'water')
c_cold = CP.PropsSI('C','H', h_cold, 'P', p, 'water')
C_cold_cells = m_cold * c_cold[:-1]
C_hot_cells = m_hot * c_hot[1:]
C_min_cells = C_cold_cells
C_r_cells = (m_cold * c_cold[:-1]) / (m_hot * c_hot[1:])
DeltaT_cells = (temp_hot[1:] - temp_cold[:-1])
Q_max = C_min_cells * DeltaT_cells
#NTU = ua_cell / C_min_cells * (1+ 1000*numpy.exp(-1.5*n))
NTU = ua_cell / C_min_cells
effect = counterflow_effectiveness(NTU, C_r_cells)
Q_cells = effect * Q_max * (10 if speed_up else 1)
Q_prime_cells = Q_cells / dx
b_cold[0] = h_cold_inlet
b_cold[1:] = Q_prime_cells * dt / (ac * rho) + h_cold[1:]
#b_cold[1:] = Q_prime_cells / (ac * rho) + h_cold[:-1]
b_hot[-1] = h_hot_inlet
b_hot[:-1] = -Q_prime_cells * dt / (ac * rho) + h_hot[:-1]
#b_hot[:-1] = -Q_prime_cells / (ac * rho) + h_hot[1:]
# Update
h_hot = a_hot_inv.dot(b_hot)
h_cold = a_cold_inv.dot(b_cold)
# Use this to speed things up
if speed_up:
h_hot = a_hot_inv_fast.dot(b_hot)
h_cold = a_cold_inv_fast.dot(b_cold)
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
#plt.plot(temp_cold,'.-')
line_hot.set_ydata(temp_hot)
line_cold.set_ydata(temp_cold)
plt.title("n = {}, effectiveness = {:f}"
.format(n, (temp_cold[-1]-temp_cold[0])/(temp_hot_inlet-temp_cold[0])))
fig.canvas.draw()
fig.canvas.flush_events()
plt.show()
Observe that the solution exhibits numerical diffusion; however, with stationary boundary conditions, it still converges to the same steady-state.
In [ ]:
Example 4¶
In [74]:
rho, ua = 1, 10000
dt, dx = 1, 1
ac_hot, m_hot = 10, 10
ac_cold, m_cold = 1, 1
cfl_cold = m_cold / (ac_cold * rho) * (dt / dx)
cfl_hot = m_hot / (ac_hot * rho) * (dt / dx)
cfl_cold, cfl_hot
Out[74]:
In [75]:
p = 10e5
temp_hot_inlet1 = 350
temp_hot_inlet2 = 450
temp_cold_inlet = 280
h_hot_inlet1 = CP.PropsSI('H','T',temp_hot_inlet1,'P',p,'water')
h_hot_inlet2 = CP.PropsSI('H','T',temp_hot_inlet2,'P',p,'water')
h_cold_inlet = CP.PropsSI('H','T',temp_cold_inlet,'P',p,'water')
J = 20
ua_cell = ua / J
In [76]:
# Stiffness matrices
a_cold = numpy.zeros((J+1, J+1))
# Assigning to diagonal doesn't work in numpy v1.9
# a_cold.diagonal(0) = m_cold / dx
# a_cold.diagonal(-1) = ac * rho / dt - m_cold / dx
a_cold[0,0] = 1
for j in range(J):
a_cold[j+1,j+1] += 1 + cfl_cold
a_cold[j+1,j] += -cfl_cold
a_hot = numpy.zeros((J+1, J+1))
a_hot[-1,-1] = 1
for j in range(J):
a_hot[j,j] += 1 + cfl_hot
a_hot[j,j+1] += -cfl_hot
a_cold_inv = numpy.linalg.inv(a_cold)
a_hot_inv = numpy.linalg.inv(a_hot)
In [77]:
n=0
h_hot = numpy.ones(J+1) * h_hot_inlet1
h_cold = numpy.ones(J+1) * h_cold_inlet
#h_hot = numpy.ones(J+1) * h_hot_inlet
#h_cold = numpy.ones(J+1) * h_cold_inlet
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
m_tank = 50
H_tank = m_tank * h_cold_inlet
load_on = False
live_plotting = False
In [78]:
%matplotlib qt5
live_plotting = True
In [79]:
steps = 1000
fig=plt.figure()
ax1 = fig.add_subplot(211)
plt.ylim(270,450)
line_hot = plt.plot(temp_hot,'r.-')[0]
line_cold = plt.plot(temp_cold,'b.-')[0]
ax2 = fig.add_subplot(212)
plt.ylim(270,370)
plt.xlim(n,n+steps)
T_tank_history = numpy.empty(steps)
T_tank_history.fill(numpy.nan)
line_tank = plt.plot(range(n,n+steps),T_tank_history,'b.-')[0]
plt.show()
fig.canvas.draw()
flushing = True
try:
fig.canvas.flush_events()
except:
flushing = False
b_cold = numpy.zeros(J+1)
b_hot = numpy.zeros(J+1)
try:
for dn in range(steps):
n += 1
h_hot_inlet = h_hot_inlet2 if ((n-1) // 20) % 2 else h_hot_inlet1
# Compute the forcing terms
c_hot = CP.PropsSI('C','H', h_hot, 'P', p, 'water')
c_cold = CP.PropsSI('C','H', h_cold, 'P', p, 'water')
C_cold_cells = m_cold * c_cold[:-1]
C_hot_cells = m_hot * c_hot[1:]
C_min_cells = C_cold_cells
C_r_cells = (m_cold * c_cold[:-1]) / (m_hot * c_hot[1:])
DeltaT_cells = (temp_hot[1:] - temp_cold[:-1])
Q_max = C_min_cells * DeltaT_cells
#NTU = ua_cell / C_min_cells * (1+ 1000*numpy.exp(-1.5*n))
NTU = ua_cell / C_min_cells
effect = counterflow_effectiveness(NTU, C_r_cells)
Q_cells = effect * Q_max
Q_prime_cells = Q_cells / dx
# Tank dynamics
H_tank += dt * m_cold * (h_cold[-1] - h_cold[0])
if load_on:
H_tank -= 500e3 * dt
h_tank = H_tank / m_tank
T_tank = CP.PropsSI('T','H',h_tank,'P',p,'water')
T_tank_history[dn] = T_tank
if load_on:
if T_tank < 320:
load_on = False
else:
if T_tank > 350:
load_on = True
b_cold[0] = h_tank
b_cold[1:] = Q_prime_cells * dt / (ac_cold * rho) + h_cold[1:]
#b_cold[1:] = Q_prime_cells / (ac * rho) + h_cold[:-1]
b_hot[-1] = h_hot_inlet
b_hot[:-1] = -Q_prime_cells * dt / (ac_hot * rho) + h_hot[:-1]
#b_hot[:-1] = -Q_prime_cells / (ac * rho) + h_hot[1:]
# Update
h_hot = a_hot_inv.dot(b_hot)
h_cold = a_cold_inv.dot(b_cold)
temp_hot = CP.PropsSI('T','H', h_hot, 'P', p, 'water')
temp_cold = CP.PropsSI('T','H', h_cold, 'P', p, 'water')
if live_plotting or divmod(dn,100)[1] == 0:
line_hot.set_ydata(temp_hot)
line_cold.set_ydata(temp_cold)
ax1.set_title("n = {}, eff = {:0.3f}, T_tank = {:0.1f}, Load = {}"
.format(n,
(temp_cold[-1]-temp_cold[0])/(temp_hot[-1]-temp_cold[0]),
T_tank,
"on" if load_on else "off"))
line_tank.set_ydata(T_tank_history)
fig.canvas.draw()
if flushing:
fig.canvas.flush_events()
except KeyboardInterrupt:
pass
line_hot.set_ydata(temp_hot)
line_cold.set_ydata(temp_cold)
ax1.set_title("n = {}, eff = {:0.3f}, T_tank = {:0.1f}, Load = {}"
.format(n,
(temp_cold[-1]-temp_cold[0])/(temp_hot[-1]-temp_cold[0]),
T_tank,
"on" if load_on else "off"))
line_tank.set_ydata(T_tank_history)
plt.show()
Interesting! There is a load that draws a fixed amount of heat from the tank, and is cycling on and off based on temperature in the tank. The time it takes to recharge and discharge the tank does not match the frequency of the intermittent heat source, so we get charge/discharge periods of different lengths:
In [124]:
#charging = numpy.where(numpy.diff(T_tank_history) > 0, 1, 0)
charging = numpy.diff(T_tank_history) > 0
switch_times = numpy.concatenate(([0], numpy.diff(charging).nonzero()[0]))
periods = numpy.diff(switch_times)
charge_periods = periods[::2]
drain_periods = periods[1::2]
In [125]:
plt.figure()
plt.hist([charge_periods,drain_periods])
plt.show()
In [10]:
%matplotlib --list
In [1]:
import matplotlib.pyplot as plt
fig=plt.figure()
fig.canvas
Out[1]:
In [ ]: