MODELING OF TIP CLEARANCE FLOWS THROUGH AN IMPROVED 3-D PRESSURE CORRECTION NAVIER-STOKES SOLVER

S source term Mechanical constraints dictate the existence of tip clearances ui curvilinear coordinates in rotating cascades, resulting to a flow leakage through this W . velocity in relative system clearance which considerably influences the efficiency and range xi Cartesian coordinates of operation of the machine. Three-dimensional Navier-Stokes I, diffusion coefficient in general transport solvers are often used for the numerical study of compressor equation and turbine stages with tip-clearance. The quality of numerical a turbulence dissipation predictions depends strongly on how accurately the blade tip ei j k alternating tensor region is modelled; in this respect the accurate modelling of tip total pressure loss coefficient region was one of the main goals of this work. In the present μeff effective viscosity paper, a 3-D Navier-Stokes solver is suitably adapted so that the μ i laminar viscosity flat tip surface of a blade and its sharp edges could be μt turbulent viscosity accurately modelled, in order to improve the precision of the p density calculation in the tip region. The adapted code solves the fully 0k, aE effective Prandtl numbers for k and c elliptic, steady, Navier-Stokes equations through a spacei shear stress marching algorithm and a pressure correction technique; the HS2 i rotational speed components type topology is retained, even in cases with thick leading edges where a special treatment is introduced herein. The analysis is applied to two different cases, a linear cascade and a Subscripts compressor rotor, and comparisons with experimental data are f,b,s,n faces of the computational cell (forward, provided, backward, south, north) LE, TE leading/trailing edge of blade P values at the first node off the wall NOMENCLATURE c , c 1 , c2 constants in the k-c turbulence model μ^^ g contravariant metric tensor Superscripts G production rate of turbulent kinetic energy * quantities before satisfying local continuity G c Coriolis/curvature influence in k-c equations corrections to satisfy local continuity

Mechanical constraints dictate the existence of tip clearances ui curvilinear coordinates in rotating cascades, resulting to a flow leakage through this W . velocity in relative system clearance which considerably influences the efficiency and range xi Cartesian coordinates of operation of the machine. Three-dimensional Navier-Stokes I, diffusion coefficient in general transport solvers are often used for the numerical study of compressor equation and turbine stages with tip-clearance. The quality of numerical a turbulence dissipation predictions depends strongly on how accurately the blade tip ei j k alternating tensor region is modelled; in this respect the accurate modelling of tip total pressure loss coefficient region was one of the main goals of this work. In the present µeff effective viscosity paper, a 3-D Navier-Stokes solver is suitably adapted so that the µ i laminar viscosity flat tip surface of a blade and its sharp edges could be µt turbulent viscosity accurately modelled, in order to improve the precision of the p density calculation in the tip region. The adapted code solves the fully 0k, aE effective Prandtl numbers for k and c elliptic, steady, Navier-Stokes equations through a space-i shear stress marching algorithm and a pressure correction technique; the H-S2 i rotational speed components type topology is retained, even in cases with thick leading edges where a special treatment is introduced herein. The analysis is applied to two different cases, a linear cascade and a loading are all results of large tip clearance, present due to mechanical constraints. Extensive experimental and computational work has been undertaken in the past in order to understand better the physics of flows with tip clearance, dictated by the need to go to more highly loaded blades. Recent experimental studies which concentrated on leakage flow are those of Inoue and Kuroumaru (1989) on an axial compressor, Yamamoto (1989) on a linear turbine cascade, Kang and Hirsch (1992) on a linear compressor cascade at design conditions. These and other efforts have identified the interaction of the leakage flow with the blade passage flow and the onset of vortex roll up and how these are affected by varying the tip clearance. A saddle point was observed in the vicinity of the leading edge that leads to the formation of a horseshoe vortex. Inviscid numerical models like the early ones of Rains (1954) and Lakshminarayana (1970) and more recently 3D Euler solvers, have captured tip leakage flow roll up and blade unloading since the leakage flow is characterized by a core of low loss fluid. This inviscid core has been observed in the experimental studies of Moore and Tilton (1987) on a turbine cascade and of Storer and Cumpsty (1990) on a compressor cascade. In the latter work it was also concluded that a 3-D Navier-Stokes code could model accurately flow patterns and losses without needing to model the exact geometry of the blade tip and thus include in the calculation the details of the flow in that region. The use of a thin blade, sometimes referred to as "pinching" of the blade tip, has widely been used since it allows the use of standard H-type grids. However it results to poor description of blade tip geometry and skewed grids in the tip gap. Kunz et al (1992) have shown that the accurate modelling of the blade tip geometry, achieved through the use of an embedded H-type grid, allows better modelling of the tip gap flow as well as of the vortex trajectory. Similar grids have been used by others like Moore and Moore (1989) and Hah and Reid (1991). Alternative grids have been constructed like the patched 0 and H type grids of Choi and Knight (1991) or the overlapping C and H grids of Watanabe et al (1991).
In this paper, an elliptic space-marching Navier-Stokes code capable of modelling tip clearance flows through the use of the "thin blade" approximation at the tip (Lapworth andElder, 1988, Tourlidakis and, was modified so as to model accurately the tips square edged shape. Unlike the grid topologies described above, the clearance gap was filled with an H-type grid which is an extension of the main H-type grid resulting to a single domain characterized by a small overhanging protrusion that covers the clearance region. This approach, in combination with the codes structure has resulted to a minimum increase in the codes complexity and cost of execution resulting from this modification. With respect to a previous work by the authors (Giannakoglou et al., 1992), improvements have also been made in the way equations are discretized in the vicinity of the leading edge where an ambiguity occurs in H-type grids. Turbulence is modeled through the k-e model, where the wall-functions approach is used in the vicinity of solid walls. The resulting code has been applied to two cases, a linear compressor cascade and a single stage compressor rotor for which experimental measurements existed in literature. Tip leakage effects were presented by examining the interaction between the tip vortex and the secondary flow within the blades passage and downstream of the trailing edge.

MATHEMATICAL MODELING
The steady-state turbulent compressible flow is governed by the Reynolds averaged Navier-Stokes equations. Governing equations are listed below, in a Cartesian coordinate system, which rotates with the turbomachinery row with a constant rotational speed 0. The continuity and momentum equations are a(PW1) =0 (1)

3(p W) a
aw, + a aWâ xe ax^^µ`fJ 43X efff axe J (2) -a -2peytn^I'x-P(nlri)D^+pQjQx, for i =1,2,3, where the last three terms account for body forces due to the rotation (Coriolis and centripetal forces) and 0 stand for the three rotational speed components. The Boussinesq hypothesis was used to express the Reynold's stresses, by introducing the effective viscosity which is the sum of the laminar and turbulent viscosities ( Peff = + pt) 'The turbulent viscosity is related to two local properties of turbulence, namely the turbulent kinetic energy (k) and the turbulence dissipation rate (c), according to the relation Modeled scalar transport equations are solved for k and c, according to the two-equation high-Reynolds k-e model, proposed by Launder and Spalding (1974 where G is the production rate of turbulent kinetic energy and Gc accounts for the influence of rotation and curvature on the turbulence structure. This term was also used by other researchers (Rhie et al, 1984, Lapworth andElder, 1988) and has the form where r and ie are the shear stress components at each node, normal to the Coriolis acceleration and the streamline curvature respectively. In the above formula, W is the local relative velocity module and r stands for the local streamline radius of curvature. Computational cost is reduced through appropriate assumptions with respect to the quantities involved (see also Lapworth and Elder, 1988 The energy equation is solved in terms of the rothalpy in an approximate form where additional source terms are ignored. Finally, for compressible flows, perfect gas relations are used to effect closure and the Sutherland's viscosity law enables prediction of laminar viscosity, while for incompressible flows a constant density is assumed.

GEOMETRICAL TRANSFORMATION
Cartesian coordinate systems are usually inconvenient for the numerical solution of governing equations and a generalized mapping from physical coordinates to appropriate computational ones is common practice. The transformation u i =u i (x ,x2 ,x3), 1=1,2,3 introduces the contravariant metric tensor gi and the Jacobian J and gJk=VuVut 3(u 1.u2,u3) where F stands for the diffusion coefficient and S summarizes the source terms. Such transport equations may be handled through various numerical schemes, with different accuracy levels and numerical stability properties (Patankar, 1980).

NUMERICAL SCHEMES AND SOLUTION ALGORITHM
The present method is based on a finite-volume, pressurecorrection SIMPLE (Patankar and Spalding, 1972) technique and a fully-elliptic algorithm, capable of handling separated and recirculating flows. For variable density flows, the linearization of the continuity equation restricts the method only to subsonic cases. A non-staggered grid arrangement is used, which is more suited to turbomachinery applications. The space-marching character of the algorithm introduces a primary flow direction, from inlet to outlet of the flow domain and an iterative process is established, based on a multi-pass treatment of successive cross-stream planes. The domain is streamwise swept, in an iterative fashion, until all variables are converged.
By assuming a pressure field, velocities on each examined cross-stream plane could be estimated by solving the momentum equation through an approximate factorization, namely the Alternating Direction Implicit (ADI) procedure (Peaceman and Rachford, 1955) where relaxation is also required. For compressible flows, the new density is calculated directly through the state equation in terms of pressure. In order to satisfy the mass conservation equation, corrections are proposed for all variables involved, in the form (8) under certain constraints, so that the already satisfied equations are not damaged. Primed quantities are used to promote local mass conservation at nodes lying on this plane, for which the discretized continuity equation results to an equation in terms of the pressure correction P. This equation contains extra entries resulting from the velocity-pressure coupling scheme (Rhie and Chow, 1983) in order to suppress pressure oscillations. Other auxiliary pressure correction steps, namely a three-dimensional and a one-dimensional in the sense of global conservation, have been implemented in order to accelerate convergence (Lapworth and Elder, 1988).
while the contravariant velocities play the role of the convector field in the transformed domain. All governing equations are typical transport ones, where the divergence of the convection and diffusion fluxes counterbalances the production or dissipation of the independent variable, which forms the source term of the equation. In this way, the transformed form of the conservation equation of any scalar 4) could be cast in a similar form, namely a (pJWi4)= a (Jrgi ±)+JS (7) aui a auK)

BOUNDARY CONDITIONS
At the inlet plane, all flow variables need to be specified. Velocity can easily be adapted to any available experimental data, while for turbulence, uniform (k) i n and (e) i n profiles are specified by imposing the inlet turbulence intensity (r u) i n and the level of inlet turbulence viscosity as a multiple of the local molecular viscosity, namely where W in is the inlet average velocity. Periodicity conditions are imposed for the now upstream and downstream of the blades and inside the tip gap.
Boundary conditions for the turbulence quantities are applied at the first node away from the wall, through the wall function technique (Launder and Spalding, 1974). This results in a considerable economy in computer memory requirements, since the number of grid nodes required to resolve the flow in the proximity of solid walls is drastically reduced. In a wall function formulation, the basic idea is to apply any boundary condition at the first node away from the wall (node denoted by P) by appropriately eliminating any contribution from the boundary node itself. The wall shear stress TW is assumed to remain unchanged between the boundary node and node P (for compressible flows what remains constant is TW /p instead of TW ) and is found from the formulae

ln(Ey, )
where the non-dimensional distance ypT of node P is given by Yy where (y") is the distance of the node P from the wall, measured along grid lines, for the sake of convenience. For the derivation of the above formulae, the assumption of "turbulence production= turbulence dissipation" is made for the logarithmic part of the boundary layer, leading also to the following expressions for k and c at node P 3i< 3n The latter equation is the Dirichlet type condition for the eequation that applies to node P directly, thus eliminating any contribution from the boundary node. A similar effect for the k-equation is obtained by using zero-Neumann condition, instead of equation (9a). For the velocity equations, the contribution of the boundary node is eliminated by using the wall shear stress definition, approximately analyzed in the three Cartesian coordinates (Lapworth and Elder, 1988).

GRID GENERATION AT THE TIP REGION
A discussion on the various griding techniques which are in use for tip clearance calculations is given in the introduction of this paper. Here, the modelling of the correct square shape of the tip was achieved by filling the clearance gap with a local Htype grid, which was "blended" with the main H-type grid covering the calculation domain. Figures la and lb show the way this extra grid is implemented in both the physical and the transformed domain, where the clearance grid (represented by the hatched region) forms a "balcony" above a notional blade. The task of grid generation is split into two parts. Firstly, the main grid (blank region of fig.1) is generated by stacking 2-D grids in the spanwise direction with appropriate clustering near the hub, the shroud and the blade tip. Then, the clearance grid is generated by linear interpolation between the grids corresponding to the pressure and suction sides of the blade at blade tip height. The main and clearance grids have the same grid density in the spanwise and through -flow directions. In the blade-to-blade direction, grid density of the clearance grid can be specified independently of the main grid.
This approach, in combination with the use of a single string for storing grid nodes and flow variables in the code, allows for a minimum increase in memory requirement. The approach is far more accurate than the use of "pinched" grids where the blade tappers to a single line at its tip, since the real squareedged tip geometry is retained and the resulting grid is nearly orthogonal. It also seems very economic, compared to other similar techniques like the embedded grids proposed by Kunz et al (1992), where a number of redundant nodes cannot in general be avoided.

THE PARTICULAR TREATMENT AT THE LEADING AND TRAILING EDGES
The use of H-type grids seems to be more appropriate in combination with a space-marching algorithm, but there are considerable inconveniences in the neighborhood of the leading and trailing edges, especially if these are thick. In a previous work by the authors (Giannakoglou et al, 1992) the same problem was handled by eliminating the real leading edge node and by forming particular cells of triangular shape between the first node before the leading edge and the first node after it, along the solid wall. This treatment, which was common for both the leading and the trailing edges, has the disadvantage of requiring grid lines that are very close to the flow direction near the leading and the trailing edges. This problem is related to the ambiguity of calculating metrics locally; Kunz and Lakshminarayana (1991) have already presented a way of circumventing the problem.
Here, a technique to alleviate this ambiguity is introduced, which is compatible with the finite-volume discretization scheme in use. For the sake of simplicity, the proposed technique will be demonstrated on a two-dimensional cascade. The extension to three-dimensions is straightforward. The discussion will be divided into two parts, one related to the satisfaction of the mass conservation equation and the other to the analysis of non-orthogonal cross-terms in the diffusion part of the transport equations.

Mass Conservation Analysis
The hatched region in figure 2 shows the finite-volume which is adjacent to the blade at its leading edge. The cell is drawn in the transformed domain and its shape deviates from a square, with its lower side coinciding with the solid boundary, so as to contribute to a fully covered computational domain.
The continuity equation, when discretized over the hatched finite volume, reads where e,w,f,b stands for the east, west, forward and backward faces of the cell, as indicated in figure 2. For any internal cell, this equation uses contravariant velocity components at the midnodes found by an interpolation scheme which is compatible with the accuracy of the whole discretization; there, of course, Dl=0=1. In the particular cell of fi ure 2, the term (JpW I ) w^t; is approximated as (JpW^) BL /2 and the contravariant component (JpW i ) at mid node B is found as: (JP W1)g= ! (JP W') A + 4(Jp W')r^= I (JP Wl)e

Non-Orthogonal Terms
The non-orthogonal terms in the diffusion part of all transport equations, integrated over the hatched computational cell of figure 2, read

+(J µg t ')^(^ jr -O,,,) -(Jµ9 ")W(OfWt b,,)
For the particular face (w) or (bw-fw) the term (Jµg 13 )w(t fw (Dbw ) is approximated through the following expression where the first term contains geometrical information, like the metric g 13 which is calculated using coordinates for grid points belonging to the solid wall or its neighborhood, while the second term represents information coming only from nodes at the "periodic boundary" region. In this way, the aforementioned ambiguity of metrics calculation is circumvented.

RESULTS AND DISCUSSION
The code has been applied to two test cases, one being the linear compressor cascade of Flot (1975) and the other the compressor rotor of Inoue et al (1985). The code has been run on an Alliant FX-80 mini supercomputer. Run-times were of the order of 0.8 msec/node/iteration for a four-processor complex.

Case 1
The first problem deals with a low speed compressor cascade for which an experimental investigation of the secondary flow behaviour was accomplished in Ecole Centrale de Lyon (Flot, 1975). The highly loaded cascade consists of NACA 65-12-A10-10 blades, with a chord length of C=0.13 m, a stagger angle of 30 deg. and a solidity of 1.25. The span/chord ratio was 2.1 and two cases were examined, one without tip clearance and another with a clearance equal to 1,1% C. The cases correspond to test problems (B) and (C) as referred to in the papers by Papailiou et al (1977) or by Pouagare and Delaney (1986); the only difference being that an exact swirl velocity profile was imposed at the inlet (lying at a distance of 2/3 C chordwise upstream of the leading edge) for case B, instead of a constant inlet flow angle mentioned in the above references. A perspective view of the hub and blade surfaces, along with a magnification of the 3-D grid used to fill the tip gap (if any) is shown in figure 3. The 2-D grid that coincides with the hub was stacked in the spanwise direction in order to generate the 3-D grid. In the tip region different clustering was used in the spanwise direction, depending on whether the case examined was with or without tip clearance. Dimensions of the main H-type grid were 25x50x48, in the pitchwise, spanwise and throughflow directions respectively. The size of the clearance grid used to fill the tip gap was 10x8x24, which results to a total number of calculation nodes of the order of 6x105 . Figure 4 shows a comparison between predicted and measured velocity distributions along the pressure and suction sides of the blade. Figure 5a shows the measured velocity vectors at the mid-span and at a plane lying at a distance of 0.5mm from the shroud. The latter field could be identified as those vectors that deviate from the "throughflow" direction and generally have smaller magnitudes. These could be directly compared to figures 5b and Sc showing predicted velocity vectors at the corresponding blade-to-blade planes. A comparison between figures 5a and 5c demonstrates that the accurate modelling of the blade tip geometry has allowed the reproduction of the characteristic patterns of clearance flows. Firstly, strong secondary flow effects can be observed in figure Sc, where flow is directed from the pressure to the suction side. This flow strongly interacts with the jet emerging from the tip clearance near the leading edge at the suction side; this jet results from the large pressure difference between the pressure and suction sides of the blade, close to the leading edge. The "collision" of the secondary flow with the leakage jet forms a stagnation line that starts at the leading edge and develops downstream halfway between the pressure and suction sides. Upon exiting the blade section, this line curves and forms a spiral, near the suction \\\ \ side. Figure 6 presents the pressure coefficient contours on the shroud for cases B (figure 6a) and C (figure 6b). This dimensionless pressure coefficient is defined as the pressure rise from the inlet static pressure divided by the dynamic pressure corresponding to the mean inlet velocity. In figure 6b, a low pressure trough may be observed at approximately 10% of the chord from the leading edge. This trough corresponds to the generation of the tip leakage vortex and has been observed by other researchers as well (for example Kang and Hirsch, 1992).  Figure 7 shows the secondary velocity vectors at a plane at x/cx =0.828 for case C. The passage vortex is placed slightly closer to the suction side than that shown by Pouagare and Delaney and the tip leakage vortex is more concentrated in the corner of the passage. The two contrarotating vortices coexist in the passage without considerable mixing. The passage and mass-averaged total pressure loss coefficient c as defined be Pouagare and Delaney (1986), will be used for the evaluation of the loss mechanisms. This loss coefficient was calculated in the trailing edge cross-plane and at a plane lying at two thirds of the chord downstream of the trailing edge and is given in the following table: The present calculation provides losses that are quite close to the ones calculated by Pouagare and Delaney, especially at the plane lying downstream of the trailing edge. In case B, it seems that the present method results to approximately 5% less losses than the other method. This could be explained by the fact that a partial separation that occurs on the suction side close to the trailing edge at mid-span was not captured by the present method and thus less losses are produced in this region. However the loss ratio between the trailing edge plane and the plane further downstream is the same in both methods. It is also to be noted that there is a difference in the inlet flow angle between the two calculations for case B; in the present study, the mid-span inlet flow angle was 56.6 deg and the imposed swirl caused a deviation in angle of approximately 3 deg near the endwalls. This angle was constant and equal to 56.2 deg for Pouagare and Delaney (1986). For case C both methods predict the same total pressure loss coefficient at the plane 66% of chord downstream the trailing edge, where all flow phenomena related to losses production are completed. At the trailing edge cross plane, losses are greater than those calculated by Pouagare and Delaney, who claim that their reduced losses are due to the blowing of the suction side separation, which of course does not appear in the present calculation. This could be also observed in figure 8, where contours of total pressure loss coefficient are presented at the trailing edge cross plane. Losses of the present calculation are more pronounced close to the endwalls where the interaction of tip vortex and the passage flow occurs.

Case 2
The tip clearance flow was studied in a low-speed axial compressor. A 12-blade isolated rotor of 449mm casing diameter and a hub-to-tip ratio equal to 0.6 was used for the taking of extensive measurements by his team (1985, 1989). The rotational speed was 1300 RPM and the mass flow rate was 22.44 Kgs 1 . Among the two rotors for which measurements are published in the aforementioned works and which were both designed for free vortex operation but differed in the blade solidities, the one (Rotor A) with the smaller solidities was examined herein. The blade profile changed with radius and belonged to the NACA 65 series. The calculation grid was thus formed by wrapping different 2-D grids, depending on the radius, on cylindrical surfaces. Different tipclearances were examined, from z=0.5 mm up to i=5 mm. Figure 9 shows a view of the grid used, depicting the hub and the blade surface grids. Again it can be observed that the grid is clustered at the blade tip, in order to accommodate the tip clearance grid.
Results of the calculation are shown in the figures that follow. Figure 10 shows the pitch angle distributions on planes at two different distances downstream of the trailing edge for the two examined tip clearances. Figures 10b and IOd are cross planes immediately after the trailing edge of the blade, while figures 10c and l0e correspond to cross planes a bit further downstream. All these figures present only the upper -15% of the blade span, plus the tip region. Measured patterns of the pitch angle distribution, for the case of the 3 mm tip clearance are shown in figure 10a (Inoue et al, 1985). Figures 10a and 10b depict the same behaviour which is asymmetric about the vortex center, that was due to the inclination of the vortex axis. In the present calculation, values of pitch angle are slightly underestimated. For the small clearance (figures 10d and 10e) the tip vortex almost vanishes while the vortex caused by the secondary flow is even more pronounced. It can also be observed that the predicted vortices diffuse as the flow moves downstream. Figures 11, 12 and 13 compare pressure coefficient contours between experiments and prediction on the shroud, for t=.5 mm, i=2 mm and t=3 mm respectively. The pressure coefficient is defined as the pressure rise from the inlet stagnation pressure divided by the dynamic pressure corresponding to the blade tip speed. It can be observed that the minimum pressure trough moves away from the leading edge with increasing tip clearance, a fact which has also been observed in the experiments. Figures 14 and 15 compare relative velocity vectors for T=3 mm at two radii lying in the clearance gap. It can be seen that the predicted flowfield has the same trends as the measured one. Finally, with respect to the general 0.9 characteristics of the flow, like axial velocity, relative peripheral velocity and exit flow angle, it can be seen from figures 16 and 17, that all predicted results agree well with measurements.
In order to compare losses predicted by the present calculation, the pressure loss coefficient defined by ( AR2 V^-QR1 V81)-Ap f p is compared with the measured values in figure 18. Locations 1 and 2 correspond to planes upstream and downstream of the rotor and Opt stands for the total pressure rise through the stage. Non-dimensionalization is done using the peripheral speed of the blade tip. It can be observed that losses are generally overpredicted by a factor of 2 in the core region. At the hub, the shear layer losses are comparable to the measured ones, while tip-clearance losses are overpredicted. With respect to the position of the pressure loss peak, this is better predicted for the .5mm clearance case than the 2 and 3 mm cases, where this peak is slightly shifted towards the endwall.

ACKNOWLEDGEMENTS
This work was funded by the EEC BRITE-EURAM Project AERO-0021-C(D), entitled "Tip Clearance Effects in Advanced Axial Flow Compressor". The authors wish to thank the project partners for permission to publish this paper. The basis for this work was code "TURBO3D" developed in Cranfield Institute of Technology.

CONCLUSIONS
The tip clearance flow in two different test cases, a linear cascade and a compressor rotor was examined using a 3-D Navier-Stokes solver. Existing software was adapted so that the flat tip surface of a blade and its sharp edges could be modelled in detail and this improved the code's prediction capabilities in the clearance region with a minimum increase in computer running costs. Improvements have been made to the