Computational Fluid Dynamics in Turbulent Flow Applications

Abstract

This chapter is intended to present to readers a general scope of the technical, theoretical, and numerical applications of computational fluid dynamics using the finite volume method, restricted to incompressible turbulent flows (Ma < 0.3). The main objective of this chapter was to provide readers of a starting point to select an adequate numerical model for the flow regime of interest. Such knowledge could be a key at the moment of extending the analysis to more complex problems, for example, the ones found in heat transfer and fluid flows, multiphase flows, and compressible flows.

1. Introduction

Computational fluid dynamics (CFD) is a scientific tool capable of producing information about the main structures of a flowing fluid. As a knowledge area, it finds its origins in the discrete solution of the fundamental equations used in fluid dynamics, such as the mass conservation equation, the momentum conservation equations (based on Newton’s second law), and the energy conservation equation (based on the first law of Thermodynamics). All the equations are solved in an approximate way using iterative algorithms to solve lineal algebraic equations systems using matrices [1, 2].

A CFD model has three sequential stages known as preprocessing, solving, and post-processing. In the first stage or preprocessing, the governing equations, properties, and boundary conditions are defined within a domain composed by small interconnected volumetric elements to model the fluid movement. During the second stage or solving process, the adequate selection of discretization schemes (temporal and spatial) for the governing equations is carried out. Also, the convergence criteria or number of iterations, relaxation factors, and pressure and velocity coupling algorithm are set up at this stage. Finally, the post-processing stage involves the analysis and interpretation of the obtained solutions using flow, temperature, and pressures fields [3].

Nowadays, there are different commercial codes developed in a friendly interface, such as ANSYS FLUENT, OPEN FOAM, CFX, X-FLOW, and COMSOL. All of them are capable of helping researchers and engineers at the three aforementioned CFD stages. Generally, the modules of these commercial softwares have been validated comparing the results obtained for specific information such as separation point, circulation length, drag and lift coefficients, and velocities, with experimental and/or analytical results. Very often, the experimental results used for comparison are obtained from carefully controlled experiments about topics such as flow around circular cylinders, flow over forward or backward facing step, flow over a venture, and jet flow, among others [4]. Respecting the programming languages commonly used to develop the CFD codes used in commercial software, the most frequently used are C++, FORTRAN, and lately MATLAB.

2. Governing equations in computational fluid dynamics

The governing equations in computational fluid dynamic mathematically express the three fundamental physical principles that describe the movement of any fluid. These equations are as follows:

• •

  Continuity equation or mass conservation principle.

• •

  Newton’s second law or momentum conservation principle.

• •

  The first law of thermodynamics or energy conservation principle.

A fluid is a substance that, due to molecular distances, does not present a defined form and adopts the form of the vessel that contains it. Therefore, it is difficult to analyze a fluid from the universal approach used for solids. In general, a fluid can be defined as a substance that deforms continually under the action of a shear stress. If a fluid is in movement, the velocity can be different at different positions within the studied domain, and their particles can rotate and deform at the same time. Considering the fluid as a continuum, there are two approaches that can be used to understand its movement, making it susceptible to be analyzed using the fundamental principles: The first one is based on the study of a volume element of infinitesimal size and fixed in the space with the fluid flowing through it (known as the Eulerian approach) and an infinitesimal fluid element that moves along a streamline with a velocity equal to the local velocity of the flow at each point (Lagrangian approach) [5].

When a differential element of a fluid moving along a streamline, and considering the motion of a particle with a defined velocity through a differential volume defined in a Cartesian space, the governing equations for a compressible viscous flow are as follows:

1. (a)

   Continuity equation:

   $\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \overrightarrow{V}\right)=0$E1.1

2. (b)

   Navier-Stokes equations:

x-direction component (î)

$\frac{\partial (\rho u)}{\partial t}+\nabla \cdot \left(\rho u\overrightarrow{V}\right)=-\frac{\partial P}{\partial x}+\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{yx}}{\partial y}+\frac{\partial {\tau }_{zx}}{\partial z}+\rho f_x$E1.2

y-direction component (Ĵ)

$\frac{\partial (\rho v)}{\partial t}+\nabla \cdot \left(\rho v\overrightarrow{V}\right)=-\frac{\partial P}{\partial y}+\frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial {\tau }_{yy}}{\partial y}+\frac{\partial {\tau }_{zy}}{\partial z}+\rho f_y$E1.3

z-direction component ($$\widehat{k}$$)

$\frac{\partial (\rho w)}{\partial t}+\nabla \cdot \left(\rho w\overrightarrow{V}\right)=-\frac{\partial P}{\partial z}+\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\tau }_{yz}}{\partial y}+\frac{\partial {\tau }_{zz}}{\partial z}+\rho f_z$E1.4

Energy equation:

$\begin{array}{l}\frac{\partial }{\partial t}\left[\rho \left(e+\frac{u^2+v^2+w^2}{2}\right)\right]+\nabla \cdot \left[\rho \left(e+\frac{u^2+v^2+w^2}{2}\right)\overrightarrow{V}\right]=\rho \dot{q}+\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\\ \frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)-\frac{\partial \left(uP\right)}{\partial x}-\frac{\partial \left(vP\right)}{\partial y}-\frac{\partial \left(wP\right)}{\partial z}+\frac{\partial \left(u{\tau }_{xx}\right)}{\partial x}+\frac{\partial \left(u{\tau }_{yx}\right)}{\partial y}+\frac{\partial \left(u{\tau }_{zx}\right)}{\partial z}+\frac{\partial \left(v{\tau }_{xy}\right)}{\partial x}+\\ \frac{\partial \left(v{\tau }_{yy}\right)}{\partial y}+\frac{\partial \left(v{\tau }_{zy}\right)}{\partial z}+\frac{\partial \left(w{\tau }_{xz}\right)}{\partial x}+\frac{\partial \left(w{\tau }_{yz}\right)}{\partial y}+\frac{\partial \left(w{\tau }_{zz}\right)}{\partial z}+\rho \overrightarrow{f}\cdot \overrightarrow{V}\end{array}$E1.5

where e is the internal energy per unit of mass, $\overrightarrow{f}$ represents the body forces that act on the centroid of the fluid element (like gravitational, electrical or magnetic forces), “k” is thermal conductivity, P is the static pressure, $\dot{q}$ is the heat transfer rate per unit of mass in the volume element (it can be heat that is a combustion product, a chemical reaction, or an electron flow), t is time, T is temperature, ρ is density, τ are viscous stresses, and $\nabla =\frac{\partial }{\partial x}\overrightarrow{i}+\frac{\partial }{\partial y}\overrightarrow{j}+\frac{\partial }{\partial z}\overrightarrow{k}$ is the DEL vector operator for Cartesian coordinates. It is necessary to point out that viscous stresses (normal and shear stresses) are related to the deformation rate of the fluid element. Shear stresses are related to the temporal deformation caused by the constant shear force acting on the fluid element. Conventionally, τ_ij represents a stress in the “j” direction, exerted on a plane perpendicular to i-axis. Normal stresses are related to the volume change rate in the fluid element (compression or tension); these stresses have in general a lower value than shear stresses.

At the end of the seventeenth century, Isaac Newton established that the shear stress in a fluid is proportional to its temporal deformation rate (velocity gradients, for example). Fluids that follow this behavior are known as Newtonian fluids, and they have many applications in fluid dynamics. For this kind of fluids, Stokes found, in 1845, the following relationships:

$\begin{array}{l}{\tau }_{xy}={\tau }_{yx}=\mu \left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right){\tau }_{xx}=\lambda \left(\nabla \cdot V\right)+2\mu \frac{\partial u}{\partial x}\\ {\tau }_{xz}={\tau }_{zx}=\mu \left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right){\tau }_{yy}=\lambda \left(\nabla \cdot V\right)+2\mu \frac{\partial v}{\partial y}\\ {\tau }_{yz}={\tau }_{zy}=\mu \left(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}\right){\tau }_{zz}=\lambda \left(\nabla \cdot V\right)+2\mu \frac{\partial w}{\partial z}\end{array}$E1.6

where μ is the molecular viscosity and λ is the second viscosity (which is not a term of common use in engineering). For gases, a good approximation can be obtained using the value λ = −2/3μ.

Non-Newtonian fluids are those that do not present a linear relationship between stresses and deformation rates due to shear stresses (velocities gradients). A non-Newtonian fluid can present viscoelastic and thixotropic properties as well as different kinds of relationships among density, pressure, and temperature in comparison with a Newtonian fluid. This kind of fluid requires a more specialized treatment to study them, which is out of the scope of this document.

The governing equations presented (1.1)–(1.5) contain seven unknown flow variables expressed in a set of five differential equations. For practical applications, additional equations can be used to “close” the system (equal number of equations and unknowns). For example, in the study of an aerodynamic phenomenon, it is in general possible to assume that a gas behaves as an “ideal gas”. For an ideal gas, the state equation is = ρRT, where R is the gas constant. This equation provides the equation system with a sixth equation. It can also be considered a seventh equation to close the equation system using a thermodynamic relationship among state variables, for example: e = e(T, P) For an ideal gas (with constant specific heats), this relationship becomese = C_vT where C_v is the constant volume-specific heat.

Finally, applying the substantial derivative definition for Cartesian coordinatesD/Dt = ∂/∂t + $u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}+w\frac{\partial }{\partial z}$ Eqs. (1.1)–(1.5) can be written as follows:

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho V\right)=0$E1.7

$\rho \frac{Du}{Dt}=-\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]+\rho f_x$E1.8

$\rho \frac{Dv}{Dt}=-\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]+\rho f_y$E1.9

$\rho \frac{Dw}{Dt}=-\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]+\rho f_z$E1.10

$\begin{array}{l}\frac{D}{Dt}\left[\rho \left(e+\frac{u^2+v^2+w^2}{2}\right)\right]=\rho \dot{q}+\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)\\ -\frac{\partial \left(uP\right)}{\partial x}-\frac{\partial \left(vP\right)}{\partial y}-\frac{\partial \left(wP\right)}{\partial z}+\frac{\partial \left(u{\tau }_{xx}\right)}{\partial x}+\frac{\partial \left(u{\tau }_{yx}\right)}{\partial y}+\frac{\partial \left(u{\tau }_{zx}\right)}{\partial z}+\\ \frac{\partial \left(v{\tau }_{xy}\right)}{\partial x}+\frac{\partial \left(v{\tau }_{yy}\right)}{\partial y}+\frac{\partial \left(v{\tau }_{zy}\right)}{\partial z}+\frac{\partial \left(w{\tau }_{xz}\right)}{\partial x}+\frac{\partial \left(w{\tau }_{yz}\right)}{\partial y}+\frac{\partial \left(w{\tau }_{zz}\right)}{\partial z}+\rho \overrightarrow{f}\cdot V\end{array}$E1.11

The governing equations presented (1.1)–(1.5) or (1.7)–(1.11) are widely known as the Navier–Stokes equations to honor the French physicist Claude Louis Navier and the English mathematician George Gabriel Stokes, who in an independent way, both obtained the equations in the first half of the nineteenth century. Originally, the Navier–Stokes equations terminology was only defined for the momentum equations. However, current CFD literature has been expanded to include the equations of mass and energy conservation. It should be noted that the governing equations can be also expressed in cylindrical and spherical coordinates, or even in generalized curvilinear coordinates.

It has to be kept in mind that these equations were originally formulated to reproduce the physics present in single-phase and non-reactive Newtonian flows. The governing equations for reactive flow systems with multiple components (like the ones found in combustion and multiphase flow problems) can be also established. However, they are more complex because they involve multiple species and phases. It is necessary to note that for this complex flow systems, it becomes a necessity to include approximation models or empirical correlations because several terms are not constant anymore and they become functions for example of temperature, pressure, location, introducing new nonlinearities in the diffusive terms.

The governing equations constitute a coupled nonlinear partial differential equations system; therefore, an analytical solution is difficult to obtain. At the present time, there is not a general solution in a closed form for this equation system. In fact, this is one of the most important reasons of the use of CFD, since it provides numerical approximations to the equations solutions, which have not been found yet using analytical methods, except for idealized cases for one-dimensional or two-dimensional flows in laminar or creeping flows.

The mathematical character of the equations has a significant impact on CFD. First of all, it is very important that the problem is adequately posed, and it means that the solution of the problem exists and it is unique. Furthermore, this solution is only affected by initial and boundary conditions. It is also important to classify the governing equations as elliptical, parabolic, or hyperbolic. Elliptical equations have to be solved simultaneously in the complete flow domain, whereas hyperbolic and parabolic equations propagate from one position to another one.

Mathematically, parabolic and hyperbolic equations are susceptible to be solved using time steps. On the other hand, for the elliptical equations, the flow variables in a given point should be solved simultaneously with the flow variables in other positions. It can be stated that the Navier-Stokes equations have, in general, a mixed nature. They are parabolic in the time domain and elliptical in the space domain [1, 3].

As it was mentioned before, the governing equations can be applied for laminar and turbulent flows by means of the use of additional terms to represent the influence of the fluctuating eddies. In a basic and rigorous way, direct numerical simulation (DNS) works solving numerically these equations to a desired accuracy degree without any additional model or correlation. However, its application is still limited because of the large quantity of resources required to solve the majority of the existent problems, which exceed the capacity of conventional computers.

For turbulent flows, there is a wide spectrum of time and length scales that have to be solved that increase the computational time in DNS simulation. However, from a scientific point of view, the information obtained using DNS is very valuable because it is considered a validation means for turbulent models (with less complexity). The traditional performance of the CFD methods for turbulent flows has been focused on predicting the average influence of the main characteristics of turbulence. Very often, such influence is approached including new transport equations to replicate the effects of generation or dissipation of Reynolds stresses , which represent the transport of momentum in the streamwise direction (x-direction) caused by the turbulent eddies in the “i” and “j” directions. Frequently, average stress models adjust to semi-empirical models, which depend on a number of constants obtained using experimental correlations. These correlations are very often obtained from experiments where the working fluid is air or water, and very frequently, the experimental setup includes a flat plate. In the literature, this kind of closing model is known as a Reynolds Average Navier–Stokes (RANS) turbulence model. This kind of model is the most frequently used for engineering applications today. Unfortunately, for complex turbulent flows there is not a “universal” model that can be applied to any kind of flow, since boundary layer effects can be very different. In fact, the description of what happens within the boundary layer is the main weakness of the RANS models. For many years, RANS models have had success reproducing turbulent zones far away from the wall for fully developed flows. However, they still have problems reproducing unsteady flows or flow separation phenomena (as a result of an adverse pressure gradient) [6, 7]. Even though it is true that for very complex turbulent flows, there exists new and more powerful RANS models, capable of giving useful information to make decisions in an industry for example. Nevertheless, in the scientific world, those models do not seem to have universal validity. Due to that, alternative approaches such as large eddy simulation (LES) and DNS, with the help of the technological evolution of computers, have gained more attention in the last three decades.

3. The ANSYS FLUENT software

The software ANSYS FLUENT is one of the most used tools to discretely solve the governing equations in fluid dynamics. It is very useful as a tool oriented to reproduce the behavior and give information for a wide range of fluid flows (compressible or incompressible, laminar or turbulent, steady or unsteady). It also possesses a wide variety of mathematical models for transport phenomena (such as heat transfer and chemical reactions) that can be applied to both simple and more complex geometries. FLUENT can be applied to laminar or turbulent flows within devices, such as heat exchangers, nozzles, complex-shaped ducts, turbomachinery (heat transfer and aerodynamic), heat engine components, external flows, flow through compressors, pumps, fans, and multiphase flows, among others.

In order to make a generalization of the governing equations to include compressibility effects and interactions among different phases, the continuity equation used by FLUENT can be described as follows:

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \overrightarrow{V}\right)=S_m$E1.12

where the source term S_m represents the mass added to the continuous phase by the disperse phase. Other source terms can be defined by the user.

With respect to the other quantities transported by the fluid, they are modeled based on an integral general equation, where variable ϕ can be substituted by different flow variables before being discretized in its matrix form. For that Eq. (1.13), represents a dissipation constant for the transported fluid, which can be, for example, the dynamic viscosity μ for the momentum transport equations in any direction, or the Eddie viscosity μ_t associated with the energy production–dissipation models for the turbulent kinetic energy, whose magnitude depends on secondary factors such as turbulent kinetic energy intensity and some semi-empirical correlations.

$\frac{\partial }{\partial t}{\displaystyle {\int }_V\rho \varphi dV}+{\displaystyle {\int }_A\rho \varphi \overrightarrow{V}\cdot dA}={\displaystyle {\int }_A{\text{Γ}}_{\varphi }\nabla \varphi \cdot dA}+{\displaystyle {\int }_V{S}_{\varphi }dV}$E1.13

Table 1 shows some equivalencies associated to the general variable.

Equation	Variable (φ)
Continuity	1
X momentum	u
Y momentum	v
Z momentum	w
Transport of turbulent kinetic energy	k
Dissipation of the turbulent kinetic energy	ϵ

Table 1.

Different variables used in the general equation.

Figure 1.

The coupled and segregated solvers comparison diagram, adapted with permission from [8].

The FLUENT solvers work based on a logical sequence to solve the discretized equations in two ways, the segregated one or the coupled one [8]. The segregated mode is recommended for flow problems where it can be expected the formation of large velocity gradients. These gradients can be produced by a boundary layer formation with or without flow separation where the compressibility effects are minimal (Mach < 0.3) and with or without temperature gradients at the boundary. Some examples of the application of the segregated solver are as follows: vehicular aerodynamics problems, internal flows, centrifugal pumps velocity fields’ analysis, cyclone separator flows, flows in refrigeration ducts, flow pattern in cooling fins, etc. Regarding the coupled solver, it is widely used for the solution of problems related to large compressibility effects in the fluid, such as supersonic flows, hypersonic flows within nozzles combustion chambers, the movement of projectiles, and other objects where the shock waves are produced due to large pressure gradients and density changes.

Numerically speaking, the biggest difference between both solver modes is the way the flow governing equations are solved. The segregated mode solve them sequentially, taking the velocity and pressure equations (or corrected pressure equations) as the main variables, the coupled one solves globally the equations for continuity, momentum, and species exchange. Figure 1 shows a diagram of both solver modes’ performances.

3.1. Main parameters and boundary conditions

Although ANSYS FLUENT is a very friendly application, its implementation requires of an adjustment for some predetermined parameters, which will guide the solution process.

Before uploading the mesh associated with the physical domain to be studied, the software requires to define the number of dimensions used to reach the solution; it can be two-dimension (2D) or 3D. To analyze uniform cross sections in bodies with infinite length, it is recommended to use 2D analysis, suppressing the third coordinate because gradients along the depth are not significant. The software offers an option to choose double precision (more significant ciphers). This option is recommended when large gradients in the flow variables are expected, especially at very close nodal position, cells with very large aspect ratio or simply when a better convergence is desirable.

Another important characteristic is the capacity to adjust the number of core processors to be used internally in parallel. This characteristics is especially helpful for problems that require many iterations to be solved (e.g.,, non-stationary models), multiphase and reactive problems, very large grid sizes, and very complex geometries such as airplanes and vehicles.

Finally, this computer program allows to scale the grid physical dimensions to required sizes using both English and international systems. It can be used to work on problems such as the flow between building zones located hundreds of meters from each other or problems related to microfluids. It also allows to select constant values to calculate non-dimensional parameters such as friction coefficient and lift and drag coefficients.

In any CFD model, boundary conditions have to represent faithfully the real conditions at the boundary. According to the mathematical characteristics of the differential equations, boundary conditions can be classified as follows:

Dirichlet boundary condition: The value of the flow variable is specified at the boundary. This kind of boundary condition is typically linked to problems involving flow inlets and isothermal walls.

Neumann boundary condition: For this kind of boundary condition, the value of the gradient normal to the flow variable is specified. These boundary conditions are often associated with symmetric boundaries and adiabatic walls.

Mixed boundary conditions: These are the conditions resulting from a combination of boundary conditions of the kind of Neumann and Dirichlet.

The boundary conditions can be classified by its physical meaning as follows: physical boundaries, such as solid walls, or artificial boundaries, such as outflow. The latter are mathematical approximations to the real behavior of the flow in certain zones. Artificial boundaries can be applied, if the computational domain constitutes part of the total flow field, to research the most interesting region and to reduce computational costs. Artificial boundary conditions require mathematical formulations that are physically significant to the real flow behavior.

Regarding ANSYS FLUENT, there are different types of boundary conditions that can be used to simulate different flow problems. Next, a general description of the most common boundary conditions used in most hydrodynamic studies is presented. The authors recommend readers to consult Ref. [8] for more detailed information.

3.1.1. Velocity inlet

The “velocity inlet” condition was conceived to prescribe a uniform velocity profile at the entrance of a computational domain (bi-dimensional or tri-dimensional). This boundary condition was created to analyze incompressible flows, such as external flows over complex geometries or internal flows in different sections of piping or machinery. It is not recommended for compressible flows. It also allows to specify the magnitude and the direction of a uniform velocity profile. Other velocity profiles such as parabolic, sinusoidal, and logarithm can be imposed using a user-defined function code (UDF). For turbulent flows, the software allows prescribing the turbulence levels required at the entrance or the domain using different methods. A hydraulic diameter and a value of turbulence intensity can be used to do so for internal flows. For example, for square ducts, the hydraulic diameter can be calculated as 4A/P, where A is the cross section area and P is the perimeter of the same area, also known as the wetted perimeter. For external flows, it is required to provide a characteristic length (for example, the diameter of a cylinder if the flow over a circular cylinder is studied). In this case, the turbulence intensity represents a percentage of the value of the fluctuating values of velocity based on the mean value. For example, if the magnitude of the fluid velocity is 5 m/s and a turbulence level of 10% is prescribed, the inlet velocity would present values from 4.5 to 5.5 m/s. Additionally, different parameters like flow temperatures or heat flux can be established at the entrance of the domain, which are important for energy studies and heat and mass transfer analysis.

3.1.2. Pressure inlet

The “pressure inlet” boundary condition is used to define the pressure characteristics at the entrance of the domain for cases where the mass flow or the velocity is unknown at the inlet of the domain. This boundary condition can be used for compressible and incompressible flows. For incompressible flows, the total pressure can be modeled as follows:

$P_{total}=P_{estatic}+\frac{1}{2}\rho V^2$E1.14

where the total pressure or “gauge total pressure” is known, and it is equivalent to the sum of the static pressure at the inlet, also known as supersonic/initial gauge pressure and a dynamic pressure, which is not necessarily known, but can be manually calculated and compared to the velocity assigned to the normal direction to the boundary condition.

$P_{total,abs}=P_{estatic}{\left[1+\frac{k-1}{2}{M}^2\right]}^{\frac{k}{k-1}}$E1.15

Regarding the compressible flows, Eq. (1.15) is used, where “k” is the specific heats ratio and M is the Mach number.

3.1.3. Outflow

The “outflow” boundary condition represents and artificial cut through the flow field, similar to the velocity inlet condition, but set up at the outlet of the flow domain. The outflow conditions are used when the characteristics looked for in the flow are developed at a fraction of the total flow field. When this is true, this boundary condition allows to avoid numerically modeling the complete domain (which can be computationally expensive). In those cases, a distance from the domain inlet is defined as the flow outlet, imposing there the outflow condition. The main difference between an inflow and an outflow condition is that there is not flow information available outside the computational domain for an outflow condition. In contrast, for the inflow condition, there is always information about the entering flow.

The flow variables in an outflow condition have to be approximated in a physically significant way in a manner that does not affect the solution of the governing equations in the computational domain. For an outflow boundary condition, the numerical effects found upstream that are generated have to be eliminated or reduced. Conventionally, in the CFD programming for laminar and steady RANS models, the fully developed conditions in the streamwise direction are expressed by ∂φ/∂x_i = 0. However, for more complex cases where recirculating flow structures exist, there can be problems due to the reentrance of mass circulating in a vortex for example. For these cases, a generalized convective boundary condition is used:

$\frac{\partial \varphi }{\partial t}+U_{conv}{\frac{\partial \varphi }{\partial x_i}|}_{outlet}=0$E1.16

The correct use of that equation warrants the vortices can get near or cross the outlet boundary without significantly perturbing the computational solution inside the domain. In the last equation, x_i represents the main flow direction, and U_conv denotes the mean convective velocity at the outlet position, approximated using the outlet mass flow value [9].

3.1.4. Pressure outlet

The “pressure outlet” boundary condition is conventionally used for the solution of coupled problems (high Mach number values and compressibility effects), where the outflow boundary condition is not convenient. It requires to specify the static pressure (gauge pressure) at the domain outlet. Such value is exclusively respected for subsonic flow cases (M < 1). When the flow is supersonic, the specified pressure will not be used, and a new value will be extrapolated based on previous values.

Figure 2.

Pressure components of the pressure outlet boundary condition, adapted with permission from [10].

To use this boundary condition correctly, it is required to establish a set of flow relations for reversing flow at the outlet during the solution process (a set of backflow relations). Such relations are important because incorrect values for these parameters can cause convergence problems.

Specifically for the coupled solver in the option density based, pressure values at the outlet faces are calculated using a splitting procedure based on the AUSM scheme developed by Liou [11].

For subsonic flow conditions at the outlet, pressure is determined using a weighted average based on the left and right states at the domain boundary (Figure 2). This average value is a mix of polynomial adjustments of fifth order based on the Mach number values at the outlet face. Then, the pressure value is finally stated P_f = f(P_c, P_e, M_n) where P_c is the pressure at the pressure value at the inside neighbor cell next the exit face f, being P_e the specified pressure and M_n the Mach number at the normal direction. If there is a supersonic flow, the pressure at the exit face is extrapolated from the inlet cell value. For incompressible flow, the exit face pressure is calculated as an average between the outlets specified pressure and the average pressure value.

$P_f=\frac{1}{2}\left(P_c+P_e\right)$E1.17

The “pressure outlet” boundary condition does not guarantee a constant pressure along the outlet of the domain. However, once the solution of the model converges, the average pressure value at the outlet will tend to reach a value close to the static pressure imposed at the exit.

3.1.5. Walls

For viscous flows, the nonslip condition has to be imposed at the walls. This can be attained prescribing u= v = w = 0 at the walls, to reproduce the formation of the dynamic boundary layer. In some cases, when solving the flow within the boundary layer is not needed (e.g., in external flows over an immerse body, there is no need to solve the boundary layer on the walls that delimitate the domain), values of shear stresses equal to zero are imposed, which is known as slip walls [12]. In addition to prescribing values for shear stresses, temperature values (Dirichlet type) or normal heat flux values (Neumann type) can be imposed and they are particularly useful for thermal analysis. It is also possible to define wall roughness values, which are especially useful in turbulent flows studies, as well as translational and rotational movement of the walls related to the same frame of reference.

Figure 3.

Rotational periodicity of a cylindrical recipient, adapted with permission from [10].

Figure 4.

Translational periodicity in a heat exchanger array.

3.1.6. Periodic boundary condition

Periodic boundary condition is useful to cut down the size of the studied domain, saving simulation time in cases where the geometry of the given problem and the flow pattern have characteristics that are repetitive and periodic along a certain characteristic length “L”. The periodic flow behavior can be found in many applications such as heat exchanger channels, flow over pipe banks, and fully developed flow in pipes and ducts. There exist two types of periodicity, such as translational and rotational periodicity (Figure 3). For rotational periodicity a constant pressure is defined along the periodicity planes, also defining a rotational axis in the fluid in studies related to turbomachinery.

Regarding translational periodicity, a finite pressure drop can be defined between two periodic Planes (Figure 4). Also, a mass flow can be defined.

3.2. Some technical suggestions for an adequate solving

To solve adequately a CFD problem, it is recommended a series of best practices whose Application can help to obtain better convergence values for CFD simulations. The first one is to verify that the skewness value of most of the mesh elements is between 0 and 0.5, being the values closest to 0 the most desirable ones. This is to minimize the numerical dissipation for the algorithms responsible of carrying out the flow balances at the cell faces. Another important suggestion to improve initial convergence in fluid, heat transfer, and multiphase flow problems is to first obtain solutions using single precision and first-order discretization models. After that, based on these solutions, change the discretization models to second order models. This should be carried out considering that second-order schemes produce lower error values. From this point, it is recommended to activate the energy equation (if needed) or any multiphase flow or reactive flow parameters (if needed) to obtain adequate convergence. It is also needed to consider that final runs have to be performed using double precision. Finally, having previously obtained convergence, relaxation factor values can be increased to obtain a better solution in each iterative process. It is necessary to remember the importance of the results validation, using for that, experimental values of the same phenomenon in a controlled environment.

3.3. Pressure–velocity coupling

In the Navier-Stokes Eqs. (1.7)–(1.10), calculated velocities are related to the momentum transport in two or three dimensions. These velocities should satisfy the mass conservation equation. However, continuity equation does not depend on the pressure as it does the momentum exchange. The convective and dissipative terms of the momentum conservation equation, which depends on the pressure, have second-order, nonlinear terms. When governing equations have to be solved, it is necessary to find the velocity and pressure fields for stationary and time-dependent problems. If the flow is incompressible, continuity equation can be used as a transport equation for density, energy conservation equation represents temperature transport, and the existing pressure gradients can be obtained using equations such as the ideal gas law or any valid state equation. For incompressible flow, however, density has a constant value and there is not a direct link among the governing equations for the mass and momentum transport equations. This situation provokes unpredictable oscillations in the solving process making it considerably difficult.

Regarding the semi-implicit method for pressure-linked equations (SIMPLE) algorithm proposed by Patankar and Spalding in 1972 [13], it has been used for many years as a useful and convenient solving strategy. This algorithm is depicted in Figure 5. Unit mass flows go through the faces of the cells, and are evaluated from an initial guessed velocity value in two or three dimensions. From these velocity and pressure fields, a momentum transport equation solution is obtained. Then, from this solution, a corrected pressure is obtained using a series of relationships obtained from the continuity equation.

Finally, such corrected pressure values will be used for the transport and continuity equations solving. The whole process has to be repeated in an iterative manner until a convergence criterion is reached. As it was mentioned before, to reach a solution, it is necessary to initialize the flow using arbitrary values for the velocity and pressure fields. Such values are often selected from meaningful flow quantities such as free stream velocity and the atmospheric pressure for external flows. In the case of an internal flow, the main velocity value(the one with larger magnitude when several inlets are present) and the internal pressure are used.

Nowadays, there are a series of improvements made to the SIMPLE algorithm. One of them is the SIMPLEC algorithm also known as SIMPLE-Consistent [14], which is one of the most popular. This algorithm works similarly to SIMPLE, being its main difference that the momentum conservation equations are modified in a way that allows the velocity correction equations are simplified. The result of this change is that for some cases, the convergence time decreases. This algorithm is optimal to solve laminar flow problems without temperature gradients. It requires less convergence time than the conventional SIMPLE algorithm. Nevertheless, its use is not recommendable to solve turbulent flow, species transport, combustion problems, or any problem involving not only fluid movements, but also a more complex physical problem.

Figure 5.

Overall stages in the SIMPLE algorithm, adapted with permission from Versteeg and Malalesekera [3].

The pressure implicit with splitting of operators (PISO) algorithm is an improved version of the SIMPLE algorithm. It was created for non-iterative calculations in the numerical solution of non-steady flows. It has been adopted for iterative solution process for steady and non-steady flows, being especially useful for the last ones. This algorithm involves a stage where velocities and pressures are arbitrarily predicted, and two correcting stages. Therefore, it has an additional correcting stage with respect to SIMPLE and SIMPLEC. In the prediction and correction stages, the obtained pressures involve functions that contain within the continuity equation [15].

Even though the PISO algorithm requires additional memory to store the values of the added correction equations, and in general, it requires relaxation factors to stabilize the calculation process, it performs in a quick and efficient way.

In relation to periodic or transient problems, although less popular than PISO and SIMPLE, the fractional step method (FSM) algorithm is a method created as a solving alternative of non-iterative nature (NITA–non-iterative time advancement). Even though its use is less popular than the SIMPLE algorithm and its modifications, due to its formulation, involving a decoupling of the mass and momentum conservation equations using mathematical operators is capable of decreasing the simulation time [16].

Regarding the pressure–velocity coupling algorithm known as COUPLED, it is recommended to solve problems involving steady-state single-phase flows. For this algorithm, the solving process of all flow components is carried out using vectors. This means using a square coefficients matrix derived from flows entering or leaving the control volume through each control surface. This matrix is then algebraically associated with a unidimensional unknown vector. The product of both matrices is then equivalent to a solution vector where boundary conditions are included. The COUPLED algorithm is an alternative to SIMPLE and its spin-offs in the segregated solver, and it is recommended for cases where the domain presents low-quality elements for non-steady problems.

Finally, before starting a solution process, it is recommended to take a careful look into the most adequate coupling schemes for the problem at hand. This can be done considering the most important flow characteristics and looking for a similar case in any classic study case where a benchmark has already been performed (e.g., jets, flow over a flat plate, internal flows).

3.4. Turbulent flows CFD

Reynolds number (Re = ρ Uℓ/μ) is above a certain limit several events take place and cause that the flow behaves in a random manner, and its velocity components fluctuate along the three spatial directions. The flow will also present an unstable nature, promoting large-scale mixing and energy dissipation at small scales. Such regime is known as turbulent regime or turbulence. To include the effect of such fluctuations in the flow, it is necessary to substitute the flow variables that depend on the velocity vector $\overrightarrow{V}=f(u,v,w)$, as well as the scalar P for the sum of an average component $\overline{i}$ and is fluctuating counterpart “i′” having then:

$\overrightarrow{V}=\overrightarrow{V}+\overrightarrow{V^{\prime }};u=\overline{u}+u^{\prime };v=\overline{v}+v^{\prime };w=\overline{w}+w^{\prime };P=\overline{P}+p^{\prime }$E1.18

Such decomposition of the different flow variables has to be included in the governing equations, provoking the appearing of three normal stresses and three shear stresses, which, due to its units, are known as the Reynolds stresses, and are Represented as follows:

${{\tau }^{\prime }}_{xx}=-\overline{\rho {u^{\prime }}^2};{{\tau }^{\prime }}_{yy}=-\overline{\rho {v^{\prime }}^2};{{\tau }^{\prime }}_{zz}=-\overline{\rho {w^{\prime }}^2}$

${{\tau }^{\prime }}_{xy}={{\tau }^{\prime }}_{yx}=-\overline{\rho u^{\prime }{v}^{\prime }};{{\tau }^{\prime }}_{xz}={{\tau }^{\prime }}_{zx}=-\overline{\rho u^{\prime }{w}^{\prime }};{{\tau }^{\prime }}_{yz}={{\tau }^{\prime }}_{zy}=-\overline{\rho v^{\prime }{w}^{\prime }}$E1.19

The vector expression for the mass conservation equation and the Navier–Stokes equations are as follows:

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \overline{V}\right)=0$E1.20

$\rho \frac{D\overline{V}}{Dt}=\rho g-\nabla \overline{p}+\nabla \cdot {\tau }_{ij}$E1.21

where:

${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right)-\overline{\rho {u^{\prime }}_1{u^{\prime }}_j}$E1.22

Laminar Turbulent

It can be noted that the mass conservation equation is practically the same for laminar and turbulent flows. However, it is necessary to point out that this equivalence is valid only for the average values with a significant number of samples, meaning that the average of the fluctuating quantities becomes zero.

Because these stresses vary in their frequencies and sizes as a function of the distance from the walls, turbulent models take into account such variations to solve the closure problem of these extended equations.

4. RANS models

In 1972, Jones and Launder [17] proposed a model with two differential equations in addition to the mass and momentum conservation equations. In that model, an equation approximated the effect of the turbulent kinetic energy production rate $k=\frac{1}{2}\left(\overline{u{\text{'}}^2}+\overline{v{\text{'}}^2}+\overline{w{\text{'}}^2}\right)$ and the other one contained the effects of the turbulent kinetic energy dissipation rate “ε”. Together, the use of both equations proved to be adequate to obtain good results in turbulent flows whose interest was focused on the study of the influence of average turbulent effects. Both equations are as follows:

$\frac{Dk}{Dt}\approx \frac{\partial }{\partial x_j}\left[\frac{{\nu }_t}{{\sigma }_K}\frac{\partial k}{\partial x_j}\right]+{\nu }_t\frac{\partial {\overline{u}}_i}{\partial x_j}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_i}{\partial x_i}\right)-ϵ$E1.23

$\frac{Dϵ}{Dt}\approx \frac{\partial }{\partial x_j}\left[\frac{{\nu }_t}{{\sigma }_{ϵ}}\frac{\partial ϵ}{\partial x_j}\right]+C_1{\nu }_t\frac{ϵ}{k}\frac{\partial {\overline{u}}_i}{\partial x_j}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_i}{\partial x_i}\right)-C_2\frac{{ϵ}^2}{k}$E1.24

where σ_K and σ_ϵ are the “effective Prandtl numbers”, which relate the eddy diffusion of K and ϵ to the eddy viscosity: ν_t/ν_K and σ_ϵ = ν_t/ν_ϵ. Eddy viscosity is modeled as follows:

${\nu }_t\approx \frac{C_{\mu }{k}^2}{ϵ}$E1.25

The five constants of this model are obtained in an experimental way for cases of boundary layer calculations in non-detached flows. They are as follows:

$C_{\mu }=0.09C_1=1.44C_2=1.92C_{\mu }=0.09{\sigma }_K=1.0{\sigma }_{ϵ}=\text{1}.\text{3}$E1.26

Unfortunately, the above presented values cannot be used for any kind of flows (they are not universal). They need to be modified to model the behavior of wakes, jets, and recirculating flows. To deal with this difficulty, different and more sophisticated models have been proposed, which try to reproduce these flow patterns. Some of them are as follows:

1. 1.

   The k – ϵ model, d on the group renormalization, Yakhot y Smith [18].

2. 2.

   The k – ω model, Wilcox [19].

3. 3.

   The k – ϵ – v² – f model, Durbin [20].

5. The LES technique

Large eddy simulation is a technique based on the assumption that the largest eddies in a flow are of sizes relative to the characteristic length of the flow, whereas the smallest ones that form on the walls are isotropic and very similar for different kinds of flow, and therefore being more susceptible to be modeled. Under this consideration, the LES technique solve the largest eddies in an explicit way and the effects of the smallest eddies on largest eddies are modeled using a subgrid scale stress model. The justification for both assumptions is that the largest eddies contain the most part of the kinetic energy in a moving fluid, they transport the majority of conservative properties and they vary the most from flow to flow. Meanwhile, the smallest eddies are considered more “universal” and less important in the total flow; therefore, when they are modeled, it is expected that the error introduced is small. Unlike RANS models, LES technique uses a filtering process instead of an averaging process. It means that the information obtained for the largest scales is instantaneous and does not represent an average information. Instead of considering the influence of the velocity fluctuations equal to zero $\overline{{u^{\prime }}_{ij}}=0$ like in the RANS models, it has a value that locally influences the solutions, it means $\overline{{u^{\prime }}_{ij}}\ne 0$. In general, the use or the LES technique has produced better results in flow simulations when RANS models usually fail, like turbulent boundary layers subjected to adverse pressure gradients, the prediction of drag coefficients for immerse bodies in periodic regimes [21] and other unsteady flows. Governing equations for the LES model in an incompressible flow are as follows:

$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$E1.27

$\frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial {\overline{u}}_i{\overline{u}}_j}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \overline{P}}{\partial x_i}-\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{\partial }{\partial x_j}\left[\nu \frac{\partial {\overline{u}}_i}{\partial x_j}\right]$E1.28

$(i=1,2,3)$

In the Eqs. (1.27) and (1.28), ${\overline{u}}_i$ represents the filtered velocity components at different Cartesian positions x_i, where $\overline{P}$ represents the filtered pressure component, “v” and “$\rho$” represent the kinematic viscosity and density, respectively. Regarding the local influence of the Reynolds stresses, one of the most used models is the classic Smagorinsky model, which is expressed as follows:

${\tau }_{ij}=\overline{u_i{u}_j}-{\overline{u}}_i{\overline{u}}_j$E1.29

${\tau }_{ij}-\frac{1}{3}{\tau }_{kk}{\delta }_{ij}=-2{\nu }_t{\overline{S}}_{ij}=-2{\left(C_s\overline{\text{Δ}}\right)}^2\Vert S\Vert {\overline{S}}_{ij}$E1.30

where the strain tensor of the fluid in each cell is $\mathrm{\Vert S\Vert }=\sqrt{2{\overline{S}}_{ij}{\overline{S}}_{ij}}$, being ${\overline{S}}_{ij}$ the filtered strain tensor, v_t the subgrid scale turbulent viscosity, “$\overline{\Delta }$”the filter characteristic length. Generally, the size of the filter is proportional to the volume “Ω” of the cell $\overline{\text{Δ}}=2{(\Omega )}^{1/3}$, and the Smagorinsky constant C_s can be varied from 0.05 to 0.1 [22].

6. RANS and LES capacities for turbulent flows analysis

Originally, RANS models were conceived to include the turbulence average effect in the Navier–Stokes equations solution. Therefore, they are widely recommended to simulate steady-state fully developed turbulent flows, being especially good reproducing off-wall phenomena. When the problem needed to solve is highly depended on the adequate prediction of pressures and forces on the wall, it is recommended to apply functions for “enhanced wall treatment” and to have at least one nodal element positioned within the viscous sublayer (y⁺ < 5), with at least ten nodal elements along the whole boundary layer to achieve acceptable results [23]. Another use of RANS models is to provide preliminary results that allow to shorten the simulation time for unsteady problems based on the same technique or the use of LES. Regarding the kind of flows that RANS models have been successful reproducing, it can be said that internal flows, flows over boundary layers, turbulence prediction in a fluid that floats due to density changes, combustion coupled problems, two-phase flows, agitated wall induced turbulence, etc. have had successful results. As for problems where RANS models perform poorly, these are periodic or oscillatory flows, boundary layers affected by large pressure gradients, flow regimes that depend on separation point, and jets.

On the other hand, the LES technique is capable of giving average and instantaneous information about the main structures of turbulent flows without the complete dependence on the average modeling of Reynolds stresses found in the traditional RANS approach. This is advantageous to capture velocities and flow patterns in three dimensional problems. Among the successes of the LES model, it can be mentioned the prediction of separation points and fluctuating forces in external flows for geometries such as circular cylinders, eolic turbines blades, as well as structures and frequencies in jets and atomization systems, etc. Unfortunately, due to the large number of temporal and spatial scales that requires modeling for an acceptable performance, computer time and requirements are still high. The walls produce high demands since it requires the use of grids with nodal elements positioned at y⁺ ≈ 1, very short temporal steps and the size of the eddy scales that the smallest grid elements are capable of capture. This has made impossible the massive application of this technique for flow problems that have many complex geometries and many scales. It results evident that for compressible flows, the Reynolds stresses are more complex than for incompressible flow, being that SGS models for complex flows those that involve multiple phases and reactive flows found in combustions are still in a stage of development and validation [24].