Investigations of Advanced Reynolds Stress Closure Models by
Using Database from DNS for Modelling Turbulent Flow

L. Yu¹, A.M. Righetto² and Jun Yu³

1. SHS - Consultoria e Projetos de Engenharia Ltda., São Carlos, SP, Brazil

2. Dept. of Hydraulics and Sanitary Engineering, EESC, University of São Paulo,

São Carlos, SP, and LARHISA-CT/Dept. Civil Engineering, University of Rio Grande do Norte, Natal, RN, Brazil

3. Student, Institute of Mathematics, Statistics & Scientific Computation, University of Campinas, SP, Brazil 

Abstract: This paper presents the numerical investigations of four turbulence Reynolds Stress Models (RSM) appeared in the recent years, which can be applicable within the near-wall sub-layer and without the help of traditional wall treatment technique (law of the wall or wall-function approximation). These models are the Launder and Shima’s model (L&S-Model, 1989), Launder and Tselepidakis’ model (L&T-Model, 1991), Launder and S-P Li’s model (L&Li-Model, 1993) and Shima’s model (S-Model, 1993), respectively. The numerical results for simulating the classical flat-plate boundary layer flow, by using these four Reynolds stress closure models (second-moment closure), were compared with the database of Direct Numerical Simulation (DNS). The authors found out that the Launder and Shima’s model is the best one among these four tested models.

1    INTRODUCTION

Numerical computation, as one of powerful tools, for modelling turbulent flows in pipes, open canals and boundary layers which are a foundation of Hydraulics has been widely adopted for many years. The conventional turbulence two-equations closure models, such as the famous k-e model, commonly used in current engineering computations, are simply based on the Boussinesque’s eddy viscosity assumption. However, the RSM based on the transport equations for the Reynolds-Stress Tensor (RST)  are high-class, refined multi-equation closure models with higher precision and stronger ability in simulation and prediction. The RSM, on the other hand, are unfortunately more complex, especially when they have to satisfy the requirements of invariance and realizability. Though RSM have great potential, they are therefore not much in use for practical applications so far, as their great complicity and difficulty for solving numerically the problems in practice. One of the main reasons which are still restricted the application is that the simplified wall treatment technique (law of the wall or wall-function approximation) has to be used in early versions of RSM. Several proposals and corresponding models to extend second-moment closure models up to the wall have appeared in the recent few years. It seems that these proposals could bring us to overcome this disadvantage. However, the ability, efficiency and realizability of these newly appeared models are not quite clear and numerical investigations and comparisons are necessary in deed.

With the help of sufficiently large and powerful computers and advanced numerical algorithms, the DNS technique through the direct solution of Navier-Stokes equations by using accurate approximations has been developed recently. Though there still has a long distance for directly adopting the DNS in turbulence computations in engineering, the DNS of simple shear flows, however, has proven to be an effective and important tool for studying turbulence structures and their near-wall effects. It is possible to examine the individual behaviours of some important terms of RSM by use of the results of DNS.

In this paper, the distributions of turbulence parameters resulted from four Reynolds stress closure models and the database resulted from DNS are compared carefully with the aim of investigating numerically the ability of the advanced turbulence Reynolds Stress Models.

2   FUNDAMENTAL EQUATIONS

For the incompressible and steady flat-plane boundary layer flow under a Cartesian co-ordinate system (x₁x₂x₃), the equations of mean movement can be expressed as follows:

                           (1)

                  (2)

where  and  (i=1,2) stand for the mean velocity and corresponding fluctuating quantity; n is the kinematic viscosity; x₁ represents the longitudinal direction x₂, the normal direction and x₃, the lateral direction, respectively.

At the level of Reynolds stress tensor closure, the second-moment  can be solved by following transport equation:

                    (3)

where ; P_ij stands for the stress generation rate by mean shear; f_ij represents the averaged product of the fluctuating pressure and stain fields (pressure-strain rate tensor); e_ij denotes the dissipation rate tensor;  and d_ij are the viscous diffusion rate and turbulence diffusion, respectively. In order to determine turbulent kinetic energy k= appeared in the terms f_ij and d_ij, three normal stress transport equations, except for the equations (1)-(3), also need to be solved simultaneously:

                      (4)

                     (5)

                           (6)

For e_ij, the assumption of local isotropy  is frequently adopted, where e, the energy dissipation rate of turbulence, needs to be obtained by solving an adequate transport equation.

In order to close the fundamental governing equations, four extra Reynolds stress transport equations (3)-(6) and one dissipation rate equation (e-equation) should be solved together with the equations (1) and (2).

3    SECOND-MOMENT CLOSURE MODELS FOR NEAR-WALL SUB-LAYER

Since P_ij and  require no approximation and d_ij usually adopts the generalised gradient diffusion hypothesis, one can easily derive the component expressions of P_ij:

, , ,  and the component expressions of total stress diffusion as follows:

          ,

          ,

The pressure-strain rate tensor  can be divided into four distinct terms for modelling near-wall sub-layer, namely the slow, rapid, slow wall-reflection and rapid wall reflection terms, respectively:

                         (7)

where the last two terms in the right member of the equation (7),  and , usually have variable and extremely complicated expresses in different models.

In L&S-Model (1989), for example, the  is defined by the following expression:

                            (7.1)

where = represents the dimensionless anisotropic part of the Reynolds stress, which is a proposal first made by Rotta as early as in 1951. The recommended form of the coefficient , where the turbulence Reynolds number  and stress flatness factor , A₂ and A₃ are two independent Reynolds stress invariants, i.e., the second invariation of a_ij:  and the third invariation of a_ij: . ,  and  can be expressed as follows:

                         (7.2)

       (7.3)

         (7.4)

where the shear production rate of turbulence energy P=P_kk/2, the near-wall damping function , x₂ stands for the distance normal to the wall and n_k represents the unit vector perpendicular to the surface, the coefficients ,  and , respectively. The turbulence energy dissipation rate e in the L&S-Model (1989) is obtained by solving following equation:

       (7.5)

where , , , empirical coefficients c_e, and  are 0.18, 1.45 and 1.9, respectively.

4    NUMERICAL METHOD

A numerical calculation procedure for momentum, heat and mass transfer in three dimensional parabolic flows (Patankar and Spalding, 1971) has been adopted to solve fundamental governing equations and second-moment transport equations. For this kind of flows, the co-ordinate in the main flow direction (x₁) becomes a “one-way” co-ordinate. It is a convenient behaviour that enables us to employ a marching integration from an upstream station to a downstream one. The fundamental governing equations of fluid movement and transport for boundary-layer flows in a Cartesian co-ordinate system can be generally written by following form:

      (8)

where  stands for the transport property such as viscosity in momentum equations. The source-sink term S_f on the right member should include all terms, which can not be expressed in the forms of convection and diffusion. By integrating equation (8) over the control volume shown in the Figure 1 by dotted lines, one can transform it into a discretized algebraic expression:

                  (9)

where A_N, A_S, A_E and A_W stand for coefficients, B^f represents source term. The finite-difference equation like (9) can be solved by successive use of the TDMA (Tri-Diagonal Matrix Algorithm) in the x₂ and x₃ directions. The foregoing procedure for calculation was based on the assumption that the longitudinal pressure gradient  is taken to be the same as the longitudinal pressure gradient prevailing in the irrotational free stream adjacent to the boundary layer. Then the solution of momentum equation for U₁ is straightforward. In the present computation, the second-moment , ,  and e, the same as the longitudinal velocity U_p, are arranged at the centre nodal point P on the x₂x₃ face. The second-moment , the same as the v_p, is arranged at the volume face between nodal points P and S.

Fig. 1    Control volume adopted to obtain finite-difference equation

5  COMPARISONS AND ANALYSES

The authors computed the incompressible, steady flat-plate boundary layer of turbulence, closed by L&S-Model (1989), L&T-Model (1991), L&Li-Model (1993) and S-Model (1993), respectively. These computations employed the genuine no-slip boundary condition at wall, instead of the frequently used law of the wall. The adopted longitudinal forward step length Dx₁ in our computations is equal to 10^-4 m and the molecular viscosity n is 1.48´10^-5 m²/s, respectively. The number of the nodal points within the computational altitude of the boundary layer equals to 175, which was determined in terms of a grid independent analysis (Yu 1999).

The DNS, carried out previously, is directly to solving the Navier-Stokes equations in the case of the flat-plate boundary layer, with the Reynolds number being 8000. The database resulted from DNS can be configured by several profiles of turbulent quantities, which are distributed across all of the boundary layer thickness, including the viscous sub-layer (0<y⁺<5), buffer layer (buffer region, 5<y⁺<40), outer zone and super-layer, respectively. The non-dimensional distance to the walls y⁺ equals to yU_w/n, where y stands for the real distance to the wall and U_w represents the fraction velocity. Figure 2 - Figure 7 present the comparisons of the data resulted from DNS, represented by the circle-dot, with those resulted from second-moment closure modelling by use of L&S-Model (1989), represented by the real line, while the number of forward steps equals to 10⁴. These figures illustrate the profiles of the dimensionless longitudinal velocity (mean velocity)  (U⁺), turbulent shear stress (Reynolds shear stress) , turbulence intensities in longitudinal direction (x or x₁-direction; longitudinal normal stress) , vertical direction (y or x₂-direction; vertical normal stress)  and lateral direction (z or x₃-direction; lateral normal stress)  as well as the dissipation rate of turbulent kinetic energy , respectively. The abscissa (in the normal direction of the flat-plate boundary layer) is adopted to illustrate y⁺.

            Fig. 2    Mean velocity               

Fig. 3    Turbulent shear stress

The results computed by L&S-Model (1989) coincide quite well with those resulted from DNS, especially in the distributions, such as the mean-velocity (see Figure 2) and three normal stresses (see Figures 4-6), but except for the dissipation rate of turbulence kinetic energy e⁺ in the range of y⁺<80 (see Figure 7). The profile of e⁺, calculated by L&S-Model (1989), seems to be unreasonable, while e⁺, which should decreases along with the increase of y⁺, is over-predicted within the range of 8<y⁺<80 with the appearance of a peak value in the buffer layer at y⁺=10. However, the corresponding comparison between the e⁺, calculated by L&S-Model (1989), and the e⁺, resulted from DNS, did not be presented in the original paper in which the L&S-Model was offered firstly (Launder and Shima 1989). The difference of the distributions between the calculated e⁺ by L&S-Model (1989), and one, resulted from DNS, may be caused by the assumption of local isotropy ( ). Also, the assumption of local isotropy may make some difference of turbulence shear stresses distributions in the range between 1<y⁺<10 (see Figure 3), while the production of turbulence calculated by L&S-Model (1989) is larger than the data resulted from DNS, though these two profiles in other ranges (y⁺<1 and y⁺>10) are quite well agreed each other.

     Fig. 4    Longitudinal turbulent intensity   

 Fig. 5    Vertical turbulence intensity

        Fig. 6    Lateral turbulent intensity

Fig. 7    Dissipation rate of turbulence kinetic energy

The numerical comparisons with the database resulted from DNS show that the results calculated by the closures of other three turbulence Reynolds Stress Models exist distinct differences (Yu 1999). For example, Figure 8 - Figure 10 present the distributions of the dimensionless Reynolds shear stresses, resulted from the closures of L&T-Model (1991), L&Li-Model (1993) and S-Model (1993) at thousandth forward step, respectively.

  Fig. 8    Turbulent shear stress calculate by  L&T-Model  

Fig. 9    Turbulent shear stress calculated by L&Li-Model

Fig. 10    Turbulent shear stress calculated by S-Model

6    CONCLUSIONS

According to the published literatures, four second-moment closure models, including L&S-Model (1989), L&T-Model (1991), L-Li-Model (1993) and S-Model (1993), have been to simulate the incompressible, steady turbulence flat-plate boundary layer flow. The foregoing procedure of Patankar and Spalding (1972) is adopted by the exercise of the genuine no-slip boundary condition at wall, instead of the frequently used law of the wall (wall function) at the region close to the rigid surface. The grid independent analysis is performed at first in order to determine the suitable and economic grid arrangement. The database resulted from DNS have been employed to compare and analyse the results calculated by the tested Reynolds stress closure models. The numerical investigations found out that L&S-Model (1989), though being provided rather earlier than other three ones, possesses good ability for modelling longitudinal velocity, various Reynolds stress components in all of the boundary layer thickness and relatively fair ability for modelling the distribution of the dissipation rate of turbulence kinetic energy, especially in the outer zone. Through these investigations, it is also to be understood that the distribution of the dissipation rate of turbulence kinetic energy in sub-layer and buffer layer of boundary layer, predicted by L&S-Model (1989), is unreasonable; the distribution of turbulence shear stress in the range between 1<y⁺<10, calculated by the same model, is unsatisfied. These limitations, in authors’ opinion, are mainly caused by the adoption of the simple assumption of local isotropy for turbulence dissipation rate e_ij in the very model.

The results calculated by L&T-Model (1991) and S-Model (1993) generally do not agree to the database resulted from DNS. The “over-smooth” velocity profiles are predicted by these two models and the “smoothness”, or say the positive velocity shift, increases along with the increase of operated forward steps. The turbulence parameters in sub-layer and buffer layer are usually predicted insufficiently, however, the corresponding turbulence quantities in outer zone are over-predicted. It is also understood that L&T-Model (1991) is more poor than S-Model (1993), as the simulation closed by L&T-Model (1991) only can run up to 5´10³ forward steps.

On the contrary, the results calculated by L&Li-Model (1993) demonstrate the “over-rough” behaviour of mean velocity profile. In sub-layer and/or buffer layer, the turbulence quantities are usually over-predicted and disagreeable to the data resulted from DNS, though their distributions basically remain unchanged with the increase of forward steps. In outer zone, a inverse changeable tendency appears, which makes the profiles of turbulence parameters, except for the dissipation rate of turbulence kinetic energy, shifted towards the direction of the increase of y⁺. L&Li-Model (1993), unfortunately the same as L&T-Model (1991), also can not run beyond 5´10³ forward steps. The fact that L&T-Model (1991), L&Li-Model (1993) and S-Model (1993) were developed after L&S-Model (1989) shows the great difficulty and complicity for improving the high-class and refined second-moment closure to model near-wall behaviours of turbulent flows. 

Acknowledgements

The support of FAPESP (Foundation for Supporting Research at São Paulo State) through Process No. 1997/13955-4 is greatly acknowledged.

References

Launder, B.E. and Shima, N.: “Second-moment closure for the near-wall sublayer: development and application”, AIAA Journal, Vol. 27 (10), (1989), pp. 1319-1325.

Launder, B.E. and Tselepidakis, D.P.: “Progress and paradoxes in modelling near-wall turbulence”, Proceedings of eighth symposium on turbulent shear flows, 29-1, (1991), Technical University of Munich, September.

Launder, B.E. and Li, S.P.: “On the prediction of flow over riblets via 2nd-moment closure”, In: Near-wall turbulent flows, (1993), pp. 739-748, (ed. R.M.C. So, C.G. Speziale and B.E. Launder), Elsevier Science Publishers B.V.

Patankar, S.V. and Spalding, D.B.: “A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows”, International Journal of Heat and Mass Transfer, Vol. 15, (1972), pp. 1787-806.

Rotta, J.: “Statistische theorie nichthomogener turbubulenz”, Zeitsch. für Physik, Vol. 129, pp. 547-572 and Vol. 131, (1951), pp. 51-77.

Shima, N.: “Reynolds stress budgets of a second-moment closure near fixed and moving walls”, Proceedings of ninth symposium on “Turbulent shear flows”, 8-3, (1993), Kyoto, Japan, August.

Yu, L.: “Numerical Analyses and Comparisons of Different Second-Moment Closures with DNS for Modelling Near-Wall Sub-layer of Turbulent Flow”, Research Report for FAPESP (1999), São Carlos, Brazil.