Incompressible Turbulent Flow Simulation Using the κ-ε Model and Upwind Schemes

In the computation of turbulent flows via turbulence modeling, the treatment of the convective terms is a key issue. In the present work, we present a numerical technique for simulating two-dimensional incompressible turbulent flows. In particular, the performance of the high Reynolds κ-ε model and a new high-order upwind scheme (adaptative QUICKEST by Kaibara et al. (2005)) is assessed for 2D confined and free-surface incompressible turbulent flows. The model equations are solved with the fractional-step projection method in primitive variables. Solutions are obtained by using an adaptation of the front tracking GENSMAC (Tomé and McKee (1994)) methodology for calculating fluid flows at high Reynolds numbers. The calculations are performed by using the 2D version of the Freeflow simulation system (Castello et al. (2000)). A specific way of implementing wall functions is also tested and assessed. The numerical procedure is tested by solving three fluid flow problems, namely, turbulent flow over a backward-facing step, turbulent boundary layer over a flat plate under zero-pressure gradients, and a turbulent free jet impinging onto a flat surface. The numerical method is then applied to solve the flow of a horizontal jet penetrating a quiescent fluid from an entry port beneath the free surface.


Introduction
With the rapid advance of computer technology, numerical modeling has become an important tool in the understanding of fluid dynamics phenomena.One of the challenging tasks is the computation of incompressible turbulent flows (mainly those with free surfaces), which can, in principle, be carried out by the direct numerical integration of instantaneous Navier-Stokes equations.Unfortunately, due to the large computational effort involved, this technique has been restricted to flows at low Reynolds numbers.Practical calculations at the present time must therefore be based on the unsteady Reynolds averaged Navier-Stokes (URANS) equations, with the high Reynolds κ-ε turbulence model.However, the performance of this modeling decisively depends on the form that the nonlinear advective terms are approximated.
A wide variety of techniques for discretizing the nonlinear advective terms has been proposed over the last 20 years, given that the combination of NVD (normalized variable diagram) [1] and TVD (total variation diminishing) [2] is one of the most popular in the CFD community.For instance, [3] proposed the VONOS (variable-order non-oscillatory scheme), an NVD scheme which emerged, according to [4,5], as an acceptable upwinding tool for simulation of free surface flows.Reference [6] proposed a third-order accurate and limited scheme named WACEB (weighted-average coefficient ensuring boundedness) TVD.Numerical results for scalar convection problems show that this scheme has the same ability of QUICK in reducing numerical diffusion without introducing spurious extrema (oscillations).However, this scheme still has convergence problems for non-Newtonian flows.As a remedy, [7] devised a high-resolution scheme called CUBISTA (convergent and universally bounded interpolation scheme for treatment of advection) TVD.The evaluation of the accuracy and convergence properties of the scheme was measured in two-dimensional cases by using linear and nonlinear problems and for Newtonian and non-Newtonian flows.
In the last years, a great effort has been made to develop high-order bounded advection schemes that combine the TVD and NVD formulations.Using this combination, recently, [8] derived an upwind scheme for unsteady flow fields (called adaptative QUICK-EST), which did a very good job in solving laminar incompressible free surface flows (see [9] for details).The main motivation for the present work is to simulate incompressible turbulent free surface flows at high Reynolds numbers.By using the standard κ-ε turbulence model and the adaptative QUICKEST scheme, the present paper describes an effective 2D finite difference methodology for the numerical solution of this class of flows.The calculations are performed by the 2D version of the Freeflow simulation system of [10].
The paper is organized as follows.First, the model equations are set out (Section 2).The initial and boundary conditions are then presented (Section 3).The adaptative QUICKEST scheme is described in Section 4. The numerical technique is given in Section 5, while the finite difference discretization is described in Section 6. Numerical results are presented and discussed in Section 7. Conclusions are presented in Section 8.

Equation models
The flow regime of interest in this paper is modeled by the time-dependent, incompressible, constant property 2D Reynolds averaged Navier-Stokes equations, mass conservation equation, and κ-ε model in the primitive variable formulation.In conservative and nondimensional forms, these equations, omitting averaging symbols, can be written ) (2.5) In the above equations, t is the time, u = u(x, y,t) and v = v(x, y,t) are, respectively, the components in the x and y directions of the local time-averaged velocity vector field of the fluid, κ = κ(x, y,t) is the local time-averaged turbulent kinetic energy of the fluctuating motion, ε = ε(x, y,t) is the turbulence dissipation rate of κ, p e = p + (2/3)(1/Re)κ is the effective scalar pressure field divided by the density, and g = (g x ,g y ) is the acceleration due to gravity.The nondimensional parameters Re = UL/ν and Fr = U/ L|g| denote the associated Reynolds and Froude numbers, respectively, in which U is a characteristic velocity scale and L is a length scale of the flow.The nondimensional turbulent viscosity ν t , turbulent shear stress production P, and turbulence time scale T t are, respectively, defined as ) ) where variables with a bar refer to their corresponding dimensional variables.

Initial and boundary conditions
Equations (2.1)-(2.5)are coupled nonlinear differential equations and, together with the eddy viscosity model (2.6), are sufficient, in principle, to solve for the five unknowns u, v, p e , κ, and ε when appropriate initial and boundary conditions are specified.In this work, a staggered grid is used where the effective pressure, the turbulent kinetic energy, and the dissipation rate are stored at the centre of a computational grid cell, while velocities are stored at the cell edges.A typical cell showing the physical locations at which these dependent variables are defined is illustrated in Figure 3.1.With this grid system, effective pressure boundary conditions are not needed.The boundary and initial conditions have been implemented as follows.
The initial conditions for the mean velocities and effective pressure are specified in the same way as in the laminar case (see [4] or [9]), that is, these variables are prescribed.It is difficult to specify initial conditions for the turbulent variables, since they must be in agreement with the physics of the problem.Thus, for the problems considered in this work, we prescribe the initial conditions for κ and ε, and hence ν t , as functions of an upstream turbulent intensity I = 8.0 × 10 −2 .This specific intensity level is used because it V. G. Ferreira et al. 5 is within the bounds of realistic values.In nondimensional form, the turbulent variables can be written as Four types of boundary conditions have been implemented, namely: inflow, outflow, free surface, and rigid-wall boundaries.At the inflow, the velocities u and v are prescribed while the values of κ and ε are estimated in such a way that they are consistent with the initial conditions (3.1).At the outflow, the streamwise gradient for each variable is required to be equal to zero.At a free surface, we consider the fluid to be moving into (or out of) a passive atmosphere (zero-pressure) and, in the absence of surface tension forces, the normal and tangential components of the stress must be continuous across the free surface; hence on such a surface, we have (see, e.g., [11]) Here, n and m are, respectively, unit normal and tangential vectors to the surface, p ext is the external (atmospheric) pressure (assumed zero in this paper), and τ = τ(p e ,ν t ,u) is the Cauchy stress-tensor given by τ = −p e I + Re −1 1 + ν t ∇u + (∇u) T , (3.4) where I denotes the identity tensor.Equations (3.2) and (3.3) are discretized by accurate local finite difference approximations on the free surface, namely, from condition (3.2) one determines the effective pressure; and from (3.3) one obtains the velocities at the free surface.Due to the complexity of the dynamics of the turbulence near the free surface, the values of the turbulent variables κ and ε at the boundary are difficult to specify.For instance, it is not known how turbulence interacts with surface tension, and therefore, it is difficult to specify the distribution of κ on an irregular moving boundary.So, as a first approximation, we assume that the free surface is locally flat and the movement of the fluid does not cause any discontinuities at the boundary.In summary, the turbulent variables at the free surface are determined by imposing The derivatives in (3.5) are approximated by first-order (either forward or backward) finite difference schemes.
The κ-ε model, as formulated in (2.4)-(2.5),cannot be applied as the calculation approaches a rigid wall.This is because the turbulent time scale in (2.5) exhibits a singular nature near a wall, since the turbulent kinetic energy κ tends to zero there [12].For this reason, the wall function is employed in the near-wall region.In this case, the fundamental equation for determining the fictitious velocities and turbulent variables near a rigid wall is the total momentum flux τ ω given by [13] 1 where u represents the mean velocity component tangential to the rigid wall, and u * is the friction velocity.The values of κ and ε in the inertial sublayer are, respectively, prescribed by the well-known relations where K = 0.41.In the viscous region close to the wall, we use the strategy of [14], that is, where y + is defined as y + = Re u * y, and l * represents the length scale proposed by [15].
Neglecting the buffer layer of the turbulent boundary layer, the critical value of y + (denoted by y + c ) in (3.8) separates the viscous sublayer from the inertial sublayer.

Wall boundary conditions.
The strategy adopted here to describe the solution of the flow near a rigid wall is the wall function which describes the asymptotic behavior of the turbulent variables near the wall.The main advantages of the wall function approach are (a) the need to extend the calculation right down to the wall is avoided, a fact which saves computing time and storage; and (b) it is not necessary to account for the viscous effects in the turbulence model.In summary, the behavior of the mean velocity profiles in the viscous and inertial sublayers is given by (see, e.g., [16,17] or [18]) ln where u + = u/u * and E = exp(KB); B is an empirical constant and is usually chosen to correspond to a hydrodynamically smooth wall.One of the central questions in the application of the wall functions (3.6)-(3.8)and (3.9)-(3.10) is the accurate determination of the friction velocity, and hence the wall shear stress.This is determined from relation (3.9) or (3.10), depending on the local Reynolds number y + .The Newton-Raphson method was used to obtain u * from (3.10), with u * = 11.60 as initial condition.This initial value was obtained from the numerical solution of the system (3.9)-(3.10)(see also [19]).To begin with, we need to know the critical Reynolds number y + c in (3.8).By neglecting the transition sub-layer, the friction velocity is estimated in the following specific way: with the tangential velocity u * known in the first grid cell adjacent to the wall, u * is updated according to the value of y + given by (3.9).If y + is less than y + c , we use (3.9); on the other hand, if it is not, we employ (3.10).The fictitious velocities are calculated by the central-difference approximation of (3.6) for a known wall shear stress.

Treatment of nonlinear advection terms
For the flow regime considered in this work, the momentum and turbulence transport equations are dominated by convection, and it is well recognized that standard discretization (i.e., QUICK, central difference, or Lax-Wendroff) for the nonlinear terms leads to oscillatory solutions.In this paper, the discretization of all advective terms in transport equations (2.1)-(2.2) and (2.4)-(2.5) is performed by using the adaptative QUICKEST scheme of [8] (see also [9]).This high-order upwind technique is derived from the normalized variable of [20] and by enforcing the total variation diminishing property of [2,21].Consequently, it satisfies the CBC of [22].The main idea in the derivation of this scheme was to combine accuracy and monotonicity, while ensuring flexibility (it depends on a free parameter).The adaptative QUICKEST scheme enjoyies the property that total variation of the variables does not increase with time, thus spurious numerical oscillations are not generated.The numerical solution can be second-or third-order accurate in the smooth parts of the solution, but only first-order near regions with large gradients.In summary, a general interfacial flow property φ f is implemented in the current Freeflow code by the functional relationship where θ = V f • δt/δx is the convective Courant number, V f is a convective velocity, and δx is the grid spacing, and [1]).The subscripts D, U, and R referring to values at the downstream, upstream, and remote-upstream locations are defined according to the sign of V f at f face (see Figure 4.1).The constants a and b in (4.1) are given by The corresponding flux limiter for the adaptative QUICKEST scheme is as follows (see details in [9]): where The implementation of this scheme will be presented later (see Section 6).

Numerical solution procedure
The governing equations (2.1)-(2.5)are solved in a partly segregated manner using an extension of the GENSMAC methodology of [23] for turbulent flow fields.A detailed description of this technique is provided in [24].Based on the predictor-corrector method (see, e.g., [25]), the numerical solution procedure is an explicit finite difference, first-or second-order accurate numerical method for calculation of free surface flows as well as confined flows.
For calculations, a uniform Cartesian staggered grid system is used, where the effective pressure, the turbulent kinetic energy, and the dissipation rate are stored at the center of a computational grid cell, while velocities are stored at the cell edges.For flows possessing free surface, this boundary generally moves, and therefore the domain of interest deforms with time.In this case, the front-tracking MAC (marker and cell) method [26] is adopted in Freeflow to determine the free surface location.In summary, the interface is represented discretely by connected Lagrangian markers to form a front which lies within and moves through an Eulerian mesh; as the front moves and deforms, interface points are added/deleted and reconnected as necessary (for details, see [10]).To advance the numerical solution in time, the projection method of [27] is employed (see also [25]).

Discretizations
The differential equations are discretized using the finite difference technique on a uniform staggered grid system.The temporal derivatives were discretized using the firstorder forward difference (Euler's method), while the spatial derivatives were approximated by standard second-order central differences with the exception of the advection terms (denoted here by CONV(•)), which are approximated by the adaptative QUICK-EST scheme.The Poisson equation in discretized scheme (see [9]) is done by using the V. G. Ferreira et al. 9 usual five-point Laplacian operator, and the corresponding symmetric-positive definite linear system is solved by the conjugate-gradient method.In summary, fluid flow equations (2.1)-(2.5)take the following discretized form.
(i) Momentum-u: where the superscript n denotes the time level, and (iv) Poisson equation for ψ: where (vi) ε-equation: The production of turbulence, the eddy viscosity, and the time scale are discretized, respectively, as follows: For the nonlinear advection terms in the momentum equations (the advection terms of κ and ε equations follow a similar procedure), the application of the adaptative QUICK-EST scheme is as follows.For simplicity, only the discretization of the nonlinear terms in u-component of the time-averaged Navier-Stokes equations will be presented.The discretization of the other nonlinear term is made in a similar way.In position (i + 1/2, j) of the mesh, this term can be approximated by the following conservative scheme: ≈ u i+1, j u i+1, j − u i, j u i, j δx where the velocities u i+1, j , u i, j , v i+1/2, j+1/2 and v i+1/2, j−1/2 are obtained by averaging.For instance, v i+1/2, j−1/2 is approximate by v i+1/2, j−1/2 ≈ 0.5(v i, j−1/2 + v i+1, j−1/2 ).The velocities u i, j and u i+1, j are calculated (the other velocities follow similar procedures) by the following.(i) When u i, j > 0 and (6.13) (ii) When u i, j < 0 and u i+1/2, j = (u i+(1/2), j − u i+(3/2), j )/(u i−(1/2), j − u i+(3/2), j ), the value of u i, j is (6.15) (iv) When u i+1, j < 0 and u i+3/2, j = (u i+(3/2), j − u i+(5/2), j )/(u i+(1/2), j − u i+(5/2), j ), the value of u i+1, j is where The Courant number is calculated in the code by analyzing the direction in which information propagates, that is, the sign of a previously calculated local normal averaged velocity u at a face of the control volume, and by computing the expression θ = u • δt/δx in each control volume.

Numerical tests
In order to validate the actual Freeflow code, incremented with the original κ-ε model and adaptative QUICKEST scheme, we report now the numerical results for a turbulent flow over a backward-facing step, the turbulent boundary layer over a flat surface, and a turbulent jet impinging onto a flat surface.Also, as an example of application, we present a turbulent planar jet penetrating into a pool.

Turbulent flow over a backward-facing step.
The turbulent flow over a backwardfacing step is a standard test case, often used for validation of turbulence models.This flow is computationally challenging, because both a primary and a secondary recirculation eddy vertex occur.The problem configuration is illustrated in Figure 7.1.By using the current Freeflow code with a fully developed parabolic velocity profile prescribed at the inlet section, we simulate this flow at Re = 1.32 × 10 5 .This is based on the maximum velocity U max = 1.0 ms −1 at that section and the height of the step h = 0.1 m.Computations were performed on three different meshes, namely, the coarse mesh (200 × 15 computational cells, δx = δ y = 0.02 m); the medium mesh (400 × 30 computational cells, δx = δ y = 0.01 m); and the fine mesh (800 × 60 computational cells, δx = δ y = 0.005 m).Table 7.1 depicts values of the reattachment length x 1 on the three meshes, using four advection schemes, including the adaptative QUICKEST.By comparing these numerical results with experimental data of [28], which obtained value of x 1 = 7.1, one can see that our numerical results underpredict the experimental reattachment point by 20%-25%.The CUBISTA scheme for the coarse mesh provides a value greater than 7.1.The good results with WACEB on coarse meshes can be attributed to the value of the y + .On the other hand, it can be seen that our numerical results are in good agreement with the numerical result of [29], which found that x 1 = 6.0.From this same table, it should also be observed that WACEB, CUBISTA, and adaptative QUICK-EST schemes provide good results, while the VONOS scheme gives much less satisfactory results.We believe that most of this difference may be attributed to the fact that the WACEB, CUBISTA, and adaptative QUICKEST schemes are TVD, whereas the VONOS scheme is not.For simple comparison, Table 7.2 shows other estimates for the reattachment length obtained by HOAB, STOIC, SMART, CLAM, and FOU schemes (see [30]).From this table and Table 7.1, one can observe that the adaptative QUICKEST scheme provided a consistent reattachment length.In addition, a convergence test of the numerical solution obtained with the adaptative QUICKEST scheme on these three meshes was made.This is illustrated in Figure 7.2, which shows how the reattachment length was estimated (the change in the sign of the u-velocity profile adjacent to the lower bounding wall) in the code.
One can see from this figure that both the velocity profile and the reattachment length tend to converge to a solution near the numerical one in the fine mesh.For illustration, Figures 7.3 and 7.4 present the pressure contours and v-component of the velocity field, on the medium mesh, using adaptative QUICKEST scheme.

Turbulent flow past a flat plate.
The interaction between the fluid and the boundary wall is of great importance in turbulent flows.Due to the strong velocity gradients occurring near the wall, a large amount of turbulence is generated.This turbulence plays a very important role in physical phenomena as reattachment of separated regions.In −0.156 0.009 0.174 0.339 0.504 0.669 0.834 0.998 this subsection, a two-dimensional turbulent boundary layer over a flat plate is simulated according to the classical zero-pressure gradient theory.This is justified by the fact that the plate is flat, and thus, the only contribution to the pressure gradient comes from the product between the dynamic viscosity and the second derivative of the longitudinal velocity with respect to the transversal coordinate which, in the present case, can be neglected.Physically, in a zero-pressure gradient turbulent boundary layer, the point of inflection is at the wall itself; there can be no flow separation [31].
This fluid flow problem has been extensively studied in the literature (see, e.g., [5]), and numerous formulae have been proposed to estimate the coefficient of skin friction (C f ).In order to solve this problem, a uniform free stream boundary condition is imposed at the inlet, and the Reynolds number, based on length and velocity scales of unity, is Re = 2.0 × 10 6 .Figure 7.5 compares the calculated dimensionless turbulent skin friction coefficient C f = 2τ w with the estimates given by Prandtl approximation, Power-law theory, and the "exact" profile of White (see [18]).These figures display C f against the local Reynolds number Re x = U 0 x/ν at the (nondimensional) time t = 6.0 calculated for the following three different-sized meshes, namely, the coarse mesh (20 × 100 computational cells, δx = δ y = 0.05 m); the medium mesh (40 × 200 computational cells, δx = δ y = 0.025 m); and the fine mesh (80 × 400 computational cells, δx = δ y = 0.0125 m).Additionally, the corresponding laminar result is also included for simple comparison.As shown in Figures 7.5(a), 7.5(b), and 7.5(c), the numerical estimates are generally satisfactory for Re x beyond 1.0 × 10 6 .It can also be observed from Figure 7.5(d) that when the coarse mesh was twice refined, there appears to be convergence of the numerical solution to a profile near the power-law theory and the "exact" White relation.On the other hand, for Re x ≤ 1.0 × 10 6 , a systematic discrepancy existed and this may be due to the uniform meshes used and/or the initial velocity profile not being sufficiently turbulent at the entrance region.the free surface boundary conditions must be specified on an arbitrarily moving boundary.This free surface flow in turbulent regime is also chosen as a representative test case because there is (see [4,32]) an approximated analytical solution for the total thickness of the fluid layer flowing on the surface (see the illustration in Figure 7.6).In summary, for a given volumetric flux Q through the inlet section of diameter L = 2a, the analytical solution is where x 4/5 , In (7.1),A = 0.239 and k = 0.260.The problem configuration is illustrated in Figure 7.6.By using three different meshes, namely, the coarse mesh (200 × 50 computational cells, δx = δ y = 0.001 m); the medium mesh (400 × 100 computational cells, δx = δ y = 0.0005 m); and the fine mesh (800 × 200 computational cells, δx = δ y = 0.00025 m), the Freeflow code, equipped with the adaptative QUICKEST advection scheme and κ-ε model, run this moving free boundary problem at Reynolds number 5.0 ×10 4 , which was based on the maximum velocity U max = 1.0 m/s and diameter of the inlet L = 0.01 m (or volumetric flux Q = ν, Re = 0.01 m 2 /s).On these three meshes, a comparison is made between the free surface height (the total thickness of the layer), obtained from our numerical solutions and from the analytical viscous solution of Watson.This is displayed in Figure 7.7 and its enlargement in Figure 7 Coarse mesh Medium mesh Fine mesh Analytical numerical results on these meshes are similar, showing, in some regions, a small difference when compared to Watson's solution.We believe that most of this difference may be attributed to insufficient grid points, being used near the rigid wall.For illustration, Figures 7.9-7.11present at the time t = 1.0 the pressure and velocities fields on the medium mesh using the adaptative QUICKEST scheme.7.4.Application: a horizontal jet penetrating a quiescent fluid.We conclude this work by presenting the numerical simulation of a horizontal jet penetrating a quiescent fluid from an entry port at depth H = 6.0 m beneath the free surface.The purpose here is to show that the actual Freeflow can simulate the largest eddies present in the flow and their nonlinear interaction with a free surface.This free surface flow problem has also been simulated by [5] using a classical upwind scheme.The geometrical configuration and parameters for this free surface fluid flow are shown in Figure 7.12.In this numerical simulation, the associated Reynolds and Froude numbers are Re = DU 0 /ν = 5.0 × 10 4 and Fr = U 0 / gD ≈ 12.77, respectively.The mesh used in this test case is 100 × 100 computational cells (δx = δ y = 0.010 m).The development of pressure and velocities distributions, together with the free surface elevation at various times, are presented in Figures 7.13 through 7.15.In this case, the interaction with the free surface occurs only at the later stages of the flow development.Initially, one can observe the growth of the instability of the boundary layers between the entering jet and the stagnant fluid, and subsequently, the formation of a pair of counter-rotating eddies.Later on, the first pair of eddies propagates towards the free surface.
The result obtained in this simulation should be interpreted as representing the 2D motion of one realization occurring at scales greater than the discretization scale.In other words, both the free surface position and the velocity fields computed here may be regarded as the deterministic motion at these larger scales.Indeed, this simulation may be 20 Mathematical Problems in Engineering thought of URANS or as VLES (very large scale simulations), as opposed to the 3D, and much more expensive, LES.Turbulent flow simulations using LES and DNS approaches have been performed by other authors, but are mostly restricted to very low (or negligible) Froude numbers.

Conclusions
In this work, a finite difference numerical technique for simulating 2D incompressible turbulent flows was described.The Freeflow simulation system coupled with the original high Reynolds κ-ε turbulence model and the high-order adaptative QUICKEST advection scheme has been applied to simulate three problems, namely, turbulent flow over a backward-facing step, the turbulent boundary layer over a flat surface (zero-pressure gradient case), and a turbulent free jet impinging onto a flat rigid wall.According to the computed results, the new version of the Freeflow code can, in fact, predict both confined and free surface turbulent flows with satisfactory accuracy.In order to illustrate the robustness and applicability of the code to compute the interactions between a free surface V. G.

Figure 7 . 2 .
Figure 7.2.Comparison on three meshes of u velocity component using the adaptative QUICKEST scheme.

Figure 7 . 5 .Figure 7 . 6 .
Figure 7.5.Comparison of the local skin friction on a flat plate for turbulent flow, showing several theoretical estimates and those obtained by the present finite difference scheme on three meshes: (a) coarse; (b) medium; (c) fine; and (d) shows comparison of the three numerical solutions.

Figure 7 . 13 .
Figure 7.13.Evolution of the pressure contours of a jet in a fluid portion.

Table 7 .
2. Other estimates for the reattachment length.

7.3. Turbulent jet impinging onto a flat surface
. A jet impinging normally onto a flat rigid surface is a good example of a free surface flow, but is difficult to simulate because 88.One can see from these figures that the Ferreira et al.21