31.3 Integral form of the Conservative Notation

31.3.1 Derivation

As stated, the notation of Eq. 31.1 is required for transferring the equation into conservative form. Conservative means that the quantity of the change of the dependent variable within the control volume changes in dependence of the fluxes of the dependent variable across the boundary. Outflowing quantity will result in a reduction of the dependent variable within the control volume and v.v.. As stated, we have derived many differential equations using such conservation approaches on infinitesimally small control volumes, within which the quantity of the dependent variable was balanced.

Eq. 31.1 is in differential notation, i.e., the changes of the dependent variable within the control volume over time are balanced by the sources, the sinks, the diffusion inside the control volume, as well as the advection into and out of the control volume. This equation is valid for every point in the computational domain. If we work on the differential form, we would obtain an FDM problem that can be solved if we take into account control volumes of infinitesimally small lateral size. However, as we have seen, FDM is computationally very expensive because it requires a large number of control volumes with very small lateral dimension in order to yield exact results. It would be desirable to “inflate” the control volumes to larger volumes. This is exactly what we do in FVM. The reason why FVM allows working with larger control volumes is due to the fact that it does not approximate the changes of the dependent variable along the control volume, but performs a conservation of the dependent variable for each control volume. This allows calculating an average value for the dependent variable within the cell. The control volume is then replaced by a single value. Obviously, the larger the control volumes, the more “averaged” the solution will become, i.e., finer features of the solution such as spikes, e.g., in concentration will not be visible anymore because they are blurred by the averaging process. The fact that FVM averages gives rise to one very important property: FVM is surprisingly stable against discontinuities and large changes in the dependent variables. These tend to be problematic for FDM because at discontinuities, a function usually has no derivative and thus, FDM cannot approximate values at these locations.

31.3.2 Integration

In order to obtain the conservative form of Eq. 31.1 we integrate this equation over the control volume resulting in

${\int }_V\frac{\partial g}{\partial t}\left(\overrightarrow{x},t\right)dV+{\int }_V\overrightarrow{\nabla }\cdot \overrightarrow{F}\left(\overrightarrow{x},t\right)dV={\int }_{V}f\left(\overrightarrow{x},t\right)dV$

Obviously, this notation is not very useful as of yet. We note two things. First, we can exchange the differential and the integral of the first term on the left-hand side. We can do this because our control volume is fixed in space and not dependent on time, i.e., it is not a function of time. This allows us to convert the partial into a regular differential

${\int }_V\frac{\partial g}{\partial t}dV=\frac{\partial }{\partial t}{\int }_{V}gdV=\frac{d}{dt}{\int }_{V}gdV$

Second, we can use Gauss’s theorem (see section 7.2.1) to convert the volume integrals into surface integrals. Doing this, we can rewrite Eq. 31.4 to

$\frac{d}{dt}{\int }_{V}gdV+{\oint }_{d\partial V}\overrightarrow{F}\cdot \overrightarrow{n}d\partial V={\int }_{V}fdV$

where we have converted the volume integrals into surface integrals. This is a major simplification: Instead of looking at the changes of the dependent quantity in the whole volume V, we simply have to consider the changes on the surface of the control volume ∂V. We can make this control volume as large as required. For finer resolution, V should be small, whereas for rougher solutions, the volume can be chosen larger.

31.3.3 Right-Hand Side Terms

If you look at the right-hand side of Eq. 31.5, you will see that we still have an integral and there seems to be no clever way to make this integral disappear, because we know nothing of the function f. In the general case, this is true. However, we often find that the function f consists mainly of three types of terms.

Constants: Sources and Sinks. Very often, we will have constants on the right-hand side that are obviously easy to integrate. In many cases, we will even have volume-specific constants with the physical unit of “quantity of the dependent variable per volume.” As an example, consider the right-hand side of the Navier–Stokes equation (see Eq. 11.40) where we find the volume forces $\overrightarrow{k}$ , which are volume-specific constants. These are simply integrated by multiplication with the volume.

Constants are physically interpreted as sources (in the case of positive values) or sinks (in the case of negative values) of the dependent variable’s quantity. As an example, a source or sink of concentration could be a chemical reaction that creates or consumes substance. A source for momentum could be forces, e.g., the volume forces in the Navier–Stokes equation.

Nabla Operators: Divergence.f will often contain divergence terms, i.e., terms which have a prepended nabla operator of the form $\overrightarrow{\nabla }$ · p. These can obviously be integrated using the Gauss’s theorem as well. As an example, consider the right-hand side of the Navier–Stokes equation where we actually have two of such terms. One is the pressure gradient $\overrightarrow{\nabla }$ p, which can be integrated using the Gauss’s theorem to

${\int }_V\overrightarrow{\nabla }\cdot pdV={\oint }_{\partial V}p.\overrightarrow{n}d\partial V$

where we have again converted the volume to a surface integral.

Laplace Operators: Diffusion. A second example, which may not be visible intuitively, is the viscous term $\eta \Delta \overrightarrow{\nu }$ of the right-hand side. Remember that the Laplace operator is nothing other than the squared nabla operator therefore

$\eta \Delta \overrightarrow{\nu }=\eta \overrightarrow{\nabla }\cdot \left(\overrightarrow{\nabla }\overrightarrow{\nu }\right)$

in which case we can use the Gauss’s theorem to find the integral

${\int }_V\overrightarrow{\nabla }\cdot \left(\eta \overrightarrow{\nabla }\overrightarrow{\nu }\right)dV={\oint }_{\partial V}\left(\eta \overrightarrow{\nabla }\overrightarrow{\nu }\right)\cdot \overrightarrow{n}d\partial V$

As you can see, even the Laplace operator can be integrated accordingly. The physical interpretation of these terms is diffusion of the independent variable’s quantity, which can obviously distribute, as well as reduce or increase the overall quantity in the control volume. As an example, a fluid is able to transport momentum via diffusion. This is what the viscous term of the right-hand side of the Navier–Stokes equation expresses. Another example is the diffusion equation (see Eq. 8.4) that has a Laplace term on the right-hand side. Obviously, the concentration of a substance c in the control volume cannot only change by advection (convection) or by sources and sinks. It can also be transported out of the control volume by means of diffusion. This is what this term accounts for.

31.3.4 Interpretation

The physical interpretation of Eq. 31.5 is simple. The total change of the dependent variable within the control volume is balanced by the advection of quantity across the boundaries of the control volume, the diffusion of the quantity across the boundaries, and the action of the sources and sinks. However, as we have seen, the right-hand side terms usually give rise either to a constant, or to a function that is again evaluated on the boundary only. This means that the action of the source, sink, and diffusion terms can be evaluated also by looking only at their contribution to advection into and out of the control volume.

As stated, the size of the control volume can be arbitrary. It can be as large as the atmosphere of earth. Obviously, such a control volume is of no use if the dynamics of earth’s atmosphere should be studied, as it would virtually be a very coarse one-voxel-solution.

Another important point is that the order of the equation has been reduced by one from the differential form (see Eq. 31.1) to the conservative or integral form (see Eq. 31.5). You can also see this if you consider the example of the viscous terms of the Navier–Stokes equation we just discussed. In the conservative form, the equation contains only first-order derivatives. As already mentioned, functions with discontinuities or very steep gradients tend to be problematic for FDM. The underlying reason is that a function with a discontinuity does not have a derivative at the location of the discontinuity. Thus, a finite difference solver would not be able to approximate a derivative at this location. In FVM, the reduction in order will avoid problems with nonexisting first-order partial differentials because they do not occur on the left-hand side. For all cases where the right-hand side of the equation in differential form contains a Laplace operator (as does the Navier–Stokes equation), there will be a first-order partial differential in the final equation. However, the value of the differential is only required on the boundary of the control volume. If the discontinuity happens not to lie on the boundary of the domain, it is of no importance. Obviously, the layout of the control volumes can be adapted to work around such critical points.

31.3.5 Volume-Averaged Conservative Form

In the following step, we perform a discretization of the computational domain resulting in a total number of N cells. Each cell has a given (known) volume given by δV^(j). We will replace the value of the dependent variable within the cell j with the volume-averaged value of the dependent variable ${\overline{g}}^{\left(j\right)}$ given by

${\overline{g}}^{\left(j\right)}=\frac{1}{\delta V\left(j\right)}{\oint }_{\delta V\left(j\right)}gdV$

in which case Eq. 31.5 can be written as

$\delta V^{\left(j\right)}\frac{d{\overline{g}}^{\left(j\right)}}{dt}+{\oint }_{\partial \delta V^{\left(j\right)}}\overrightarrow{F}\cdot \overrightarrow{n}d\partial V^{\left(j\right)}={\int }_{\partial \delta V^{\left(j\right)}}fdV$

Please note that Eq. 31.5 and Eq. 31.6 are exact. There was no simplification made during the derivation.

31.5.1 Determining the Constants c_i

Similar to the Galerkin method, we require a set of N functions that allow us to determine c_i. The Galerkin method used the constraint that the integral of the residual across the domain must be zero. Here, we use a similar approach. We state that the integral of the approximation function within each cell must equal the average value $\overline{g}$ within the cell. If this is not the case, the approximation function is inexact, as it does not yield the same average value. We therefore find

${\int }_{\delta V^{\left(j\right)}}{\tilde{g}}^{\left(j\right)}dV\stackrel{!}{=}\delta V^{\left(j\right)}{\overline{g}}^{\left(j\right)}$

Obviously this must be true for each cell. So if our polynomial has a total of N terms, we have a total of N unknown coefficients c_i, and therefore we require a total of N cells. For each cell j in the computational domain we can thus state

$\begin{array}{ll}{\int }_{\delta V^{\left(j\right)}}{\tilde{g}}^{\left(j\right)}dV\hfill & \stackrel{!}{=}\delta V^{\left(j\right)}{\overline{g}}^{\left(j\right)}\hfill \\ {\displaystyle \sum _{i=1}^N{c}_i{\int }_{\delta V^{\left(j\right)}}{\psi }_i\left(\overrightarrow{x}\right)dV}\hfill & \stackrel{!}{=}\delta V^{\left(j\right)}{\overline{g}}^{\left(j\right)}\hfill \end{array}$

which gives us a total of N independent equations to determine the coefficients c_i.

31.5.2 Usage of $\tilde{\mathit{g}}$

One may ask the question why we go through all the work of determining the coefficients c_i, i.e., why do we need the approximation function $\tilde{g}$ at all? The reason for this is that having the average values ${\overline{g}}^{\left(j\right)}$ for each of the cells j will give us the distribution of our dependent variable “inside” of the cell. As stated, we expect this value to be constant over the whole volume of the cell. We usually refer to these values as being cell-centered, which means that we know the value of the dependent variable at the centers of each of the cells (see Fig. 31.2). However,

Eq. 31.6 requires us to evaluate the function on the boundaries between two cells. Obviously, if we assume g to be piecewise-linear, the value at the boundary will be discontinuous depending on the direction from which we approach the discontinuity. It will be either the value of the first piecewise-linear function (if approaching from the left) or the value of the second piecewise-function (if approaching from the right). Usually, this problem can be resolved by stating that the value at the boundary is the average of the two values. In fact, we will use this simplification in an example in just a moment.

However, if the right-hand side function f in Eq. 31.6 contains a diffusion term, we end up with a first derivative in Eq. 31.6 that we need to evaluate at the discontinuity at the boundary of the cells. However, we cannot determine the derivative of the function at the boundary.

These problems are circumvented if we approximate g by the function $\tilde{\mathit{g}}$ , which we can design in a way that it has the values ${\overline{g}}^{\left(j\right)}$ at each cells’ center and is continuous on the boundaries. Obviously, the type of function we expect only depends on the basis function ψ_i. In the following we will discuss the most commonly encountered basis functions for the one-dimensional case. In general, the approximation function used is

${\tilde{g}}^{\left(j\right)}\approx {\displaystyle \sum _{i=0}^N{c}_i}{\left(\frac{x-x^{\left(j\right)}}{\partial x^{\left(j\right)}}\right)}^i$

which uses the simple polynomials $\left(\frac{x-x^{\left(j\right)}}{\delta x^{\left(j\right)}}\right)$ as basis functions ψ_i. Please note that the polynomials are shifted to x − x^(j), i.e., into the center of the cell. Additionally, the shifted term is normalized to the width of the interval. We often use the notation

${\tilde{g}}^{\left(j\right)}\approx {\displaystyle \sum _{i=0}^N{c}_i}{\left({\xi }^{\left(j\right)}\right)}^i$

With

$\begin{array}{l}{\xi }^{\left(j\right)}=\frac{x-x^{\left(j\right)}}{\delta x^{\left(j\right)}}\\ d{\xi }^{\left(j\right)}=\frac{\partial x}{\delta x^{\left(j\right)}}\to dx=\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\end{array}$

where $-\frac{1}{2}\le \xi \le \frac{1}{2}$ . From Eq. 31.9 we can derive the following equations for calculating c_i

${\displaystyle \sum _{i=0}^N{c}_i}{\int }_{-\frac{1}{2}}^{\frac{1}{2}}{\left({\xi }^{\left(j\right)}\right)}^i\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\stackrel{!}{=}\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}$

which is a system of N linearly independent equations (one for each cell j).

31.5.3 Piecewise-Constant Profile

The piecewise-constant profile (see Fig. 31.2a) is the simplest profile. It assumes that the value ${\tilde{g}}^{\left(j\right)}$ is constant within a cell. Therefore, N = 0 in Eq. 31.13 in which case we find

$\begin{array}{ll}{c}_0{\int }_{-\frac{1}{2}}^{\frac{1}{2}}{\left({\xi }^{\left(j\right)}\right)}^0\delta x^{\left(j\right)}d{\xi }_j\hfill & \stackrel{!}{=}\delta x^{\left(j\right)}{g}^{-\left(j\right)}\hfill \\ c_0{x}^{\left(j\right)}\hfill & \stackrel{!}{=}\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}\hfill \\ c_0\hfill & \stackrel{!}{=}{\overline{g}}^{\left(j\right)}\hfill \end{array}$

and therefore

$\tilde{g}\left(x\right)={\overline{g}}^{\left(j\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}{x}^{\left(j\right)}\le x<x^{\left(j+1\right)}$

As we have already discussed, the value at the boundary of the cells is also of importance. For piecewiseconstant profiles, the value at the interface between cell j and cell j + 1 is defined as the average of the two values

${\tilde{g}}^{\left(j\to j+1\right)}=\frac{{\tilde{g}}^{\left(j\right)}+{\tilde{g}}^{\left(j+1\right)}}{2}$

As stated, this profile will be unsuitable if we have a diffusion term in our equation which requires the evaluation of the derivative, i.e., the flux term at the boundary.

31.5.4 Piecewise-Linear Profile

The second commonly used profile is the piecewise-linear profile (see Fig. 31.2b). It assumes the value of ${\tilde{g}}^{\left(j\right)}$ to change linearly between the interpolation points that are located at the cells’ centers. We obtain this approximation by using N = 1 in Eq. 31.10 from which we find

${\tilde{g}}^{\left(j\right)}=c_0+c_1{\xi }^{\left(j\right)}\stackrel{!}{=}\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}$

Obviously, we require N = 2 cells in order to solve for the unknowns c₀ and c₁.

Stencils. This is a very common mechanism used in the reconstruction of polynomials from fixed (given values). Fitting algorithms have to solve similar tasks. Given a finite number of points, they need to fit a given curve, such that a minimal error is obtained. For polynomials, methods based on stencils are very common for this task. In addition to the cell j, we can choose the cells j − 1 or j + 1 to contribute. The choice is arbitrary. Usually, the selection pattern of neighboring cells is referred to as a reconstruction stencil and there are certain patterns which lend themselves better than others. The stencil is selected due its ability to achieve stable solutions which do not oscillate. Therefore, this type of reconstruction is commonly referred to as essentially non-oscillatory (ENO) reconstruction. In this type of reconstruction, only one stencil is selected. The second type of reconstruction schemes is weighted essentially nonoscillatory (WENO) reconstruction. Here, all potential stencils are taken into account and weighted in order to obtain the most suitable solution.

A Simple ENO Reconstruction. In the following we will use a simple ENO scheme which is sufficient for linear reconstruction. WENO reconstruction becomes more important for higher number of terms, i.e., for higher values of N. Here, we chose cell j + 1 as second cell for our stencil. We then state that in order to be a good approximation, the integral of $\tilde{\mathit{g}}$ must

• amount to the average value ${\overline{g}}^{\left(j\right)}$ in cell j, i.e., for $-\frac{1}{2}\le \xi \le \frac{1}{2}$ and

• amount to the average value ${\overline{g}}^{\left(j+1\right)}$ in cell j+ 1, i.e., for $\frac{1}{2}\le \xi \le \frac{3}{2}$ and

From this we derive the following equations

$\begin{array}{l}{\int }_{-\frac{1}{2}}^{\frac{1}{2}}\left(c_0+c_1{\xi }^{\left(j\right)}\delta x^{\left(j\right)}\right)d{\xi }^{\left(j\right)}=\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}\\ {\int }_{\frac{1}{2}}^{\frac{3}{2}}\left(c_0+c_1{\xi }^{\left(j\right)}\delta x^{\left(j\right)}\right)d{\xi }^{\left(j\right)}=\delta x^{\left(j+1\right)}{\overline{g}}^{\left(j+1\right)}\end{array}$

or

$\begin{array}{l}{\int }_{-\frac{1}{2}}^{\frac{1}{2}}\left(c_0+c_1{\xi }^{\left(j\right)}\right)d{\xi }^{\left(j\right)}={\overline{g}}^{\left(j\right)}\\ {\int }_{-\frac{1}{2}}^{\frac{3}{2}}\left(c_0+c_1{\xi }^{\left(j\right)}\right)d{\xi }^{\left(j\right)}=\frac{\delta x^{\left(j+1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j+1\right)}\end{array}$

After performing the integrations we obtain

$\begin{array}{ll}{c}_0\hfill & ={\overline{g}}^{\left(J\right)}\hfill \\ c_0+c_1\hfill & =\frac{\delta x^{\left(j+1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j+1\right)}\hfill \end{array}$

from which we derive

$\begin{array}{ll}{c}_0\hfill & ={\overline{g}}^{\left(j\right)}\hfill \\ c_1\hfill & =\frac{\delta x^{\left(j+1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j+1\right)}-{\overline{g}}^{\left(j\right)}\hfill \\ \tilde{g}\left({\xi }^{\left(j\right)}\right)\hfill & ={\overline{g}}^{\left(j\right)}\left(1-{\xi }^{\left(j\right)}\right)+\frac{\delta x^{\left(j+1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j+1\right)}{\xi }^{\left(j\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}-\frac{1}{2}\le {\xi }^{\left(j\right)}\frac{3}{2}\hfill \end{array}$

Obviously, if the grid is regular and x^(j) = x^(j+ 1) we obtain

$\tilde{g}\left({\xi }^{\left(j\right)}\right)={\overline{g}}^{\left(j\right)}\left(1-{\xi }^{\left(j\right)}\right)+{\xi }^{\left(j\right)}{\overline{g}}^{\left(j+1\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}-\frac{1}{2}\le {\xi }^{\left(j\right)}\frac{3}{2}$

This function can be evaluated at the cell’s edges by simply setting $\xi =-\frac{1}{2}$ for the left edge and $\xi =\frac{1}{2}$ for the right edge. Compared to the piecewise-constant function, we can form the derivative of Eq. 31.15 which is

$\frac{d{\tilde{g}}^{\left(j\right)}}{d{\xi }^{\left(j\right)}}={\overline{g}}^{\left(j+1\right)}-{\overline{g}}^{\left(j\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}-\frac{1}{2}\le {\xi }^{\left(j\right)}\le \frac{3}{2}$

which is simply the difference between the two cells’ average values. As you can see, this function can be differentiated without problems and the fluxes can be calculated at the cells’ boundaries.

Conversion to the Differential$\frac{\mathit{d}}{\mathit{d}\mathit{x}}\text{.}$ If you consider Eq. 31.16 closely, you will see that the derivative has the same unit as the average value, which obviously is incorrect. It should have the unit “dependent value unit per length.” The reason for this is that Eq. 31.16 is the derivative with respect to ξ^(j). If we require the derivative with respect to x, we can see from Eq. 31.12 that we have to divide Eq. 31.16 by dξ^(j). This also results in the correct unit.

31.5.5 Piecewise-Parabolic Profile

The third commonly used approximation is piecewise-parabolic with N = 2 in Eq. 31.10 from which we find

${\tilde{g}}^{\left(j\right)}=c_0+c_1{\xi }^{\left(j\right)}+c_2{\left({\xi }^{\left(j\right)}\right)}^2\stackrel{!}{=}\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}$

where we again have the problem that we need to apply a stencil in order to determine the three unknowns. We take into account the left- and right-hand side neighbor cells j − 1 and j + 1 and state that the approximation ${\tilde{g}}^{\left(j\right)}$ should be able to recover the average values ${\overline{g}}^{\left(j-1\right)}$ , ${\overline{g}}^{\left(j\right)}$ and ${\overline{g}}^{\left(j+1\right)}$ correctly. In this case we find from integrating Eq. 31.17

$\begin{array}{ll}{\displaystyle {\int }_{-\frac{3}{2}}^{\frac{1}{2}}{c}_0+c_1{\xi }^{\left(j-1\right)}+{\left(c_2{\xi }^{\left(j-1\right)}\right)}^2}\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\hfill & =\delta x^{\left(j-1\right)}{\overline{g}}^{\left(j-1\right)}\hfill \\ {\displaystyle {\int }_{-\frac{1}{2}}^{\frac{1}{2}}{c}_0+c_1{\xi }^{\left(j\right)}+{\left(c_2{\xi }^{(j}\right)}^2}\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\hfill & =\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}\hfill \\ {\displaystyle {\int }_{\frac{1}{2}}^{\frac{3}{2}}{c}_0+c_1{\xi }^{\left(j+1\right)}+{\left(c_2{\xi }^{\left(j+1\right)}\right)}^2}\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\hfill & =\delta x^{\left(j+1\right)}{\overline{g}}^{\left(j+1\right)}\hfill \end{array}$

for which we find after performing the integration

$\begin{array}{ll}{c}_0-c_1+\frac{13}{12}{c}_1\hfill & =\frac{\delta x^{\left(j-1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j-1\right)}\hfill \\ c_0+\frac{1}{12}{c}_2\hfill & ={\overline{g}}^{\left(j\right)}\hfill \\ c_0+c_1+\frac{13}{12}{c}_1\hfill & =\frac{\delta x^{\left(j+1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j+1\right)}\hfill \end{array}$

which can be solved to yield

$\begin{array}{l}{c}_0=\frac{1}{24}\frac{26\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}-6\delta x^{\left(j-1\right)}{\overline{g}}^{\left(j-1\right)}-\delta x^{\left(j+1\right)}{\overline{g}}^{\left(j+1\right)}}{\delta x^{\left(j\right)}}\\ c_1=-\frac{1}{2}\frac{\delta x^{\left(j-1\right)}{\overline{g}}^{\left(j-1\right)}-\delta x^{\left(j+1\right)}{\overline{g}}^{\left(j+1\right)}}{\delta x^{\left(j\right)}}\\ c_2=-\frac{1}{2}\frac{2\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}-6\delta x^{\left(j-1\right)}{\overline{g}}^{\left(j-1\right)}-\delta x^{\left(j+1\right)}{\overline{g}}^{\left(j+1\right)}}{\delta x^{\left(j\right)}}\end{array}$

In case of a regular grid δx^(j− 1) = δx^(j)= δx^(j+ 1) which simplifies the constants to

$\begin{array}{l}{c}_0=\frac{-{\overline{g}}^{\left(j+1\right)}+26{\overline{g}}^{\left(j\right)}-{\overline{g}}^{\left(j-1\right)}}{24}\\ c_1=\frac{{\overline{g}}^{\left(j+1-{\overline{g}}^{\left(j-1\right)}\right)}}{2}\\ c_2=\frac{{\overline{g}}^{\left(j+1\right)}-2{\overline{g}}^{\left(j\right)}+{\overline{g}}^{\left(j-1\right)}}{2}\end{array}$

in which case our approximation function from Eq. 31.17 becomes

$\begin{array}{l}\tilde{g}\left({\xi }^{\left(j\right)}\right)=\frac{-{\overline{g}}^{\left(j+1\right)}+26{\overline{g}}^{\left(j\right)}-{\overline{g}}^{\left(j-1\right)}}{24}+\frac{{\overline{g}}^{\left(j+1\right)}-{\overline{g}}^{\left(j-1\right)}}{2}{\xi }^{\left(j\right)}+\frac{{\overline{g}}^{\left(j+1\right)}-2{\overline{g}}^{\left(j\right)}+{\overline{g}}^{\left(j-1\right)}}{2}{\left({\xi }^{\left(j\right)}\right)}^2,\\ \begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}\frac{-3}{2}\le {\xi }^{\left(j\right)}\le \frac{3}{2}\end{array}$

Again, this function can be evaluated at the boundaries, i.e., for ${\xi }^{\left(j\right)}=-\frac{3}{2}$ (left boundary of cell j− 1) ${\xi }^{\left(j\right)}=-\frac{1}{2}$ (right boundary of cell j− 1; left boundary of cell j), ${\xi }^{\left(j\right)}=\frac{1}{2}$ (right boundary of cell j; left boundary of cell j + 1), and ${\xi }^{\left(j\right)}=-\frac{3}{2}$ (right boundary of cell j + 1). Similar to the piecewise-linear function Eq. 31.18 can be differentiated finding

$\frac{d\tilde{g}}{d{\xi }^{\left(j\right)}}=\frac{{\overline{g}}^{\left(j+1\right)}-{\overline{g}}^{\left(j-1\right)}}{2}+{\overline{g}}^{\left(j+1\right)}-2{\overline{g}}^{\left(j\right)}+{\overline{g}}^{{}^{\left(j-1\right)}}{\xi }^{\left(j\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}-\frac{3}{2}\le {\xi }^{\left(j\right)}\le \frac{3}{2}$

which can be evaluated at the boundaries for calculating the fluxes.

31.5.6 WENO Schemes

We already briefly touched on the problems of the selection of a suitable stencil. Stencils tend to form oscillatory results, especially near discontinuities in the dependent function. These oscillations may lead to instability and (ultimately) in the inability to produce a useful function approximation. In order to prevent this from happening, stencils are weighted by their susceptibility to oscillation and combined in order to reduce these problems. If we return to the case of the piecewise-linear approximation, we found the stencil for the cells j and j + 1 to be (see Eq. 31.15)

${\tilde{g}}_{j,j+1}\left({\xi }^{\left(j\right)}\right)={\overline{g}}^{\left(j\right)}\left(1-{\xi }^{\left(j\right)}\right)+{\xi }^{\left(j\right)}{\overline{g}}^{\left(j+1\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}-\frac{1}{2}\le {\xi }^{\left(j\right)}\le \frac{3}{2}$

We stated that we could also have taken a stencil using the cells j − 1 and j. In this case, we would have found

$\begin{array}{ll}{\int }_{-\frac{3}{2}}^{-\frac{1}{2}}{c}_0+c_1{\xi }^{\left(j\right)}\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\hfill & =\delta x^{\left(j-1\right)}{\overline{g}}^{\left(j-1\right)}\hfill \\ {\int }_{-\frac{1}{2}}^{\frac{1}{2}}{c}_0+c_1{\xi }^{\left(j\right)}\delta x^{\left(j\right)}d{\xi }^{\left(j\right)}\hfill & =\delta x^{\left(j\right)}{\overline{g}}^{\left(j\right)}\hfill \end{array}$

which we can solve to

$\begin{array}{ll}{c}_0-c_1\hfill & =\frac{\delta x^{\left(j-1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j-1\right)}\hfill \\ c_0\hfill & ={\overline{g}}^{\left(j\right)}\hfill \end{array}$

from which we obtain

$c_0={\overline{g}}^{\left(j\right)}{c}_1={\overline{g}}^{\left(j\right)}-\frac{\delta x^{\left(j-1\right)}}{\delta x^{\left(j\right)}}{\overline{g}}^{\left(j-1\right)}$

which results in the approximation

$\tilde{g}\left({\xi }^{\left(j\right)}\right)={\overline{g}}^{\left(j\right)}\left(1+{\xi }^{\left(j\right)}\right)-{\xi }^{\left(j\right)}\frac{\delta x^{\left(j-1\right)}}{\delta x^{{}^{\left(j\right)}}}{\overline{g}}^{\left(j-1\right)}$

In case of a regular grid with $\delta x^{\left(j-1\right)}=\delta x^{\left(j\right)}$ we find this second stencil to be

${\tilde{g}}_{j-1,j}\left({\xi }^{\left(j\right)}\right)={\overline{g}}^{\left(j\right)}\left(1+{\xi }^{\left(j\right)}\right)-{\xi }^{\left(j\right)}{\overline{g}}^{\left(j-1\right)},\begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}\frac{-3}{2}\le {\xi }^{\left(j\right)}\le \frac{-1}{2}$

Both of these functions ${\tilde{g}}_{j,j+1}$ and ${\tilde{g}}_{j-1,j}$ produced from these stencils are said to be second-order stencil because they take two cells into consideration.

Transition to Third-Order WENO Stencils. We will now create a WENO stencil by using a weighted sum of the two functions ${\tilde{g}}_{j,j+1}$ and ${\tilde{g}}_{j-1,j}$ . The ideal weighting constants are $\frac{2}{3}$ for ${\tilde{g}}_{j,j+1}$ and $\frac{1}{3}$ for ${\tilde{g}}_{j-1,j}$ . These numbers are not random, as we will see in a moment. In doing to, we create a new approximation function that takes the cells j − 1, j, j + 1 into account and that is given by

$\begin{array}{l}{\tilde{g}}_{j-1,j,j+1}\left({\xi }^{\left(j\right)}\right)=\frac{2}{3}\left({\overline{g}}^{\left(j\right)}\left(1-{\xi }^{\left(j\right)}\right)+{\xi }^{\left(j\right)}{\overline{g}}^{\left(j+1\right)}\right)+\frac{1}{3}\left({\overline{g}}^{\left(j\right)}\left(1+{\xi }^{\left(j\right)}\right)-{\xi }^{\left(j\right)}{\overline{g}}^{\left(j-1\right)}\right)\\ ={\overline{g}}^{\left(j\right)}\left(1-\frac{1}{3}{\xi }^{\left(j\right)}\right)+{\xi }^{\left(j\right)}\left(\frac{2}{3}{\overline{g}}^{\left(j+1\right)}-\frac{1}{3}{\xi }^{\left(j\right)}{\overline{g}}^{\left(j-1\right)}\right)\text{'}\\ \begin{array}{llllll}\hfill & \hfill & \hfill & \hfill & \hfill & \hfill \end{array}\frac{-3}{2}\le {\xi }^{\left(j\right)}\le \frac{-3}{2}\end{array}$

If we evaluate this function for ${\xi }^{\left(j\right)}=\frac{1}{2}$ we find

${\tilde{g}}_{j-1,j,j+1}\left(\frac{1}{2}\right)=-\frac{1}{6}{\overline{g}}^{\left(j-1\right)}+\frac{5}{6}{\overline{g}}^{\left(j\right)}+\frac{1}{3}{\overline{g}}^{\left(j+1\right)}$

We have already derived a scheme that uses three cells when deriving the piecewise-parabolic approximation (see Eq. 31.18). If we evaluate this function for ${\xi }^{\left(\mathrm{j}\right)}=\frac{1}{2}$ we find

$\begin{array}{l}\tilde{g}\left(\frac{1}{2}\right)=\frac{-{\overline{g}}^{\left(j+1\right)}+26{\overline{g}}^{\left(j\right)}-{\overline{g}}^{\left(j-1\right)}}{24}+\frac{{\overline{g}}^{\left(j+1\right)}-{\overline{g}}^{\left(j-1\right)}}{4}+\frac{{\overline{g}}^{\left(j+1\right)}-2{\overline{g}}^{\left(j\right)}+{\overline{g}}^{\left(j-1\right)}}{8}\\ ={\overline{g}}^{\left(j-1\right)}\left(-\frac{1}{24}-\frac{1}{4}+\frac{1}{8}\right)+{\overline{g}}^{\left(j\right)}\left(\frac{26}{24}-\frac{1}{4}\right)+{\overline{g}}^{\left(j+1\right)}\left(-\frac{1}{24}-\frac{1}{4}+\frac{1}{8}\right)\\ =-\frac{1}{6}{\overline{g}}^{\left(j-1\right)}+\frac{5}{6}{\overline{g}}^{\left(j\right)}-\frac{1}{3}{\overline{g}}^{\left(j+1\right)}\end{array}$

As you can see, Eq. 31.23 and Eq. 31.24 yield the same result. By adequate weighting with the (seemingly random weighting factors $\frac{1}{3}$ and $\frac{2}{3}$ ) we have obtained the same result as obtained using the stencil from the piecewise-parabolic scheme, which is of third-order. As one would expect, the higher the order of the scheme, the more stable it will be.

Fifth-Order WENO Stencils: WENO-5. If we return to the derivation of piecewise-parabolic approximation in section 31.5.5, we could have chosen one of three stencils

• cells j − 2, j − 1,and j

• cells j − 1,j, and j + l

• cells j, j + 1, and j + 2

We derived the second stencil of this list which is given by Eq. 31.18. If this stencil is evaluated for $\xi =\frac{1}{2}$ we find

$\begin{array}{l}{\tilde{g}}_{j-1,j,j+1}\left(\frac{1}{2}\right)=\frac{-{\overline{g}}^{\left(j+1\right)}+26{\overline{g}}^{\left(j\right)}-{\overline{g}}^{\left(j-1\right)}}{24}+\frac{{\overline{g}}^{\left(j+1\right)}-{\overline{g}}^{\left(j-1\right)}}{4}+\frac{{\overline{g}}^{\left(j+1\right)}-2{\overline{g}}^{\left(j\right)}-{\overline{g}}^{\left(j-1\right)}}{8}\\ ={\overline{g}}^{\left(j-1\right)}\left(-\frac{1}{24}-\frac{1}{4}+\frac{1}{8}\right)+{\overline{g}}^{\left(j-1\right)}\left(\frac{26}{24}-\frac{1}{4}\right)+{\overline{g}}^{\left(j+1\right)}\left(-\frac{1}{24}+\frac{1}{4}+\frac{1}{8}\right)\\ =\frac{-{\overline{g}}^{\left(j-1\right)}+5{\overline{g}}^{\left(j\right)}+2{\overline{g}}^{\left(j+1\right)}}{6}\end{array}$

We will not derive the two other stencils, as they are calculated likewise. When evaluating them for $\xi =\frac{1}{2}$ we find

${\tilde{g}}_{j-2,j-1,j}\left(\frac{1}{2}\right)=\frac{2{\overline{g}}^{\left(j-2\right)}-7{\overline{g}}^{\left(j-1\right)}+11{\overline{g}}^{\left(j\right)}}{6}$

${\tilde{g}}_{j,j+1,j+2}\left(\frac{1}{2}\right)=\frac{2{\overline{g}}^{\left(j\right)}+5{\overline{g}}^{\left(j+1\right)}+11{\overline{g}}^{\left(j+2\right)}}{6}$

Each of the three stencils used to calculate Eq. 31.25, Eq. 31.26, and Eq. 31.27 is third-order. Adding them up forms a fifth-order scheme, which is usually referred to as WENO-5. The weighting can be determined by performing a piecewise-polynomial approximation for the cells j – 2, j − 1, j, j + 1, and j + 2. The scheme is given by

${\tilde{g}}_{j-2,j-1,j,j+1,j+2}\left(\frac{1}{2}\right)=\frac{1}{10}{\tilde{g}}_{j-2,j-1,j}\left(\frac{1}{2}\right)+\frac{6}{10}{\tilde{g}}_{j-1,j,j+1}\left(\frac{1}{2}\right)+\frac{3}{10}{\tilde{g}}_{j,j+1,j+2}\left(\frac{1}{2}\right)$

WENO-5 is a very stable scheme. As you can see, these “taking more cells into account” resemble the procedure taken for numerical schemes for solving ODEs. As an example, take the Adams-Bashforth schemes (see section 27.2.6). Depending on the order of the method, they rely on more data points “before” and “after” the current location j. The more values are taken into account, the better and more stable the numerical scheme will be. This is the same with stencil methods. Obviously, the disadvantage is the fact that higher-order stencils require knowledge about more points to the left and the right of the current location. As we approach the boundaries, these points may fall outside of the domain.

31.5.7 Summary

Stencil methods are required whenever an approximation function must be constructed. As we have seen, if a polynomial of order N = 0 is required, we only need to take cell j into account. These methods are first-order. When the order of the polynomial is N = 1, we obtain a second-order stencil requiring two cells, either j − 1 and j or j and j + 1. For a second-order polynomial, we obtain a third-order stencil which can be based either on the cells j – 2, j – 1, j or cells j – 1, j, j + 1 or cells j, j + 1, j + 1. As we have also seen, it is possible to construct higher-order stencils by appropriately weighting stencils of lower order.

As stated, the approximation function (and its derivative) are required for evaluating the function and the gradient on the boundary of cell j. As you have seen, we usually only evaluate the value for $\xi =\frac{1}{2}$ , which is the right-hand boundary of the current cell.

31.7 Two-Dimensional Problems of First Order in Time and Space

31.7.1 Introduction

The example we just derived was a one-dimensional example containing a first-order partial derivative with respect to time, and a second-order partial derivative with respect to space. In the next step, we want to solve a problem where we have more than one dimension. Returning to the general form of Eq. 31.1, we assume that the vector function F is now not only dependent on the spatial coordinate x₁ = x, but also on the spatial coordinate x₂ = y (see Eq. 31.2). In this case we find

$\begin{array}{ll}\frac{\partial g}{\partial t}\left(x,y,t\right)+\overrightarrow{\nabla }\cdot \overrightarrow{F}\left(x,x,t\right)\hfill & =0\hfill \\ \frac{\partial g}{\partial t}\left(x,y,t\right)+\frac{\partial F_x}{\partial x}\left(x,y,t\right)+\frac{\partial F_y}{\partial y}\left(x,y,t\right)\hfill & =0\hfill \end{array}$

We restrict ourselves to the homogeneous case here. As we have seen, right-hand side terms are treated likewise and they have been covered in section 31.3.3. The domain on which we perform FVM is a plane and shown in Fig. 31.6. The respective cell is indexed with j along the x-directiou and k along the y-direction. As usual for FVM, the surfaces of the cells may not be identical. The domain of the cell (j, k) is given by Ω with the counter-clockwise circumference $\mathit{\partial }$ Ω.

31.7.2 Integration

The concept is identical to the general derivation of FVM as outlined in section 31.3. We first perform the integration of Eq. 31.39 on the domain Ω resulting in

$\begin{array}{ll}{\int }_{\Omega }\frac{\partial g}{\partial t}d\Omega +{\int }_{\Omega }\frac{\partial F_x}{\partial x}d\Omega +{\int }_{\Omega }\frac{\partial F_y}{\partial y}d\Omega \hfill & =0\hfill \\ \frac{d}{dt}{\int }_{\Omega }gd\Omega +{\int }_{\Omega }\left(\frac{\partial F_x}{\partial x}+\frac{\partial F_y}{\partial y}\right)d\Omega \hfill & =0\hfill \end{array}$

where we have rewritten the left-hand side partial differential with respect to time to an ordinary differential as before

31.7.3 Green’s Theorem

We now apply Green’s theorem (see section 7.2.3) that helps in simplifying volume and surface integrals that originate from vector functions (see Eq. 7.13). It states that

${\int }_{\Omega }\left(\frac{\partial F_y}{\partial x}-\frac{\partial F_x}{\partial y}\right)d\Omega ={\oint }_{\partial \Omega }\overrightarrow{F}\cdot d\overrightarrow{s}$

where the vector $d\overrightarrow{s}$ is the infinitesimal vector on the circumference given by

$d\overrightarrow{s}=\left(\begin{array}{c}\hfill dx\hfill \\ \hfill \mathit{dy}\hfill \end{array}\right)$

in which case we can rewrite Eq. 31.41 to

${\int }_{\Omega }\left(\frac{\partial F_y}{\partial x}-\frac{\partial F_x}{\partial y}\right)d\Omega ={\oint }_{\partial \Omega }{F}_{x}dx+F_y\mathit{dy}$

Using Green’s theorem, we have converted the surface integral into an integral along the contour. Please note that the contours of the cells are usually straight lines, which are not necessarily parallel to one of the coordinate axis. When the cells are parallel to the coordinate system, we can make a couple of additional simplifications.

31.7.4 Linearization

Using Eq. 31.42 we can rewrite the integral of Eq. 31.40 as follows

$\begin{array}{l}{\int }_{\Omega }\frac{\partial F_x}{\partial x}+\frac{\partial F_y}{\partial y}d\Omega ={\oint }_{\partial \Omega }{F}_x\mathit{dy}-F_{y}dx={\oint }_{1\to 2\to 3\to 4\to 1}{F}_x\mathit{dy}-F_{y}dx\\ \begin{array}{cc}\hfill \hfill & \hfill \hfill \end{array}={\oint }_{1\to 2}{F}_x\mathit{dy}-F_{y}dx+{\oint }_{2\to 3}{F}_x\mathit{dy}-F_{y}dx-\left({\oint }_{3\to 4}{F}_x\mathit{dy}-F_{y}dx+{\oint }_{4\to 1}{F}_x\mathit{dy}-F_{y}dx\right)\end{array}$

where we have rewritten the line integral on the circumference based on the points of Fig. 31.6. Please note that for all sections of the circumference where $\overrightarrow{r}$ is oriented in the negative direction of the coordinate system, the contribution to the line integral is negative. We now need to determine the individual integrals. For this we need an approach for the functions F_x and F_y along the circumference. As for the one-dimensional case, we assume linear changes of the functions within the cell, which allows us to linearize the differentials. For cell (j, k) we therefore find

$\begin{array}{l}d{x}_{1\to 2}^{\left(j,k\right)}\approx \Delta x_{1\to 2}^{\left(j,k\right)}={x_2}^{\left(j,k\right)}-{x_1}^{\left(j,k\right)}\\ d{x}_{2\to 3}^{\left(j,k\right)}\approx \Delta x_{2\to 3}^{\left(j,k\right)}={x_3}^{\left(j,k\right)}-{x_2}^{\left(j,k\right)}\\ d{x}_{3\to 4}^{\left(j,k\right)}\approx \Delta x_{3\to 4}^{\left(j,k\right)}=x_4^{\left(j,k\right)}-x_3^{\left(j,k\right)}\\ d{x}_{4\to 1}^{\left(j,k\right)}\approx \Delta x_{4\to 1}^{\left(j,k\right)}=x_1^{\left(j,k\right)}-x_4^{\left(j,k\right)}\\ d{y}_{1\to 2}^{\left(j,k\right)}\approx \Delta y_{1\to 2}^{\left(j,k\right)}=y_2^{\left(j,k\right)}-y_1^{\left(j,k\right)}\\ d{y}_{2\to 3}^{\left(j,k\right)}\approx \Delta y_{2\to 3}^{\left(j,k\right)}=y_3^{\left(j,k\right)}-y_2^{\left(j,k\right)}\\ d{y}_{3\to 4}^{\left(j,k\right)}\approx \Delta y_{3\to 4}^{\left(j,k\right)}=y_4^{\left(j,k\right)}-y_3^{\left(j,k\right)}\\ d{y}_{4\to 1}^{\left(j,k\right)}\approx \Delta y_{4\to 1}^{\left(j,k\right)}=y_1^{\left(j,k\right)}-y_4^{\left(j,k\right)}\end{array}$