Pergamon PII:
Computers & Fluids Vol. 26, No. 4, pp. 417-423, 1997 0 1997 El&w Science Ltd. All rights reserved Printed in Great Britain SOO45-7930(96)0004&5 00457930/97 $17.00 + 0.00
TECHNICAL
NOTE
NUMERICAL PREDICTION OF SEPARATION AND REATTACHMENT OF TURBULENT FLOW IN AXISYMMETRIC DIFFUSER I). XV*, M. A. LESCHZINER2,
B. C. KHOO’ and C. SHU’
‘Centre for Computational Mechanics and Mechanical & Production Engineering Department, National University of Singapore, Singapore, 119260, Singapore :‘Mechanical Engineering Department, UMIST, Manchester M60 lQD, U.K.
(Received
8 April 1996; in revised form
12 August
1996)
Abstract-The turbulent flow inside an axisymmetric diffuser with a curved-surface centre body is numerically simulated using high Reynolds number k-6 + wall function model, high Reynolds number k-c+ one-equation model, and low Reynolds number k-c turbulence model. For the separation and reattachment positions, the comparisons made between numerical results and experimental measurements show that the high Reynolds k-6 + one-equation turbulence model is superior to other models in the present calculation. 0 1997 Elsevier Science Ltd.
1. INTRODUCTION
Separation and reattachment are very common features in most flows of engineering interest, and phenomena associated with separation often dominate the flow fields, therefore drastically affecting the efficiency of engineering devices. A classical example of this, in external flows, is the dramatic increase in form drag associated with separation. However, many engineering devices operate at their high efficienlcy close to separation. It is this fact that has led to intensive research in this area. The principal ‘objective of the present research effort is to develop and validate the various turbulence mode.ls that will successfully predict the effects of flow curvature on the associated turbulence structure of a boundary layer. One important outcome of the present work is the ability to predict the separation and reattachment of turbulent flows inside an axisymmetric diffuser with a curved-surface centre body. Such findings can be extended to a wide range of flows such as those over wings of aer’oplanes and airfoils of cascades where flow separation and attachment occur over surfaces with curvature. 2. MATHEMATICAL
FORMULATIONS
2.1. Governing equation In tensor notation, the incompressible, two-dimensional, time-averaged Navier-Stokes of motion and continuity of fluid are given below
qfg
= -
equations
j&v+ $A)+-&[@+v,,(~ +tp)]
2.2. Turbulence models Equations (1) (and (2) by themselves are not closed because of the turbulent kinetic energy term (k) and turbulent viscosity term (v,). To enable the closure of the time-averaged Navier-Stokes *Corresponding
author. 417
418
Technical
Note
equations, one has to find effective and suitable equations relating k and v,. The various models presented are listed below. 2.2.1. One-equation model. In the standard one-equation model, the kinetic energy k can be obtained from the modelled transport equation
u,g=q(fg +z)+ here U, is obtained
&[(v+$A$]-&
from following: v, = C,k’121,,
(4)
F = CDk3i2/1D
(5)
f, = C,y[l.O - exp( - A&,.)]
(6)
ID = C,y[ 1.O - exp( - A&)]
(7)
R> = k”‘y/v
(8)
where C,, = 0.09, C, = KC’;‘.~, CD = 1.0, A, = 0.014, A, = 0.20, K = 0.42. y is the distance to the wall. 2.2.2. Two-equation k-c model. The two-equation k-6 model can be further divided into high Reynolds number k-t model and low Reynolds number k-c model. 2.2.2.1.
High Reynolds
Number
k-c Model
In the region away from the wall (i.e. y+ > 60) the standard high Reynolds number version of the k-c model, originally devised by Launder and Spalding [l], was adopted. The turbulent viscosity is now taken as
C, takes the same value as in the above one-equation model. The turbulent kinetic energy is still governed by equation (3) as in the one-equation model, except that the rate of dissipation in that equation is now obtained from a modelled transport equation
(10) here C, = 0.09, C,, = 1.44, C,, = 1.92, crh = 1.0, cr, = 1.3. The reason for this present form of the t equation used in the study is that it has been tested extensively in computations for shear-free flows and wall flows, and has proven to be successful [ 1, 21. 2.2.2.2.
Low Reynolds
Number
k-c Model
In the near wall region, a more general low Reynolds number turbulent model is required, in which the turbulent viscosity can be determined purely through the transport equations. In the present study, Launder-Sharma’s low Reynolds number k-c model is adopted, which is given as v, = C,'k=/E .
(11)
Here i? = E -
2V(&,6C/dX,)2
.
It assumes the value zero at the wall while it differs only slightly from 6 in the fully turbulent The following transport equations are provided for both k and 8: uig
4.33
+ $3,
&[(v+
-$+-2v(g)’
(12)
region.
(13)
Technical Note
uig
=C$J,(~
+ $)2
+ &[(v+
419
$$]-C/$
+2vv,($)2
(14)
where C,’ = C, exp[ - 3.4/(1 + R,/50)2]
(15)
CE2’= C,,[l - 0.3 exp( - R:)]
(16)
R, = k2jEv .
(17)
In the fully turbulent region away from the wall, the second-order derivative term, a’U,/ax,ax,, becomes negligibly small and both the coefficients C,’ and Cc; take on the constants recommended for the high Reynolds number flows. Consequently, the model automatically reverts to the high Reynolds number version of k-c model in the high Reynolds number flow region. 3. NUMERICAL
SOLUTION
The governing equations may be summarized as follows:
(18) where 4 is a general dependent variable which may stand for Vi, k and c; Fb is the diffusive coefficient, and Sqbrepresents all the source terms which cannot be expressed as either convection or diffusion. With control volume methods, equation (18) is integrated over a control volume surrounding P, as shown in Fig. 1, and yields
Applying Gauss’ theorem, the volume integral on the left-hand side may be expressed in terms of surface integral; equation (19) can thus be rewritten as F,-F,+F,-F,=s,*Vol
(20)
where the subscr:ipts e, w, n and s indicate the east, west, north and south faces of the control volume, respectively. With the approximations for convection, diffusion and the source terms, equation (20) may be rearranged, such as
(A, - SP)~,=
CA,&
+ S, .
(21)
Equation (21) is then solved by a ‘line-by-line’ iterative method. The present study adopts the SIMPLE [3] algorithm to calculate the pressure field. Three types of boundary conditions are imposed in the present study. At the solid wall, the velocity is due to the no-slip condition, zero, and the turbulent intensity is either set to be zero
Table 1. Grid Indetwtdent
-ax+s
1
ss
Fig. 1. Control volume.
test
Mesh
K/H
XRIH
80 x 40 100 x 50 150 x 60 17.5 x 70 200 x 70 250 x 70 300 x 70
1.459 I ,411 1.363 1.332 1.332 1.293 1.301
6.619 6.215 6.616 6.658 6.598 6.583 6.512
420
Technical
0
5
Note
10 Fig. 2. Geometry
15
20
and mesh (200 x 70).
or obtained from a wall function according to the turbulence models. The inlet velocity is obtained from experiment; the kinetic turbulence energy is set to kinlet= CkV&,
(22)
with C, = 0.01-0.001. The turbulence energy dissipation rate, L, is determined as follows: a = k’:‘/C;
3141
where I= min(0.41y,0.096). Here, y is distance to the wall and 6 is boundary-layer the outlet boundary, all governing parameters are determined by
a4jac = 0.
(23) thickness. At (24)
The iterative solution may be generally considered to have converged when the sum of the normalized absolute residuals across all nodes is less than a prescribed small value 6
In our computations,
we set 6 to be 10 - 4. 4.
RESULTS
AND DISCUSSION
The numerical solutions are carried out for an axisymmetric diffuser, where experimental data are available in Refs [4, 51. 4.I. Grid independent
The grid independence of the results is examined before the main calculations are carried out. The case using high Reynolds k-t turbulence model and wall function is repeatedly calculated with seven kinds of grids from coarse to fine. The separation and reattachment locations of streamlines denoted by X,/H and X,/H, respectively, are computed and compared for different grid sizes. From Table 1, a convergence to grid independence is fairly obvious. It may be said that the result from a 200 x 70 grid is sufficiently accurate for prediction purposes because the separation and reattachment locations for the 200 x 70 grid differ from the 300 x 70 grid by less than 0.03H and O.O8H, respectively. The former uses about 45% less CPU time than the latter. From the above, further computations pertaining to different turbulence models are carried out for the typical grid size of 200 x 7 (Fig. 2). It may be noted that more grid lines are intensively distributed near to the solid wall and back-facing curved-surface region to reflect the expected rapid change or much sharper gradient of flow field. 4.2. Results of high Reynolds
k-c f wall function model
In the first calculation, the high Reynolds k-c turbulence model is used to deal with the turbulent flow in the main flow region, and the wall function method to treat the flow near the wall where the assumption used by the high Reynolds turbulence model is not valid. The geometry of the grids is shown in Fig. 2. The upper and lower edges correspond to the solid
Technical
t_L
I
I
I
421
Note
I
I
I
I
I
I
0.10
0.00 Fig. 3. Streamline
distribution,
high Reynolds
number
k-c + wall function
model.
Table 2. Separation and reattachment positions Model
&IH
Test High Re k-a + wall function High Re k-c + one-equation Low Re k-c
1.78 1.332
XRIH
1.638 1.332
7.5C-8.00 6.598 7.991 7.170
wall; the left is the flow entry and the right is exit. The curved back-facing step, H, 16.87 mm in height, is used to normalize all the lengths. The inlet conditions are taken from the measurements. The under-relaxation iteration technique is adopted to increase numerical stability with relaxation factors a, = CI,= c(k= a, = 0.5 and ap = 0.7. The criteria of convergence is 6 = 0.0001. There are 200 x 70 grid nodes in the calculations. The grid systems are obtained by firstly solving the Poisson equation, then forcing the grid lines to the upper and lower walls. The final convergent result is obtained after 3270 iterations. Figure 3 shows the locally enlarged streamline distributions around the backward curved-surface. It can be clearly observed that there is a recirculation zone developed just behind the curved surface. The separation starts at the curved surface and reattachment occurs downstream. The separation and reattachment positions are listed in Table 2. 4.3. Results of high Reynolds k-e + one-equation model This is sometimes called the ‘zonal modelling’. The computational
0.10
0.00 Fig. 4. Streamline
domain is divided into two
distribution,
high Reynolds
number
k-c + one-equation
model.
422
Technical
Note
2.5
0
10
5 Fig. 5. Streamlines
distribution,
low Reynolds
number
k-t model.
regions. The high Reynolds number k-t model is used in the main flow region, and the one-equation turbulence model (Wolfshtein Model) is adopted in the region near the wall. The same number of 5 and n lines as in Section 4.1 is used for this model. On the upper wall, the wall function method is still being used instead of Wolfshtein’s model, since the main interest is in the flow region around and near the central body, and not along the upper wall. The boundary conditions and the under-relaxation factors used are the same as before. The final convergent result is obtained after 16,485 iterations. Figure 4 shows the locally enlarged streamline distributions around the backward curved-surface. The separation and reattachment positions are also listed in Table 2. 4.4. Results of low Reynolds
k-c model
In this section, the low Reynolds number Launder-Sharma turbulence model is used to calculate the turbulent flow. In order to reduce CPU time, the high Reynolds number k-c and wall function model is used to deal with the flow near the upper wall region as our main interest is flow around the centre body. In order to accommodate the sharp increase of gradients of variables near the wall, one has to distribute many and denser grid lines near the wall. It should be noted that the y+ value of the near-wall node must be less than 1.O. Another thing worth pointing out is that more attention needs to be paid to the grid distribution near the wall. The second-order derivatives for both the k and E equations require that grid lines inside the boundary layer be very smooth. It is precisely due to the second-order derivatives in both the k and t equations that much more CPU time is used by this model. The under relaxation factors for u, v, p, k and 6 are 0.2, 0.2, 0.1, 0.2 and 0.2, respectively. The final convergent result is obtained after more than 46,000 iterations. Figure 5 shows the locally enlarged streamline distributions around the backward curved-surface. The separation and reattachment positions are also listed in Table 2. According to the figures and tables, one may conclude that the high Reynolds number k-c + Wolfshtein’s one-equation turbulence model gives the best prediction for reattachment position. The other two calculations give smaller reattachment lengths than the measured ones. All the calculations predict earlier separation positions than the experiment in the present study. 5. CONCLUSION Turbulent flows inside an axisymmetric diffuser are numerically studied in the present study by using three different turbulent models, high Reynolds number k-c + wall function model, high Reynolds number k-c + one-equation model and low Reynolds number k-c model. From numerical calculations and comparisons made between numerical solutions and experimental measurements, the following conclusions can be drawn: 1. High Reynolds number k-c + wall function model uses the least CPU time.
Technical
Note
423
2. High Reynokls number k-6 + one-equation model gives the best numerical prediction for reattachment and uses reasonable CPU time. 3. High Reynolds number k-t + wall function model and low Reynolds number k-6 Launder-Sharma’s model under-predict the reattachment position. 4. Low Reynolds number k-6 Launder-Sharma’s model uses more CPU time but does not improve much the numerical predictions. 5. All the models predict an earlier separation positions than the experiment. REFERENCES 1. Launder, B. E. and Spalding, D. B., The numerical computation of turbulent flow. Comp. Mech. in Appl. Mech. Engng, 1974, 3, 269-289. 2. Launder, B. E., Nlorse, A. and Spalding, D. B., The prediction of free shear stress flows-a comparison of the performance of six turbulent models. Proc. Langley Free Shear Flows Conf., 1973, NASA SP320. 3. Patankar, S. V. and Spalding, D. B., A calculation procedure for heat, mass, and momentum transfer in three-dimensional Hows. Znt. J. Heat Muss Transfer, 1972, 15, 178771806. 4. Xu, D., Turbulent flow pass curved surfaces, part 1: experimental study. Technical Report, Mechanical Engineering Department, UMIIGT, 1995. 5. Xu, D., Turbulent flow pass curved surfaces, part 2: computational study. Technical Report, Mechanical Engineering Department, UMI!;T, 1995.