Refined Turbulence Modeling for Swirl Velocity in Turbomachinery Seals

A generalized new form of the rotation-sensitive source term coefficient previously proposed by Bardina and colleagues as an extension of the standard kturbulence model was developed. The proposal made by Bardina and colleagues focused on rotating flows without significant turbulence generation, and the result was a negative-valued constant coefficient. The new functional form developed here for the coefficient has global as well as local dependence. The new model predictions of laser Doppler anemometry measurements of swirling flows in labyrinth seals were compared with the swirl distribution measurements and with the standard kmodel (i.e., no rotation source term) predictions. It was found that for the labyrinth seal cases for which detailed measurements are available, the standard kmodel gives unsatisfactory predictions, whereas the new model gives significantly improved predictions.

the rotation-dependent viscosity, and so are of greatest interest here.It is the primary purpose of present study to develop a modification of the k-ε turbulence model based on an improved rotation-dependent form of the dissipation rate equation.The improved swirl velocity will be obtained indirectly through the changed values of the eddy viscosity.
In this study, the high-Reynolds number (Re) k-ε model (Launder and Spalding, 1974) was used as the basic turbulence model.In the current development of the modified turbulence model, the focus has been confined to a labyrinth seal with enough cavities that a streamwise-periodic flow condition can be assumed.The numerical technique based on the finite volume method (Patankar, 1980) was used together with the SIMPLEC scheme (van Doormaal and Raithby, 1984) to couple the pressure and velocity fields.

PREVIOUS WORK
A series of measurements was made using three-dimensional laser Doppler anemometry (LDA) for turbulent swirl flows in labyrinth seals (Morrison et al., 1988a(Morrison et al., , 1988b(Morrison et al., , 1996)).The measurements provide three-dimensional time-averaged mean velocity distributions in the axial and radial directions, and they can be used as model validation data for comparison with the present numerical predictions.Jacquin (1990) and colleagues experimentally studied the effect of rotation imposed on a simply decaying homogeneous turbulent flow in a simple cylindrical configuration.According to their results, the existence of rotation reduces the dissipation (or decay rate) of turbulence.This result seems attributable to the role of rotation in inhibiting the transfer of turbulent kinetic energy from larger to smaller eddies.Another significant finding of their work was that there exist anisotropic characteristics in the turbulence length scales in which the axial length scale is generally larger than that of the radial and circumferential directions except far downstream where all the length scales have a similar magnitude.Very similar experimental results could be found in the experimental work done by Wigeland (1978).
In a computational effort, Cambon and Jacquin (1989) investigated the anisotropic effects of solid-body rotation imposed on a homogeneous, simply decaying turbulent flow.Their interest was focused on the anisotropic behavior of the turbulent kinetic energy spectrum incorporated with the anisotropic turbulence length scale using the eddy-damped quasi normal markovian (EDQNM) approach.A similar study was made by Teissedre and Dang (1987) using a direct numerical simulation (DNS).The common result obtained in both studies was that the magnitude of turbulence energy transfer is not the same in all directions (anisotropy).It was also found that the total turbulence energy transfer, integrated over the wave number space, from larger to smaller eddies, decreased in comparison to the case without rotation.In other words, the total amount of dissipation, which is believed to occur in the smallest eddy sizes, decreased.The results obtained by both the EDQNM and the DNS approaches are in good agreement with that of an earlier work by Bardina and colleagues (1985) using DNS.From the numerical results, Bardina and colleagues suggested an additional source term in the dissipation rate to account for the rotational effect.The new term was expressed as the product of the local vorticity and the dissipation rate.Their modeling was, however, based on solidbody rotation swirling flow assuming no turbulence generation.Therefore, it is not proper to implement their original model without an appropriate modification for the effect of turbulence generation on the dissipation mechanism.

OBJECTIVE
The purpose of the present study was the development of a new turbulence model extension which can be applied to the complex swirling flow in labyrinth seals.The basic framework of the model is the k-ε turbulence model (Launder and Spalding, 1974); a new source term was added to the dissipation rate equation to account for the rotation effect (Bardina et al., 1985).
In the development of the present turbulence model extension, the flow in a labyrinth seal was chosen as a model case so that a comparison with three-dimensional LDA measurements would allow the evaluation of some relevant constants.Then the new model was tested by comparing predictions with other measurements for various flow conditions, such as these with positive, negative, and zero inlet swirl.

GOVERNING EQUATIONS
According to the eddy viscosity concept, the Reynoldsaveraged momentum equations become, in tensor form, The Boussinesq eddy viscosity (Schlichting, 1979) assumes that the Reynolds stresses are proportional to the local gradients of mean velocity, as where δ i j is Kroneker delta and ν t is the eddy viscosity defined as ν t = c µ k 2 /ε.The transport equations for turbulent kinetic energy k and its dissipation rate ε are modeled as follows (Launder and Spalding, 1974): where the empirical constants are given as c 1 = 1.44, c 2 = 1.92, σ k = 1.0, and σ ε = 1.3.The swirling flow in labyrinth seals shows obvious anisotropic turbulent characteristics.As was discussed earlier, the primary observation made from experiments (Jacquin et al., 1990;Wigeland, 1978) and from computations (Bardina et al., 1985;Teissdre and Dang, 1987) for simply decaying homogeneous, rotating flows provides evidence that there exist scales of different lengths between the axial and the other directions, thus indicating the anisotropic characteristics of swirling flows.
In the present study an accounting for the anisotropic length scale was incorporated and is described as follows.The new turbulence source term introduced in the dissipation rate equation will alter ε values and, thus, the Boussinesq eddy viscosity given above.The altered eddy viscosity will have more influence on the swirl momentum equation than on the axial momentum equation.The reason for this different influence of the altered eddy viscosity can easily be seen if one considers the order of magnitude of each term in the momentum equations.Specifically, much more dependence on eddy viscosity is found in the swirl momentum equation than in other momentum equations.Bardina and colleagues (1985) suggested that a new negative source term be added to the dissipation rate equation of the standard k-ε model to represent the effect of rotation in a swirling flow: where = [ i j ji ] 1/2 , i j = rotation tensor, and c 3 = −0.15.
The negative constant value of c 3 = −0.15 was suggested by Bardina and colleagues only for the simply decaying swirling flow where no generation of turbulent kinetic energy occurs.Even though the dominant effect of rotation is to decrease the dissipation rate of turbulence for simply decaying turbulence, an increment of the dissipation rate was also observed immediately downstream of the rotating grids where the generation of turbulence is not negligible (Wigeland, 1978).Moreover, in labyrinth seals the generation terms are usually very strong in the cavity regions near the dividing streamline and the rotor, and their interference with the dissipation mechanism is no longer negligible.Therefore, it is more reasonable to view the role of flow rotation in general applications as either a negative or a positive source term in the dissipation rate equation.In other words, the sign of c 3 can be positive as well as negative.A detailed description for the present evaluation of c 3 is given subsequently.

TURBULENCE MODEL DEVELOPMENT
From the preliminary computation of labyrinth seal flows implemented with the constant coefficient c 3 that was suggested by Bardina and colleagues it was found that there is no single value of the coefficient c 3 that results in good agreement with the mean swirl velocity measurements.Specifically, some selected values of constant coefficient c 3 produced good agreement for the inlet region and a significantly worse result for the exit region.For other values, the reverse was found, or the agreement was poor for the entire seal region.As a result of the preliminary examination, two important facts were found: (1) the new source term in the ε-equation gives significant influence on the swirl velocity and (2) |c 3 | should be 0.2 or less for the best agreement with measurements.The final chosen grid had 204 × 25 nonuniformly spaced lines after a grid independence test (Table 1) comparing the rms variation among several different grids for the preliminary computation.It will be used for the remaining modeling and test cases.
In the definition, the superscripts f and c stand for the fine and coarse grids, and U and W represent the average velocities at the fine grid.Note that the u-rms is the area-weighted value.
The coefficient c 3 in Equation ( 5) is herein assumed to be a function that should depend on the local flow properties such as wall shear stress and swirl velocity, and on the global flow conditions such as the Reynolds number and Taylor number.The function c 3 is defined in terms of an independent variable N * to where θ = Re/Ta, and N * and g(N * ) are to be determined in order to account for the streamwise development of the circumferential (swirl) velocity and its net driving (circumferential shear) force.
In the above expression, c B is the constant term corresponding to the Bardina constant coefficient c 3 = −0.15.Therefore, the term c 3 (θ )g(N * ) is the modification made to the original model to account for the effect of the generation of turbulence as well as for the effect of the non-solid-body rotation through the dependence on θ and N * .Considering the fact that most turbulence generation is produced by the axial velocity at walls and sharp corners of any cavities (Rhode and Guidry, 1995), c 3 (θ )g(N * ) should vanish when θ = 0 because this situation corresponds to the solid-body-rotation that Bardina and colleagues (1985) considered.It is intended in the formulation that c 3 (θ ) represent the global flow characteristics in terms of the relative significance of the fluid rotation to the axial flow, while g(N * ) describes its gradual change in the downstream direction-that is, the increasing effect of fluid rotation in the downstream direction.
The following expression was used as an intermediate step: where a and b are undetermined constants, x is the axial coordinate, L is the total length of the seal, and n is a positive number to be determined from the modeling case.The function c 3 (θ ) was assumed to be linear (a + bθ ), partly because there are not enough experimental cases to allow the constants of a more sophisticated function to be properly evaluated.Specifically, one must (of course) avoid using the test case experiments to evaluate constants in a proposed turbulence model modification.
Regarding the function g(N * ), the solutions of a series of exploratory Computational Fluid Dynamics (CFD) computations were studied, and the general algebraic form of a preliminary expression for g(N * ) was found.This form involving x/L is shown in brackets in Equation ( 7).Notice that the constant a in Equation (7) must be zero in order for c 3 to equal c B for Bardina's case of negligible generation and solid-body rotation.To evaluate the constant b in Equation ( 7), the tentative (best estimate) expression for c 3 (i.e., Equation ( 7)) was used with the CFD code in a series of exploratory computations.By comparing with measured swirl velocity in the modeling experiment case (described in the next paragraph), it was found that b = 0.068 and n = 4.0 give the best agreement with the LDA measurements of swirl.
Available LDA swirl velocity measurements in a labyrinth seal were used as the modeling case (Morrison et al., 1988a).The geometry of the seal with seven cavities is shown in Figure 1, and the global flow conditions are Re = 28,000 and Ta = 7000.The flow in the model case, as well as in the three test cases in

FIGURE 1
The geometry and scales of the stator and rotor for computation.a, 79.002 mm; b, 82.05 mm; c, 83.32 mm; d, 1.524 mm; e, 0.762 mm; total seal length L, 33.5 mm.
Table 2, was computed using the finite volume method (Patankar, 1980) with the SIMPLEC algorithm (van Doormaal and Raithby, 1984).The well-known log-law of the wall was applied as the boundary condition at the rotor and stator walls.The near-wall grid points were carefully arranged so that most of the y + values   Morrison and Johnson (1996).
were within the range of 30 to 90 to make sure that realistic solutions would be obtained (Virr et al., 1994).
The next step in developing the turbulence model modification was to define and evaluate the flow variable N * that gives a more generalized streamwise development of swirl veolcity while maintaining a similar streamwise dependence as the function in brackets (from CFD solutions) in Equation ( 7).The following definition of the normalized swirl friction coefficient N * was proposed in the present study as (tooth) from Morrison et al., 1988a.)where the swirl friction coefficient N is defined as and the subscript j signifies the j th tooth cavity; J gives the total number of tooth cavities-that is, J = 7 for the present seal geometry.Observe that N j was chosen to account for the effect of the gradual streamwise variation of swirl on c 3 for the j th streamwise segment of the flowfield length L. Note that the swirl friction coefficient N j is the ratio of the net driving (circumferential shear) force to the resulting change of circumferential (swirl) momentum for the j th streamwise segment.The choice of the variable N * gives a dimensionless form which is very similar to the parameter proposed by Rhode et al. (1986) in a parametric study of the sealing performance of labyrinth seals.
After evaluating the N * values from a CFD solution of the modeling case using the best available estimate of c 3 (i.e.Equation ( 7)), it was found that the following relation between N * and x was appropriate for the modeling case: As a result, by eliminating x/L in Equations ( 7) and ( 10), the final form of the modified source term was determined to be c 3 (θ; N * ) = −0.15+ 0.068θ{1 − (1 − N * ) 4/5 }. [11] Figure 2 shows the mean swirl velocity compared with the corresponding measurement for the modeled case.The thick line at the bottom signifies the geometry of the labyrinth seal considered.A reasonably good prediction from the present modeling (Equation ( 11)) can be observed when compared with the large discrepancy from measurement of the k-ε model.Better agreement with measurement can be observed downstream where the effects of rotation and turbulence generation are fully present in the flow field.Note that no actual improvement is found by the present modeling in the tooth tip area.
Figures 3 through 6 show the radial distributions of the axial and swirl velocities at the axial center of the third and seventh teeth.Two interesting observations can be made.The first observation is that there is no improvement (and even no change) in the distribution of the axial velocity between the original and the modified k-ε models.The reason is that the pressure gradient and the axial inertia terms are dominant in the axial momentum equation, so that the local value of the eddy viscosity from the modification of the dissipation rate equation makes only a small  8) and ( 9)).

TEST CASES
Here, three other flows (not used in developing the model extension) in labyrinth seals are considered to test the general applicability of the presently modified turbulence model to other flow conditions.The computed results are compared with the measurements made by a three-dimensional LDA system.Because of the lack of alternative detailed measurements, the geometry of the labyrinth seal considered here is the same as that shown in Figure 1.The flow conditions are summarized in Table 2.
The law of the wall was applied within the range of y + = 30 ∼ 90 at the nearest grid points from the seal surface for all tested cases.The computing procedure was arranged to switch from the standard k-ε model to the present new model after 300 initial iterations so that a favorable distribution of the shear function N * could be obtained before it acts as the independent variable in the new model.The inlet boundary conditions for the axial and swirl velocity were obtained from the experimental measurements, and the radial velocity was assumed to be zero.Note that even without the inlet preswirl device (in the modeling case and test case 1), a positive preswirl can be induced in the inlet region by turbulent momentum diffusion.The outlet boundary conditions were assumed to be fully developed for all variables.
For test case 1, c 3 (θ ) = 0.1 was obtained from Equation (11) because θ is 1.46 for the given Re and Ta values.The significance of the generation of turbulence is relatively weaker than for the modeling case, where θ equals 4.0.Figure 7 shows the bulk swirl velocities obtained from the present model and the original k-ε model.It is rather unfortunate to find that the original k-ε model predicts reasonably well for this particular operating condition; thus, there is only a slight possibility for improvement here.The fact that the standard k-ε model works well for this case can be explained by an argument made in the literature.Specifically, as the intensity of rotation (swirl) increases, the structure of turbulence tends to return to the two-dimensional state where the anisotropy of length scale vanishes.Therefore, as the anisotropy disappears, the assumption of a single isotropic length scale made in the derivation of the standard k-ε model becomes more reasonable.Figures 8 and 9 show the radial distributions of the swirl velocities at the selected tooth tips and cavities for this case.As was discussed earlier, there is no difference between the CFD results of the two models for the axial velocity, so they are not discussed here.Also, no practical improvement is found by the present model in the distribution of the swirl velocity.As a result, it is apparent that the new turbulence model does not make a significant improvement over the standard k-ε model for flows with such low θ values because the standard k-ε model provides reasonably good results for swirl velocity in such cases.
Figure 10 shows the bulk swirl velocity distribution for case 2, which has negative preswirl and a θ value of 3.6.The standard k-ε model significantly overpredicts the swirl near the inlet, and underpredicts it near the outlet.The results for the second through the fifth cavities seem reasonable.The results from the new model improve the predictions at both the inlet and the outlet, whereas they are very similar to the results from the standard k-ε model in the middle cavities.The radial distribution of the swirl velocities is shown in Figure 11 for the seventh tooth-cavity as an example.Improvements are clearly shown near the rotor surface, and no practical change is found near the stator surface.For test case 3, a flow with positive preswirl was considered.In this case, θ is the same as that of case 2. Similar to case 2, the new model shows improvement near the inlet and outlet, while no significant improvements are shown in the middle cavities, as indicated by Figure 12. Figure 13 shows the radial distribution of the swirl velocities for the seventh tooth cavity.A similar conclusion can be made for test case 2. One needs to note here that the θ values of test cases 2 and 3 are larger than that of test case 1, where improvements were not as significant.

SUMMARY
A new turbulence model extension has been proposed and tested.The foundation of the new extended model is the standard k-ε turbulence model with the modification for the rotation (swirl) proposed by Bardina and colleagues.In the new, extended model, consideration of the additional effect of the generation of turbulence, which has not been considered before in such a rotation-sensitive model, was made to cover a wide range of swirling axial flows such as those of turbomachinery labyrinth seals.Unlike the suggestion made by Bardina and colleagues, the coefficient of the rotation source term is modeled as a function to be either positive or negative instead of being a negative constant value.The refined effect of rotation on swirl is obtained in an automatic way during the iterative solution procedure by using a dimensionless friction coefficient.
Compared with the standard k-ε turbulence model, the new model gives generally improved results for both the bulk and the local distribution of swirl velocity.More important, swirl velocity improvements are more evident near the rotor surface, which is where the rotordynamic characteristics of seals are determined.

TABLE 1
Results of the Grid Independence Test

TABLE 2
Specification of the Test Cases *  Morrison et al. (1988b).* * c seal tooth clearance c 1 , c 2 empirical coefficients of the k-e turbulence model c 3 coefficient of Bardina and colleagues' source, i.e.,