STUDY ON GRID DEPENDENCE OF NUMERICAL SIMULATION
 FOR CIRCULAR PIPE FLOW USING LOW- REYNOLDS
 NUMBER k-ε MODEL

Yo Kumagai¹ and Hitoshi Tanaka²

1 Ph.D. Student, Department of Civil Engineering, Tohoku University, Aoba, Sendai 980-8579, Japan.  

Phone: (+81)-22-225-2897, fax: (+81)-22-222-0453

E-mail: w910735@tohoku-epco.co.jp

2 Professor, Department of Civil Engineering, Tohoku University, Aoba, Sendai

980-8579, Japan.  Phone: (+81)-22-217-7451, fax: (+81)-22-217-7451

E-mail: tanaka@tsunami2.civil.tohoku.ac.jp

Abstract: In a numerical analysis of turbulent flow with a boundary fitted coordinate system, dependence of computation accuracy on grid system should be fully considered. Up to now, many computations have been carried out for various kinds of hydraulic structures using curvilinear coordinate, though, there have been no subjective and quantitative criterion for the choice of grid system which yields sufficient accuracy. In this study, an analysis for a flow in a circular pipe is made and compared with direct numerical simulation (DNS) data, in order to evaluate the dependence of computation accuracy on grid topologies. A grid distortion parameter  is newly introduced in the present study and the relationship between  and the accuracy of corresponding computation is discussed. It is found that a grid system which satisfies >0.7 in the region of y⁺<20 produces a considerable error when a low Reynolds k- turbulence model is applied.

Keywords: k- turbulence model, grid topology, circular pipe flow, DNS, computational hydraulics

1  INTRODUCTION

There has been recently an increase in application of computational hydraulics to the practical hydraulic structure design because computer price is getting lower from year to year. Generally, hydraulic structure varies with the width, height and the shape. When calculation using structured grid is carried out, an error originating from grid is generated, so that it is difficult to anticipate an adequate accuracy of calculated results. Boundary fitted coordinate system and unstructured grid system are used for hydraulic structure design, however, there is no quantitative criterion for the choice of grid system yielding sufficient accuracy. This study deals with a turbulence analysis of flow in a circular pipe, being a fundamental shape of flow. A relationship between grid distortion and the error originating from it is examined, by comparing numerical analysis results with DNS data.

2   NUMERICAL COMPUTATION METHOD

Basic equations

The basic equations are the continuity equations and the Reynolds-averaged Navier-Stokes equations:

,                      (1)

                         (2)

where ρ is the fluid density,  (i=1,2,3) is the time-averaged-velocities, p is the pressure, ν is the kinetic viscosity, ν_t is the eddy viscosity,  is the turbulent stresses calculated with the k-εmodel, and k and εare turbulent kinetic energy and its dissipation rate respectively. The turbulence model adopted in the present study is the low Reynolds number k-ε model proposed by Abe, Nagano and Kondo^{1), 2)}, as shown in Eqs.(3), (4) and (5).

                     (3)

                        (4)

                                   (5)

where σ_k and σ_ε are the model constants, the first term of the right hand side of Eq.(3) is turbulent diffusion TD, the second term is production PR, and the third term is dissipation DS. The five constants and the model functions f_μ and f₂ in this model are as follows:

C_t =0.09       σ_k=1.4       σ_ε=1.4       C₁=1.5       C₂=1.9        (6)

                     (7)

                      (8)

Dimensionless length y*, y⁺, turbulent Reynolds number R_t , and Kolmogorov velocity scale  U_ε are defined as follows:

                     (9)  

                       (10)

                      (11)  

 U_ε = (νε)^1/4ν                     (12)

where U_τ =(τ_w /ρ)^1/2 is the wall friction velocity andτ_w is shear stress. The wall boundary condition is a no slip one at y⁺<1, and a log-law at y⁺>11.6. The wall boundary condition for the intermediate region between y⁺<1 and y⁺>11.6 is given so that the no slip condition and log-law regions are bridged. The above-mentioned condition is summarized as follows, Eqs.(13) - (16), where  (=11.6) is the intersection point of linear-law and log-law.

                            (13)

                                  (14)

                     (15)

                                   (16)

In the numerical computation, the finite volume method and the SIMPLE algorithm are adopted, and the convective acceleration terms are discretized by mean of MUSCL.

Dns data for model verification

Eggels et al.³⁾ carried out DNS computation for turbulent pipe flow in a cylindrical geometry, in which the length of computational domain L is L=5D where D is the pipe diameter. The computation has been performed with 96×128×256 elements equally spaced in the r, θ, z direction, respectively. The Reynolds numbers is 5,300 based on the pipe diameter and the averaged velocity U_b. The value y⁺ is 0 on the wall, and 180 at the center of cylindrical flow. From their DNS computation, the mean flow properties as well as turbulence statistics are obtained such as axial time-averaged-velocity U, turbulence energy k, and the terms in the energy budget of turbulence kinetic energy k, i.e., turbulent diffusion TD, production PR and dissipation DS.

Grid system

In the present study, calculation using the low Reynolds number k-εmodel is carried out in the computational domain with the length of L=D, as indicated in Fig.1. Calculation cases and their grid system are shown in Table 1 and Fig.2, in which Case 1 through 4 have 41 or 48 elements in y and z directions. In the present study, calculation using the low Reynolds number k-g model is carried out in the computational domain with the length of L=D, as indicated in Fig.1. Calculation cases and their grid system are shown in Table 1 and Fig.2, in which Case 1 through 4 have 41 or 48 elements in y and z directions. The characteristics of each grid system are summarized as follows: (1) In Case 1, the distortion g is relatively smaller in the computational domain as compared with other three cases, (2) In Case2, the maximum of is higher than Case 1 and is located in the fully turbulent region, (3) In Case 3, the peak of g is located between log-law and linear-low region, and (4) In Case 4, the peak of g is located within the viscous sublayer. More than three elements are placed within 0 <y⁺<11.6 in each case. Computation is carried out only in the upper half region of the cylindrical pipe, using the periodical boundary condition. Here, we define the magnitude of grid distortion parameter g as,

γ=cos|(e_x , e_y)|+cos|(e_y , e_z)|+cos|(e_z , e_x)|                     (17)

Fig. 1  Coordinate system

                                   Table 1   Calculation Conditions

Fig. 2  Grid system and cross-sectional distribution of g

where e_x, e_y, e_z is the unit vectors in x-, y- and z-directions. In the present study, there is no distortion in x-direction and γ can be then simplified as Eq.(18).

γ= cos|(e_y , e_z)|                              (18)

3  RESULTS AND DISCUSSIONS

Comparison of computation with dns

U, k and γ in DNS data is completely axisymmetric, i.e., independent of the angle θ, though, the result of this computation is not axisymmetric due to distortion of the grid system shown in Fig.2. Fig.3 shows the vertical profile U from DNS data along with the present computation for Case 3 at θ=0 deg. and 45 deg. These show good agreement with each other at the angle of 0 deg., but small difference can be observed between computational results at 0 deg. and 45 deg. in the range of 0＜y⁺＜10. Fig.4 shows a comparison of k between DNS and the numerical simulation for Case 3. The difference of the numerical results at two angles is more in the range of 0＜y⁺＜100, being bigger difference as compared with the mean velocity U shown Fig.3. It is seen that the maximum of k at 0 deg. is 20% higher than that at 45 deg. Figs.5 and 6 illustrate vertical profile of TD, PR and DS defined in Eq.(3). As seen in Fig.5, the difference between 45 deg. and 0 deg. attains to 4 times. Fig.6 shows similar behavior of difference dependent on the angle in the range of 0＜y⁺＜20. The grid system for Case 3 shows remarkable distortion the near-wall region as seen in Fig.2.  Thus, based on the results shown in Figs.3 through 6, it is expected that the error in the computation is dependent on γ as well as on y⁺.

Fig. 3   Distribution of U/U_b

Fig. 4　Distribution of  k

Fig. 5  TD terms in k equation

Fig. 6  PR and DS terms in k equation

Evaluation of computation accuracy

In order to make a quantitative investigation of the error in the computation associated with the grid distortion γ and y⁺ , the following quantities are defined.

,   , ,

,                          (19)

where the subscripts (k–ε) and (DNS) denote the model employed and max means maximum value in the computation domain. It is noted that DS is normalized by the minimum value,  because DS takes negative value at the range of 0<y⁺ <180.

The contours of error of U, k, TD, PR and DS calculated by Eq.(19) is plotted in Fig.7, where the abscissa is the dimensionless distance y⁺ and the ordinate is the grid distortion γ. The gray area in the figures denotes that there in no grid point corresponding to this condition. In Case 2, for example, the grids with high distortion are located at y⁺ ≒90, whereas they are located in the vicinity of the wall in Case 3 and Case4. Therefore, the contours of error can be drawn in the white area. In some figures, the area close to the origin is exaggerated so that detailed contours can be understood easily.

In Case 1, the maximum of γ is about 0.5 and it appears sufficiently outside of the viscous sublayer. In this case, the error of U, k, TD, PR and DS calculated by Eq.(19) is negligibly small as seen in Fig.7. In Case 2, the maximum of γ is located far outside of the viscous sublayer, though, its maximum value attains to about 1.0, being much bigger than Case 1. Although the error of TD in as high as 60%, it does not affect to other quantities such as k and U.

In Case 3 and Case 4, in contrast, the maximum of γ appears much closer to the wall as compared with Case 1 and Case 2, and the error of U, k, TD, PR and DS calculated by Eq.(19) is remarkably higher. It should be here noted that TD shows an appreciable error beyond 200%. It is thus concluded that a grid system for the low Reynolds number k-ε turbulence model should not satisfy γ＞0.7 in the range of  y⁺＜20 to avoid a considerable error in the computation.

Fig. 7　Error of U, k, TD, PR and DS (%)

4  CONCLUSIONS

(1) Error of TD appears when γ＞0.7 is satisfied, but that this does not cause the error of  U and k as far as the maximum of γ is located far outside the viscous sublayer. As the peak of γ approaches the near-wall region, the error with higher magnitude will be  generated not only in TD but also in PR, and DS, resulting in remarkable error of k and U .

(2) The grid system for an application of the low Reynolds k-ε turbulence model should not take the range of y⁺＜20 and γ＞0.7 to achieve high accuracy of the computation.  

References

[1]  Abe, K., Nagano, Y. and Kondoh, T.: An improved k-ε model for prediction of turbulent flows with separation and reattachment, JSME No.58, Vol.554, pp.3003-3010, 1992.

[2]  Abe, K., Nagano, Y. and Kondoh, T.: Numerical prediction of separating and reattaching flows with a modified low-Reynolds-number k-ε model, J. Wind Engineering, No. 52, pp.213-218, 1992.

[3]  Eggels, J.G.M., Unger, F., Weiss, M.H., Westerweel, J., Adrian, R.J., Friedrich, R. and Nieuwstadt, F.T.M.: Fully developed turbulent pipe flow: A comparison between direct numerical simulation and experiment, J. Fluid Mech., Vol.268, pp.175-209, 1994.