Finite Difference Method (FDM)

This section describes the formulation and methodology of finite difference method to solve the governing equations on a computational domain.

The finite difference method is the oldest method for the numerical solution of partial differential equations. It is also the easiest to formulate and program for problems which have a simple geometry.

When calculating derivatives finite difference replaces the infinitesimal limiting process with a finite quantity. The derivative of a function at point expressed as: (1)

\frac{\partial f}{\partial x} = lim Δ x \to 0 \frac{f (x + Δ x) - f (x)}{Δ x}

$\frac{\partial f}{\partial x} = \lim_{Δ x \to 0} \frac{f (x + Δ x) - f (x)}{Δ x}$

is replaced by: (2)

\frac{\partial f}{\partial x} = \frac{f (x + Δ x) - f (x)}{Δ x} + O (Δ x)

$\frac{\partial f}{\partial x} = \frac{f (x + Δ x) - f (x)}{Δ x} + O (Δ x)$

The preceding expression is commonly referred to as the forward difference approximation. This derivative can have more refined approximations using a number of approaches such as truncated Taylor series expansions and polynomial fitting.

The term

O (Δ x)

$O (Δ x)$ gives an indication of the magnitude of the error as a function of the mesh spacing and is therefore termed as the order of magnitude of the finite difference method. In the above formulation, the finite difference approximation employed is first-order accurate. A second-order approximation would have the order of magnitude expressed as

O (Δ x^{2})

$O (Δ x^{2})$ . Most of the finite difference methods used in practice are second-order accurate. An example of a second-order method is the central difference approximation of the first derivative. It is expressed as: (3)

\frac{\partial f}{\partial x} = \frac{f (x + Δ x) - f (x - Δ x)}{2 Δ x} + O (Δ x^{2})

$\frac{\partial f}{\partial x} = \frac{f (x + Δ x) - f (x - Δ x)}{2 Δ x} + O (Δ x^{2})$

The basic methodology in the simulation process using finite difference approach consists of the following steps.

The finite difference formulation generally employs a structured grid. The most commonly used indices are $i$ $i$ , $j$ $j$ and $k$ $k$ for the grid lines at $x = x_{i}$ $x = x_{i}$ , $y = y_{j}$ $y = y_{j}$ and $z = z_{k}$ $z = z_{k}$ , respectively. The function value at such a grid point is expressed as $f_{i, j, k} \equiv f_{i j k} \equiv f (x_{i}, y_{j}, z_{k})$ $f_{i, j, k} \equiv f_{i j k} \equiv f (x_{i}, y_{j}, z_{k})$ .

Consider the one dimension advection-diffusion equation for a scalar quantity

φ

$φ$ governed by (4)

\frac{d (ρ u φ)}{d x} = \frac{d}{d x} (ϵ \frac{d φ}{d x})

$\frac{d (ρ u φ)}{d x} = \frac{d}{d x} (ϵ \frac{d φ}{d x})$

where $u$ $u$ is the specified velocity, $ρ$ $ρ$ is the density of the fluid and $ϵ$ $ϵ$ is the diffusivity.

If the domain is defined by the boundary $x = 0$ $x = 0$ and $x = L$ $x = L$ and the boundary conditions by $φ (0) = φ_{0}$ $φ (0) = φ_{0}$ and $φ (L) = φ_{L}$ $φ (L) = φ_{L}$ , the domain can be discretised (non-uniformly) by a total of $N + 1$ $N + 1$ grid points for the finite difference solution of the problem.

The diffusion term can be approximated using the central difference (both the inner and outer derivative) as: (5)

\frac{d}{d x} {(ϵ \frac{d φ}{d x})}_{i} \approx \frac{{(ϵ \frac{d φ}{d x})}_{i + \frac{1}{2}} - {(ϵ \frac{d φ}{d x})}_{i - \frac{1}{2}}}{\frac{1}{2} (x_{i + 1} - x_{i - 1})}

$\frac{d}{d x} {(ϵ \frac{d φ}{d x})}_{i} \approx \frac{{(ϵ \frac{d φ}{d x})}_{i + \frac{1}{2}} - {(ϵ \frac{d φ}{d x})}_{i - \frac{1}{2}}}{\frac{1}{2} (x_{i + 1} - x_{i - 1})}$

(6)

\frac{d}{d x} {(ϵ \frac{d φ}{d x})}_{i} \approx \frac{(ϵ_{i + \frac{1}{2}} \frac{φ_{i + 1} - φ_{i}}{x_{i + 1} - x_{i}}) - (ϵ_{i - \frac{1}{2}} \frac{φ_{i} - φ_{i - 1}}{x_{i} - x_{i - 1}})}{\frac{1}{2} (x_{i + 1} - x_{i - 1})}

$\frac{d}{d x} {(ϵ \frac{d φ}{d x})}_{i} \approx \frac{(ϵ_{i + \frac{1}{2}} \frac{φ_{i + 1} - φ_{i}}{x_{i + 1} - x_{i}}) - (ϵ_{i - \frac{1}{2}} \frac{φ_{i} - φ_{i - 1}}{x_{i} - x_{i - 1}})}{\frac{1}{2} (x_{i + 1} - x_{i - 1})}$

The convection term can also be approximated using the central difference as: (7)

\frac{d}{d x} {(ρ u φ)}_{i} \approx \frac{{(ρ u φ)}_{i + 1} - {(ρ u φ)}_{i - 1}}{x_{i + 1} - x_{i - 1}}

$\frac{d}{d x} {(ρ u φ)}_{i} \approx \frac{{(ρ u φ)}_{i + 1} - {(ρ u φ)}_{i - 1}}{x_{i + 1} - x_{i - 1}}$

If the grid is uniform and the values of density, velocity and coefficient of diffusion are constant throughout the domain the above equations reduce to: (8)

ρ u \frac{φ_{i + 1} - φ_{i - 1}}{2 Δ x} = ϵ \frac{φ_{i + 1} - 2 φ_{i} + φ_{i - 1}}{Δ x^{2}}

$ρ u \frac{φ_{i + 1} - φ_{i - 1}}{2 Δ x} = ϵ \frac{φ_{i + 1} - 2 φ_{i} + φ_{i - 1}}{Δ x^{2}}$

The resulting set of linear algebraic equations can be solved to get the values of at the grid points.

Finite difference methods are advantageous for the numerical solution of partial differential equations because of their simplicity, efficiency and low computational cost. Their major drawback is their geometrical inflexibility, such as application on an unstructured grid or moving boundaries. Their formulation increases in complexity as the complexity of the domain increases.

The restrictions resulting from the above mentioned geometrical complexities can be alleviated by the use of methods such as grid transformation and immersed boundary techniques.