1.3. Advection

In the first chapter, it was shown that for incompressible flow, constant diffusion coefficient (\(D\)) and accounting for mass conservation, the equation of transport can be expressed as:

\[ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = D \nabla^2 \phi +F_{\phi} \]

A pure advection problem can be used to illustrate the finite volume method and to isolate the problems that may arise while treating this term. Consider the idealized problem without diffusion and source terms, a quantity \(\phi(\mathbf{x},t)\) traveling in a fluid with velocity field \(\mathbf{u}(\mathbf{x},t)\):

\[ \frac{\partial \phi(\mathbf{x},t)}{\partial t} + \mathbf{u}(\mathbf{x},t) \cdot \nabla \phi(\mathbf{x},t) = 0 \]

This relationship describes phenomenon that can be thought of as a “pulse” of quantity \(\phi\) traveling through the fluid, for example a single wave (\(\phi\) being the wave amplitude), or a pollutant (\(\phi\) as concentration). In the example presented here, we will consider the pulse being carried (advected) through a fluid with constant velocity \(\mathbf{c}\), where the components in each direction are:

\[ \mathbf{c} = c_u \,\hat{i} + c_v \,\hat{j} + c_w \,\hat{k} \]

As the fluid velocity is constant in time and space, we can re-write the conservation equation as follows:

\[ \frac{\partial \phi(\mathbf{x},t)}{\partial t} + \mathbf{c} \cdot \nabla \phi(\mathbf{x},t) = 0 \]

Application in 1D

The finite volume method will be applied in 1D for simplicity following the steps:

1. Integrate PDE over volume of interest \(V\)

2. Apply Gauss’s Theorem on the bounding surface \(S\)

3. Discretize the volume and integral equations over finite volumes (in space): meshing

4. Apply a numerical scheme to discretize equations in time, then solve

Integral Form

Consider only the \(x\)-direction to represent a 1D problem:

\[ \frac{\partial \phi}{\partial t} + c \, \frac{\partial \phi}{\partial x} = 0 \]

where the subscript \(u\) (denoting the \(u\)-component of \(\mathbf{c}\)) has been dropped. The situation represents a “pulse” traveling along the \(x\)-axis with speed \(c\), and is illustrated in Fig. 1.6 for initial conditions such that a pulse with width \((b-a)\) has amplitude \(\phi_0\) between the interval \((a,b)\) in \(x\).

Fig. 1.6 Propagation of wave in 1D with constant speed \(c\).

Integrate the linear convection equation in \(x\) over an arbitrary volume \(V\) with surface \(S\):

\[ \int_{V}\frac{\partial \phi}{\partial t}dV + \int_{V}c\frac{\partial \phi}{\partial x} \, dV = 0 \]

Applying Gauss’s theorem

Recognizing \(\frac{\partial \phi}{\partial x}\) as the 1D divergence \(\nabla \cdot u\) and applying Gauss’s Theorem with \(\mathbf{n}\) being the normal surface vector:

\[ \int_{V}\frac{\partial \phi}{\partial t}dV + \int_{S}c\, \phi \cdot \mathbf{n}\, dS = 0 \]

Finite Volume Mesh

The \(x\)-domain will be divided into six finite volumes, shown in Fig. 1.7, numbered 1 through 6 from left to right, such that the volume number refers to the \(x\)-coordinate of the geometric center. For example, volume 3 is located with center at \(x_3\). The quantity \(\phi\) will only being calculated at the center of the finite volume, \(\phi(x_i)\).

Fig. 1.7 Schematic of 6 finite volumes aligned with \(\hat{i}\) (\(x\)-direction).

Discretization

The integral equations may now be discretized over each finite volumes with the “unknown” at the center, \(\phi_i\), and with dimensions \(\Delta x\), \(\Delta y\) and \(\Delta z\). .

For the first term: velocity \(c\) is constant over the domain, the integral is applied over the finite volume \(i\) and a forward Euler scheme is applied with time-step \(\Delta t\):

\[ \int_{V}\frac{\partial \phi}{\partial t}dV = \frac{d\phi_i}{dt}\int{dV} = \frac{d\phi_i}{dt} \Delta x\Delta y \Delta z \]

\[ = \frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta t} \Delta x\Delta y \Delta z \]

Note that even though it is a 1D problem, the volume is still represented with the three dimension \(\Delta\) to represent \(\Delta V\)

Fig. 1.8 6 finite volumes showing discretization around \(\phi_{i}\) (\(x\)-direction) and the Western and Eastern boundary .

For the second term: the general surface integral can be represented as the sum of integrals over the active volume surfaces (east and west) as visualized in the figure above:

\[ \int_{S} c \, \phi \cdot \mathbf{n} \, dS = \int_{S_w}c \, \phi_{w} \cdot \mathbf{n}_w \, dS + \int_{S_e}c \, \phi_{e} \cdot \mathbf{n}_e \, dS \]

In the previous equation the west and east faces are explicitly written as subscript in the normal vector \(\mathbf{n}\) to emphasize that it is the normal to those surfaces pointing outwards the volume, so the normal vector of the west surface is \(-i\) (pointing outwards) and the normal vector of the east surface is \(i\), i.e. one integral has a negative sign and the other a positive one:

\[ \int_{S} c \, \phi \cdot \mathbf{n} \, dS = -\int_{S_w}c \, \phi_{w} \, dS + \int_{S_e}c \, \phi_{e} \, dS \]

For a 1D problem there wont be any variations of \(\phi\) along the surface, if there were we can already approximate that value with a numerical integration technique (Left Riemman, midpoint, etc). Recalling that the midpoint evaluation is second order accurate and convenient, we obtain:

\[ \int_{S} c \, \phi \cdot \mathbf{n} \, dS = -c \, \phi_{w} \, \Delta y \Delta z + c \phi_{e} \, \Delta y \Delta z = c\left(\phi_{e}-\phi_{w}\right)\Delta y \Delta z \]

where \(\phi_{e,w}\) are evaluated at the midpoint of the east and west surfaces. As those values are not yet known at the surface (they are known at the volume centroid), an approximation is needed. A linear interpolation can be used, thus only two points per surface value are needed:

\[ \phi_{e} = \frac{\phi_{i+1} + \phi_i}{2} \qquad \phi_{w} = \frac{\phi_{i-1} + \phi_i}{2} \quad \rightarrow \quad \phi_{e} - \phi_{w} = \frac{\phi_{i+1} + \phi_i - \phi_{i-1} - \phi_i}{2} \]

Putting both integral terms together:

\[ \frac{d\phi_i}{dt} \Delta x\Delta y \Delta z + c\left(\frac{\phi_{i+1}-\phi_{i-1}}{2}\right)\Delta y \Delta z = 0 \]

This is an ODE! We can use one of the methods of the previous quarter to solve it. For example, using Explicit Euler and dividing by \(\Delta x \Delta y \Delta z\):

\[ \frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta t} + \frac{c}{2\Delta x}\left(\phi_{i+1}^n-\phi_{i-1}^n\right) = 0 \]

Rearranging, the algebraic equation for the unknown term at the next time step:

\[ \phi_{i}^{n+1} = \phi_{i}^{n} - \frac{c \Delta t}{2\Delta x}\left(\phi_{i+1}^n-\phi_{i-1}^n\right) \]

Turns out that the 1D problem final algebraic expression is the same as for the finite difference method using central differences to approximate the spatial derivative! The differences become evident once non-rectangular volumes are used.

Unfortunately, the use of Explicit Euler in this case is useless as it will make the problem unstable for any time step.

Stability

As shown with ODE’s, the existance of a stability criteria only applies for explicit schemes. Unfortunately, using Explicit Euler in this case is useless as it will make the problem unstable for any time step (to be tested in the workshop). The physical reason why this combination worked fine with the diffusion equation is that diffusivity occurs in all directions (elliptic, parabolic problem depending on unsteadiness), but for advection (hyperbolic equation) information mainly comes from upstream (the direction of information propagation), so it is not logical to have equal weights from contiguous volumes to define the flux of \(\phi\)! It should be more relevant the quantity coming from upstream.

This can be circumvented in two ways. One is by using central differences in time:

\[ \phi_{i}^{n+1} = \phi_{i}^{n-1} - \frac{c \Delta t}{\Delta x}\left(\phi_{i+1}^n-\phi_{i-1}^n\right) \]

Indirectly, there is more relevance from the upstream direction by using \(\phi^{n-1}\) rather than \(\phi^{n}\).

The other and more easily appreciated is by applying an upwind scheme. This, if \(c\) is positive:

\[ \phi_{e} = \phi_{i} \qquad \phi_{w} = \phi_{i-1}\qquad \]

if \(c\) is negative :

\[ \phi_{e} = \phi_{i+1} \qquad \phi_{w} = \phi_{i} \quad \]

.

The fist option scheme is named leapfrog because the value \(\phi^n_i\) is skipped. The resulting stability criteria, here shown without proof, is the Courant-Friedrich-Lewy criteria (CFL):

\[\frac{|u|\Delta t}{\Delta x}<1\]

The physical meaning of this criteria is that time step should be less than the time that takes the information to cross to the next volume. Thus, if a very fine spatial resolution is used, the time-set must be very small to maintain the CFL criteria.

Solution as System of Equations

We will now assemble a system of equations for the entire problem and formulate it as a matrix. This is instructive not only for understanding the solution and how the individual finite volumes relate to each other geometrically, as well as to the boundary conditions, but it also facilitates a comparison with other numerical methods (e.g., finite difference and finite element methods).

Observe that the general forward Euler solution above for a non-boundary volume can be expressed in vector/matrix form as:

\[\begin{split} \phi_{i}^{n+1} = \phi_{i}^{n-1} - \frac{c \Delta t}{\Delta x} \begin{bmatrix} -\frac{1}{2} & 0 & \frac{1}{2} \end{bmatrix} \begin{bmatrix} \phi_{i-1} \\ \phi_{i} \\ \phi_{i+1} \end{bmatrix}^{n} \end{split}\]

Recognizing that the vector of coefficients \(\begin{bmatrix} -\frac{1}{2} & 0 & \frac{1}{2} \end{bmatrix}\) can be written \(\begin{bmatrix} a_{i,i-1} & a_{i,i} & a_{i,i+1} \end{bmatrix}\) and is the \(i^{\textrm{th}}\) equation in the system allows generalizing to the global matrix formulation:

\[ \begin{bmatrix} \: \phi \: \end{bmatrix}^{n+1} = \begin{bmatrix} \: \phi \: \end{bmatrix}^{n-1} + \begin{bmatrix} \: A \: \end{bmatrix}^{n} \begin{bmatrix} \: \phi \: \end{bmatrix}^{n} \]

Where the coefficients of \(A\), \(a_{ij}\), are defined by the finite volume discretization scheme:

\[\begin{split} \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \\ \phi_{4} \\ \phi_{5} \\ \phi_{6} \end{bmatrix}^{n+1} = \quad \begin{bmatrix} \:I\: \end{bmatrix} \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \\ \phi_{4} \\ \phi_{5} \\ \phi_{6} \end{bmatrix}^{n-1} + \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} & a_{15} & a_{16} \\ a_{21} & a_{22} & a_{23} & a_{24} & a_{25} & a_{26} \\ a_{31} & a_{32} & a_{33} & a_{34} & a_{35} & a_{36} \\ a_{41} & a_{42} & a_{43} & a_{44} & a_{45} & a_{46} \\ a_{51} & a_{52} & a_{53} & a_{54} & a_{55} & a_{56} \\ a_{61} & a_{62} & a_{63} & a_{64} & a_{65} & a_{66} \\ \end{bmatrix}^{n} \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \\ \phi_{4} \\ \phi_{5} \\ \phi_{6} \end{bmatrix}^{n} \end{split}\]

Applying and extending the pattern found above:

\[\begin{split} \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \\ \phi_{4} \\ \phi_{5} \\ \phi_{6} \end{bmatrix}^{n+1} = \quad \begin{bmatrix} \:I\: \end{bmatrix} \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \\ \phi_{4} \\ \phi_{5} \\ \phi_{6} \end{bmatrix}^{n-1} - \frac{c\Delta t}{\Delta x} \begin{bmatrix} 0 & \frac{1}{2} & 0 & 0 & 0 & 0 \\ -\frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 & 0 \\ 0 & -\frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & -\frac{1}{2} & 0 & \frac{1}{2} & 0 \\ 0 & 0 & 0 & -\frac{1}{2} & 0 & \frac{1}{2} \\ 0 & 0 & 0 & 0 & -\frac{1}{2} & 0 \\ \end{bmatrix} \begin{bmatrix} \phi_{1} \\ \phi_{2} \\ \phi_{3} \\ \phi_{4} \\ \phi_{5} \\ \phi_{6} \end{bmatrix}^{n} \end{split}\]

The initialization of the problem is delicate as only after the second time step the algebraic equation is fully satisfied without additional requirements.

The coefficient matrix \(A\) only has non-zero elements along the upper and lower diagonal due to the regular geometry incorporated in the 1D discretization. In the derivation of the system of equations, the boundary conditions of end volumes were not presented. Actually, the advection equation has only one boundary condition due to the first spatial derivative! It is defined at the boundary from where \(\phi\) is being advected. The formulation of the coefficients in the matrix implies that the values of the \(0^{\textrm{th}}\) and \(7\textrm{th}\) volume (\(\phi_{0}\) and \(\phi_7\)) were taken as 0, the latter is an artifice to vanish the quantity \(\phi\). Rigorously speaking, it should not be a boundary condition.

Solution Techniques

This system can be solved without the need for building the matrix as each value of the next time step can be computed from known values of previous steps. So, it can be solved with a for loop as well!

Note that the matrix in the above example contains a lot of zeros: 10 out of 36 elements are non-zero (28%). This is called a sparse matrix, and the proportion of non-zero entries grows with the size of the problem. For example, \(N=1000\) finite volumes requires \(2N-2\) non-zero entries (1998), which is only 0.1% of the matrix! Since matrix calculations are computationally expensive, we would never carry out the matrix calculations as represented here. Instead, various vectorization approaches are implemented in all numerical analysis software packages. Typically, a mapping of volumes and coordinates are stored in a vector format that retains the relationship between volumes and the interpolation points (e.g., neighbor volumes). Then calculations per volume are carried out using this vectorized information. This can be easily implemented in any programming language, for example, with a for loop over each volume. However, note that as geometry of the domain of interest and the discretization scheme changes (e.g., non-square volumes!) more advanced vectorization techniques are required.

Diffusion

Lets look at the diffusion equation.

\[ \frac{\partial \phi}{\partial t} - D \nabla^2 \phi = F_{\phi} \]

Previously, we analyzed the advection term and derived it for the finite volume method. Now, let’s examine the diffusion equation, which is a second-order differential equation. For simplicity, assume that \(F_{\phi}=0\). We will follow a similar approach as before::

Integrate the diffusion equation over an arbitrary volume:

\[ \int_{V} \frac{\partial \phi}{\partial t}dV = \int_{V}D \frac{\partial^2 \phi}{\partial x^2} \, dV \]

Applying Gauss’s theorem

Recognizing \(\frac{\partial \phi}{\partial x}\) as the 1D divergence \(\nabla \cdot u\) and applying Gauss’s Theorem with \(\mathbf{n}\) being the normal surface vector:

\[ \int_{V} \frac{\partial \phi}{\partial t}dV = \int_{S}D \frac{\partial \phi}{\partial x} \cdot \mathbf{n}\, dS \]

We first focus on the diffusion term \( \int_{S}D \frac{\partial \phi}{\partial x} \cdot \mathbf{n}\, dS\) :

\[ \int_{S} D \frac{\partial \phi}{\partial x} \cdot \mathbf{n}\, dS = \int_{S_w} D \, \frac{\partial \phi_w}{\partial x} \cdot \mathbf{n}_w \, dS + \int_{S_e} D \, \frac{\partial \phi_e}{\partial x} \cdot \mathbf{n}_e \, dS \]

Here, the west and east faces are explicitly indicated with subscripts \(w\) and \(e\) in the normal vector \(\mathbf{n}\) to show that it points outward from the volume. For example, the normal vector on the west surface is \( -\hat{i} \) (pointing outwards), while on the east surface it is \( \hat{i} \), leading to a negative sign for one integral and a positive sign for the other:

\[ \int_{S} D \frac{\partial \phi}{\partial x} \cdot \mathbf{n}\, dS = -D \int_{S_w} \, \frac{\partial \phi_w}{\partial x} \, dS + D \int_{S_e} \, \frac{\partial \phi_e}{\partial x} \, dS \]

In Q1, we learned that we can approximate the first-order derivative forward and backward over the grid as::

\[ \frac{\partial \phi_w}{\partial x} \approx \frac{\phi_i - \phi_{i-1}}{\Delta x} \]

\[ \frac{\partial \phi_e}{\partial x} \approx \frac{\phi_{i+1} - \phi_{i}}{\Delta x} \]

Substituting these approximations, we get:

\[ \int_{S} D \frac{\partial \phi}{\partial x} \cdot \mathbf{n}\, dS = - D \frac{\phi_i - \phi_{i-1}}{\Delta x} \Delta y \Delta z + D \frac{\phi_{i+1} - \phi_{i}}{\Delta x} \Delta y \Delta z \]

Substituting this expression back together with the first term, we get:

\[ \frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta t} \Delta x\Delta y \Delta z = - D \frac{\phi_i^{n} - \phi_{i-1}^{n}}{\Delta x} \Delta y \Delta z + D \frac{\phi_{i+1}^{n} - \phi_{i}^{n}}{\Delta x} \Delta y \Delta z \]

This simplifies to:

\[ \frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta t} \Delta x\Delta y \Delta z = D \frac{ \phi_{i-1}^{n}- 2\phi_i^{n} + \phi_{i+1}^{n} }{\Delta x} \Delta y \Delta z \]

Finally , dividing by \(\Delta x \Delta y \Delta z\):

\[ \frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\Delta t} = D \frac{ \phi_{i-1}^{n}- 2\phi_i^{n} + \phi_{i+1}^{n} }{\Delta x^2} \]

Stability

This equation can be rewritten to isolate the amplification factor:

\[ \phi_i^{n+1} = (1-2*\Delta t D/\Delta x^2)\phi_i^n + \Delta t D \frac{\phi_{i-1}^n + \phi_{i+1}^n}{\Delta x^2} \]

where the amplification factor is multiplying \(\phi_i^n\) and should be lower than 1 to ensure stability:

\[ -1 < 1-2D \frac{\Delta t}{\Delta x^2} < 1 \]

thus

\[ D \frac{\Delta t}{\Delta x^2} <= 1/2 \]

Attribution

This chapter is written by Robert Lanzafame and Jaime Arriaga Garcia. {ref}Find out more here <finite_volume_method_credit>`.