**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).

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

$$\frac{\text{dm}}{\text{dt}} = \sum\_{\text{in}} \dot{m} - \sum\_{\text{out}} \dot{m} \tag{1}$$

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

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

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[\rho u + \frac{\partial(\rho u)}{\partial \mathbf{x}} \Delta x\right] \Delta y \Delta z + \left[\rho v + \frac{\partial(\rho v)}{\partial \mathbf{y}} \Delta y\right] \Delta x \Delta z + \left[\rho w + \frac{\partial(\rho w)}{\partial \mathbf{z}} \Delta x\right] \Delta x \Delta y \tag{3}$$

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

$$\begin{split} \frac{\partial(\rho \Delta x \Delta y \Delta z)}{\partial t} &= \rho u \Delta y \Delta z + \rho v \Delta x \Delta z + \rho w \Delta x \Delta y \\ &\quad - \left[ \rho u + \frac{\partial(\rho u)}{\partial x} \Delta x \right] \Delta y \Delta z \\ &\quad - \left[ \rho v + \frac{\partial(\rho v)}{\partial y} \Delta y \right] \Delta x \Delta z \\ &\quad - \left[ \rho w + \frac{\partial(\rho w)}{\partial x} \Delta z \right] \Delta x \Delta y \end{split} \tag{4}$$

**Figure 1.** *The schematic of mass conservation in a control volume.*

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(\rho u)}{\partial x} + \frac{\partial(\rho v)}{\partial y} + \frac{\partial(\rho w)}{\partial z} = \mathbf{0} \tag{5}$$

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) = \mathbf{0} \tag{6}$$

#### **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 subchamber, there are 8 moles of N2 and 2 moles of CO, while in the right chamber there are 2 moles of N2 and 8 moles of CO. Once the valve is open, there is no convective flow due to the equal pressure on both sides, but N2 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]:

$$J\_i = D \frac{\mathbf{d}(\rho \mathbf{X}\_i)}{\mathbf{d}\mathbf{x}} \tag{7}$$

where *Ji* 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 *Xi* is the mass fraction of component *I* in the fluid.

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

**Figure 2.** *Demonstration of diffusive flow.*

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

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

where *m*<sup>i</sup> is the mass of the *i*th component in a control volume. The *i*th 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 *Xi*, the flow of component *I* into the control volume together with the total convective flow is *XiρuΔyΔz*. The convective flow of *i*th 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 {*Xiρu +* [*∂*(*Xiρu*)/*∂x*]*Δx*}*ΔyΔz*. Likewise, the *i*th component entering and exiting along y and z directions are *XiρvΔxΔz*, *XiρwΔyΔx*, {*Xiρv +* [*∂*(*Xiρv*)/*∂y*]*Δy*}*ΔxΔz*, {*Xiρw +* [*∂*(*Xiρw*)/*∂z*]*Δz*}*ΔyΔx*, respectively.

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*[*∂*(*ρXi*)/*∂x*]*ΔyΔz*. Then, the diffusive flow exiting from the opposite surface is *D*{*∂*(*ρXi*)/*∂x* + [*∂*<sup>2</sup> (*ρXi*)/ *∂x*<sup>2</sup> ]*Δx*}*ΔyΔz*. Likewise, along the y direction, the diffusive flows entering and exiting are *D*[*∂*(*ρXi*)/*∂y*]*ΔxΔz* and *D*{*∂*(*ρXi*)/*∂y* + [*∂*<sup>2</sup> (*ρXi*)/*∂y* 2 ] *Δy*}*ΔxΔz*. Similarly, along the z direction, the diffusive flow entering and exiting are *D*[*∂*(*ρXi*)/*∂z*]*ΔxΔy* and *D*{*∂*(*ρXi*)/ *∂z* + [*∂*<sup>2</sup> (*ρXi*)/*∂z*<sup>2</sup> ] *Δz*}*ΔxΔy*.

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{split} \frac{\partial}{\partial t} (\rho X\_i \Delta x \Delta y \Delta z) &= D \frac{\partial^2 (\rho X\_i)}{\partial x^2} \Delta x \Delta y \Delta z + D \frac{\partial^2 (\rho X\_i)}{\partial y^2} \Delta y \Delta x \Delta z + D \frac{\partial^2 (\rho X\_i)}{\partial z^2} \Delta z \Delta x \Delta y \\ &- \frac{\partial (\rho u X\_i)}{\partial x} \Delta x \Delta y \Delta z - \frac{\partial (\rho v X\_i)}{\partial y} \Delta y \Delta x \Delta z - \frac{\partial (\rho w X\_i)}{\partial z} \Delta z \Delta x \Delta y \end{split} \tag{9}$$

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

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

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

$$\frac{\partial(\rho \mathbf{X}\_i)}{\partial t} + \frac{\partial(\rho u \mathbf{X}\_i)}{\partial \mathbf{x}} + \frac{\partial(\rho v \mathbf{X}\_i)}{\partial \mathbf{y}} + \frac{\partial(\rho w \mathbf{X}\_i)}{\partial \mathbf{z}} = D \left[ \frac{\partial^2(\rho \mathbf{X}\_i)}{\partial \mathbf{x}^2} + \frac{\partial^2(\rho \mathbf{X}\_i)}{\partial \mathbf{y}^2} + \frac{\partial^2(\rho \mathbf{X}\_i)}{\partial \mathbf{z}^2} \right] \tag{10}$$

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}\tag{11}
$$

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\_{\chi} = \frac{\mathrm{D}u}{\mathrm{D}t} \tag{12}$$

$$a\_{\gamma} = \frac{\mathbf{D}v}{\mathbf{D}t} \tag{13}$$

$$a\_x = \frac{\text{D}w}{\text{D}t} \tag{14}$$

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

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} \tag{15}$$

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.

$$\mathfrak{a}\_{\mathfrak{x}} = \frac{\sum F\_{\mathfrak{x}}}{m} \tag{16}$$

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.

The summation of all the surface forces is

$$
\sum F\_{\text{surf}} = \frac{\partial \sigma\_{\text{xx}}}{\partial \mathbf{x}} \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} + \frac{\partial \tau\_{\text{yx}}}{\partial \mathbf{y}} \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} + \frac{\partial \tau\_{\text{xx}}}{\partial \mathbf{z}} \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} \tag{17}
$$

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\_{\text{xx}} = -p + 2\mu \frac{\partial u}{\partial \mathbf{x}} \tag{18}
$$

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

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

$$
\tau\_{\rm xy} = \tau\_{\rm yx} = \mu \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial \mathbf{x}} \right) \tag{19}
$$

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

$$\begin{split} \sum F\_{\text{surf}} &= \frac{\partial}{\partial \mathbf{x}} \left( -p + 2\mu \frac{\partial u}{\partial \mathbf{x}} \right) \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} + \mu \frac{\partial}{\partial \mathbf{y}} \left( \frac{\partial u}{\partial \mathbf{y}} + \frac{\partial v}{\partial \mathbf{x}} \right) \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} \\ &+ \mu \frac{\partial}{\partial \mathbf{z}} \left( \frac{\partial u}{\partial \mathbf{z}} + \frac{\partial v}{\partial \mathbf{x}} \right) \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} \end{split} \tag{20}$$

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

$$
\sum F\_{\text{body}} = \rho \Delta x \Delta y \Delta \text{zg}\_x \tag{21}
$$

Then, the total force in the x direction is

$$\begin{split} \sum F\_{\mathbf{x}} &= \frac{\partial}{\partial \mathbf{x}} \left( -p + 2\mu \frac{\partial u}{\partial \mathbf{x}} \right) \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} + \mu \frac{\partial}{\partial \mathbf{y}} \left( \frac{\partial u}{\partial \mathbf{y}} + \frac{\partial v}{\partial \mathbf{x}} \right) \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} \\ &+ \mu \frac{\partial}{\partial \mathbf{z}} \left( \frac{\partial u}{\partial \mathbf{z}} + \frac{\partial w}{\partial \mathbf{x}} \right) \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} + \rho \Delta \mathbf{x} \Delta \mathbf{y} \Delta \mathbf{z} \mathbf{g}\_{\mathbf{x}} \end{split} \tag{22}$$

Eqs. (22) and (16) yield

$$\mathbf{a}\_{\mathbf{x}} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial \mathbf{y}^2} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial \mathbf{z}^2} + \frac{\mu}{\rho} \frac{\partial}{\partial \mathbf{x}} \left( \frac{\partial u}{\partial \mathbf{x}} + \frac{\partial v}{\partial \mathbf{y}} + \frac{\partial w}{\partial \mathbf{z}} \right) + \mathbf{g}\_{\mathbf{x}} \tag{23}$$

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

$$
\frac{\partial u}{\partial \mathbf{x}} + \frac{\partial v}{\partial \mathbf{y}} + \frac{\partial w}{\partial \mathbf{z}} = \mathbf{0} \tag{24}
$$

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

$$\mathfrak{a}\_{\mathfrak{x}} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathfrak{x}} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial \mathfrak{x}^2} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial \mathfrak{y}^2} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial \mathfrak{z}^2} + \mathfrak{g}\_{\mathfrak{x}} \tag{25}$$

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) + \mathbf{g}\_x \tag{26}$$

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:

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

$$\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) + \mathbf{g}\_y \tag{27}$$

$$\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 x} + \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) + \mathbf{g}\_z \tag{28}$$

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{DE}{Dt} \Delta x \Delta y \Delta z = \sum \dot{Q} + \sum \dot{W} \tag{29}
$$

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.

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.

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

$$\begin{split} \rho \frac{DE}{Dt} &= -\frac{\partial q\_x}{\partial \mathbf{x}} - \frac{\partial q\_y}{\partial y} - \frac{\partial q\_x}{\partial \mathbf{z}} + \frac{\partial (u \sigma\_{\mathbf{x}x})}{\partial \mathbf{x}} + \frac{\partial (v \sigma\_{\mathbf{yy}})}{\partial \mathbf{y}} \frac{\partial (w \sigma\_{\mathbf{z}x})}{\partial \mathbf{z}} \\ &+ \frac{\partial (u \tau\_{\mathbf{yx}})}{\partial \mathbf{y}} + \frac{\partial (u \tau\_{\mathbf{xx}})}{\partial \mathbf{z}} + \frac{\partial (v \tau\_{\mathbf{xy}})}{\partial \mathbf{x}} + \frac{\partial (v \tau\_{\mathbf{xy}})}{\partial \mathbf{z}} + \frac{\partial (w \tau\_{\mathbf{xx}})}{\partial \mathbf{x}} + \frac{\partial (w \tau\_{\mathbf{yx}})}{\partial \mathbf{y}} \end{split} \tag{30}$$

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

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

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

$$q\_{\chi} = -\lambda \frac{\partial T}{\partial \chi} \tag{31}$$

$$q\_y = -\lambda \frac{\partial T}{\partial y} \tag{32}$$

$$q\_x = -\lambda \frac{\partial T}{\partial \mathbf{z}}\tag{33}$$

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

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

$$
\rho \frac{DE}{Dt} = \lambda \frac{\partial}{\partial \mathbf{x}} \left( \frac{\partial T}{\partial \mathbf{x}} \right) + \lambda \frac{\partial}{\partial \mathbf{y}} \left( \frac{\partial T}{\partial \mathbf{y}} \right) + \lambda \frac{\partial}{\partial \mathbf{z}} \left( \frac{\partial T}{\partial \mathbf{z}} \right) + \Phi \tag{34}
$$

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{DE}{Dt} = \lambda \frac{\partial}{\partial \mathbf{x}} \left( \frac{\partial T}{\partial \mathbf{x}} \right) + \lambda \frac{\partial}{\partial \mathbf{y}} \left( \frac{\partial T}{\partial \mathbf{y}} \right) + \lambda \frac{\partial}{\partial \mathbf{z}} \left( \frac{\partial T}{\partial \mathbf{z}} \right) \tag{35}
$$

If the energy *E* only refers to internal energy, the change of energy D*E* could be represented by *C*pD*T*. Then, Eq. (35) could be further rearranged to Eq. (36)

$$
\rho C\_{\rm p} \frac{DT}{Dt} = \lambda \frac{\partial}{\partial \mathbf{x}} \left( \frac{\partial T}{\partial \mathbf{x}} \right) + \lambda \frac{\partial}{\partial \mathbf{y}} \left( \frac{\partial T}{\partial \mathbf{y}} \right) + \lambda \frac{\partial}{\partial \mathbf{z}} \left( \frac{\partial T}{\partial \mathbf{z}} \right) \tag{36}
$$

Expanding the total derivative of *T*, we get

$$\frac{\partial T}{\partial t} + u \frac{\partial T}{\partial \mathbf{x}} + v \frac{\partial T}{\partial \mathbf{y}} + w \frac{\partial T}{\partial \mathbf{z}} = \frac{\lambda}{\rho \mathbf{C}\_p} \left( \frac{\partial^2 T}{\partial \mathbf{x}^2} + \frac{\partial^2 T}{\partial \mathbf{y}^2} + \frac{\partial^2 T}{\partial \mathbf{z}^2} \right) \tag{37}$$

Eq. (37) is the final expression of the energy equation. The first term on the lefthand 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 (\rho u)}{\partial \mathbf{x}} + \frac{\partial (\rho v)}{\partial \mathbf{y}} + \frac{\partial (\rho w)}{\partial \mathbf{z}} = \mathbf{0} \tag{38}$$

$$\frac{\partial(\rho \mathbf{X}\_i)}{\partial t} + \frac{\partial(\rho u \mathbf{X}\_i)}{\partial \mathbf{x}} + \frac{\partial(\rho v \mathbf{X}\_i)}{\partial \mathbf{y}} + \frac{\partial(\rho w \mathbf{X}\_i)}{\partial \mathbf{z}} = D \left[ \frac{\partial^2(\rho \mathbf{X}\_i)}{\partial \mathbf{x}^2} + \frac{\partial^2(\rho \mathbf{X}\_i)}{\partial \mathbf{y}^2} + \frac{\partial^2(\rho \mathbf{X}\_i)}{\partial \mathbf{z}^2} \right] \tag{39}$$

$$\frac{\partial(\rho u)}{\partial t} + \frac{\partial(\rho uu)}{\partial x} + \frac{\partial(\rho vu)}{\partial y} + \frac{\partial(\rho uu)}{\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) + \mathbf{g}\_x \tag{40}$$

$$\frac{\partial(\rho v)}{\partial t} + \frac{\partial(\rho uv)}{\partial x} + \frac{\partial(\rho vv)}{\partial y} + \frac{\partial(\rho uv)}{\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) + \mathbf{g}\_y \tag{41}$$

$$\frac{\partial(\rho w)}{\partial t} + \frac{\partial(\rho uw)}{\partial x} + \frac{\partial(\rho vw)}{\partial y} + \frac{\partial(\rho ww)}{\partial z} = -\frac{\partial p}{\partial x} + \mu \left( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2} \right) + \mathbf{g}\_x \tag{42}$$

$$\frac{\partial(\rho T)}{\partial t} + \frac{\partial(\rho uT)}{\partial \mathbf{x}} + \frac{\partial(\rho vT)}{\partial \mathbf{y}} + \frac{\partial(\rho wT)}{\partial \mathbf{z}} = \frac{\lambda}{C\_p} \left( \frac{\partial^2 T}{\partial \mathbf{x}^2} + \frac{\partial^2 T}{\partial \mathbf{y}^2} + \frac{\partial^2 T}{\partial \mathbf{z}^2} \right) \tag{43}$$

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

$$\frac{\partial(\rho\phi)}{\partial t} + \frac{\partial(\rho u\phi)}{\partial \mathbf{x}} + \frac{\partial(\rho v\phi)}{\partial \mathbf{y}} + \frac{\partial(\rho w\phi)}{\partial \mathbf{z}} = \Gamma \left( \frac{\partial^2 \phi}{\partial \mathbf{x}^2} + \frac{\partial^2 \phi}{\partial \mathbf{y}^2} + \frac{\partial^2 \phi}{\partial \mathbf{z}^2} \right) \tag{44}$$

where *Γ* is the diffusion coefficient. If *ϕ =* 1, Eq. (44) is the continuity equation. If *ϕ = Xi*, 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

$$\frac{\partial(\rho r \Delta r \Delta \theta \Delta z)}{\partial t} = u\_r (r \Delta \theta \Delta z) - \left(u\_r + \frac{\partial u\_r}{\partial r} \Delta r\right) (r + \Delta r) \Delta \theta \Delta z$$

$$+ u\_\theta (\Delta r \Delta z) - \left(u\_\theta + \frac{\partial u\_\theta}{\partial \theta} \Delta \theta\right) (\Delta r \Delta z) \tag{45}$$

$$+ u\_x \left[ \left(r + \frac{\Delta r}{2}\right) \Delta r \Delta \theta \right] - \left(u\_x + \frac{\partial u\_x}{\partial \mathbf{z}} \Delta \mathbf{z} \right) \left[ \left(r + \frac{\Delta r}{2}\right) \Delta r \Delta \theta \right]$$

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 (ru\_r)}{\partial r} + \frac{1}{r} \frac{\partial u\_\theta}{\partial \theta} + \frac{\partial u\_x}{\partial z} = \mathbf{0} \tag{46}$$

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

$$
\rho \stackrel{\textstyle \cdot}{\mathbf{\hat{a}}} = \frac{\stackrel{\textstyle \cdot}{\mathbf{\hat{F}}}}{V} = \frac{\stackrel{\textstyle \cdot}{\mathbf{\overline{F\_{surf}}}} + \stackrel{\textstyle \cdot}{\mathbf{\overline{F\_{body}}}}}{V} = \frac{\stackrel{\textstyle \cdot}{\mathbf{\overline{F\_{surf}}}}}{V} + \rho \stackrel{\textstyle \cdot}{\mathbf{\hat{g}}} \tag{47}
$$

which is then expressed by velocity vector as

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

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

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

In cylindrical coordinate, the velocity vector is given as

$$
\overrightarrow{\mathbf{u}}(r,\theta,z,t) = u\_r(r,\theta,z,t)\overrightarrow{\mathbf{e}}\_r^\flat + u\_\theta(r,\theta,z,t)\overrightarrow{\mathbf{e}}\_\theta^\flat + u\_z(r,\theta,z,t)\overrightarrow{\mathbf{e}}\_z^\flat \tag{49}
$$

And the acceleration is then given by

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

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<sup>θ</sup>* and *u*<sup>z</sup> is.

$$\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\_x}{\partial t} \overrightarrow{\mathbf{e}\_x} \tag{51}$$

$$
\mu\_r = \frac{\text{d}r}{\text{d}t} \tag{52}
$$

$$
\mu\_{\theta} = \frac{r \mathbf{d}\theta}{\mathbf{d}t} \tag{53}
$$

$$
\mu\_x = \frac{\text{d}z}{\text{d}t} \tag{54}
$$

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

$$\frac{\mathbf{D}\overrightarrow{\mathbf{u}}}{\mathbf{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} \tag{55}$$

**Figure 9.** *The schematic of ∂***u***/∂t in cylindrical coordinate.*

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

*∂***u** ! *<sup>∂</sup><sup>r</sup>* <sup>¼</sup> *<sup>∂</sup> ur***e***<sup>r</sup>* ! <sup>þ</sup> *<sup>u</sup>θ***e***<sup>θ</sup>* ! <sup>þ</sup> *uz***e***<sup>z</sup>* ! *<sup>∂</sup><sup>r</sup>* <sup>¼</sup> *<sup>∂</sup>ur ∂r* **e***r* ! <sup>þ</sup> *ur ∂***e***<sup>r</sup>* ! *∂r* þ *∂u<sup>θ</sup> ∂r* **e***θ* ! <sup>þ</sup> *<sup>u</sup><sup>θ</sup> ∂***e***<sup>θ</sup>* ! *∂r* þ *∂uz ∂r* **e***z* ! <sup>þ</sup> *uz ∂***e***<sup>z</sup>* ! *∂r* (56) *∂***u** ! *<sup>∂</sup><sup>θ</sup>* <sup>¼</sup> *<sup>∂</sup> ur***e***<sup>r</sup>* ! <sup>þ</sup> *<sup>u</sup>θ***e***<sup>θ</sup>* ! <sup>þ</sup> *uz***e***<sup>z</sup>* ! *<sup>∂</sup><sup>θ</sup>* <sup>¼</sup> *<sup>∂</sup>ur <sup>∂</sup><sup>θ</sup>* **<sup>e</sup>***<sup>r</sup>* ! <sup>þ</sup> *ur ∂***e***r* ! *∂θ* þ *∂u<sup>θ</sup> <sup>∂</sup><sup>θ</sup>* **<sup>e</sup>***<sup>θ</sup>* ! <sup>þ</sup> *<sup>u</sup><sup>θ</sup> ∂***e***<sup>θ</sup>* ! *∂θ* þ *∂uz <sup>∂</sup><sup>θ</sup>* **<sup>e</sup>***<sup>z</sup>* ! <sup>þ</sup> *uz ∂***e***<sup>z</sup>* ! *∂θ* (57) *∂***u** ! ! <sup>þ</sup> *<sup>u</sup>θ***e***<sup>θ</sup>* ! <sup>þ</sup> *uz***e***<sup>z</sup>* ! *∂***e***r* ! *∂u<sup>θ</sup> ∂***e***<sup>θ</sup>* ! *∂uz ∂***e***<sup>z</sup>* !

$$\frac{d\overrightarrow{\mathbf{u}}}{d\mathbf{z}} = \frac{\partial \left( u\_r \overrightarrow{\mathbf{e}\_r} + u\_\theta \overrightarrow{\mathbf{e}\_\theta} + u\_z \overrightarrow{\mathbf{e}\_z} \right)}{d\mathbf{z}} = \frac{\partial u\_r}{d\mathbf{z}} \overrightarrow{\mathbf{e}\_r} + u\_r \frac{\partial \overrightarrow{\mathbf{e}\_r}}{d\mathbf{z}} + \frac{\partial u\_\theta}{d\mathbf{z}} \overrightarrow{\mathbf{e}\_\theta} + u\_\theta \frac{\partial \overrightarrow{\mathbf{e}\_\theta}}{d\mathbf{z}} + \frac{\partial u\_z}{d\mathbf{z}} \overrightarrow{\mathbf{e}\_z} + u\_z \frac{\partial \overrightarrow{\mathbf{e}\_z}}{d\mathbf{z}} \tag{58}$$

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^{\rightarrow}}}{\partial r} = \frac{\partial \overrightarrow{\mathbf{e}\_\theta^{\rightarrow}}}{\partial r} = \frac{\partial \overrightarrow{\mathbf{e}\_x^{\rightarrow}}}{\partial r} = \mathbf{0} \tag{59}$$

$$\frac{\partial \overrightarrow{\mathbf{e}\_r^{\rightarrow}}}{\partial \theta} = \overrightarrow{\mathbf{e}\_\theta^{\rightarrow}} \tag{60}$$

$$\frac{\partial \overrightarrow{\mathbf{e}\_{\theta}^{\rightarrow}}}{\partial \theta} = -\overrightarrow{\mathbf{e}\_{r}^{\rightarrow}}\tag{61}$$

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\_x}{\partial r} \overrightarrow{\mathbf{e}\_x} \tag{62}$$

$$\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\_x}{\partial \theta} \overrightarrow{\mathbf{e}\_x} \tag{63}$$

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

$$\frac{\partial \overrightarrow{\mathbf{u}}}{\partial \mathbf{z}} = \frac{\partial u\_r}{\partial \mathbf{z}} \overrightarrow{\mathbf{e}\_r} + \frac{\partial u\_\theta}{\partial \mathbf{z}} \overrightarrow{\mathbf{e}\_\theta} + \frac{\partial u\_z}{\partial \mathbf{z}} \overrightarrow{\mathbf{e}\_z} \tag{64}$$

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

D**u** ! <sup>D</sup>*<sup>t</sup>* <sup>¼</sup> *<sup>∂</sup>***<sup>u</sup>** ! *∂t* þ *∂***u** ! *∂r* d*r* d*t* þ *∂***u** ! *∂θ* d*θ* d*t* þ *∂***u** ! *∂z* d*z* d*t* <sup>¼</sup> *<sup>∂</sup>ur ∂t* **e***r* ! <sup>þ</sup> *∂u<sup>θ</sup> ∂t* **e***θ* ! <sup>þ</sup> *∂uz ∂t* **e***z* ! þ*ur ∂ur ∂r* **e***r* ! <sup>þ</sup> *∂u<sup>θ</sup> ∂r* **e***θ* ! <sup>þ</sup> *∂uz ∂r* **e***z* ! þ *uθ r ∂ur <sup>∂</sup><sup>θ</sup>* � *<sup>u</sup><sup>θ</sup>* **<sup>e</sup>***<sup>r</sup>* ! <sup>þ</sup> *ur* <sup>þ</sup> *∂u<sup>θ</sup> ∂θ* **<sup>e</sup>***<sup>θ</sup>* ! <sup>þ</sup> *∂uz <sup>∂</sup><sup>θ</sup>* **<sup>e</sup>***<sup>z</sup>* ! þ*uz ∂ur ∂z* **e***r* ! <sup>þ</sup> *∂u<sup>θ</sup> ∂z* **e***θ* ! <sup>þ</sup> *∂uz ∂z* **e***z* ! <sup>¼</sup> *<sup>∂</sup>ur ∂t* þ *ur ∂ur ∂r* þ *uθ r ∂ur <sup>∂</sup><sup>θ</sup>* � *<sup>u</sup>*<sup>2</sup> *θ r* þ *uz ∂ur ∂z* **<sup>e</sup>***<sup>r</sup>* ! þ *∂u<sup>θ</sup> ∂t* þ *ur ∂u<sup>θ</sup> ∂r* þ *uθ r ∂u<sup>θ</sup> ∂θ* þ *uθur r* þ *uz ∂u<sup>θ</sup> ∂z* **<sup>e</sup>***<sup>θ</sup>* ! þ *∂uz ∂t* þ *ur ∂uz ∂r* þ *uθ r ∂uz <sup>∂</sup><sup>θ</sup>* <sup>þ</sup> *uz ∂uz ∂z* **<sup>e</sup>***<sup>z</sup>* ! (65)

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

$$\mathcal{a}\_r = \frac{\mathcal{D}u\_r}{\mathcal{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} \tag{66}$$

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

$$\mathbf{a}\_r = \frac{\mathbf{1}}{\rho} \frac{F\_r}{\delta \mathbf{V}} + \mathbf{g}\_r \tag{67}$$

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

$$\begin{aligned} F\_r &= \sigma\_{rr}(r+\Delta r)\Delta\theta\Delta z + \frac{\partial\sigma\_{rr}}{\partial r}\Delta r(r+\Delta r)\Delta\theta\Delta z - \sigma\_{rr}(r\Delta\theta\Delta z) \\ &+ \sigma\_{\theta r}(\Delta r\Delta z)\cos\frac{\Delta\theta}{2} + \frac{\partial\sigma\_{\theta r}}{\partial \theta}\Delta\theta(\Delta r\Delta z)\cos\frac{\Delta\theta}{2} - \sigma\_{\theta r}(\Delta r\Delta z)\cos\frac{\Delta\theta}{2} \\ &+ \sigma\_{\text{zr}}\left[\left(r+\frac{\Delta r}{2}\right)\Delta r\Delta\theta\right] + \frac{\partial\sigma\_{\text{zr}}}{\partial z}\Delta z \left[\left(r+\frac{\Delta r}{2}\right)\Delta r\Delta\theta\right] - \sigma\_{\text{r}xz}\left[\left(r+\frac{\Delta r}{2}\right)\Delta r\Delta\theta\right] \\ &- \sigma\_{\theta\theta}(\Delta r\Delta z)\sin\frac{\Delta\theta}{2} - \frac{\partial\sigma\_{\theta\theta}}{\partial\theta}\Delta\theta(\Delta r\Delta z)\sin\frac{\Delta\theta}{2} - \sigma\_{\theta\theta}(\Delta r\Delta z)\sin\frac{\Delta\theta}{2} \end{aligned} \tag{68}$$

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.

#### **Figure 10.**

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

$$\begin{aligned} F\_r &= \sigma\_{rr} \Delta r \Delta \theta \Delta z + \frac{\partial \sigma\_{rr}}{\partial r} \Delta r (r + \Delta r) \Delta \theta \Delta z \\ &+ \frac{\partial \sigma\_{\theta r}}{\partial \theta} \Delta \theta (\Delta r \Delta z) \cos \frac{\Delta \theta}{2} \\ &+ \frac{\partial \sigma\_{\theta r}}{\partial z} \Delta z \left[ \left( r + \frac{\Delta r}{2} \right) \Delta r \Delta \theta \right] \\ &- 2 \sigma\_{\theta \theta} (\Delta r \Delta z) \frac{\Delta \theta}{2} - \frac{\partial \sigma\_{\theta \theta}}{\partial \theta} \Delta \theta (\Delta r \Delta z) \frac{\Delta \theta}{2} \end{aligned} \tag{69}$$

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

$$\frac{F\_r}{\Delta V} = \frac{\sigma\_{rr}}{r} + \frac{\partial \sigma\_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma\_{\theta r}}{\partial \theta} + \frac{\partial \sigma\_{xr}}{\partial z} - \frac{\sigma\_{\theta\theta}}{r} \tag{70}$$

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\_x \frac{\partial u\_r}{\partial z} = \frac{1}{\rho} \left( \frac{\sigma\_{rr}}{r} + \frac{\partial \sigma\_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma\_{\theta r}}{\partial \theta} + \frac{\partial \sigma\_{xr}}{\partial z} - \frac{\sigma\_{\theta \theta}}{r} \right) + \mathbf{g}\_r \tag{71}$$

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

$$
\sigma\_{rr} = -p + 2\mu \frac{\partial u\_r}{\partial r} \tag{72}
$$

$$
\sigma\_{\partial\theta} = -p + 2\mu \left( \frac{1}{r} \frac{\partial u\_{\theta}}{\partial \theta} + \frac{u\_r}{r} \right) \tag{73}
$$

$$
\sigma\_{xx} = -p + 2\mu \frac{\partial u\_x}{\partial \mathbf{z}}\tag{74}
$$

$$
\sigma\_{r\theta} = \sigma\_{\theta r} = \mu \left( \frac{1}{r} \frac{\partial u\_r}{\partial \theta} + \frac{\partial u\_\theta}{\partial r} - \frac{u\_\theta}{r} \right) \tag{75}
$$

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

$$
\sigma\_{rz} = \sigma\_{xr} = \mu \left( \frac{\partial u\_r}{\partial z} + \frac{\partial u\_x}{\partial r} \right) \tag{76}
$$

$$
\sigma\_{\partial \mathbf{z}} = \sigma\_{x\theta} = \mu \left( \frac{\mathbf{1}}{r} \frac{\partial u\_x}{\partial \theta} + \frac{\partial u\_\theta}{\partial \mathbf{z}} \right) \tag{77}
$$

Combining Eqs. (70) and (71) yields

$$\begin{split} &\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\_x \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] + \mathbf{g}\_r \end{split} \tag{78}$$

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{split} \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\_x \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] + \mathbf{g}\_{\theta} \end{split} \tag{79}$$
 
$$\begin{split} \frac{\partial u\_x}{\partial t} + u\_r \frac{\partial u\_x}{\partial r} + \frac{u\_{\theta}}{r} \frac{\partial u\_x}{\partial \theta} + u\_x \frac{\partial u\_x}{\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\_x}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 u\_x}{\partial \theta^2} + \frac{\partial^2 u\_x}{\partial z^2} \right] + \mathbf{g}\_x \end{split} \tag{80}$$

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

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

$$
\rho \frac{DE}{Dt} r \Delta r \Delta \theta \Delta z = -\left(q\_r + \frac{\partial q\_r}{\partial r} \Delta r\right) (r + \Delta r) \Delta \theta \Delta z + q\_r (r \Delta \theta \Delta z)
$$

$$
$$

$$
$$

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

$$
\rho \frac{DE}{Dt} = -\frac{1}{r} \frac{\partial \left( r q\_r \right)}{\partial r} - \frac{1}{r} \frac{\partial q\_\theta}{\partial \theta} - \frac{\partial q\_x}{\partial x} \tag{82}
$$

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

$$\rho C\_{\text{p}} \left( \frac{\partial T}{\partial t} + u\_r \frac{\partial T}{\partial r} + \frac{u\_\theta}{r} \frac{\partial T}{\partial \theta} + u\_x \frac{\partial T}{\partial x} \right) = -\frac{1}{r} \frac{\partial (rq\_r)}{\partial r} - \frac{1}{r} \frac{\partial q\_\theta}{\partial \theta} - \frac{\partial q\_x}{\partial x} \tag{83}$$

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} \tag{84}$$

$$q\_{\theta} = -\frac{\lambda}{r} \frac{\partial T}{\partial \theta} \tag{85}$$

$$q\_x = -\lambda \frac{\partial T}{\partial x} \tag{86}$$

$$\begin{aligned} \rho \mathbf{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\_x \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{aligned} \tag{87}$$
