Open access peer-reviewed chapter

A Numerical Procedure for 2D Fluid Flow Simulation in Unstructured Meshes

By António M. G. Lopes

Submitted: October 23rd 2015Reviewed: March 14th 2016Published: August 24th 2016

DOI: 10.5772/63077

Downloaded: 1165


The present work addresses the numerical simulation of fluid flow for 2D problems. The physical principles and numerical models implemented in the software package EasyCFD are presented in a synthetic and clear way. The 2D form of the Navier-Stokes equations is considered, using the eddy-viscosity concept to take into account turbulence effects upon the mean flow field. The k-ε and the k-ω Shear Stress Transport (SST) turbulence models allow for the calculation of the turbulent viscosity. The numerical model is based on a control volume approach, using the SIMPLEC algorithm on an unstructured quadrilateral mesh. The mesh arrangement is a non-staggered type. The coordinate transformation, integration discretization and solution method for the governing equations are fully described. As an example of application, the airflow around a NACA 0012 airfoil is calculated and the results for the aerodynamic coefficients are compared with available experimental data.


  • CFD
  • unstructured mesh
  • fluid flow
  • 2D simulation

1. Introduction

The software EasyCFD [1] is a 2D simulation tool aimed at an initiation in the field of computational fluid dynamics. The main guiding lines on its development were the simplicity and the intuitiveness of utilization, in a self-contained package. The physical domain is discretized with quadrilateral unstructured meshes, allowing the simulation to deal with complex geometries of any configuration virtually. The present work synthesizes the main physical principles and numerical models implemented for the solution of 2D fluid flow problems, including heat transfer, in arbitrarily shaped geometries.


2. Basic transport equations

For the description of the transport equations, the xand zcoordinates of the Cartesian system will be taken as the independent variables, to which correspond, for velocity, the uand the wcomponents. The coordinate transformation for the generalized mesh will be presented later on the present work.

The Navier-Stokes equations describe momentum conservation and, for a 2D situation, may be stated as follows:


where p[N/m2] stands for pressure and Irepresents the buoyancy forces. The diffusion coefficient includes the contributions of the dynamic and turbulent viscosity, i.e., Γ=μ+μt[N s/m2]

In turn, the conservation of mass law, or continuity equation, is


The energy conservation equation is obtained considering the transport equation for the enthalpy ϕ=cp T, where cp[J/kg K]is the specific heat and T[K] is its temperature. For a fluid medium, the corresponding equation is


with the source term STrepresenting the heat generation rate per unit volume [W/m3] and Γ=(μ/Pr+μt/Prt)cp, with Pr and Prt as the laminar and the turbulent Prandtl number. For solid regions, heat transfer is governed by conduction:


where λ[W/mK]is the material thermal conductivity. In a multi-component fluid flow situation, different fluids are present, sharing the same velocity, pressure, temperature and turbulence quantities. The concentration of one of the components (the secondary fluid) is computed with a normal transport equation for a scalar:


where ϕc2is the mass fraction, or concentration, of the secondary fluid. The diffusion coefficient is given by


where Dϕc[m2/s] is the kinematic diffusivity that characterizes the fluids pair and Sct is the turbulent Schmidt number. The concentration of the secondary fluid plays an active role in the flow field since fluid properties depend on the concentration of each component. The mixture properties (such as viscosity, specific heat, etc.) are determined by a weighted average of the properties of the components. Thus, for a generic property X, we have


where the constraint ϕc1+ϕc2=1holds. Exception is made for density, where a simple analysis may prove that the mixture density ρis given by


3. Turbulence modelling

Two of the most popular turbulence models are presented next.

3.1. The k-εturbulence model

The standard formulation of this turbulence model is described in, e.g., [2]. The turbulent viscosity is given by


The turbulence kinetic energy, k[m2/s2], and its dissipation rate, ε[m2/s3], are computed with the following transport equations:


The term Pkrepresents the production rate of kas the results of the velocity gradients:


while the term GT accounts for the production or destruction of kand εdue to the thermal gradients:


with g[m/s2] being the gravity acceleration and β[K−1] being the fluid dilatation coefficient. The remaining constants are


In the proximity of a wall, the previous equations should be modified to account for the viscous effects that become predominant. Wall functions ensure the connection between the viscous sub-layer and the inertia layer, at a location established by the y+ value:


where uτrepresents the friction velocity and y[m] is the distance to the wall. The utilization of wall laws does not require mesh nodes in the viscous sub-layer, which is a major advantage from the computational standpoint. The momentum flux per unit area along the direction normal to the wall is given by the wall shear stress, which is computed differently depending on the location of the wall neighbour node relatively to the transition between the viscous and the inertia sub-layers. Denoting by vthe generic velocity component parallel to the wall (which may, actually, be uor w) and by ythe distance to the wall, we have


where v0is the wall velocity, E*=9.793for smooth walls and χ=0.4187is the von Karman constant. In the wall neighbourhood, the production term given by Eq. (12) is computed assuming a Couette flow:


where, once again, vis the generic velocity component parallel to the wall and yis the generic distance to the wall. Due to its significant variations near the wall, εis averaged for the term ρεin Eq. (10):


For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. For the dissipation rate, Eq. (11) is not employed in the node adjacent to the wall. Instead, the dissipation rate is given by


As for momentum, the energy flux is computed differently depending on the y+ values. In this case, we have


where T0 is the wall temperature and P* is the so-called Jayatillaka function:


3.2. The k-ωSST turbulence model

The k-ωSST model [3] represents a combination of the k-εand the k-ωmodels. According to Menter et al. [4], the k-ωmodel is more accurate near the wall but presents a high sensitivity to the ωvalues in the free-stream region, where the k-εmodel shows a better behaviour. The k-ωSST model represents a blend of the two, through a weighting factor computed based on the nearest wall distance. The governing equations are


where ωis the frequency of dissipation of turbulent kinetic energy [s-1]. The production of turbulent kinetic energy is limited to prevent the build-up of turbulence in stagnant regions:


The weighting function F1 is given by




where yrepresents the distance to the neighbour wall and vis the laminar dynamic viscosity. F1 is zero away from the wall (k-εmodel) and changes to unit inside the boundary layer (k-ωmodel) with a smooth transition based on y. The turbulent viscosity is computed as


where Sis the invariant measure of the strain rate given by




The constants are computed as a blend of the k-εand k-ωmodels by the following generic equation:


The constants are α1= 5/9; β1= 3/40; σk1 = 0.85; σω1 = 0.5; α2= 0.44; β2 = 0.0828; σk2 = 1; σω2 = 0.856; β*= 0.09.

The near wall treatment for momentum and turbulence equations follows the proposal described in [5]. The basic principle behind the automatic wall functions is to switch from a low-Reynolds number formulation to a wall function based on the grid nodes proximity to the wall. According to these authors, the automatic wall treatment avoids the deterioration of results typical of low-Reynolds models when applied on too coarse meshes.

The known solutions for ωin the viscous (linear) and in the logarithmic near wall region are


The imposed value for ωat the first node close to a wall is


For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. In turn, for the momentum equations, a similar reasoning applies, with expressions for the shear velocity in the viscous and in the logarithmic region:


with U1 being the fluid velocity relative to the wall velocity. The wall shear stress is computed as follows:


4. Numerical method

4.1. Transformation of coordinates

Discretizationand integration of the transport equations described previously are performed using a non-orthogonal generalized mesh as shown in Figure 1. The independent Cartesian coordinates (x, z), describing the physical domain, are thus replaced by a boundary-fitted coordinate system (ξ, ζ), defined by the mesh lines which may have, locally, any orientation and inclination. The computational domain (as opposed to the physical domain) is the space considered in terms of the boundary-fitted coordinate system(ξ, ζ) as depicted in Figure 1. In the computational domain, the mesh spacing is considered, for convenience, as unitary, i.e., Δξ=Δζ=1, and the mesh lines are always horizontal (ξlines) or vertical (ζlines). The mesh arrangement is of the collocated type (as opposed to the staggered mesh) with the two velocity components and scalar quantities (temperature, pressure, turbulence kinetic energy and its dissipation rate or frequency, as well as concentrations) located at the control volume (CV) centre.

Figure 1.

The two domains. (a) Physical domain, showing the mesh lines; (b) Computational domain.

Transformation of the original equations is accomplished by replacing the independent variables, using the chain rule, which states that, generically:


The Jacobian of the transformation


represents the ratio between the physical size of the control volume and its computational size (unitary). The derivatives ξx,ξz,ζxζzare the contravariant metrics of the transformation. They are computed from the covariant metrics, xξ,zξ,xζ,zζas follows:


To obtain the strong conservative form in the boundary-fitted coordinate system (ξ, ζ), the transport equations are transformed through the application of the chain rule (35) and multiplied by the Jacobian of the transformation. The metric identity


is then used to recast some terms. The result, for a generic variable ϕ, is


The terms g11, g13and g33are the contravariant metric relations given by


The non-orthogonal term g13is null if the mesh is locally orthogonal. Uand W, in Eq. (39), are the contravariant velocities. The terms JρUand JρWrepresent mass fluxes through the control volume faces along the computational directions ξand ζ, respectively, and are computed as follows:


Note that, in this case, the sub-index for the fluxes (such as in Fξ) represents the flux direction (and not a derivative, such as for the metrics). Eq. (39) may be rewritten as


where the cross-derivatives were incorporated into the source term Scross. A similar procedure is applied to obtain the generalized form of remaining equations, leading to the following form for the Navier-Stokes and continuity equations:


4.2. Integration

The integration and solution method for the transport equations are entirely based on the methodology in [6], with some of the suggestions described in [7] and incorporating the necessary modifications for the generalized mesh approach. The general Eq. (43) may be written as


The integration of the previous equation in its control volume leads to


where the source term has been written as a linear combination involving ϕP. Fand Drepresent the advective fluxes and the diffusive coefficients, respectively:


In the previous expressions, the subscripts indicate the location relative to the CVcentre, in the computational domain, with the uppercase denoting neighbour nodes and the lowercase denoting neighbour faces: E,e: East; O,o: West; T,t: Top; B,b: Bottom. For simplicity, in Eq. (48), the subscripts ξand ζ  were dropped from the fluxes F(in fact, fluxes across the eastern and western faces are always along the ξdirection and fluxes across the top and bottom faces are always along the ζ  direction).

For the solution of the equations, it is necessary to evaluate the values of ϕin the CVfaces (i.e., ϕe, ϕo, ϕt, ϕb). These values are computed as a function of both ϕPand the values in the neighbour nodes ϕE, ϕO, ϕT, ϕBaccording to the adopted advection scheme (to be described later on). Eq. (48) may, then, be written in the following standard form:


or, in a more compact manner,


with “nb” indicating that the sum is to be performed for all the neighbouring locations.

It is necessary to compute the face values ϕe, ϕo, ϕt, ϕbas a function of ϕP, ϕE, ϕO, ϕT, ϕBin order to obtain the coefficients aE, aO, aT, aBof Eq. (51). It is known that a simple arithmetic average is not a physically plausible solution since, due to the presence of a flow, the property ϕ, being advected, tends to assume a value closer to the upwind value. Several advection schemes may be adopted, being the simplest one the upwindscheme. According to this scheme, the property ϕin the CVface takes the upwind value, i.e.,for example, Fe>0ϕe=ϕP; Fe<0ϕe=ϕE. The complete mathematical formulation for the upwind scheme is


Due to its first-order character, the upwind scheme is often not used due to associated numerical diffusion. The Quick scheme is third-order accurate. The deferred correction version of Hayase [8] combines the first-order upwind scheme with a third-order correction, bquick, added to the source term as follows:


The Quick scheme, although third-order accurate, presents oscillations that may lead to some unrealistic behaviour. Total variation diminishing (TVD) schemes, also implemented in the present code, were developed to provide second-order accurate solutions that are free or nearlyfree from oscillations. For further information on this, please refer to, e.g., [9].

4.3. Pressure-velocity coupling

EasyCFD adopts the SIMPLEC algorithm (Semi-Implicit Method for Pressure-Linked Equations-Consistent) [7], which is based on the original formulation SIMPLE[6]. Due to the non-staggered mesh arrangement (collocated mesh), the Rie-Chow interpolation procedure [10], with the modifications proposed in [11] and [12], is implemented. Let us consider the umomentum conservation Eq. (44). After integration, the evaluation of this equation leads to


During the iterative process, velocities u* are computed from the available velocity field umand pressure field p* obtained at the previous step as follows:


where Eis the under-relaxation factor [7]. During the iterative process, unless convergence is reached, the velocity field u*obtained from the solution of the momentum equations does not satisfy the continuity requirement unless the correct pressure field pis employed (instead of an incorrect pressure field p*). This means that, if the correct pressure field was employed, a mass-conservative velocity field would be obtained:


Thus, the pressure and velocity fields p* and u* should be corrected by a certain amount p' and u' :


Subtracting Eq. (56) from Eq. (55) and taking into account Eq. (57), one obtains


The keystone of the SIMPLECalgorithm consists on the subtraction of the term nbanbuP'on both sides of the equation and subsequent dropping, leading to


or, for simplicity:




The equations for the velocity correction are obtained through the pressure correction field, recurring to the previous equation:


leading to


4.4. The pressure correction equation

As previously stated, the objective of the pressure correction is to produce a pressure field such that the solution of the momentum equations is a mass-conservative velocity field. Consequently, the equations for solving the p’field must be obtained from the continuity equation. The discretized form of this equation is obtained directly from the integration of Eq. (2):


As one may see, velocities are, now, needed at the control volume faces. Taking Eq. (63), for the uvelocity component, at the “e“ and “o” faces and the same equation for the wcomponent at the “t” and “b” faces, substituting into Eq. (64) and rearranging the terms leads to


The terms gijare the covariant metric relations defined as


The derivatives are evaluated as


and, for the cross-derivatives,


Introducing the discretization expressed by Eqs. (68) and 67 into Eq. (65) allows us to obtain the pressure correction equation:




The ⟨ ⟩symbol denotes a linear interpolation from the control volume centre (where the momentum equations are defined) to the control volume faces (where the fluxes for the continuity equations are needed). The pressure correction field p’obtained from the solution of Eq. (69) is employed for correcting pressure and velocity corrections are obtained from Eq. (63). Note that, since these equations are defined at the control volume centre, p’derivatives are evaluated as


One may note that pressure values at the control volume centre are not included in the evaluation of these derivatives (the same applies for the pressure derivatives in the momentum equations). This may lead to the well-known checkerboard pattern for the pressure field. To avoid this effect, the Chie-Row interpolation method proposes that the mass fluxes, to be evaluated at the control volume interfaces for all the transport equations (Fe, Fo, Ftand Fb, in Eq. (48)), be corrected using the pressure correction p’field (instead of being computed from the corrected control volume centre velocities). The correction equations for fluxes are obtained from the correction equations for velocities:


The “starred” fluxes F*at the control volume interfaces are obtained by interpolating the momentum equation. The keystone of the method is that the pressure gradient is not interpolated, but, instead, is obtained directly from the pressure at contiguous control volume centres. Considering Eq. (55), evaluated in terms of fluxes, one may write


The terms F^e,oand F^t,brepresent the fluxes


where the revised velocities ûand ŵare obtained from the momentum equations as follows:


Note that the source term bincludes all the contributions (e.g., transient term, cross-derivatives and buoyancy for the wvelocity component), except the under-relaxation and pressure gradient.

4.5. Solution of the equations

The solution of the equations previously described is carried out with an iterative procedure. For accelerating the convergence rate, two relaxation factors (described next) are applied.

The solution of the equation is sub- or over-relaxed in the following manner:


Values of flower than unity lead to sub-relaxation, while values greater than unity over-relax the solution process. In the present code, the value f =1.1is employed for the transport equations, while, for the pressure correction equation, the Point Successive Over Relaxation (PSOR) method [13] combined with the additive correction multigrid method [14] is employed.

The whole flow field calculation is considered to be converged when all the normalized residuals are lower that a predefined value Rconv:


The total normalized residual for the transport equations of ϕ(ϕ=u,w,T,k,ε,ω) is determined as follows:


where (ϕmaxϕmin)quantifies the amplitude of the variable ϕin the calculation domain.

5. Examples of application

To exemplify the numerical method described previously, two test cases are presented next along with a comparison of results with published data.

5.1. Flow past an airfoil

A calculation was performed to compute the aerodynamic coefficients of a NACA 0012 airfoil operating at a Reynolds number of 6×106. The obtained values for the drag and lift coefficients are compared with existing experimental data. The first step is to define the calculation domain, which should be large enough in order to avoid numerical blockage effects. Figure 2 represents the domain, for a 1 m airfoil cord length. Lateral boundaries are assigned a free slip condition and a uniform velocity profile with 5% turbulence intensity is imposed at the inlet. A mass conservative condition is applied at the outlet boundary. After a mesh independency study, a total of approximately 250,000 mesh nodes were employed, with three mesh refinement regions. Particular care was taken near the airfoil surface, were y+values ranging from, typically, 0.1 to 6, with an average value of 1.7 all around, were obtained. Figure 3 depicts the mesh employed.

For the present simulations, both the first-order hybrid [6] and the Quick advection schemes were employed, along with the k-εand k-ωSST turbulence models. Experimental data are reported by Abbott and Von Doenhoff [15] and Ladson [16].

Figure 2.

Domain dimensions. Airfoil cord is 1 m.

Figure 3.

Non-structured quadrilateral mesh.

Results for the dependence of the lift coefficient with the airfoil angle of attack αare shown in Figure 4. As expected, the lift coefficient presents a linear dependence with the angle of attack αup to the onset of separation, which occurs at α=16. The two advection schemes give similar results up to separation, after which the lift drop in the stall region is more pronounced with the hybrid scheme. Separation is completely established at α=18and for α20the flow becomes unsteady. Both turbulence models perform very well in the linear region. After separation, data are more spread. The difficulty to obtain reliable experimental data in this circumstance is commonly recognized as the flow is unsteady and presents a 3D behaviour. Comparing the two advection schemes, the Quick scheme performs better, particularly after separation.

Figure 4.

Lift coefficient vs. angle of attack. Influence of (a) turbulence modeland (b) advection scheme.

Figure 5 depicts the relation between lift and drag coefficients. In this case, the k-ωSST turbulence model performs substantially better than the k-εmodel. This is certainly due to the fact that the friction component plays an important role in drag, and thus correctly resolving the boundary layer in the very proximity of the wall is crucial for the drag computation. It is also interesting to note that the advection scheme plays a very important role, with the higher order scheme Quick showing much better than the hybrid scheme, when results are compared with the experimental data.

Figure 5.

Drag coefficient vs. lift coefficient. Influence of (a) turbulence model and (b) advection scheme.

5.2. Natural convection inside a cavity

The natural convection flow in a cavity is a classical test for numerical methods in fluid dynamics. The cavity is a square shape (cf. Figure 6) with adiabatic horizontal walls. A constant temperature is imposed in each the vertical wall.

Figure 6.

Problem definition for the natural convection inside a cavity.

The problem is governed by the following non-dimensional parameters:


where the thermal diffusivity is α=λ/(ρcp), νis the kinematic viscosity, λis the thermal conductivity, ρis the density and cpis the specific heat:


where βis the thermal expansion coefficient, gis the gravity acceleration and ΔT=ThTcis the temperature difference between the vertical walls. The transition between laminar and turbulent flow takes place approximately for Ra=107. In the present work, simulations were conducted for Ra=105(laminar regime). Air is the operating fluid, for which Pr = 0.71. The hybrid advection scheme was used and the Boussinesq approximation was adopted. Computations were performed in a non-uniform grid, with 82 × 82 = 6400 nodes. Reference results are reported in [17] and [18] for several Rayleigh numbers in laminar regime, comparing solutions given by several authors. Results for laminar and turbulent flow are also presented in [19].

Figure 7(a) and (b) displays isothermal lines generated using a constant value spacing between the minimum and the maximum verified within the domain. Figure 8(a) and (b) shows the flow streamlines. The flow, in the steady-state situation, is characterized by a large vortex filling the cavity, rotating in the clockwise direction. Two small vortices rotating in the same direction are located near the cavity centre. For this case, the minimum and the maximum streamline values used in the visualization do not correspond to the total amplitude of the stream function within the domain. These values were, instead, adjusted in EasyCFD to correspond to those employed in [19]. The agreement between the calculations and those reported in the literature is very good. Vahl Davis and Jones [18] present normalized maximum values for the uand for wcomponents of velocity


occurring in the vertical and horizontal symmetry lines, respectively. Table 1 shows the results obtained with EasyCFD, the reference values in [17] and the range of variation for the 37contributions reported in [18]. This range does not include the minimum and maximum reported values since these clearly fall outside the general trend of the remaining contributions.

Figure 7.

Isothermal lines. Ra = 105. (a) EasyCFD and (b) Dixit and Babu [19].

Figure 8.

Streamlines. Ra = 105. (a) EasyCFD and (b) Dixit and Babu [19].

EasyCFDVahl Davis [17]Variation range [18]

Table 1.

Normalized maximum values for the uand wvelocity components.


6. Concluding remarks

The numerical simulation of fluid flow for 2D problems was addressed. The physical principles and numerical models here presented correspond to the implementation in the software package EasyCFD. Transformation of the original equations to cope with a non-orthogonal generalized mesh is described in detail, along with the coupling of momentum and continuity with an adapted SIMPLEC algorithm for non-staggered meshes. Although not addressed in the present chapter, this software package was developed entirely based on a graphical interface, aiming at an easy and intuitive utilization. With a fast learning curve, this package is very suitable for learning the principles and application methods in computational fluid dynamics and has a great value both as a didactic and an applied tool. Although, at first, the restriction to 2D situations may seem very limitative, a great number of practical situations may be addressed with this approach.

© 2016 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

António M. G. Lopes (August 24th 2016). A Numerical Procedure for 2D Fluid Flow Simulation in Unstructured Meshes, Numerical Simulation - From Brain Imaging to Turbulent Flows, Ricardo Lopez-Ruiz, IntechOpen, DOI: 10.5772/63077. Available from:

chapter statistics

1165total chapter downloads

1Crossref citations

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

On a New Numerical Approach on Micropolar Fluid, Heat and Mass Transfer Over an Unsteady Stretching Sheet Through Porous Media in the Presence of a Heat Source/Sink and Chemical Reaction

By Stanford Shateyi, Fazle Mabood and Gerald Tendayi Marewo

Related Book

First chapter

Some Commonly Used Speech Feature Extraction Algorithms

By Sabur Ajibola Alim and Nahrul Khair Alang Rashid

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us