The Basic Theory of CFD Governing Equations and the Numerical Solution Methods for Reactive Flows

Abstract

The universal principles of fluid motion are the conservation of mass, momentum and energy. This chapter will introduce the CFD governing equations and describe how the continuity equation, component equation, Navier-Stokes equation and energy equation were derived from the principles above. With the expanding application of CFD simulation technology, some processes such as fluid-involved reactions, adsorption and permeation, which break the conservation of mass, momentum and energy for fluid phase, should be coupled to CFD model. In view of this, this chapter provided the theories about source terms for the mass equation, momentum equation and thermal energy equation. The technology for solving these governing equations remained a challenge for a long period due to the complexity. Thanks to the development of numerical methods, such as the finite difference method and the finite volume method, these equations can be solved and provide reasonable numerical results of flows, heat transfer and reactions. This chapter also demonstrates the basics of these two major numerical techniques.

1. Introduction

Using Computational Fluid Dynamics (CFD) simulation is of great importance in modeling reactive flows because it allows for the prediction and analysis of complex fluid dynamics and chemical reactions in a controlled and cost-effective environment [1]. In fluid reactors, the reaction rate depends on many factors such as concentrations of reactant fluid and product fluid, temperature and turbulence, and these key factors vary with time and space coordinate [2]. Knowing the evolution and distribution of the heat and mass transfer factors is important in reactor design and operating condition optimization. CFD simulations can be used to model and simulate the behavior of reactive flows, such as combustion [3], pyrolysis [4, 5, 6], gasification [7] and some other chemical reactions-driven fluid flow [8], and get detailed information on heat and mass transfer inside reactors.

There are many advantages of CFD over other traditional research methods. For example, CFD simulations can predict the behavior of combustion, which can be used to optimize combustion systems and reduce emissions. CFD simulations can also be used to optimize the design of reaction systems, such as gas combustion in turbines, boilers, and engines, pyrolysis or gasification in thermal reactors by predicting the performance of different designs and optimizing the operating conditions for maximum efficiency. Safety analysis can be performed by CFD simulations as well, by predicting the behavior of the system under different conditions and identifying potential risks. A distinct advantage of CFD simulations is that it can provide a cost-effective means of analyzing and optimizing reactive flow systems, as they can be used to simulate and analyze the behavior of the system under a wide range of conditions without the need for expensive physical testing [9].

The results of a CFD simulation are obtained through numerical methods that solve the fluid flow governing equations using a finite difference or a finite volume method. The numerical methods used in CFD simulations involve dividing the fluid domain into a grid or mesh of discrete cells and then solving the equations governing fluid flow and heat transfer at each of these grid points or cells. The numerical methods used in CFD simulations are complex and require high-performance computing resources. The accuracy of the results depends on the quality of the mesh, the numerical methods used and the convergence criteria [10]. Therefore, CFD simulations require a deep understanding of numerical methods with careful validation and verification to ensure that the results are accurate and reliable.

To better understand CFD simulations, learning the basics of CFD theory is essential for a few reasons. CFD basic theory provides a fundamental understanding of the governing equations and numerical methods used in CFD simulations. This understanding is essential for developing accurate and reliable CFD simulations. A thorough understanding of CFD theory allows for the development of more efficient and accurate simulations. This is because it enables the user to select appropriate numerical methods, grid sizes and convergence criteria. CFD simulations can be complex, and errors can arise due to various reasons. Understanding CFD theory allows the user to identify the source of the error and develop solutions to fix it. CFD theory is constantly evolving, and researchers are continually developing new methods and algorithms to improve the accuracy and efficiency of CFD simulations. A strong understanding of CFD theory is essential for contributing to the advancement of the field.

2. Governing equations

As fluid is one of the existing states of matter, it follows the laws of motion, such as mass conservation, momentum conservation and energy conservation. These laws govern the motion and force of fluid, and thus are called governing equations. These equations are generally expressed in differentiation forms.

2.1 Continuity equation

If our object is an infinitely small grid, continuity indicates that the mass of fluid entering this small grid equals the mass of fluid flowing out of this grid plus the mass of fluid accumulating in this grid. Since the continuity equation is derived from this mass conservation concept, it can also be regarded as a mass conservation equation which is expressed as Eq. (1).

$$\frac{\mathrm{d}m}{\mathrm{d}t}=\sum _{\text{in}}\dot{m}-\sum _{\text{out}}\dot{m}$$E1

where m is the mass of this small grid and is defined as m = ρ∆x∆y∆z. The mass entering and exiting the small grid in the Cartesian coordinate is demonstrated in Figure 1. In the actual scenario, the flow direction does not have to be the same as Figure 1 shows. Any reversed flow direction can be considered if a minus sign is in the front of a term. According to Figure 1, the mass entering this grid with the flow is

Figure 1.

The schematic of mass conservation in a control volume.

$$\sum _{\text{in}}\dot{m}=\mathit{\rho u}\mathrm{\Delta }y\mathrm{\Delta }z+\mathit{\rho v}\mathrm{\Delta }x\mathrm{\Delta }z+\mathit{\rho w}\mathrm{\Delta }x\mathrm{\Delta }y$$E2

Based on the Taylor series expansion for the ρu, ρv, and ρw terms by omitting the higher order terms, the mass exiting out of this grid is

$$\sum _{\text{out}}\dot{m}=\left[\mathit{\rho u}+\frac{\partial \left(\mathit{\rho u}\right)}{\partial x}\mathrm{\Delta }x\right]\mathrm{\Delta }y\mathrm{\Delta }z+\left[\mathit{\rho v}+\frac{\partial \left(\mathit{\rho v}\right)}{\partial y}\mathrm{\Delta }y\right]\mathrm{\Delta }x\mathrm{\Delta }z+\left[\mathit{\rho w}+\frac{\partial \left(\mathit{\rho w}\right)}{\partial z}\mathrm{\Delta }z\right]\mathrm{\Delta }x\mathrm{\Delta }y$$E3

Substituting Eqs. (2) and (3) into Eq. (1), one can obtain Eq. (4).

$$\begin{array}{ll}\frac{\partial \left(\rho \mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\right)}{\partial t}& =\mathit{\rho u}\mathrm{\Delta }y\mathrm{\Delta }z+\mathit{\rho v}\mathrm{\Delta }x\mathrm{\Delta }z+\mathit{\rho w}\mathrm{\Delta }x\mathrm{\Delta }y\\ & -\left[\mathit{\rho u}+\frac{\partial \left(\mathit{\rho u}\right)}{\partial x}\mathrm{\Delta }x\right]\mathrm{\Delta }y\mathrm{\Delta }z\\ & -\left[\mathit{\rho v}+\frac{\partial \left(\mathit{\rho v}\right)}{\partial y}\mathrm{\Delta }y\right]\mathrm{\Delta }x\mathrm{\Delta }z\\ & -\left[\mathit{\rho w}+\frac{\partial \left(\mathit{\rho w}\right)}{\partial z}\mathrm{\Delta }z\right]\mathrm{\Delta }x\mathrm{\Delta }y\end{array}$$E4

After canceling the identical terms with opposite signs, dividing ∆x∆y∆z on both sides, and taking the limit for ∆x → 0, ∆y → 0 and ∆z → 0, we have Eq. (5).

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\mathit{\rho u}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho v}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho w}\right)}{\partial z}=0$$E5

Then the continuity equation is derived. Sometimes it is also written in a concise vector form as Eq. (6).

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \overrightarrow{u}\right)=0$$E6

2.2 Component balance equation

Even though the continuity equation governs the total mass balance of fluid during its flowing process, for a multi-component flow, the mass balance of each component must also be ensured. The mass transfer of a component is more complex than the total flow. In addition to convective flow, a component could also diffuse under the driving force of the concentration gradient from the high concentration region to the low concentration region.

Figure 2 is a demonstration of diffusive flow. Initially, a chamber is separated by a valve into two sub-chambers with equal volume and pressure. In the left sub-chamber, there are 8 moles of N₂ and 2 moles of CO, while in the right chamber there are 2 moles of N₂ and 8 moles of CO. Once the valve is open, there is no convective flow due to the equal pressure on both sides, but N₂ would diffuse from left to right due to the concentration difference between the left and the right sides. In the meantime, CO would diffuse the other way around owing to the same reason. During this process, there is no convective flow at all, but the mass transfer of each component really occurred. Therefore, when considering the mass balance of a component in a multi-component flow, the diffusive flow or diffusion must be included. Generally, the diffusion follows Fick’s law [11]:

Figure 2.

Demonstration of diffusive flow.

$$J_i=D\frac{\mathrm{d}\left(\rho X_i\right)}{\mathrm{d}x}$$E7

where J_i denotes the diffusive flux of component I, D is called diffusivity which is a function of temperature and pressure, and could be evaluated based on Fuller Eq. [12], ρ is the flow density and X_i is the mass fraction of component I in the fluid.

In multi-component flow, the ith component follows the mass balance of Eq. (8)

$$\frac{\mathrm{d}{m}_i}{\mathrm{d}t}=\sum _{\text{in}}{\dot{m}}_i-\sum _{\text{out}}{\dot{m}}_i$$E8

where m_i is the mass of the ith component in a control volume. The ith component entering and exiting the control volume in terms of convective flow in Cartesian coordinate is demonstrated in Figure 3. Taking the x direction as an example, the total flow into this control volume from the upwind y-o-z surface is ρu∆y∆z (the left blue arrow). If the mass fraction of component I is X_i, the flow of component I into the control volume together with the total convective flow is X_iρu∆y∆z. The convective flow of ith component exiting from the downwind y-o-z surface (the right blue arrow) is then expressed by Taylor series expansion by omitting higher-order terms {X_iρu + [∂(X_iρu)/∂x]∆x}∆y∆z. Likewise, the ith component entering and exiting along y and z directions are X_iρv∆x∆z, X_iρw∆y∆x, {X_iρv + [∂(X_iρv)/∂y]∆y}∆x∆z, {X_iρw + [∂(X_iρw)/∂z]∆z}∆y∆x, respectively.

Figure 3.

The schematic of the convective flow of I component entering and exiting a control volume.

In addition to the convective flow depicted in Figure 3, there is also diffusive flow driven by the concentration gradient (Figure 4). According to Fick’s law, the diffusive flow of component I entering from the upwind surface of y-o-z is D[∂(ρX_i)/∂x]∆y∆z. Then, the diffusive flow exiting from the opposite surface is D{∂(ρX_i)/∂x + [∂²(ρX_i)/∂x²]∆x}∆y∆z. Likewise, along the y direction, the diffusive flows entering and exiting are D[∂(ρX_i)/∂y]∆x∆z and D{∂(ρX_i)/∂y + [∂²(ρX_i)/∂y²] ∆y}∆x∆z. Similarly, along the z direction, the diffusive flow entering and exiting are D[∂(ρX_i)/∂z]∆x∆y and D{∂(ρX_i)/∂z + [∂²(ρX_i)/∂z²] ∆z}∆x∆y.

Figure 4.

The schematic of the diffusive flow of I component entering and exiting a control volume.

Substituting all the convective flow and diffusive flow expressions into Eq. (8), and canceling out the identical terms with opposite signs leads to Eq. (9).

$$\begin{array}{ll}\frac{\partial }{\partial t}\left(\rho X_i\mathit{\text{ΔxΔyΔz}}\right)& =D\frac{{\partial }^2\left(\rho X_i\right)}{\partial x^2}\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z+D\frac{{\partial }^2\left(\rho X_i\right)}{\partial y^2}\mathrm{\Delta }y\mathrm{\Delta }x\mathrm{\Delta }z+D\frac{{\partial }^2\left(\rho X_i\right)}{\partial z^2}\mathrm{\Delta }z\mathrm{\Delta }x\mathrm{\Delta }y\\ & -\frac{\partial \left(\rho {\mathit{uX}}_i\right)}{\partial x}\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z-\frac{\partial \left(\rho {\mathit{vX}}_i\right)}{\partial y}\mathrm{\Delta }y\mathrm{\Delta }x\mathrm{\Delta }z-\frac{\partial \left(\rho {\mathit{wX}}_i\right)}{\partial z}\mathrm{\Delta }z\mathrm{\Delta }x\mathrm{\Delta }y\end{array}$$E9

After dividing ∆x∆y∆z for both sides, Eq. (9) could be rearranged to Eq. (10)

$$\frac{\partial \left(\rho X_i\right)}{\partial t}+\frac{\partial \left(\rho {\mathit{uX}}_i\right)}{\partial x}+\frac{\partial \left(\rho {\mathit{vX}}_i\right)}{\partial y}+\frac{\partial \left(\rho {\mathit{wX}}_i\right)}{\partial z}=D\left[\frac{{\partial }^2\left(\rho X_i\right)}{\partial x^2}+\frac{{\partial }^2\left(\rho X_i\right)}{\partial y^2}+\frac{{\partial }^2\left(\rho X_i\right)}{\partial z^2}\right]$$E10

which is the final expression of the component balance equation. Component equation is especially important in governing the flow with chemical reactions, because the reasonable prediction of reaction rates depends on accurate information about component distribution in a reactor.

2.3 Momentum equation

Newton’s second law is the principle to express the correlation between force and motion of the fluid. Namely, the force acting on an infinitesimal control volume equals the product of its mass and acceleration, which is given as (Eq. (11)).

$$\overrightarrow{F}=m\overrightarrow{a}$$E11

Actually, both force and acceleration are vectors and should be expressed by the three components in the x, y and z directions. The acceleration of fluid motion is defined as the total differential with respect to time (Eq. (12)).

$$a_x=\frac{\mathrm{D}u}{\mathrm{D}t}$$E12

$$a_y=\frac{\mathrm{D}v}{\mathrm{D}t}$$E13

$$a_z=\frac{\mathrm{D}w}{\mathrm{D}t}$$E14

Taking the x direction as an example, the acceleration in the x direction is

$$a_x=\frac{\mathrm{D}u}{\mathrm{D}t}=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}$$E15

According to Newton’s second law, the acceleration in the x direction should equal the quotient of the forces along the x direction divided by mass.

$$a_x=\frac{\sum F_x}{m}$$E16

Basically, there are two types of forces. Namely, surface forces such as normal stress, shear stress, and body forces such as gravity force and magnetic force. Figure 5 depicts all the surface forces acting on a control volume along the x direction. Firstly, there is a pair of normal stresses imposed on the left and right surfaces that are indicated as blue arrows. This surface force is caused by internal pressure and velocity gradient in x direction. Secondly, there is a pair of shear stresses on the top and bottom surfaces that are marked as yellow arrows. This surface force is caused by the velocity gradient. Thirdly, a pair of shear stresses are also acting on the front and back surfaces shown as red arrows.

Figure 5.

The schematic of the surface forces along x direction on a control volume.

The summation of all the surface forces is

$$\sum F_{\text{surf}}=\frac{\partial {\sigma }_{\mathit{xx}}}{\partial x}\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z+\frac{\partial {\tau }_{\mathit{yx}}}{\partial y}\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z+\frac{\partial {\tau }_{\mathit{zx}}}{\partial z}\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z$$E17

In fluid mechanics, the normal stress is related to the pressure and velocity gradient in the fluid, and is given as Eq. (18) [13]

$${\sigma }_{\mathit{xx}}=-p+2\mu \frac{\partial u}{\partial x}$$E18

The shear stress follows Newton’s law of viscosity as depicted in Eq. (16) [13]

$${\tau }_{\mathit{xy}}={\tau }_{\mathit{yx}}=\mu \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$E19

Bringing Eqs. (18) and (19) to Eq. (17) yields Eq. (20) for the total surface forces in the x direction.

$$\begin{array}{ll}\sum F_{\text{surf}}& =\frac{\partial }{\partial x}\left(-p+2\mu \frac{\partial u}{\partial x}\right)\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z+\mu \frac{\partial }{\partial y}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\\ & +\mu \frac{\partial }{\partial z}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\end{array}$$E20

In most cases, the body force only includes gravity force.

$$\sum F_{\text{body}}=\rho \mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }{\mathit{zg}}_x$$E21

Then, the total force in the x direction is

$$\begin{array}{ll}\sum F_x& =\frac{\partial }{\partial x}\left(-p+2\mu \frac{\partial u}{\partial x}\right)\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z+\mu \frac{\partial }{\partial y}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\\ & +\mu \frac{\partial }{\partial z}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z+\rho \mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }{\mathit{zg}}_x\end{array}$$E22

Eqs. (22) and (16) yield

$$a_x=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial x^2}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial y^2}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial z^2}+\frac{\mu }{\rho }\frac{\partial }{\partial x}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)+g_x$$E23

The term inside the parenthesis is the expression of the continuity equation which equals zero

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$E24

Then, Eq. (23) could be further simplified to

$$a_x=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial x^2}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial y^2}+\frac{\mu }{\rho }\frac{{\partial }^{2}u}{\partial z^2}+g_x$$E25

Combining Eqs. (15) and (25), we have

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)+g_x$$E26

Finally, Eq. (26) is the formula of the momentum equation of fluid in the x direction. Likewise, in the y direction and z direction, the momentum equations of fluid flow could also be derived in a similar way as follows:

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}+\frac{{\partial }^{2}v}{\partial z^2}\right)+g_y$$E27

$$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}+\frac{{\partial }^{2}w}{\partial z^2}\right)+g_z$$E28

Since Eq. (26) was progressively developed by Claude-Louis Navier and George Gabriel Stokes, this equation is generally called the Navier-Stokes Equation.

2.4 Energy equation

The heat transfer of fluid should follow the law of energy conservation. The total energy of fluid includes internal energy, kinetic energy and gravitational energy. Practically, it is not convenient to consider the total energy of fluid during the calculation. A better way to include kinetic energy and gravitational energy is by adding source terms in the energy equation. In a CFD model, the main energy form considered is thermal energy, so the thermal energy balance is usually considered for the energy equation. The thermal energy accumulated in an infinitesimal control volume equals the heat entering the control volume in terms of work, convection, conduction and radiation minus the heat going out. As a tradition, radiation is always not taken into account in the energy equation. The rate of energy change in a control volume should equal the net rate of heat added plus the net rate of work done. The heat added and work done could be either positive or negative depending on the direction of energy flow. The energy conservation is given as

$$\rho \frac{\mathit{DE}}{\mathit{Dt}}\mathit{\text{ΔxΔyΔz}}=\sum \dot{Q}+\sum \dot{W}$$E29

where E is the energy contained in unit mass. For most fluid flow cases with heat exchange, reactions or adsorptions, temperature is the key parameter, so the energy of interest is usually thermal energy.

Since we have analyzed all the x-direction surface forces on a control volume when we derived the Navier-Stokes equation, the rate of heat generation due to work done by each normal force and shear force is shown in Figure 6. Likewise, there are also same number of normal forces and shear forces along y-direction and z-direction, and they generate heat as well.

Figure 6.

The schematic of the work done by surface forces on a control volume.

Heat conduction caused by temperature gradient is also happening in heat transfer. The thermal energy entering and exiting the control volume in x-direction is described in Figure 7. Besides this pair of heat fluxes, in the pairs of x-o-z surfaces and x-o-y surfaces the conductive heat fluxes also exist.

Figure 7.

The schematic of the conductive heat flow in a control volume.

Bring all the work heat and conductive heat into Eq. (29) and dividing ∆x∆y∆z, we have

$$\begin{array}{ll}\rho \frac{\mathit{DE}}{\mathit{Dt}}& =-\frac{\partial q_x}{\partial x}-\frac{\partial q_y}{\partial y}-\frac{\partial q_z}{\partial z}+\frac{\partial \left(u{\sigma }_{\mathit{xx}}\right)}{\partial x}+\frac{\partial \left(v{\sigma }_{\mathit{yy}}\right)}{\partial y}\frac{\partial \left(w{\sigma }_{\mathit{zz}}\right)}{\partial z}\\ & +\frac{\partial \left(u{\tau }_{\mathit{yx}}\right)}{\partial y}+\frac{\partial \left(u{\tau }_{\mathit{zx}}\right)}{\partial z}+\frac{\partial \left(v{\tau }_{\mathit{xy}}\right)}{\partial x}+\frac{\partial \left(v{\tau }_{\mathit{zy}}\right)}{\partial z}+\frac{\partial \left(w{\tau }_{\mathit{xz}}\right)}{\partial x}+\frac{\partial \left(w{\tau }_{\mathit{yz}}\right)}{\partial y}\end{array}$$E30

The heat flux q normally follows Fourier’s law of heat conduction which is given by [14]:

$$q_x=-\lambda \frac{\partial T}{\partial x}$$E31

$$q_y=-\lambda \frac{\partial T}{\partial y}$$E32

$$q_z=-\lambda \frac{\partial T}{\partial z}$$E33

Then, the energy balance equation (Eq. (30)) becomes Eq. (34)

$$\rho \frac{\mathit{DE}}{\mathit{Dt}}=\lambda \frac{\partial }{\partial x}\left(\frac{\partial T}{\partial x}\right)+\lambda \frac{\partial }{\partial y}\left(\frac{\partial T}{\partial y}\right)+\lambda \frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)+\mathrm{\Phi }$$E34

where Φ is the dissipation function which includes all terms relating to the energy due to deformation work done on the control volume. In most practical cases, dissipation function Φ can be neglected, because the term is significantly small compared with heat conduction. Then, Eq. (34) is simplified as

$$\rho \frac{\mathit{DE}}{\mathit{Dt}}=\lambda \frac{\partial }{\partial x}\left(\frac{\partial T}{\partial x}\right)+\lambda \frac{\partial }{\partial y}\left(\frac{\partial T}{\partial y}\right)+\lambda \frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)$$E35

If the energy E only refers to internal energy, the change of energy DE could be represented by C_pDT. Then, Eq. (35) could be further rearranged to Eq. (36)

$$\rho C_{\mathrm{p}}\frac{\mathit{DT}}{\mathit{Dt}}=\lambda \frac{\partial }{\partial x}\left(\frac{\partial T}{\partial x}\right)+\lambda \frac{\partial }{\partial y}\left(\frac{\partial T}{\partial y}\right)+\lambda \frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)$$E36

Expanding the total derivative of T, we get

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}=\frac{\lambda }{\rho C_p}\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)$$E37

Eq. (37) is the final expression of the energy equation. The first term on the left-hand side is called the transient term, the second to the fourth terms are called advection terms, and the terms on the right-hand sides are diffusion terms. Energy equation related temperature to time and space coordinate. With the energy equation, obtaining the temperature evolution and distribution in a fluid field becomes possible.

2.5 Summary of the governing equations

Previously, we have derived the continuity equation, component equation, momentum equation and energy equation based on the conservation principles.

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\mathit{\rho u}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho v}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho w}\right)}{\partial z}=0$$E38

$$\frac{\partial \left(\rho X_i\right)}{\partial t}+\frac{\partial \left(\rho {\mathit{uX}}_i\right)}{\partial x}+\frac{\partial \left(\rho {\mathit{vX}}_i\right)}{\partial y}+\frac{\partial \left(\rho {\mathit{wX}}_i\right)}{\partial z}=D\left[\frac{{\partial }^2\left(\rho X_i\right)}{\partial x^2}+\frac{{\partial }^2\left(\rho X_i\right)}{\partial y^2}+\frac{{\partial }^2\left(\rho X_i\right)}{\partial z^2}\right]$$E39

$$\frac{\partial \left(\mathit{\rho u}\right)}{\partial t}+\frac{\partial \left(\mathit{\rho uu}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho vu}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho wu}\right)}{\partial z}=-\frac{\partial p}{\partial x}+\mu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)+g_x$$E40

$$\frac{\partial \left(\mathit{\rho v}\right)}{\partial t}+\frac{\partial \left(\mathit{\rho uv}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho vv}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho wv}\right)}{\partial z}=-\frac{\partial p}{\partial y}+\mu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}+\frac{{\partial }^{2}v}{\partial z^2}\right)+g_y$$E41

$$\frac{\partial \left(\mathit{\rho w}\right)}{\partial t}+\frac{\partial \left(\mathit{\rho uw}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho vw}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho ww}\right)}{\partial z}=-\frac{\partial p}{\partial z}+\mu \left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}+\frac{{\partial }^{2}w}{\partial z^2}\right)+g_z$$E42

$$\frac{\partial \left(\mathit{\rho T}\right)}{\partial t}+\frac{\partial \left(\mathit{\rho uT}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho vT}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho wT}\right)}{\partial z}=\frac{\lambda }{C_p}\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)$$E43

If we introduce a universal variable ϕ, a generic form of the governing equations could be written as follows,

$$\frac{\partial \left(\mathit{\rho \varphi }\right)}{\partial t}+\frac{\partial \left(\mathit{\rho u\varphi }\right)}{\partial x}+\frac{\partial \left(\mathit{\rho v\varphi }\right)}{\partial y}+\frac{\partial \left(\mathit{\rho w\varphi }\right)}{\partial z}=\mathrm{\Gamma }\left(\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}\right)$$E44

where Γ is the diffusion coefficient. If ϕ = 1, Eq. (44) is the continuity equation. If ϕ = X_i, and Γ = D, Eq. (44) transforms to the component equation. If ϕ = u, v or w, Γ = μ, and the pressure term and gravity term are represented by source terms, Eq. (44) becomes the Navier-Stokes equation. If ϕ = T, and Γ = λ/C_p, Eq. (44) turns out to be the energy equation.

2.6 Governing equations in cylindrical coordinate

In some common cases, fluid flows through pipes or reactors which are mostly in cylindrical shape, so it is more convenient to express the flow in a cylindrical coordinate. Writing CFD code based on the governing equations expressed in cylindrical coordinate is also more suitable than those in Cartesian coordinate [10].

The continuity equation could be derived from analyzing the mass balance in the fluid microelement as shown in Figure 8. The blue arrows represent the flow in radial direction. The red color arrows represent the flow in the tangential direction, and the yellow arrows stand for the flow in the axial direction. Therefore, the conservation of mass in the control volume in Figure 8 is given as

Figure 8.

The schematic of mass flux of a control volume in cylindrical coordinate.

$$\begin{array}{ll}\frac{\partial \left(\mathit{\rho r}\mathrm{\Delta }r\mathrm{\Delta }\theta \mathrm{\Delta }z\right)}{\partial t}& =u_r\left(r\mathrm{\Delta }\theta \mathrm{\Delta }z\right)-\left(u_r+\frac{\partial u_r}{\partial r}\mathrm{\Delta }r\right)\left(r+\mathrm{\Delta }r\right)\mathrm{\Delta }\theta \mathrm{\Delta }z\\ & +u_{\theta }\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)-\left(u_{\theta }+\frac{\partial u_{\theta }}{\partial \theta }\mathrm{\Delta }\theta \right)\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)\\ & +u_z\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]-\left(u_z+\frac{\partial u_z}{\partial z}\mathrm{\Delta }z\right)\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]\end{array}$$E45

Dividing r∆r∆θ∆z on both sides, and taking the limit for ∆r → 0, ∆θ → 0 and ∆z → 0, Eq. (45) can be simplified to Eq. (46), which is the continuity equation in cylindrical coordinate.

$$\frac{\partial \rho }{\partial t}+\frac{1}{r}\frac{\partial \left({\mathit{ru}}_r\right)}{\partial r}+\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{\partial u_z}{\partial z}=0$$E46

The momentum equation in cylindrical coordinate is also derived from Newton’s second law.

$$\rho \overrightarrow{\mathbf{a}}=\frac{\overrightarrow{\mathbf{F}}}{V}=\frac{\overrightarrow{{\mathbf{F}}_{\mathbf{\text{surf}}}}+\overrightarrow{{\mathbf{F}}_{\mathbf{\text{body}}}}}{V}=\frac{\overrightarrow{{\mathbf{F}}_{\mathbf{\text{surf}}}}}{V}+\rho \overrightarrow{\mathbf{g}}$$E47

which is then expressed by velocity vector as

$$\frac{\mathrm{D}\overrightarrow{\mathbf{u}}}{\mathrm{D}t}=\frac{1}{\rho }\frac{\overrightarrow{{\mathbf{F}}_{\mathbf{\text{surf}}}}}{V}+\overrightarrow{\mathbf{g}}$$E48

In cylindrical coordinate, the velocity vector is given as

$$\overrightarrow{\mathbf{u}}\left(r,\theta ,z,t\right)=u_r\left(r,\theta ,z,t\right)\overrightarrow{{\mathbf{e}}_r}+u_{\theta }\left(r,\theta ,z,t\right)\overrightarrow{{\mathbf{e}}_{\theta }}+u_z\left(r,\theta ,z,t\right)\overrightarrow{{\mathbf{e}}_z}$$E49

And the acceleration is then given by

$$\overrightarrow{\mathbf{a}}=\frac{\mathrm{D}\overrightarrow{\mathbf{u}}}{\mathrm{D}t}=\frac{\partial \overrightarrow{\mathbf{u}}}{\partial t}+\frac{\partial \overrightarrow{\mathbf{u}}}{\partial r}\frac{\mathrm{d}r}{\mathrm{d}t}+\frac{\partial \overrightarrow{\mathbf{u}}}{\partial \theta }\frac{\mathrm{d}\theta }{\mathrm{d}t}+\frac{\partial \overrightarrow{\mathbf{u}}}{\partial z}\frac{\mathrm{d}z}{\mathrm{d}t}$$E50

The partial derivative of velocity vector with respect to time can be expressed as the summation of three vectors in radial ®, tangential (θ) and axial (z) directions (Figure 9) and the definition of u_r, u_θ and u_z is.

Figure 9.

The schematic of ∂u/∂t in cylindrical coordinate.

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial t}=\frac{\partial u_r}{\partial t}\overrightarrow{{\mathbf{e}}_r}+\frac{\partial u_{\theta }}{\partial t}\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial t}\overrightarrow{{\mathbf{e}}_z}$$E51

$$u_r=\frac{\mathrm{d}r}{\mathrm{d}t}$$E52

$$u_{\theta }=\frac{r\mathrm{d}\theta }{\mathrm{d}t}$$E53

$$u_z=\frac{\mathrm{d}z}{\mathrm{d}t}$$E54

With the definitions in Eq. (52), Eq. (50) is then transformed to

$$\frac{\mathrm{D}\overrightarrow{\mathbf{u}}}{\mathrm{D}t}=\frac{\partial \overrightarrow{\mathbf{u}}}{\partial t}+u_r\frac{\partial \overrightarrow{\mathbf{u}}}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial \overrightarrow{\mathbf{u}}}{\partial \theta }+u_z\frac{\partial \overrightarrow{\mathbf{u}}}{\partial z}$$E55

The partial derivatives of velocity vector with respect to r, θ and z are

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial r}=\frac{\partial \left(u_r\overrightarrow{{\mathbf{e}}_r}+u_{\theta }\overrightarrow{{\mathbf{e}}_{\theta }}+u_z\overrightarrow{{\mathbf{e}}_z}\right)}{\partial r}=\frac{\partial u_r}{\partial r}\overrightarrow{{\mathbf{e}}_r}+u_r\frac{\partial \overrightarrow{{\mathbf{e}}_r}}{\partial r}+\frac{\partial u_{\theta }}{\partial r}\overrightarrow{{\mathbf{e}}_{\theta }}+u_{\theta }\frac{\partial \overrightarrow{{\mathbf{e}}_{\theta }}}{\partial r}+\frac{\partial u_z}{\partial r}\overrightarrow{{\mathbf{e}}_z}+u_z\frac{\partial \overrightarrow{{\mathbf{e}}_z}}{\partial r}$$E56

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial \theta }=\frac{\partial \left(u_r\overrightarrow{{\mathbf{e}}_r}+u_{\theta }\overrightarrow{{\mathbf{e}}_{\theta }}+u_z\overrightarrow{{\mathbf{e}}_z}\right)}{\partial \theta }=\frac{\partial u_r}{\partial \theta }\overrightarrow{{\mathbf{e}}_r}+u_r\frac{\partial \overrightarrow{{\mathbf{e}}_r}}{\partial \theta }+\frac{\partial u_{\theta }}{\partial \theta }\overrightarrow{{\mathbf{e}}_{\theta }}+u_{\theta }\frac{\partial \overrightarrow{{\mathbf{e}}_{\theta }}}{\partial \theta }+\frac{\partial u_z}{\partial \theta }\overrightarrow{{\mathbf{e}}_z}+u_z\frac{\partial \overrightarrow{{\mathbf{e}}_z}}{\partial \theta }$$E57

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial z}=\frac{\partial \left(u_r\overrightarrow{{\mathbf{e}}_r}+u_{\theta }\overrightarrow{{\mathbf{e}}_{\theta }}+u_z\overrightarrow{{\mathbf{e}}_z}\right)}{\partial z}=\frac{\partial u_r}{\partial z}\overrightarrow{{\mathbf{e}}_r}+u_r\frac{\partial \overrightarrow{{\mathbf{e}}_r}}{\partial z}+\frac{\partial u_{\theta }}{\partial z}\overrightarrow{{\mathbf{e}}_{\theta }}+u_{\theta }\frac{\partial \overrightarrow{{\mathbf{e}}_{\theta }}}{\partial z}+\frac{\partial u_z}{\partial z}\overrightarrow{{\mathbf{e}}_z}+u_z\frac{\partial \overrightarrow{{\mathbf{e}}_z}}{\partial z}$$E58

Particularly, in cylindrical coordinate the unit vectors in radial (r), tangential (θ) and axial (z) directions have the following correlations

$$\frac{\partial \overrightarrow{{\mathbf{e}}_r}}{\partial r}=\frac{\partial \overrightarrow{{\mathbf{e}}_{\theta }}}{\partial r}=\frac{\partial \overrightarrow{{\mathbf{e}}_z}}{\partial r}=0$$E59

$$\frac{\partial \overrightarrow{{\mathbf{e}}_r}}{\partial \theta }=\overrightarrow{{\mathbf{e}}_{\theta }}$$E60

$$\frac{\partial \overrightarrow{{\mathbf{e}}_{\theta }}}{\partial \theta }=-\overrightarrow{{\mathbf{e}}_r}$$E61

Employing the correlations in Eq. (59), Eq. (56) would become

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial r}=\frac{\partial u_r}{\partial r}\overrightarrow{{\mathbf{e}}_r}+\frac{\partial u_{\theta }}{\partial r}\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial r}\overrightarrow{{\mathbf{e}}_z}$$E62

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial \theta }=\left(\frac{\partial u_r}{\partial \theta }-u_{\theta }\right)\overrightarrow{{\mathbf{e}}_r}+\left(u_r+\frac{\partial u_{\theta }}{\partial \theta }\right)\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial \theta }\overrightarrow{{\mathbf{e}}_z}$$E63

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial z}=\frac{\partial u_r}{\partial z}\overrightarrow{{\mathbf{e}}_r}+\frac{\partial u_{\theta }}{\partial z}\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial z}\overrightarrow{{\mathbf{e}}_z}$$E64

Bringing Eqs. (51) and (62) to Eq. (55), we can obtain

$$\begin{array}{ll}\frac{\mathrm{D}\overrightarrow{\mathbf{u}}}{\mathrm{D}t}& =\frac{\partial \overrightarrow{\mathbf{u}}}{\partial t}+\frac{\partial \overrightarrow{\mathbf{u}}}{\partial r}\frac{\mathrm{d}r}{\mathrm{d}t}+\frac{\partial \overrightarrow{\mathbf{u}}}{\partial \theta }\frac{\mathrm{d}\theta }{\mathrm{d}t}+\frac{\partial \overrightarrow{\mathbf{u}}}{\partial z}\frac{\mathrm{d}z}{\mathrm{d}t}\\ & =\frac{\partial u_r}{\partial t}\overrightarrow{{\mathbf{e}}_r}+\frac{\partial u_{\theta }}{\partial t}\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial t}\overrightarrow{{\mathbf{e}}_z}\\ & +u_r\left[\frac{\partial u_r}{\partial r}\overrightarrow{{\mathbf{e}}_r}+\frac{\partial u_{\theta }}{\partial r}\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial r}\overrightarrow{{\mathbf{e}}_z}\right]\\ & +\frac{u_{\theta }}{r}\left[\left(\frac{\partial u_r}{\partial \theta }-u_{\theta }\right)\overrightarrow{{\mathbf{e}}_r}+\left(u_r+\frac{\partial u_{\theta }}{\partial \theta }\right)\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial \theta }\overrightarrow{{\mathbf{e}}_z}\right]\\ & +u_z\left[\frac{\partial u_r}{\partial z}\overrightarrow{{\mathbf{e}}_r}+\frac{\partial u_{\theta }}{\partial z}\overrightarrow{{\mathbf{e}}_{\theta }}+\frac{\partial u_z}{\partial z}\overrightarrow{{\mathbf{e}}_z}\right]\\ & =\left(\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_r}{\partial \theta }-\frac{u_{\theta }^2}{r}+u_z\frac{\partial u_r}{\partial z}\right)\overrightarrow{{\mathbf{e}}_r}\\ & +\left(\frac{\partial u_{\theta }}{\partial t}+u_r\frac{\partial u_{\theta }}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_{\theta }{u}_r}{r}+u_z\frac{\partial u_{\theta }}{\partial z}\right)\overrightarrow{{\mathbf{e}}_{\theta }}\\ & +\left(\frac{\partial u_z}{\partial t}+u_r\frac{\partial u_z}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_z}{\partial \theta }+u_z\frac{\partial u_z}{\partial z}\right)\overrightarrow{{\mathbf{e}}_z}\end{array}$$E65

Now, we took the radial direction as an example to derive its momentum equation. The acceleration in the radial direction is defined as

$$a_r=\frac{\mathrm{D}{u}_r}{\mathrm{D}t}=\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_r}{\partial \theta }-\frac{u_{\theta }^2}{r}$$E66

According to Newton’s law, the acceleration also equals to

$$a_r=\frac{1}{\rho }\frac{F_r}{\mathit{\delta V}}+g_r$$E67

The surface tensions in the radial direction are described in Fig. 10. The summation of all the surface tensions along the radial direction is

Figure 10.

The schematic of the surface tensions along radial direction on a control volume in cylindrical coordinate.

$$\begin{array}{ll}{F}_r& ={\sigma }_{\mathit{rr}}\left(r+\mathrm{\Delta }r\right)\mathrm{\Delta }\theta \mathrm{\Delta }z+\frac{\partial {\sigma }_{\mathit{rr}}}{\partial r}\mathrm{\Delta }r\left(r+\mathrm{\Delta }r\right)\mathrm{\Delta }\theta \mathrm{\Delta }z-{\sigma }_{\mathit{rr}}\left(r\mathrm{\Delta }\theta \mathrm{\Delta }z\right)\\ & +{\sigma }_{\mathit{\theta r}}\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)cos\frac{\mathrm{\Delta }\theta }{2}+\frac{\partial {\sigma }_{\mathit{\theta r}}}{\partial \theta }\mathrm{\Delta }\theta \left(\mathrm{\Delta }r\mathrm{\Delta }z\right)cos\frac{\mathrm{\Delta }\theta }{2}-{\sigma }_{\mathit{\theta r}}\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)cos\frac{\mathrm{\Delta }\theta }{2}\\ & +{\sigma }_{\mathit{zr}}\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]+\frac{\partial {\sigma }_{\mathit{zr}}}{\partial z}\mathrm{\Delta }z\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]-{\sigma }_{\mathit{\text{rzrz}}}\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]\\ & -{\sigma }_{\mathit{\theta \theta }}\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)sin\frac{\mathrm{\Delta }\theta }{2}-\frac{\partial {\sigma }_{\mathit{\theta \theta }}}{\partial \theta }\mathrm{\Delta }\theta \left(\mathrm{\Delta }r\mathrm{\Delta }z\right)sin\frac{\mathrm{\Delta }\theta }{2}-{\sigma }_{\mathit{\theta \theta }}\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)sin\frac{\mathrm{\Delta }\theta }{2}\end{array}$$E68

Taking the limit for ∆r → 0, ∆θ → 0 and ∆z → 0, Eq. (68) is then achieved to Eq. (69) which is the continuity equation in cylindrical coordinate.

$$\begin{array}{ll}{F}_r& ={\sigma }_{\mathit{rr}}\mathrm{\Delta }r\mathrm{\Delta }\theta \mathrm{\Delta }z+\frac{\partial {\sigma }_{\mathit{rr}}}{\partial r}\mathrm{\Delta }r\left(r+\mathrm{\Delta }r\right)\mathrm{\Delta }\theta \mathrm{\Delta }z\\ & +\frac{\partial {\sigma }_{\mathit{\theta r}}}{\partial \theta }\mathrm{\Delta }\theta \left(\mathrm{\Delta }r\mathrm{\Delta }z\right)cos\frac{\mathrm{\Delta }\theta }{2}\\ & +\frac{\partial {\sigma }_{\mathit{zr}}}{\partial z}\mathrm{\Delta }z\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]\\ & -2{\sigma }_{\mathit{\theta \theta }}\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)\frac{\mathrm{\Delta }\theta }{2}-\frac{\partial {\sigma }_{\mathit{\theta \theta }}}{\partial \theta }\mathrm{\Delta }\theta \left(\mathrm{\Delta }r\mathrm{\Delta }z\right)\frac{\mathrm{\Delta }\theta }{2}\end{array}$$E69

Bringing ∆V = r∆r∆θ∆z into Eq. (69) gives

$$\frac{F_r}{\mathrm{\Delta }V}=\frac{{\sigma }_{\mathit{rr}}}{r}+\frac{\partial {\sigma }_{\mathit{rr}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\mathit{\theta r}}}{\partial \theta }+\frac{\partial {\sigma }_{\mathit{zr}}}{\partial z}-\frac{{\sigma }_{\mathit{\theta \theta }}}{r}$$E70

Combining Eqs. (66), (67) and (70), we can obtain

$$\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_r}{\partial \theta }-\frac{u_{\theta }^2}{r}+u_z\frac{\partial u_r}{\partial z}=\frac{1}{\rho }\left(\frac{{\sigma }_{\mathit{rr}}}{r}+\frac{\partial {\sigma }_{\mathit{rr}}}{\partial r}+\frac{1}{r}\frac{\partial {\sigma }_{\mathit{\theta r}}}{\partial \theta }+\frac{\partial {\sigma }_{\mathit{zr}}}{\partial z}-\frac{{\sigma }_{\mathit{\theta \theta }}}{r}\right)+g_r$$E71

In fluid mechanics, the definitions of the surface tensions in Eq. (70) are given as follows:

$${\sigma }_{\mathit{rr}}=-p+2\mu \frac{\partial u_r}{\partial r}$$E72

$${\sigma }_{\mathit{\theta \theta }}=-p+2\mu \left(\frac{1}{r}\frac{\partial u_{\theta }}{\partial \theta }+\frac{u_r}{r}\right)$$E73

$${\sigma }_{\mathit{zz}}=-p+2\mu \frac{\partial u_z}{\partial z}$$E74

$${\sigma }_{\mathit{r\theta }}={\sigma }_{\mathit{\theta r}}=\mu \left(\frac{1}{r}\frac{\partial u_r}{\partial \theta }+\frac{\partial u_{\theta }}{\partial r}-\frac{u_{\theta }}{r}\right)$$E75

$${\sigma }_{\mathit{rz}}={\sigma }_{\mathit{zr}}=\mu \left(\frac{\partial u_r}{\partial z}+\frac{\partial u_z}{\partial r}\right)$$E76

$${\sigma }_{\mathit{\theta z}}={\sigma }_{\mathit{z\theta }}=\mu \left(\frac{1}{r}\frac{\partial u_z}{\partial \theta }+\frac{\partial u_{\theta }}{\partial z}\right)$$E77

Combining Eqs. (70) and (71) yields

$$\begin{array}{l}\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_r}{\partial \theta }-\frac{u_{\theta }^2}{r}+u_z\frac{\partial u_r}{\partial z}\\ =-\frac{1}{\rho }\frac{\partial p}{\partial r}+\frac{\mu }{\rho }\left[-\frac{u_r}{r^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_r}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2{u}_r}{\partial {\theta }^2}+\frac{{\partial }^2{u}_r}{\partial z^2}-\frac{2}{r^2}\frac{\partial u_{\theta }}{\partial \theta }\right]+g_r\end{array}$$E78

Finally, Eq. (78) is the momentum equation of fluid in the radial direction. Likewise, in the tangential and axial directions, the momentum equations of fluid flow could also be derived in a similar way as follows:

$$\begin{array}{l}\frac{\partial u_{\theta }}{\partial t}+u_r\frac{\partial u_{\theta }}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_{\theta }}{\partial \theta }-\frac{u_r{u}_{\theta }}{r}+u_z\frac{\partial u_{\theta }}{\partial z}\\ =-\frac{1}{\rho }\frac{\partial p}{\partial \theta }+\frac{\mu }{\rho }\left[-\frac{u_{\theta }}{r^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{\theta }}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2{u}_{\theta }}{\partial {\theta }^2}+\frac{{\partial }^2{u}_{\theta }}{\partial z^2}+\frac{2}{r^2}\frac{\partial u_r}{\partial \theta }\right]+g_{\theta }\end{array}$$E79

$$\begin{array}{l}\frac{\partial u_z}{\partial t}+u_r\frac{\partial u_z}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_z}{\partial \theta }+u_z\frac{\partial u_z}{\partial z}\\ =-\frac{1}{\rho }\frac{\partial p}{\partial z}+\frac{\mu }{\rho }\left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_z}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2{u}_z}{\partial {\theta }^2}+\frac{{\partial }^2{u}_z}{\partial z^2}\right]+g_z\end{array}$$E80

The energy equation could be derived from analyzing the energy balance in the fluid control volume as shown in Figure 11.

Figure 11.

The schematic of the conductive heat fluxes in a control volume in cylindrical coordinate.

According to Figure 11 for cylindrical coordinate, the heat energy balance is expressed as the energy accumulated in the control volume equals the heat into this control volume minus the heat out of the control volume

$$\begin{array}{ll}\rho \frac{\mathit{DE}}{\mathit{Dt}}r\mathrm{\Delta }r\mathrm{\Delta }\theta \mathrm{\Delta }z& =-\left(q_r+\frac{\partial q_r}{\partial r}\mathrm{\Delta }r\right)\left(r+\mathrm{\Delta }r\right)\mathrm{\Delta }\theta \mathrm{\Delta }z+q_r\left(r\mathrm{\Delta }\theta \mathrm{\Delta }z\right)\\ & -\left(q_{\theta }+\frac{\partial q_{\theta }}{\partial \theta }\mathrm{\Delta }\theta \right)\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)cos\frac{\mathrm{\Delta }\theta }{2}+q_{\theta }\left(\mathrm{\Delta }r\mathrm{\Delta }z\right)cos\frac{\mathrm{\Delta }\theta }{2}\\ & -\left(q_z+\frac{\partial q_z}{\partial z}\mathrm{\Delta }z\right)\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]+q_z\left[\left(r+\frac{\mathrm{\Delta }r}{2}\right)\mathrm{\Delta }r\mathrm{\Delta }\theta \right]\end{array}$$E81

Which can be further transformed to Eq. (82) by canceling out identical terms with opposite signs.

$$\rho \frac{\mathit{DE}}{\mathit{Dt}}=-\frac{1}{r}\frac{\partial \left({\mathit{rq}}_r\right)}{\partial r}-\frac{1}{r}\frac{\partial q_{\theta }}{\partial \theta }-\frac{\partial q_z}{\partial z}$$E82

If the change of thermal energy is all induced by temperature change, then Eq. (82) is expressed as:

$$\rho C_{\mathrm{p}}\left(\frac{\partial T}{\partial t}+u_r\frac{\partial T}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial T}{\partial \theta }+u_z\frac{\partial T}{\partial z}\right)=-\frac{1}{r}\frac{\partial \left({\mathit{rq}}_r\right)}{\partial r}-\frac{1}{r}\frac{\partial q_{\theta }}{\partial \theta }-\frac{\partial q_z}{\partial z}$$E83

By applying Fourier’s law of heat conduction (Eq. 84), Eq. (89) is then derived as the energy equation in cylindrical coordinate.

$$q_r=-\lambda \frac{\partial T}{\partial r}$$E84

$$q_{\theta }=-\frac{\lambda }{r}\frac{\partial T}{\partial \theta }$$E85

$$q_z=-\lambda \frac{\partial T}{\partial z}$$E86

$$\begin{array}{l}\rho C_{\mathrm{p}}\left(\frac{\partial T}{\partial t}+u_r\frac{\partial T}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial T}{\partial \theta }+u_z\frac{\partial T}{\partial z}\right)\\ =\lambda \left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)+\frac{1}{r^2}\frac{\partial }{\partial \theta }\left(\frac{\partial T}{\partial \theta }\right)+\frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)\right]\end{array}$$E87

3. Source and sinks

The expansion of CFD applications to the fields such as chemical engineering, environmental engineering, thermal engineering and pharmaceutical engineering was based on the development of multiscale and multiphase models which benefited from the coupling of source or sink terms to the conservation equations. A sink term is actually a negative source term, so both source terms and sink terms could be called by a joint name—source terms.

3.1 Source term for the continuity equation

Processes such as solid gasification or combustion, and gas adsorption by solid sorbents would break the mass conservation of the fluid phase. If there are reactions are generating additional fluid or remove the existing fluid in a control volume, the mass balance in this control volume must be disrupted by these reactions. The rate of fluid entering a control volume minus the rate of fluid flowing out equals the rate of fluid accumulating in this control volume and the net gain of fluid generated inside the control volume. If in the control volume there is a reaction generating fluid at a rate of r, Eq. (1) for mass balance will be altered to

$$\frac{\mathrm{d}m}{\mathrm{d}t}=r\mathrm{d}x\mathrm{d}y\mathrm{d}z+\sum _{\text{in}}\dot{m}-\sum _{\text{out}}\dot{m}$$E88

Eq. (45) could also include the case where there is a depletion of fluid by using a negative value for r. Therefore, the continuity equation is adjusted by a source term as follows.

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\mathit{\rho u}\right)}{\partial x}+\frac{\partial \left(\mathit{\rho v}\right)}{\partial y}+\frac{\partial \left(\mathit{\rho w}\right)}{\partial z}=S_{\text{continuity}}$$E89

where the source term is the reaction rate S_continuity = r (kg m⁻³ s⁻¹) if the reaction is a body reaction. For a surface reaction, the source term is the multiplication of the surface reaction rate (kg m⁻² s⁻¹) with the reaction area in unit volume (A_react/V_body, m² m⁻³) [10].

3.2 Source term for component equation

Similar to the imbalance for total mass conservation, the conservation of each component in a multi-component flow with reactions could also be broken. The components participating as reactants would be deemed as sinks and the components generated as products would be regarded as sources. If a component is generated at a rate of r_i, then Eq. (7) is adjusted accordingly as follows

$$\frac{\mathrm{d}{m}_i}{\mathrm{d}t}=r_i\mathrm{d}x\mathrm{d}y\mathrm{d}z+\sum _{\text{in}}{\dot{m}}_i-\sum _{\text{out}}{\dot{m}}_i$$E90

After a similar derivation from Eq. (7) to Eq. (9), the component equation should include a source term as expressed by Eq. (91)

$$\begin{array}{c}\frac{\partial \left(\rho X_i\right)}{\partial t}+\frac{\partial \left(\rho {\mathit{uX}}_i\right)}{\partial x}+\frac{\partial \left(\rho {\mathit{vX}}_i\right)}{\partial y}+\frac{\partial \left(\rho {\mathit{wX}}_i\right)}{\partial z}\\ =D\left[\frac{{\partial }^2\left(\rho X_i\right)}{\partial x^2}+\frac{{\partial }^2\left(\rho X_i\right)}{\partial y^2}+\frac{{\partial }^2\left(\rho X_i\right)}{\partial z^2}\right]+S_i\end{array}$$E91

where S_i is the source term of the component equation, and it equals reaction rate r_i if it is a body reaction. If it is a surface reaction, the source term of Eq. (91) equals the multiplication of the surface reaction rate (kg m⁻² s⁻¹) with the reaction area in unit volume (A_react/V_body, m² m⁻³).

3.3 Source term for the momentum equation

When it comes to cases with a multiphase flow such as gas-solid flow in fluidized bed, momentum conservation still holds, but interphase interactions need to be taken into account [15]. There might be extra forces between different phases, which can affect the momentum transfer between different phases. Therefore, in multiphase flow, appropriate interphase force terms need to be introduced to correct the momentum equation.

When the fluid phase and solid phase interact with each other, a reactive force that causes flow resistance is generated, known as drag force. If the drag force is expressed as F_D,i, Eq. (13) can be expanded as follows

$$\frac{\mathrm{D}u}{\mathrm{D}t}=\frac{\sum F_x}{m}=\frac{\sum F_{\text{surf}}+\sum F_{\text{body}}+\sum F_{\mathrm{D},i}}{m}$$E92

The drag force F_D,I is calculated by:

$${\mathit{F}}_{\mathrm{D},i}=\frac{{\beta }_{\mathrm{fs}}}{{\theta }_{\mathrm{f}}}\left({\mathit{u}}_{\mathrm{s}}-{\mathit{u}}_{\mathrm{f}}\right)$$E93

where u_s is the solid phase velocity, u_f is the fluid phase velocity, θ_f is the volume fraction of the fluid phase and β_fs is the drag coefficient. Many researchers have proposed the models for drag coefficients. The most widely used is the Gidaspow model which combines Ergun and Wen-Yu equations to accurately simulate the gas-solid multiphase flow:

$${\beta }_{\mathrm{p},\mathrm{f}}=\left\{\begin{array}{cc}150\frac{{\theta }_{\mathrm{s}}^2{\mu }_{\mathrm{f}}}{{\theta }_{\mathrm{f}}^2{d}_{\mathrm{s}}}+1.75\frac{{\theta }_{\mathrm{s}}{\rho }_{\mathrm{f}}}{{\theta }_{\mathrm{f}}{d}_{\mathrm{s}}}\left|{\mathit{u}}_{\mathrm{s}}-{\mathit{u}}_{\mathrm{f}}\right|& {\theta }_{\mathrm{f}}<0.8\\ \frac{3}{4}{C}_{\mathrm{D}}\frac{{\theta }_{\mathrm{s}}{\rho }_{\mathrm{f}}}{d_{\mathrm{s}}}\left|{\mathit{u}}_{\mathrm{s}}-{\mathit{u}}_{\mathrm{f}}\right|{{\theta }_{\mathrm{f}}}^{-2.65}& {\theta }_{\mathrm{f}}\ge 0.8\end{array}\right.$$E94

$$C_{\mathrm{D}}=\left\{\begin{array}{l}\frac{24\left(1+0.15{\mathit{Re}}_{\mathrm{p}}^{0.687}\right)}{{\mathit{Re}}_{\mathrm{p}}}{\mathit{Re}}_{\mathrm{p}}\le 1000\\ 0.44{\mathit{Re}}_{\mathrm{p}}>1000\end{array}\right.$$E95

θ_s is the particle volume fraction, θ_f is the fluid volume fraction, μ_f is the fluid phase viscosity, ρ_f is the fluid density, d_s is the particle diameter of the solid phase and C_D is the drag coefficient. Particle Reynolds number Re_p is described by

$${\mathit{Re}}_{\mathrm{p}}={\theta }_{\mathrm{f}}{\rho }_{\mathrm{f}}{d}_{\mathrm{s}}\left|{\mathit{u}}_{\mathrm{s}}-{\mathit{u}}_{\mathrm{f}}\right|/{\mathit{u}}_{\mathrm{f}}$$E96

Therefore, Eq. (40), the momentum equation in the x direction, is adjusted by a source term as follows

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)+g_x+S_{\mathrm{mom},x}$$E97

where S_mom is the source term for momentum which can be expressed as:

$$S_{\mathrm{mom},x}=\sum F_{\mathrm{D},i,x}/\left(\rho \mathrm{d}x\mathrm{d}y\mathrm{d}z\right)$$E98

Similarly, the momentum equation in the y direction (Eq. 41) and the momentum equation in the z direction (Eq. 42) are adjusted as follows

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}+\frac{{\partial }^{2}v}{\partial z^2}\right)+g_y+S_{\mathrm{mom},y}$$E99

$$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\frac{\mu }{\rho }\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}+\frac{{\partial }^{2}w}{\partial z^2}\right)+g_z+S_{\mathrm{mom},z}$$E100

3.4 Source term for the energy equation

As discussed in the energy equation section, the energy equation ensures the thermal energy conservation. In some cases, chemical reactions occur in the fluid domain, and these reactions are either exothermic or endothermic. If the reaction is exothermic, it could release heat to the surrounding fluid and increase the temperature. If the reaction is endothermic, it consumes the heat from the surrounding fluid and decreases the temperature. Even though total energy is conserved, the thermal energy is not. In order to consider the imbalance of thermal energy in fluid, the source term of heat must be included in the energy equation. By considering reaction heat, the energy balance is

$$\rho \frac{\mathit{DE}}{\mathit{Dt}}\mathit{\text{ΔxΔyΔz}}=\sum \dot{Q}+r\mathrm{\Delta }\mathit{\text{hΔxΔyΔz}}$$E101

where ∆h (J kg⁻¹) is the heat generated when unit mass reactant is converted to products. Bringing Fourier’s law of heat conduction and dividing both sides by ∆x∆y∆z, the energy equation with source term S_energ = r∆h is

$$\rho \frac{\mathit{DE}}{\mathit{Dt}}=\lambda \frac{\partial }{\partial x}\left(\frac{\partial T}{\partial x}\right)+\lambda \frac{\partial }{\partial y}\left(\frac{\partial T}{\partial y}\right)+\lambda \frac{\partial }{\partial z}\left(\frac{\partial T}{\partial z}\right)+r\mathrm{\Delta }h$$E102

With the inclusion of source terms, processes with chemical reactions or multiphase interactions could be accurately simulated in CFD models.

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:

Figure 12.

One-dimensional grids for the finite-difference method.

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

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

$${\varphi }_{i-1}={\varphi }_i-{\left(\frac{\partial \varphi }{\partial x}\right)}_i\mathrm{\Delta }x+{\left(\frac{{\partial }^2\varphi }{\partial x^2}\right)}_i\frac{\mathrm{\Delta }{x}^2}{2}-{\left(\frac{{\partial }^3\varphi }{\partial x^3}\right)}_i\frac{\mathrm{\Delta }{x}^3}{6}+\left(\text{Higher Order Terms}\right)$$E104

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

$${\left(\frac{\partial \varphi }{\partial x}\right)}_i=\frac{{\varphi }_{i+1}-{\varphi }_{i-1}}{2\mathrm{\Delta }x}+\frac{\left(\mathrm{\Delta }{x}^3/3\right)\left({\partial }^3\varphi /\partial x^3\right)}{2\mathrm{\Delta }x}$$E105

Since ∆x is approaching zero, the second term is very small compared to the first term. By omitting the ∆x³ term, the central difference is achieved

$${\left(\frac{\partial \varphi }{\partial x}\right)}_i=\frac{{\varphi }_{i+1}-{\varphi }_{i-1}}{2\mathrm{\Delta }x}$$E106

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

$${\left(\frac{\partial \varphi }{\partial x}\right)}_i=\frac{{\varphi }_{i+1}-{\varphi }_i}{\mathrm{\Delta }x}+O\left(\mathrm{\Delta }x\right)$$E107

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

$${\left(\frac{\partial \varphi }{\partial x}\right)}_i=\frac{{\varphi }_i-{\varphi }_{i-1}}{\mathrm{\Delta }x}+O\left(\mathrm{\Delta }x\right)$$E108

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 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

$${\varphi }_{i+1}+{\varphi }_{i-1}=2{\varphi }_i+{\left(\frac{{\partial }^2\varphi }{\partial x^2}\right)}_i\mathrm{\Delta }{x}^2+{\left(\frac{{\partial }^4\varphi }{\partial x^4}\right)}_i\frac{\mathrm{\Delta }{x}^4}{24}$$E109

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

$$\left(\frac{{\partial }^2\varphi }{\partial x^2}\right)=\frac{{\varphi }_{i+1}-2{\varphi }_i+{\varphi }_{i-1}}{\mathrm{\Delta }{x}^2}$$E110

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:

Figure 14.

Three-dimensional grids for the finite-difference method.

$${\left(\frac{\partial \varphi }{\partial y}\right)}_j=\frac{{\varphi }_{j+1}-{\varphi }_{j-1}}{2\mathrm{\Delta }y}$$E111

$${\left(\frac{\partial \varphi }{\partial z}\right)}_k=\frac{{\varphi }_{k+1}-{\varphi }_{k-1}}{2\mathrm{\Delta }z}$$E112

$$\left(\frac{{\partial }^2\varphi }{\partial y^2}\right)=\frac{{\varphi }_{j+1}-2{\varphi }_j+{\varphi }_{j-1}}{\mathrm{\Delta }{y}^2}$$E113

$$\left(\frac{{\partial }^2\varphi }{\partial z^2}\right)=\frac{{\varphi }_{k+1}-2{\varphi }_k+{\varphi }_{k-1}}{\mathrm{\Delta }{z}^2}$$E114

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 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 \varphi }{\partial t}\right)=\frac{{{\varphi }_{i,j,k}}^n-{{\varphi }_{i,j,k}}^{n-1}}{\mathrm{\Delta }t}$$E115

$$\left(\frac{\partial \varphi }{\partial t}\right)=\frac{{{\varphi }_{i,j,k}}^{n+1}-{{\varphi }_{i,j,k}}^n}{\mathrm{\Delta }t}$$E116

$$\left(\frac{\partial \varphi }{\partial t}\right)=\frac{{{\varphi }_{i,j,k}}^{n+1}-{{\varphi }_{i,j,k}}^{n‐1}}{2\mathrm{\Delta }t}$$E117

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.

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{array}{ll}\underset{\mathrm{\Omega }}{\mathrm{\iiint }}\frac{\partial {\varphi }_x}{\partial x}+\frac{\partial {\varphi }_y}{\partial y}+\frac{\partial {\varphi }_z}{\partial z}\mathrm{d}V& ={∯}_{\mathrm{\Sigma }}\overrightarrow{\varphi }\left(x,y,z\right)\mathrm{d}S\\ & ={∯}_{\mathrm{\Sigma }}{\varphi }_x\mathrm{d}y\mathrm{d}z+{∯}_{\mathrm{\Sigma }}{\varphi }_y\mathrm{d}x\mathrm{d}z+{∯}_{\mathrm{\Sigma }}{\varphi }_z\mathrm{d}x\mathrm{d}y\end{array}$$E118

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

$$\underset{\mathrm{\Omega }}{\iiint }\frac{\partial \varphi }{\partial z}\mathrm{d}V=\underset{\mathrm{\Sigma }}{\iint }\varphi \left(x,y,z\right)\mathrm{d}x\mathrm{d}y$$E119

$$\underset{\mathrm{\Omega }}{\iiint }\frac{\partial \varphi }{\partial x}\mathrm{d}V=\underset{\mathrm{\Sigma }}{\iint }\varphi \left(x,y,z\right)\mathrm{d}y\mathrm{d}z$$E120

$$\underset{\mathrm{\Omega }}{\iiint }\frac{\partial \varphi }{\partial y}\mathrm{d}V=\underset{\mathrm{\Sigma }}{\iint }\varphi \left(x,y,z\right)\mathrm{d}x\mathrm{d}z$$E121

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

$$\underset{\mathrm{\Omega }}{\iiint }\frac{\partial \varphi }{\partial x}\mathrm{d}V=\frac{\partial \varphi }{\partial x}\underset{\mathrm{\Omega }}{\iiint }\mathrm{d}V=\frac{\partial \varphi }{\partial x}\mathrm{\Delta }V$$E122

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

$$\underset{\mathrm{\Sigma }}{\iint }\varphi \mathrm{d}y\mathrm{d}z=\varphi A^x$$E123

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

Figure 16.

The representation of a tetrahedral grid.

$$\begin{array}{ll}\frac{\partial \varphi }{\partial x}& =\frac{1}{\mathrm{\Delta }V}\underset{\mathrm{\Delta }V}{\int }\frac{\partial \varphi }{\partial x}\mathrm{d}V=\frac{1}{\mathrm{\Delta }V}\underset{A}{\int }\varphi \mathrm{d}{A}^x\\ & =\frac{1}{\mathrm{\Delta }V}\left({\varphi }_1{A}_1^x+{\varphi }_2{A}_2^x+{\varphi }_3{A}_3^x+{\varphi }_4{A}_4^x\right)\\ & =\frac{1}{\mathrm{\Delta }V}\sum _{i=1}^N{\varphi }_i{A}_i^x\end{array}$$E124

where A_i^x 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{array}{ll}\frac{\partial \varphi }{\partial y}& =\frac{1}{\mathrm{\Delta }V}\underset{\mathrm{\Delta }V}{\int }\frac{\partial \varphi }{\partial y}\mathrm{d}V=\frac{1}{\mathrm{\Delta }V}\underset{A}{\int }\varphi \mathrm{d}{A}^y\\ & =\frac{1}{\mathrm{\Delta }V}\left({\varphi }_1{A}_1^y+{\varphi }_2{A}_2^y+{\varphi }_3{A}_3^y+{\varphi }_4{A}_4^y\right)\\ & =\frac{1}{\mathrm{\Delta }V}\sum _{i=1}^N{\varphi }_i{A}_i^y\end{array}$$E125

$$\begin{array}{ll}\frac{\partial \varphi }{\partial z}& =\frac{1}{\mathrm{\Delta }V}\underset{\mathrm{\Delta }V}{\int }\frac{\partial \varphi }{\partial z}\mathrm{d}V=\frac{1}{\mathrm{\Delta }V}\underset{A}{\int }\varphi \mathrm{d}{A}^z\\ & =\frac{1}{\mathrm{\Delta }V}\left({\varphi }_1{A}_1^z+{\varphi }_2{A}_2^z+{\varphi }_3{A}_3^z+{\varphi }_4{A}_4^z\right)\\ & =\frac{1}{\mathrm{\Delta }V}\sum _{i=1}^N{\varphi }_i{A}_i^z\end{array}$$E126

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

$$\begin{array}{ll}\frac{{\partial }^2\varphi }{\partial x^2}& =\frac{1}{\mathrm{\Delta }V}\underset{\mathrm{\Delta }V}{\int }\frac{{\partial }^2\varphi }{\partial x^2}\mathrm{d}V=\frac{1}{\mathrm{\Delta }V}\underset{A}{\int }\frac{\partial \varphi }{\partial x}\mathrm{d}{A}^x\\ & =\frac{1}{\mathrm{\Delta }V}\left[{\left(\frac{\partial \varphi }{\partial x}\right)}_1{A}_1^x+{\left(\frac{\partial \varphi }{\partial x}\right)}_2{A}_2^x+{\left(\frac{\partial \varphi }{\partial x}\right)}_3{A}_3^x+{\left(\frac{\partial \varphi }{\partial x}\right)}_4{A}_4^x\right]\\ & =\frac{1}{\mathrm{\Delta }V}\sum _{i=1}^N{\left(\frac{\partial \varphi }{\partial x}\right)}_i{A}_i^x\end{array}$$E127

$$\frac{{\partial }^2\varphi }{\partial y^2}=\frac{1}{\mathrm{\Delta }V}\sum _{i=1}^N{\left(\frac{\partial \varphi }{\partial y}\right)}_i{A}_i^y$$E128

$$\frac{{\partial }^2\varphi }{\partial z^2}=\frac{1}{\mathrm{\Delta }V}\sum _{i=1}^N{\left(\frac{\partial \varphi }{\partial z}\right)}_i{A}_i^z$$E129

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.

5. Conclusion

Understanding the basic theory of CFD governing equations and numerical solutions is the key to developing CFD simulation codes as well as using CFD packages. Basically, all flows can be universally described by the continuity equation, component equation, momentum equation and energy equation, because fluid follows mass conservation, momentum conservation and energy conservation. With source terms, CFD applications were broadened to a wider field like chemical engineering, environmental engineering or thermal engineering. However, these equations are almost impossible to be solved analytically, thus we have to resort to the numerical method to get the approximate results. This involves using the finite difference method and finite volume method. The numerical methods offered a feasible way to get the solutions from the governing equations, which facilitated the analysis of complex phenomena associated with fluid flow, heat transfer and chemical reactions.

Acknowledgments

Authors appreciate the financial support from teaching reform fund for postgraduate students in Dalian University of Technology (No. 20), Liao Ning Revitalization Talents Program (grant number: XLYC2007179) and the Australian Research Council (grant number DE210101549).