**4. Numerical methods to solve governing equations**

The governing equations discussed above are almost impossible to solve analytically. The only feasible pathway to solve these partial differential equations is employing numerical methods. This is the reason most CFD software need to discretize the fluid domain to generate a mesh or grid. The core idea of numerical methods is transforming differential to difference quotient. The commonly used numerical methods in most CFD software include the finite difference method and the finite volume method.

#### **4.1 Finite difference method**

The finite difference method is the first method for solving partial differential equations. In finite difference method, the Taylor series is employed to generate approximations to the partial derivatives of the CFD governing equations. According to the Taylor series, the value of a general variable ϕ at node i + 1 (**Figure 12**) could be expressed as:

$$\phi\_{i+1} = \phi\_i + \left(\frac{\partial \phi}{\partial x}\right)\_i \Delta x + \left(\frac{\partial^2 \phi}{\partial x^2}\right)\_i \frac{\Delta x^2}{2} + \left(\frac{\partial^3 \phi}{\partial x^3}\right)\_i \frac{\Delta x^3}{6} + \left(\text{Higher Order Terms}\right) \tag{103}$$

And the value of a general variable φ at node *i*-1 is:

$$\phi\_{i-1} = \phi\_i - \left(\frac{\partial \phi}{\partial \mathbf{x}}\right)\_i \Delta \mathbf{x} + \left(\frac{\partial^2 \phi}{\partial \mathbf{x}^2}\right)\_i \frac{\Delta \mathbf{x}^2}{2} - \left(\frac{\partial^3 \phi}{\partial \mathbf{x}^3}\right)\_i \frac{\Delta \mathbf{x}^3}{6} + \left(\text{Higher Order Terms}\right) \tag{104}$$

If we subtract Eq. (104) from Eq. (103), we get

$$
\left(\frac{\partial \phi}{\partial \mathbf{x}}\right)\_i = \frac{\phi\_{i+1} - \phi\_{i-1}}{2\Delta \mathbf{x}} + \frac{(\Delta \mathbf{x}^3/\mathfrak{Z})\left(\partial^3 \phi/\partial \mathbf{x}^3\right)}{2\Delta \mathbf{x}}\tag{105}
$$

Since Δ*x* is approaching zero, the second term is very small compared to the first term. By omitting the Δ*x*<sup>3</sup> term, the central difference is achieved

$$
\left(\frac{\partial \phi}{\partial \mathbf{x}}\right)\_i = \frac{\phi\_{i+1} - \phi\_{i-1}}{2\Delta \mathbf{x}}\tag{106}
$$

From Eq. (103) and assuming Δ*x* is approaching zero, we can get the forward difference

$$\left(\frac{\partial \phi}{\partial \mathbf{x}}\right)\_i = \frac{\phi\_{i+1} - \phi\_i}{\Delta \mathbf{x}} + O(\Delta \mathbf{x}) \tag{107}$$

From Eq. (104) and assuming Δ*x* is approaching to zero, we can get the backward difference

$$\left(\frac{\partial\phi}{\partial\mathbf{x}}\right)\_i = \frac{\phi\_i - \phi\_{i-1}}{\Delta\mathbf{x}} + O(\Delta\mathbf{x}) \tag{108}$$

Eqs. (106)-(108) successfully approximated the first-order derivative by difference quotient. First-order derivative is actually a slope of a variable. **Figure 13** geometrically represents the slopes by backward difference, central difference and forward difference. It is clear that the slope by central difference (blue line) is closer to the real slope (black line), so the central difference could represent the exact partial differential with a higher accuracy.

**Figure 12.** *One-dimensional grids for the finite-difference method.*

*The Basic Theory of CFD Governing Equations and the Numerical Solution Methods… DOI: http://dx.doi.org/10.5772/intechopen.113253*

**Figure 13.** *Finite-difference representation of the derivative for ∂*ϕ*/∂*x*.*

In CFD governing equations, there are second-order partial derivatives in diffusion terms, so the second-order derivatives should also be approximated by finite difference method. Adding Eqs. (103) and (104) together yields

$$
\Delta \phi\_{i+1} + \phi\_{i-1} = 2\phi\_i + \left(\frac{\partial^2 \phi}{\partial \mathbf{x}^2}\right)\_i \Delta \mathbf{x}^2 + \left(\frac{\partial^4 \phi}{\partial \mathbf{x}^4}\right)\_i \frac{\Delta \mathbf{x}^4}{24} \tag{109}
$$

By omitting higher-order terms, the difference quotient form in finite difference method for second-order derivative is

$$\left(\frac{\partial^2 \phi}{\partial \mathbf{x}^2}\right) = \frac{\phi\_{i+1} - 2\phi\_i + \phi\_{i-1}}{\Delta \mathbf{x}^2} \tag{110}$$

In a 3-D fluid domain such as **Figure 14**, Taylor series expansion also applies for y and z directions. The central differences for first- and second-order derivatives are:

$$
\left(\frac{\partial \phi}{\partial \mathbf{y}}\right)\_j = \frac{\phi\_{j+1} - \phi\_{j-1}}{2\Delta \mathbf{y}}\tag{111}
$$

$$
\left(\frac{\partial \phi}{\partial z}\right)\_k = \frac{\phi\_{k+1} - \phi\_{k-1}}{2\Delta z} \tag{112}
$$

$$\left(\frac{\partial^2 \phi}{\partial \mathbf{y}^2}\right) = \frac{\phi\_{j+1} - 2\phi\_j + \phi\_{j-1}}{\Delta \mathbf{y}^2} \tag{113}$$

$$\left(\frac{\partial^2 \phi}{\partial \mathbf{z}^2}\right) = \frac{\phi\_{k+1} - 2\phi\_k + \phi\_{k-1}}{\Delta \mathbf{z}^2} \tag{114}$$

For time term derivatives, only the first-order derivative needs to be approximated by the finite difference method. In most cases, the backward difference is used (Eq. (115)). This is because the variable value at the previous time step has already been calculated, and then, it is available to calculate the current time step in explicit expression. Since the computation performs in chronological order, the variable

**Figure 14.** *Three-dimensional grids for the finite-difference method.*

values at future steps are not available. Forward difference and central difference involve the variable at future steps (Eqs. (116) and (117)), so it is not recommended to be used in explicit solver. It can only be solved in an implicit solver.

$$\left(\frac{\partial\phi}{\partial t}\right) = \frac{\phi\_{i,j,k}{}^n - \phi\_{i,j,k}{}^{n-1}}{\Delta t} \tag{115}$$

$$\left(\frac{\partial\phi}{\partial t}\right) = \frac{\phi\_{i,j,k}{}^{n+1} - \phi\_{i,j,k}{}^{n}}{\Delta t} \tag{116}$$

$$\left(\frac{\partial\phi}{\partial t}\right) = \frac{\phi\_{i,j,k}{}^{n+1} - \phi\_{i,j,k}{}^{n\cdot 1}}{2\Delta t} \tag{117}$$

As we know it is easy to transform difference quotient to linear equation. With the finite element method, all the CFD partial differential equations can be converted to linear equations. By applying the techniques of solving linear equations, the numerical solutions of fluid field could be obtained.

#### **4.2 Finite volume method**

For structured mesh (**Figure 14**), it is ideal to use the finite difference method introduced above. However, most CFD simulations have to use an unstructured mesh (**Figure 15**) due to the complexity of geometries, so it is difficult to define the node sequence *i*, *j* and *k*. Therefore, it is almost impossible to use the finite difference method. In the case of unstructured mesh, we can resort to another popular numerical method which is called the finite volume method.

*The Basic Theory of CFD Governing Equations and the Numerical Solution Methods… DOI: http://dx.doi.org/10.5772/intechopen.113253*

**Figure 15.** *Examples of unstructured mesh.*

Finite volume method is derived from Gauss's divergence theorem. In divergence theorem, the integration of a variable divergence over a body equals the integration of the variable over the surface which encloses the body.

$$\begin{aligned} \iiint\limits\_{\Omega} \frac{\partial \phi\_x}{\partial \mathbf{x}} + \frac{\partial \phi\_y}{\partial \mathbf{y}} + \frac{\partial \phi\_x}{\partial \mathbf{z}} \mathrm{d}V &= \mathfrak{F}\_{\Sigma} \overrightarrow{\phi}(\mathbf{x}, y, \mathbf{z}) \mathrm{d}\mathbf{S} \\ &= \mathfrak{F}\_{\Sigma} \phi\_x \mathrm{d}y \mathrm{d}z + \mathfrak{F}\_{\Sigma} \phi\_y \mathrm{d}\mathbf{x} \mathrm{d}z + \mathfrak{F}\_{\Sigma} \phi\_x \mathrm{d}\mathbf{x} \mathrm{d}y \end{aligned} \tag{118}$$

It can also be written in x, y, z directions as follows

$$\iiint\_{\Omega} \frac{\partial \phi}{\partial z} \, \mathrm{d}V = \iint\_{\Sigma} \phi(\mathbf{x}, y, z) \, \mathrm{d}x \, \mathrm{d}y \tag{119}$$

$$\iiint\_{\varOmega} \frac{\partial \phi}{\partial \mathbf{x}} \mathrm{d}V = \iint\_{\varXi} \phi(\mathbf{x}, \mathbf{y}, \mathbf{z}) \mathrm{d}\mathbf{y} \mathrm{d}\mathbf{z} \tag{120}$$

$$\iiint\_{\Omega} \frac{\partial \phi}{\partial \mathbf{y}} \mathrm{d}V = \iint\_{\Sigma} \phi(\mathbf{x}, \mathbf{y}, \mathbf{z}) \mathrm{d}\mathbf{x} \mathrm{d}\mathbf{z} \tag{121}$$

In CFD, we assume ∂ϕ/∂*x* only has variation between different grids. Inside a single grid, it has only a constant value of ∂ϕ/∂*x*. The variation inside a single grid is very small and could be neglected. Therefore, for a grid, the left side of Gauss's divergence theorem could be simplified as

$$\iiint\_{\varOmega} \frac{\partial \phi}{\partial \mathbf{x}} \, \mathrm{d}V = \frac{\partial \phi}{\partial \mathbf{x}} \iiint\_{\varOmega} \mathrm{d}V = \frac{\partial \phi}{\partial \mathbf{x}} \, \mathrm{d}V \tag{122}$$

Similarly, the variation of ϕ on a grid surface could also be neglected if we assume a grid surface is too small to have variation of ϕ. The right side of Gauss's divergence theorem could be simplified as well

*Computational Fluid Dynamics – Recent Advances, New Perspectives and Applications*

$$\iint\_{\Sigma} \phi \mathbf{d}y \mathbf{d}z = \phi \mathbf{A}^{\mathbf{x}} \tag{123}$$

The divergence theorem correlated the ∂ϕ/∂*x* to ϕ, and we can utilize this correlation to get the approximation of the first-order derivative in a tetrahedral grid (**Figure 16**) as follows

$$\begin{split} \frac{\partial \phi}{\partial \mathbf{x}} &= \frac{1}{\Delta V} \int\_{\Delta V} \frac{\partial \phi}{\partial \mathbf{x}} \mathrm{d}V = \frac{1}{\Delta V} \int\_{A} \phi \mathrm{d}A^{x} \\ &= \frac{1}{\Delta V} \left( \phi\_{1} A\_{1}^{x} + \phi\_{2} A\_{2}^{x} + \phi\_{3} A\_{3}^{x} + \phi\_{4} A\_{4}^{x} \right) \\ &= \frac{1}{\Delta V} \sum\_{i=1}^{N} \phi\_{i} A\_{i}^{x} \end{split} \tag{124}$$

where *Ai <sup>x</sup>* is the projected area to the y-o-z surface.

The first-order derivative with respect to y and z can be similarly obtained as follows

$$\begin{split} \frac{\partial \phi}{\partial \mathbf{y}} &= \frac{1}{\Delta V} \int\_{\Delta V} \frac{\partial \phi}{\partial \mathbf{y}} \mathrm{d}V = \frac{1}{\Delta V} \int\_{A} \phi \mathrm{d}A^{\mathbf{y}} \\ &= \frac{1}{\Delta V} \left( \phi\_1 A\_1^{\mathbf{y}} + \phi\_2 A\_2^{\mathbf{y}} + \phi\_3 A\_3^{\mathbf{y}} + \phi\_4 A\_4^{\mathbf{y}} \right) \\ &= \frac{1}{\Delta V} \sum\_{i=1}^{N} \phi\_i A\_i^{\mathbf{y}} \end{split} \tag{125}$$

**Figure 16.** *The representation of a tetrahedral grid.*

*The Basic Theory of CFD Governing Equations and the Numerical Solution Methods… DOI: http://dx.doi.org/10.5772/intechopen.113253*

$$\begin{aligned} \frac{\partial \phi}{\partial z} &= \frac{1}{\Delta V} \int\_{\Delta V} \frac{\partial \phi}{\partial z} \mathrm{d}V = \frac{1}{\Delta V} \int\_{A} \phi \mathrm{d}A^{z} \\\\ &= \frac{1}{\Delta V} \left( \phi\_{1} A\_{1}^{z} + \phi\_{2} A\_{2}^{z} + \phi\_{3} A\_{3}^{z} + \phi\_{4} A\_{4}^{z} \right) \\\\ &= \frac{1}{\Delta V} \sum\_{i=1}^{N} \phi\_{i} A\_{i}^{z} \end{aligned} \tag{126}$$

In the same way, the second-order derivative can be approximated by the firstorder derivative

$$\begin{split} \frac{\partial^2 \phi}{\partial \mathbf{x}^2} &= \frac{1}{\Delta V} \int\_{\Delta V} \frac{\partial^2 \phi}{\partial \mathbf{x}^2} \mathbf{d}V = \frac{1}{\Delta V} \int\_A \frac{\partial \phi}{\partial \mathbf{x}} \mathbf{d}A^x \\ &= \frac{1}{\Delta V} \left[ \left( \frac{\partial \phi}{\partial \mathbf{x}} \right)\_1 A\_1^x + \left( \frac{\partial \phi}{\partial \mathbf{x}} \right)\_2 A\_2^x + \left( \frac{\partial \phi}{\partial \mathbf{x}} \right)\_3 A\_3^x + \left( \frac{\partial \phi}{\partial \mathbf{x}} \right)\_4 A\_4^x \right] \\ &= \frac{1}{\Delta V} \sum\_{i=1}^N \left( \frac{\partial \phi}{\partial \mathbf{x}} \right)\_i A\_i^x \end{split} \tag{127}$$

$$\frac{\partial^2 \phi}{\partial \mathbf{y}^2} = \frac{1}{\Delta V} \sum\_{i=1}^N \left( \frac{\partial \phi}{\partial \mathbf{y}} \right)\_i A\_i^{\mathbf{y}} \tag{128}$$

$$\frac{\partial^2 \phi}{\partial x^2} = \frac{1}{\Delta V} \sum\_{i=1}^N \left(\frac{\partial \phi}{\partial x}\right)\_i A\_i^x \tag{129}$$

From the equations above, it is clear that the finite volume method has better conservativeness of the CFD governing equations. Meanwhile it is more adaptive to complicated geometries. In most CFD commercial software, the finite volume method is the mainstream technique to solve the CFD governing equations.
