Prerequisites:

In this section, we will introduce the finite volume method. To fully understand this section, it is desirable to have a basic knowledge of numerical methods and the ability to apply the finite difference method to solve ordinary and partial differential equations. If you lack these prerequisites & speak spanish, we invite you to take the free MOOC Computational Momentum, Heat, and Mass Transfer.

Introduction to the Finite Volume Method

The finite volume method is a numerical method for solving partial differential equations that describe transport phenomena. The starting point of this method is the integral form of the conservation equations. This differentiates the method from the finite difference or finite element methods, where the differential form of the conservation equations is solved. The method involves dividing the spatial domain into a set of finite computational control volumes, over which the conservation equation of the physical quantity being transported is applied. The computational variables are stored at the centre of each volume, and the values of the variables at the faces of the volume are obtained by interpolation. Then, a system of algebraic equations is obtained and solved to get the approximate solution in each control volume. The finite volume method is very flexible and can be adapted to complex geometries and varied boundary conditions. Additionally, the method ensures the global and local conservation of the physical quantity being transported. The finite volume method has been successfully applied to problems of heat transfer by conduction, convection, and radiation, both in homogeneous and heterogeneous media, and in different thermal regimes.

The Integral Form of the Conservation Equation for a Scalar Quantity

Before introducing the method, it is convenient to recall the integral form of the conservation equation for a scalar quantity, such as energy or mass, in a control volume $V$:

$$\frac{\partial}{\partial t} \int_V \rho \phi dV + \int_S \rho \phi \mathbf{v} \cdot \mathbf{n} dS = \sum_{k} f_{\phi} $$

Where the first term on the left represents the accumulation of the quantity $\phi$ in the control volume and the second term represents the net transport of $\phi$ by convection. The term on the right represents the sum of all transport mechanisms other than convection, plus all the generation and consumption terms of the quantity $\phi$. Unlike the finite difference or finite element methods, the finite volume method is a conservative method. This means that when performing a global property balance over the entire computational domain, the quantity is conserved exactly. This property is very useful in Process and Chemical Engineering, particularly in compressible flow problems where density, temperature, and pressure profiles are coupled.

Conservation of Heat and Mass

In the case of heat transfer or mass transfer, there is always a molecular transport term as a summand on the right side of the equation. For mass transfer, $\phi = c$ and the flux of species $A$ is governed by Fick's Law for mass transfer,

$$f_{c} = - \mathfrak{D}_{AB} \nabla c_A$$

where $\mathfrak{D}_{AB}$ is the diffusivity of species A in medium B in a binary mixture, and $c_A$ is the concentration of species A. Recall that the flux of a property is the amount of property flowing through a specific point in space, per unit area per unit time.

For energy conservation, $\phi = \hat{e}$ and Fourier's Law governs heat transfer by thermal conduction

$$ f_{e} = - k \nabla T $$

The Finite Volume Method

1. The first step is to define control volumes (CV) on a suitable grid and assign the computational node at its centre. Alternatively, the locations of the nodes are defined and control volumes are constructed around them. This results in the faces of the control volumes being between two nodes. In either case, it is important that the CVs do not overlap to ensure global conservation. In this course, the first method will be favoured. As an introduction, it will be assumed that the velocity vector $\mathbf{v}$ and thermophysical properties are known. To obtain a system of algebraic equations from the integral form of the conservation equations, the integrals must be approximated using various quadrature formulas.

Quadrature Formulas for Approximating Surface Integrals

The finite volume method is applicable to 1-D, 2-D, and 3-D systems. The surface of the control volume consists of two (1-D), four (2-D), or six (3-D) flat faces. The net flux through a control volume is the sum of the integrals over all the faces:

$$ \int_S f dS = \sum_k \int_{S_k} f dS $$

where $f$ is the component of the convective flux ($\rho \phi \mathbf{v} \cdot n$) or diffusive flux ($\Gamma \nabla \phi \cdot \mathbf{n}$) in the direction normal to the control volume face. For illustration, the method will be applied to a 2-D domain. Each node has four neighbours: to the right or east (E), to the left or west (W), forward or north (N), and backward or south (S), as shown in Figure 1.

Figure 1: Illustration of a control volume (CV) in a 2-D Cartesian grid. Reprinted from Ferziger et al. Computational methods for fluid dynamics. Springer, 2019.

Note that strictly speaking, for each face, the value of $f$ at each point on the face is required. This information is not directly available, as the grid only stores the values of the variable $\phi$ at the face centres. Therefore, a two-level approximation is necessary:

• The integral is approximated in terms of the values of the variables at one or more points on each face.
• The values of $f$ over the entire face are approximated to the nodal or CV centre values.

Midpoint Rule

The simplest approximation is the midpoint rule, where the integral is approximated by the product of the integrand evaluated at the centre of the face by the area of the face. Without loss of generality, the discretization formulas will be exemplified for the right face of the node, $e$. These formulas are valid for all faces of the CV. The flux leaving through the east face, called $e$, by the mechanism $f$ is

$$F_e = \int_{S_e} f dS = \bar{f_e} S_e \approx f_e S_e $$

Where $S_e$ is the area of the east face of the CV and $f_e$ is the value of $f$ at the centre of the face, which in turn is an approximation of the average of $f$ over face $e$. This integral approximation is second-order accurate, provided that $f_e$ is known. Since this is not the case, $f_e$ must be approximated, and to maintain the second-order accuracy of the integral, the approximation of $f_e$ must also be second-order.

Trapezoidal Rule

The value of $f$ at the north-east (ne) and south-east (se) corners of the face of interest is required. This method is second-order accurate.

$$ F_e = \int_{S_e} f dS = \frac{S_e}{2} (f_{ne} + f_{se})$$

Simpson's Rule

The value of $f$ at three points is required: the north-east (ne) and south-east (se) corners and the centre of the face of interest. This method is fourth-order accurate. To maintain the same order of accuracy for the system solution, the values of $f$ must be approximated with an approximation of at least the same order of accuracy as Simpson's rule. This can be achieved by fitting cubic polynomials.

$$ F_e = \int_{S_e} f dS = \frac{S_e}{6} (f_{ne} + 4f_e + f_{se})$$

References:

Ferziger, Joel H., Milovan Perić, and Robert L. Street. Computational methods for fluid dynamics. Springer, 2019. Chapters 1 and 4.