Finite Volume Method for Heat Transfer Problem

In this section, we will apply the finite volume method to a one-dimensional initial and boundary value problem describing heat transfer by conduction and convection in a fluid and its cooling by natural convection.

We consider a circular pipe of length $L$ (1 m) and internal diameter $d_i$ (0.0508 m) transporting a cryogenic liquid, in this case, liquid nitrogen (LN2) at 70 K and atmospheric pressure (101325 Pa). The pipe has a very thin wall, so the external and internal diameters are approximately equal. The fluid enters the pipe at $x=0$ with an inlet temperature $T_{in}$ (K) of 70 K and heats up by natural convection as it flows along the pipe since the ambient temperature $T_{\text{air}}$ (298.15 K) is higher than the fluid temperature. The velocity profile of the fluid inside the pipe reaches a steady-state profile. For simplicity, we will assume a plug flow profile, meaning it is independent of the axial and radial coordinates and has a fixed value of $v_0$ (0.1 m $\textrm{s}^{-1}$).

The thermophysical properties of the fluid and problem parameters:

• Natural convection heat transfer coefficient, h = 25 $\text{ W m}^{-2}~\text{K}^{-1}$
• Fluid density, $\rho$ = 838.645 $\text{kg m}^{-3}$
• Viscosity, $\mu=$ 0.00022 Pa s
• Thermal conductivity, k = 0.16 $\text{W m}^{-1}~\text{K}^{-1}$
• Specific heat at constant pressure, $\hat{c}_p = 2014.03 \text{ J kg}^{-1}~\text{K}^{-1}$

Differential Equation and Boundary and Initial Conditions

First, we write the differential equation describing the change in internal energy (or heat energy equation) per unit volume:

$$ \begin{alignat}{5} &\frac{\partial (\rho \hat{\mathcal{U}})}{\partial t} + \left(\rho\mathbf{v}\cdot\nabla \hat{\mathcal{U}}\right) = && && && &&\qquad\text{(Rate of change of internal energy)}\ & && -\left(\nabla\cdot \mathbf{q}\right) && && &&\qquad\text{(Net transport of internal energy by conduction)}\ & && && -p\left(\nabla\cdot \mathbf{v}\right) && &&\qquad\text{(Reversible work done on fluid)}\ & && && && -\left(\mathbf{\tau}:\nabla\mathbf{v}\right) &&\qquad\text{(Irreversible work done by viscous dissipation)} \end{alignat} $$

Using the thermodynamic definition of internal energy, $\hat{\mathcal{U}} = \hat{\mathcal{H}} - p\hat{V}$, and considering that enthalpy can be expressed as a function of temperature and pressure, $d\hat{\mathcal{H}} = \hat{C}_p dT +\left[\hat{V} - T\left(\frac{\partial \hat{V}}{\partial T}\right)_p\right]dp$, we can express the internal energy change equation in terms of temperature:

$$ \rho \hat{C}_p \frac{\partial T}{\partial t} + \rho\hat{C}_p\mathbf{v}\cdot\nabla T = -\left(\nabla\cdot \mathbf{q}\right) - \left(\frac{\partial \ln{\rho}}{\partial\ln{T}}\right)_p\left(\frac{\partial p}{\partial t} +\mathbf{v}\cdot\nabla p\right)-\left(\mathbf{\tau}:\nabla\mathbf{v}\right) $$

This equation is useful in Chemical Engineering as it allows us to obtain the change in energy in terms of the system's temperature, an observable and measurable parameter in various applications.

Here we consider two assumptions: (1) the viscous dissipation term is negligible since it is generally only important in fluids with enormous velocity gradients, and (2) we consider a fluid with constant density, thus eliminating the third term on the right side of the equation $\left(\left(\frac{\partial \ln{\rho}}{\partial\ln{T}}\right)_p=0\right)$, resulting in:

$$ \rho \hat{C}_p \frac{\partial T}{\partial t} + \rho\hat{C}_p\mathbf{v}\cdot\ \nabla T = -\left(\nabla\cdot \mathbf{q}\right) $$

Finally, applying this equation to our one-dimensional problem, using Fourier's law for heat conduction, a plug flow velocity profile with a fixed velocity $v_0$, and considering a natural convection source term as a homogeneous source term ($S$), we have:

$$ \rho \hat{C}_p \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial x^2} - \rho \hat{C}_p v_0 \frac{\partial T}{\partial x} + S $$

This equation is subject to the initial condition:

$$ T(x,t=0) = T_{in} $$

And the boundary conditions:

At $x=0$ (BC1), we set the inlet temperature of nitrogen (Dirichlet):

$$ T(x=0,t) = T_{in} $$

At $x = L$ (BC2), we use an outlet condition by setting a zero temperature gradient:

$$ \frac{d T}{d x}(x=L,t) = 0 $$

We have written the natural convection heat transfer as a homogeneous source term, $S$ (W $\textrm{m}^{-3}$), but it actually occurs on the pipe surface. To find the volumetric source term, we perform a differential balance on a differential cylinder element of thickness $\Delta x$:

$$ S = \frac{\text{heat flux} \times \text{differential transfer area}}{\text{volume of differential element}} = \frac{-h(T-T_\text{aire})(\pi d_i \Delta x)}{\pi \frac{d_i^2}{4} \Delta x} = \frac{-4 h(T-T_\text{aire})}{d_i} $$

Grid

We generate a one-dimensional (1-D) regular grid as shown in Figure 2, creating $N$ points along the pipe representing the centers of $N$ control volumes. The distance between points is $\Delta x = \delta x$, which is also the width of each control volume. The distance between the center of the initial and final control volumes and the pipe edge is $\Delta x/2$ (half the width of a control volume).

In Figure 2, we highlight a control volume of interest, whose center we designate with the letter $P$, representing any internal control volume in our grid. We also designate the centers of the neighboring control volumes with the letters $W$ and $E$ (west and east), and define the boundaries or planes, located at intermediate points between $WP$ and $PE$, identified with the lowercase letters $w$ and $e$, respectively.

Discretization

We discretize the differential equation using the finite volume method. We begin by integrating the differential equation in space over the control volume with central point $P$:

$$ \begin{alignat}{5} &\int_{V_P}\frac{\partial}{\partial t} \rho C_p T d V && && && &&\qquad\text{(Energy accumulation in volume P)}\ & &&+ \int_{V_P} \rho C_p v_0 \frac{\partial T}{\partial x} d V && && &&\qquad\text{(Net energy transport by convection)}\ & && && - \int_{V_P} k \frac{\partial^2 T}{\partial x^2} d V && &&\qquad\text{(Net energy transport by conduction)}\ & && && && = \int_{V_P} S d V &&\qquad\text{(Net energy production in volume P)} \end{alignat} $$

A fundamental step of the finite volume method is to convert the volumetric integrals appearing in our approximation into surface integrals in the control volumes. For this, we use the divergence theorem or Gauss's theorem, which states that the outward flux of a vector field through a closed surface is equal to the volume integral of the divergence over the region within the same surface. The divergence theorem in its general form for multiple dimensions is defined as:

$$ \int_{V_P} \nabla \cdot \mathbf{\Phi} d V = \int_{A_P} \mathbf{\Phi} \cdot d \mathbf{A} $$

where $\mathbf{\Phi}$ is a vector field, and $A_P$ represents the surface of the control volume $V_P$.

For the integration of the time derivative term, we assume the value of $T_P$ in the volume remains uniform over the interval $\Delta t = t_{n+1} - t_n$, resulting in the following discretized term:

$$ \rho C_p V_P \frac{T_P^{n+1} - T_P^n}{\Delta t} $$

In our one-dimensional problem, the spatial domain comprises the control volumes $i-1$ and $i+1$, represented by the letters $E$ and $W$, respectively. These control volumes and the grid form a line. The number of control volumes, or nodes, is $N$, and the coordinate $x_i$ is the position of the central point $P$ of the $i^{\text{th}}$ control volume. Additionally, $T_i^n$ represents the average value of the temperature in the control volume $i$ at time $t_n$.

The energy accumulation term is given by:

$$ \rho \hat{c}_p V_P \frac{T_P^{n+1} - T_P^n}{\Delta t} $$

where $V_P = A_P \Delta x$ is the volume of the control volume.

The convection term is given by:

$$ \rho C_p v_0 \int_{V_P} \frac{\partial T}{\partial x} dV = \rho C_p v_0 A_P \left(T_e - T_w\right) $$

The conduction term is given by:

$$ -k \int_{V_P} \frac{\partial^2 T}{\partial x^2} dV = -k A_P \left(\frac{T_E - T_P}{\Delta x} - \frac{T_P - T_W}{\Delta x}\right) $$

The source term is given by:

$$ \int_{V_P} S dV = S V_P = \frac{-4 h A_P (T_P - T_\text{aire})}{d_i} \Delta x $$

Implementation in Python

Below is the Python implementation of the finite volume method for the given heat transfer problem:

import numpy as np
import matplotlib.pyplot as plt

# Define problem parameters
L = 1.0  # Length of the pipe (m)
d_i = 0.0508  # Internal diameter of the pipe (m)
T_in = 70  # Inlet temperature (K)
T_air = 298.15  # Ambient temperature (K)
h = 25  # Natural convection heat transfer coefficient (W/m^2.K)
rho = 838.645  # Density (kg/m^3)
C_p = 2014.03  # Specific heat at constant pressure (J/kg.K)
k = 0.16  # Thermal conductivity (W/m.K)
v_0 = 0.1  # Fluid velocity (m/s)
dx = 0.01  # Spatial step (m)
dt = 1  # Time step (s)
time_steps = 3000  # Number of time steps

# Define number of grid points
N = int(L/dx)

# Initialize temperature array
T = np.full(N, T_in)

# Define source term
S = lambda T: -4 * h * (T - T_air) / d_i

# Time integration loop
for n in range(time_steps):
    T_new = T.copy()
    for i in range(1, N-1):
        # Convection term
        conv_term = -rho * C_p * v_0 * (T[i+1] - T[i-1]) / (2 * dx)
        # Conduction term
        cond_term = k * (T[i+1] - 2*T[i] + T[i-1]) / dx**2
        # Source term
        source_term = S(T[i])
        # Update temperature
        T_new[i] = T[i] + dt * (conv_term + cond_term + source_term) / (rho * C_p)
    # Apply boundary conditions
    T_new[0] = T_in
    T_new[-1] = T[-2]
    # Update temperature array
    T = T_new.copy()

# Plot results
x = np.linspace(0, L, N)
plt.plot(x, T, label='Temperature')
plt.xlabel('x (m)')
plt.ylabel('Temperature (K)')
plt.title('Temperature Distribution along the Pipe')
plt.legend()
plt.show()