Numerical simulation of turbulent flow through a straight square duct using a near wall linear k – ε model

Ahmed Rechia1, Hassan Naji2,*, Gilmar Mompean2, and Abdelatif EI Marjani3 1Faculté des Sciences et Technologies de Tanger, B.P 416 Tanger, Maroc ahmed.rechia@polytech-lille.fr 2Université des Sciences et Technologies de Lille, Polytech-Lille, LML UMR CNRS 8107, Boulevard P. Langevin, 59655 Villeneuve d’ Ascq, France hassan.naji@polytech-lille.fr gilmar.mompean@polytech-lille.fr 3Ecole Mohammadia d’Ingénieurs, Laboratoire de TurboMachines Avenue Ibn Sina B.P 765 Agdal, Rabat, Maroc marjani@emi.ac.ma


INTRODUCTION
The prediction of turbulent flows involving secondary motions and using numerical simulations of the Reynolds Average Navier-Stokes (RANS) equations has great practical value in fluid mechanics and many applications can be found in centrifugal machinery design.In these simulations, the aim is to obtain accurate values of pressure and velocity for the flow field.Among the various RANS models for turbulence, the two-equation turbulence k − ε models have found their broad applications of feasibility in the majority of engineering practice for the predictions of the complex turbulent flows.By far the most popular turbulence model is the standard k − ε model which is based on the linear stress-strain relation initially proposed by Boussinesq.This type of closure has been revealed robust and efficient with respect to CPU time than more highorder models [1,2].
In this work, the numerical investigation of a low Reynolds number turbulent flow through a straight square duct is carried out in order to describe its turbulent characteristics.In spite of appearances, this flow is difficult to carry out numerically and/or experimentally.This configuration has been frequently chosen by many authors [3,4,5,6,7] since it is a relatively simple geometry which provides a good case test to improve turbulence models.Note that this flow is anisotropic and involves a secondary flow in the cross-stream plane, which is absent in the case of a plane channel.This secondary flow convects mean flow momentum from the centre of the duct to the corners.This causes a bulging of the streamwise velocity contours towards the corners.By studying the origin of the secondary flow for the square duct, Mompean [8] and Huser et al. [9], using DNS data of this problem, concluded that the streamwise vorticity is responsible for the generation of the secondary flow.The equation of mean vorticity in the streamwise direction, , reads as: (1) where V and W are the mean spanwise and normal velocities, v 2 and w 2 are the turbulent normal stresses, vwis the cross-stream correlation and v is the Kinematic viscosity.
From this formulation, the three terms on the right side of the above equation can be interpreted as the vorticity production (which is linked to the anisotropy of the cross-stream normal stresses), the production term due to the cross-stream correlation vwand the vorticity viscous dissipation.As the Reynolds number of the flow increases, the viscous dissipation term tends to be negligible and other terms become relevant.According to these authors, it seems therefore possible that the streamwise vorticity has its origin in the inequality of the normal Reynolds stresses combined with a presence of lateral wall gradient.In other words, the vorticity production term, , plays a crucial role in the generation of secondary flows.
It is well known that traditional two-equation models are incapable of describing such flows.While, traditional two-equation models, the standard k − ε model, fail to predict the near wall behaviour and anisotropy of the normal Reynolds stresses, it is shown in this paper that these standard k − ε eddy viscosity models can be improved via the introduction of damping functions of Van Driest type to produce the anisotropy.Previous studies have been carried out and provided interesting results, see Gavrilakis [10], Huser & Biringen [11] and Xu and Pollard [4] for a low Reynolds number turbulent flow.In this work, the predicted results are compared with the DNS data of Gavrilakis [10], the LES prediction of Xu [4] and the experimental data of Nishino [7] at a Reynolds number 4800 based on the duct hydraulic diameter and on the bulk velocity.As the turbulent characteristics are in good agreement with the DNS and experiment, this provides the motivation for the study of this flow from further models as non-equilibrium models [12] and nonlinear turbulence models [13,14].
The paper is organised as follows: the turbulence model k − ε and the introduction of damping functions are briefly reported within section 2; the finite volume numerical method is presented in section 3.In section 4, the main results of this work are presented, discussed and finally some conclusions are drawn.

GOVERNING EQUATIONS
To predict the turbulent flow, the statistical approach is used applying the Reynolds decomposition, which consists in splitting velocity and pressure into an average and a fluctuating part.As in all studies of Reynolds stress modelling, any flow variable φ can be decomposed as follows: ( where φ -is the mean turbulent value and ϕ is the fluctuating component.
The equations governing the mean velocity U i and the mean pressure P are obtained from the RANS equations for an incompressible flow: ( where ρ is the fluid density and is the Reynolds stress tensor.In order to achieve closure of this system of equations, a Reynolds stress model that lies τ ij to the global history of the mean velocity field must be supplemented.A variety of expressions have all independently established that, to the lowest order, τ ij can be represented by the standard eddy viscosity form [15,16]: (5) where k is the turbulent kinetic energy, v t is the eddy viscosity, is the mean strain rate tensor and δ ij is the Kronecker tensor.
Note that in this approach, the problem was reduced from deriving equations describing the evolution of τ ij to deriving equations for k with an appropriate eddy viscosity v t .In this study, the eddy viscosity is expressed as: (6) where ε is the dissipation rate of the turbulent kinetic energy k and C µ is a model constant.
In the two-equation k -ε context, the turbulent kinetic energy transport equation is derived from the contraction of the Reynolds stress transport equation, and is given by: (7) where the right-hand side represents the transport of k by the combined effects of turbulent transport and viscous diffusion , the turbulent production and the isotropic turbulent dissipation ε.
As for the transport equation of the dissipation rate ε, a rather general expression from which many of the forms can be derived.The modelling of this equation is based on analogies with the k equation and on phenomenological considerations.In its most commonly used form, this equation can be written as [17]: (8) The Eqns. ( 7)-( 8) contain five constants C µ , C ε1 , C ε2 , σ k and σ ε whose values are usually obtained from calibrations with homogeneous shear flows and from the decay rate of homogeneous and isotropic turbulence.The most commonly used values of these constants are: ; ; ; ; . This k -ε model is most commonly used in turbulent shear flow computations with a great deal of success.In the two-equation k -ε models, Eqns.( 7) and ( 8) in their high-Reynolds number form, do not provide the correct asymptotic behaviour in the wall region.In order to handle near wall effects, damping functions f µ are associated with the eddy viscosity definition (see for example [18,19]): with (10) where the constants 'a' and 'b' can take various values for each element of the Reynolds stress tensor τ ij .Previous works have been done using a priori test and interesting results have been obtained, see references [20] and [21].It is this assumption that we adopted in this paper.These corrections allow to incorporate viscous molecular diffusion in the vicinity of solid boundaries.The variables y + and z + are the wall co-ordinates based on the friction velocity and the cinematic viscosity v, y + = u τ y/v and z + = u τ z/v.

NUMERICAL METHOD
The method used to solve the set of equations ( 3)-( 4)-( 7)-( 8) is the classical finite volume method (FVM).The conservation equations are integrated over a control volume, and the Gauss theorem is used to transform the volume integrals into surface integrals.

SPACE DISCRETIZATION
To describe the numerical algorithm, the RANS equations can be written in the form: (11) where can be any of the above variables and the overbar represents the ensemble average.This is the so-called transport equation for property .It clearly highlights φ the various transport processes.For a Cartesian system , the expressions of the various terms of the equation ( 11) can be given by (see Table 1).The key step of the finite volume method is the integration of the transport equation ( 11) over a control volume of boundary γ (see Figure 1).Then, applying the Gauss divergence theorem, the resulting equation can be written as follows: The equation (12) gives the assessment of any scalar flow property φ in the control volume Ω k .This equation is the base of the space discretization of FVM.There is no need for all variables to share the same grid; a staggered mesh may turn out to be advantageous.This arrangement is shown in Figures 2 and 3. To solve the differential equation, a finite volume 322 Numerical simulation of turbulent flow through a straight square duct using a near wall linear k -ε model numerical method has been employed.The pressure, the turbulent kinetic energy, k, the dissipation, ε, and the normal Reynolds stress components are treated in the center of the control volumes; the velocities are computed in the center of the faces and the cross components of the Reynolds tensor are attached to nodes located at the mid-edges.The biggest advantage of the staggered arrangement is the strong coupling between the velocities and pressure.This helps the avoid some types of the convergence problems and oscillations in pressure and velocity fields.
The integrals volume of the temporal and source term of the equation ( 12) are approximated by adopting an Eulerian scheme and using the theorem of the average.The advection terms for the momentum equation arising from the finite volume integration are approximated by the quadratic upstream interpolation scheme for the convective kinetics (QUICK) scheme of Leonard [21].This scheme is second order accurate in space.More details about the numerical procedure are given in [22,23].

TEMPORAL DISCRETIZATION
The stationary solution was obtained by a time-marching algorithm.For a first order time scheme, the system (3)-( 4)-( 7)-( 8) is written as: (15) (16) where the superscript n is the previous time step and n + 1 is the next one.
The convective, diffusive, production and dissipation terms of the different transport equations are treated by an explicit Euler scheme.The advection terms in the k -ε equation are discretized in space using the first order upwind scheme.The diffusion terms are discretized with a second-order cell-centred scheme.The pressure is treated by an implicit scheme, where the decoupling procedure for the pressure is derived from the Marker and Cell (MAC) algorithm proposed by Harlow and Welch [24].The method of solution consists in substituting the equation for the velocity at the new time n + 1 within the discretized equation for mass conservation at the same time level.With this method, the pressure at the next time step n +1 is obtained by the resolution of the discretized Poisson equation.Then, the velocities are calculated from the momentum equations.Since the linear system for the pressure equation is symmetric and positive definite, it can be solved by the Cholesky procedure or by a preconditioned conjugate method.Note that the principle of mass flux continuity is imposed indirectly via the solution of this equation.Convergence was declared when the maximum normalized sum of absolute residual ε ~n for each variable φ over all the computational domain was less than 10 -4 .This maximum error is for each variable determined as: where N is the number of control volumes. )

COMPUTATIONAL DOMAIN AND BOUNDARY CONDITIONS
The geometrical configuration the square duct with the reference axes is shown in figure 4. The x-axis designates the streamwise direction.The normal direction is parallel to the y-axis and the spanwise direction is parallel to the z-axis.The cross-section is divided into four quadrants.Because of symmetry, only quadrant of the square need to be calculated with symmetry conditions at the center (see figure 5).
For the first time-step, at the inlet of the duct, a constant profile was given to U -, k and ε.The secondary velocities were initialised as nil (V -= 0 and W -= 0) all over the domain.The k and ε inlet values were obtained from the DNS data using the r.m.s velocities ( ) and the eddy viscosity, that was about four times the molecular viscosity.At the outlet, a homogeneous Neumann boundary condition was used for all variables.In the next time-step, the calculated outlet values were used for the inlet condition.The same procedure is used for the following time-steps up to the convergence.
The calculations were carried out using 41×41 grid points, regularly spaced in the crosssection and five grid points in the streamwise direction.The outcome with this grid was found to be satisfactory.The grid convergence was checked using 21×21 grid points in the cross-section, the maximum difference observed between the two calculations were less than 1% in the streamwise velocity near the corner (z/h = 0.1).A test case with 20 grid points in the streamwise direction and 21×21 in the cross-section was made in order to check the grid independence in the streamwise direction, and no relevant difference was observed because the maximum difference was about 2% for the spanwise velocity at z/h = 0.16.
The boundary condition values for k and ε, at the first grid point near the wall, was calculated taking into account the fact that this point was in the viscous sub-layer (always ).Also, due to the use of a staggered grid, the value of k and ε are not defined at the wall.In this paper, we consider the following boundary conditions for the equation of k and ε, which have been used by Patel et al. [19]:

or , and or
When this is done, it is necessary to modify the model itself near the wall.It is argued that the effects that need to be modelled are due to the low Reynolds number near the wall and a number of low Reynolds number modifications of the k -ε model have been proposed; see Patel et al. [19] and Wilcox [25] for a review of these modifications.
A condition for symmetry, homogeneous Neumann, was used for all the variables along the wall bisectors of the square duct.

RESULTS
This section shows the results for a turbulent flow in a square duct quadrant using the linear k -ε model with and without damping functions f µ .Direct Numerical Simulation (DNS) of Gavrilakis [10] is used to compare some turbulence characteristics.The developed turbulent flow through a straight square duct has been simulated.The DNS has been carried out for the Reynolds number based on the hydraulic diameter of the duct (2h) and the mean flow velocity U m is 4800.The Reynolds number based on the friction velocity u τ , is Re τ = u τ 2h/v = 320.The velocity ratio U -0 /U m for this configuration was 1.33, U -0 being the mean centreline velocity.The maximum Kolmogorov scale is 1.5v/u τ .In the following, the presentation of the characteristics of the flow are confined to one duct quadrant.
Figure 6 shows the contours of streamwise velocity for one quadrant, obtained with the standard k -ε model, the k -ε -f µ model and the DNS.By the examination of this figure, we remark that the flow is symmetric according the corner bisector.The predictions are in good agreement with DNS data.The model k -ε is unable to predict the effect of vortex induced by the secondary flow on the contours.
The profiles of the mean streamwise velocity (U -) for several section are shown in Figure 7(a)-(c) respectively, for y/h = 0.1, 0.5 and 1.0.In these figures, we show the comparisons between the linear k -ε model, the k -ε -f µ model, DNS of Gavrilakis [10] and the experimental data of Cheesewright et al. [27].At the section y/h = 0.1, Figure 7(a) shows a strong distortion on the mean velocity generated by the secondary flow.On the other hand, the linear k -ε model used was unable to predict this distortion on the mean streamwise velocity   in this region.After y/h = 0.5, the agreement between the results of the different models and experiments starts to improve, and is good at the wall bisector (see Figure 7c, y/h = 1.0).
The comparisons between the two models, DNS of Gavrilakis and experimental data [26] for the spanwise velocity are presented in figure 8(a) and 8(b), respectively, for two sections, z / h 0.1 and z /h = 0.8.
There is a good qualitative agreement between the various datasets of DNS, measurements and k -ε -f µ.At z /h = 0.8, we note that the spanwise velocity is overall well predicted by V ( ) ) is underpredicted at the position of strong velocity, e.g. at y/h Ϸ 0.2 (see Figure 8(a)).It is also noted that the difference between the k -ε -f µ model and the DNS data, in the region of strong velocity, decreases as z /h increases.The linear k -ε does not give any secondary flow (V -= 0) for all the values of z /h.
The figure 9 depicts the secondary velocity vectors in the quadrant obtained from DNS, k -ε -f µ and k -ε models.The DNS results are shown on the left side of the diagonal, and the k -ε -f µ or the k -ε results model on the right side (cf. Figure 9(a) and 9(b) respectively).As can be seen, the k -ε -f µ model is able to capture reasonably the vortices  structure.In other words, it is apparent that the cross-flow plane contains a vortex along the y and z -walls and the secondary vector field seems to be similar to that obtained by Gavrilakis [26].
In Figure 10(c), we can observed clearly that the spanwise contours obtained by the k -ε model are not satisfactory.On the other hand, the spanwise contours obtained by the k -ε -f µ model are compared favourably with those of the DNS data (see Figures 10(a In order to well predict the profiles of the Reynolds stresses, we have introduced the damping functions f µ (see Eqn. (10)) whose the main role is to describe better the pattern of the flow in the near-wall regions.According to our simulations, the best values of constants 'a' and 'b' which allow the reproduction of the Reynolds stresses and are given in Table 2.
by about 70% by the k -ε model.However, when the k -ε -f µ model is used, this difference can be reduced up to 35% (cf. Figure 11d).Figure 12), the results from the k -ε -f µ model are in overall agreement with the DNS at y/h = 0.1, 0.3, 0.7, whereas at y/h = 1.0, some substantial differences remain in particular for (cf. Figure 12d).As for the k -ε model, it overpredicts this quantity for starting from the section y/h = 0.1.Considering the normal stress (see Figure 13), the k -ε -f µ model agrees well with the DNS, especially in the region z + ≥ 40 at y/h = 0.7 and y/h = 1.0.However, at y/h = 0.3 (see Figure 13b), the two models give completely different   -+ and the primary shear stress -uw -+ profiles along the wall bisector of the square duct at centreline (y/h = 1.0).These values are compared with the DNS of Gavrilakis [10], the LES of Xu and Pollard [4] and the experimental data from Nishino and Kasagi [7].In global sense, the predicted results by the k -ε -ƒµ model disclose that this model can gives an overall agreement.It is seen that the discrepancy is rather weak the secondary normal stresses vv -+ , ww -+ and for the shear stress -uw -+ indicating that our methodology may be considered as sufficient.
The distribution of the turbulent kinetic energy from the DNS and the k -ε -f µ model are compared in Figure 16.The differences among our results and the DNS are not substantial.The k -ε -f µ model yields better results than the results of the k -ε model (see Fig. 16c).The figure 16b shows the results of the k -ε -f µ model indicating a qualitatively good agreement of k with the profile obtained by DNS.

CONCLUSION
The purpose of this paper is to describe a study of the turbulent flow in a straight square duct which involves a secondary flow.The spatial discretization of the RANS equations is performed by a finite volume method with a staggered variable arrangement.The linear k -ε model is employed for the prediction of the considered flow.This model has been modified by using different damping functions.The present corrections were validated -uw + by comparing the calculated velocity profiles and the turbulent quantities with the available direct numerical simulation data given in Ref. [10].These comparisons yielded the following conclusions: • Mean velocity profiles are, in general, in good agreement with the DNS data in particular when one approaches the center of the duct (y/h = 1.0 or z/h = 1.0).

•
Although some discrepancies exist for the turbulent quantities, the k -ε -f µ model leads to similar results.

•
The evaluation of the considered model with DNS data reveals that some anisotropy close to the wall can be captured.

•
The k -ε -f µ model tested in this work yields considerably better predictions than those obtained by the standard k -ε model.

Figure 2
Figure 2 Arrangement for a staggered grid.

Figure 3
Figure 3 Control volume for a staggered and notation used.

Figure 5
Figure 5 Quadrant of the square.

326
Numerical simulation of turbulent flow through a straight square duct using a near wall linear k -ε model

330
Numerical simulation of turbulent flow through a straight square duct using a near wall linear k -ε model

Figure 11 Figure 12
Figure 11 The normal Reynolds stress profiles along the wall bisector at y/ h = 0.1, 0.3, 0.7 and 1.0.Comparison of the predictions of the k -ε and k -ε -f µ models with DNS.uu + Figures 12, 13 and 14 respectively.It seems clear that the k -ε model cannot correctly reproduce all these stresses whereas the k -ε -f µ model gives very good results in particular for y/h = 0.7 and y/h = 1.0.When considering the secondary normal stress (cf.

Figure 13
Figure 13 The normal Reynolds stress along the wall bisector at y/h = 0.1, 0.3, 0.7 and 1.0.Comparison of the predictions of the k -ε and k -ε -f µ models with DNS.ww +

Figure 14
Figure 14 shows the profile of the primary shear stress obtained by the k -ε and k -ε -f µ models.At y/h = 0.7 and y/h = 1.0, the differences among results of the k -ε -f µ models are not substantial.The results of y/h = 0.7 from the k -ε model, indicate that this model underpredicts slightly this component for , whereas at y/h = 1.0, they indicate that the model underpredicts even more this quantity for .It is also noted that the k -ε -f µ model predicts the right position of the peak of the primary shear stress at z + ≈ 25 (cf.Figure 14).Figure 15 displays the normal Reynolds stresses uu

Figure 15
displays the normal Reynolds stresses uu

Figure 14
Figure 14 The dimensionless Reynolds stress along the wall bisector at y/h = 0.7 and 1.0.Comparison with DNS.− + uw

334Figure 15 Figure 16
Figure 15 The dimensionless Reynolds stresses and along the wall bisector at y/h = 1.0.Comparison of the predictions of the k -ε and k -ε -f µ models with DNS, LES and experiment.

Table 2
Values of constants 'a' and 'b' used in the f µ function