Shape Optimization of a Labyrinth Seal Applying the Simulated Annealing Method

In the present work, an environment for the shape optimization of a labyrinth seal is described. A program for a parameterized, automated grid generation is coupled with a commercial Computational Fluid Dynamics (CFD) flow solver and an optimization algorithm. Standard optimization strategies, like gradient-based methods, mostly are trapped to local optima. Therefore, the simulated annealing method is applied. It allows the finding of global minima or maxima of arbitrary functions. The presented optimization method is used to minimize the leakage through a three-finned, stepped labyrinth seal. For the optimization, the step position and the step height are chosen to be variable. The characteristics of the flow fields of selected seal configurations are compared and discussed against the background of the leakage behavior of the seal.


INTRODUCTION
Labyrinth seals represent an elementary component of gas turbines.Their leakage losses strongly affect the overall efficiency of the engine.Labyrinth seals have been designed in the past in a variety of different shapes with different leakage mass flow rates.The differing sealing efficiencies prove the existence of an optimal seal geometry.Against this background this article was initialized to attain an optimization software, which is able to determine a seal with a minimal leakage mass flow rate.
In the past, CFD has been successfully applied to investigate selected labyrinth seal configurations (e.g., Wittig et al., 1987;Rhode et al., 1994;Schramm et al., 2000;Rhode and Adams, 2001).The current available computer performance allows conducting CFD calculations for a high number of different seal geometries within a tolerable computing time.This allows the use of CFD for automated shape optimization of labyrinth seals, like it was frequently done to determine optimal shaped airfoils and diffusers.(Cholaseuk et al., 1999;Shahpar, 2000).
As a manual approach, Rhode et al. (1994) systematically varied the knife pitch (KP), fin width (FW), step height (SH), and clearance (C) of a stepped labyrinth seal.In all cases the step was positioned at the same relative axial position at SD/KP = 0.6.For the evaluation of the seals, Rhode conducted CFD calculations.The determined optimal seal was finally tested experimentally.This approach led to a distinct leakage reduction in comparison to a baseline design.Aboulaich et al. (2000) developed mathematical proof of the existence of an optimal shaped domain for their single cavity straight labyrinth seal.Using a finite element code, they found an optimal shape for incompressible, laminar flow conditions with minimal leakage.Since labyrinth seal flow is highly turbulent and in many cases compressible, the applicability of the results is limited.
Currently, there is no general tool available which allows optimizing labyrinth seals using a parameterized grid generator coupled with a Navier-Stokes solver.The objective of this article is to demonstrate the applicability of a shape optimization software to find the seal geometry with the highest sealing efficiency for given boundary conditions.The position and height of the labyrinth step were chosen to be the variable design parameters for this study.

SEAL GEOMETRY
For the investigations, a three-finned, stepped labyrinth seal was chosen, as it typically can be found in the secondary air system of gas turbines.Figure 1 defines the seal geometry for the investigations in the present study.The distance of the step to the fin (SD) and the step height (SH) are the variable design parameters and therefore not fixed.Schramm et al. (2000) presented numerical and experimental results for the same seal geometry,

FIGURE 1
Definition of the seal dimensions.
where the step was positioned in the middle between two fins (SD/KP = 0.5) and the step height amounted to SH/C = 3.33.In order to allow Laser Doppler Velocimetry (LDV)measurements, Schramm et al. scaled up the geometry of a real engine sized seal by a factor of four.The numerically predicted and measured discharge coefficients showed a relative deviation of less than 3.6%.
Since the seal geometry in this study was planar, no rotational effects were considered.A correction of the discharge behavior from non rotating geometries for rotational effects can be obtained by applying the results from Waschka et al. (1990Waschka et al. ( , 1993)).Waschka et al. found, that rotation reduces the discharge coefficient, if the circumferential velocity of the rotor exceeds the axial velocity of the flow.The influence of the scaling of the seal on the leakage characteristics was investigated by Wittig et al. (1982) and Willenborg (2001).
For the optimization, the step distance (SD) and step height (SH) were varied.
Step height values of 1.3 • C < SH < 5.0 • C were admitted.The axial position of the step varied between 0.2 • KP < SD < 0.85 • KP.Besides the actual labyrinth seal, the discretized flow region included an inflow region with a length of 6 fin heights (FH) and an outflow region with a length of 12 fin heights.

LEAKAGE CHARACTERISTICS OF LABYRINTH SEALS
To quantify the sealing efficiency, a dimensionless discharge coefficient is defined: This coefficient describes the ratio of the leakage mass flow rate to an ideal mass flow rate.The ideal mass flow is a theoretical mass flow for an isentropic expansion through a single nozzle: In this formula, π represents the ratio of the inlet and outlet pressure (π = p 0 / p n ), A is the cross-sectional area of the labyrinth gap and, T 0 is the temperature at the inlet.

FLOW SOLVER
For the numerical simulation of the labyrinth seal flow, the commercial finite-volume-code, TASCflow3D, was used.This software solves the compressible, time-averaged Navier-Stokes equations on a non orthogonal, body-fitted, structured grid.A second-order discretization scheme, the so-called Linear Profile Scheme, combined with Physical Advection Correction terms was applied (TASCflow User Documentation, 1996).The turbulence characteristics of the flow were modeled by the standardk-ε equations.The logarithmic wall function was used to describe the near-wall velocity.As shown in earlier publications (Wittig et al., 1987;Benim and Arnal, 1994;Rhode et al., 1994), the use of the standard-k-ε turbulence model in combination with the logarithmic law is capable of simulating the labyrinth flow.

PRELIMINARY STUDY
To assure grid independent results for the CFD calculations, a preliminary CFD study was conducted.For this study, the seal geometry displayed in Figure 1 with SD/KP = 0.5 and SH/C = 3.33 was discretized with several grid resolutions.The geometry of this seal is the same seal as presented in Schramm et al. (2000) therefore experimental values were available to evaluate the precision of the leakage predictions.At the inflow boundary condition, a constant inlet velocity of 4 m/s was set with air as a working fluid.At the outlet a pressure of 1•10 5 Pa was specified.With these boundary conditions, a pressure ratio of approximately π = 1.08 was attained.
Due to the fact that TASCflow3D solely supports threedimensional flow problems, three grid planes lateral to the flow direction were discretized.For this reason the number of grid points raises by a factor of three.The grid point number includes, besides the actual labyrinth, also the grid points in the inflow and outflow region.
The results of a variation in the grid resolution are shown in Figure 2.For the coarsest discretization with 4740 grid points, the leakage mass flow respective to the discharge coefficient C D deviated from the experimental value by more than 2%.With increasing grid resolution, the computed leakage mass flow rates matched the experimental results with an accuracy of better than 2%.Besides the error introduced by the turbulence-model, two effects influenced the precision of the results.One is the discretization error, which affected the prediction accuracy for coarse grid.The second effect is the influence of the wall boundary layer model, which particularly occurred for very fine grids.For very fine grids, the wall nearest the grid point left the region of validity of the logarithmic law and fell below values of y+ = 30.Therefore, TASCflow uses a polynomial approach

FIGURE 2
Influence of the grid resolution on the computed discharge coefficient.
to predict the flow conditions in the region between the low logarithmic law region and the viscous sublayer.This approach reduces, but does not eliminate the introduced error by small wall distances.The computational effort (Figure 3) showed an approximately linear progression with increasing grid point number.The calculations were stopped after the maximum residuum fell below a value of 5 • 10 −5 .Except for the coarsest grid, all calculations in the preliminary study agreed with the measured value within a tolerance of less than 2%.Due to this similar high precision of the predictions for differing grid resolutions, the computation time was the second important aspect in determining a suited grid resolution, since during the optimization process hundreds of CFD calculations had to be conducted.For the planned optimization, a grid resolution of about 10.000 grid points was chosen.A section of such a grid is shown in Figure 4.With this resolution, a CPU-time of approximately 25 minutes on a 375 MHz power3-II node (IBM RS/6000 SP-SMP) was necessary for the computation.With the selected grid resolution, values of y + >30 can be realized for most of the wall nearest grid points.Nevertheless, in regions of stagnation points, too low values of y+ cannot be avoided.

FIGURE 3
Influence of the grid resolution on the CPU-time.

FIGURE 4
Detail of a computational grid (10182 grid points).

OPTIMZATION ENVIRONMENT
The optimization environment consists of five major components: 1.A parameterized grid generation software.2. A preprocessing software (TASCbob3d) for the setup of the boundary and initial conditions (TASCtool). 3. The flow-solver (TASCflow3D).4. A postprocessing tool (TASCtool) for the evaluation of the current seal efficiency. 5.The optimization algorithm, which varies the geometry parameters in order to obtain an optimal seal.
Figure 5 shows an outline of the optimization procedure.Starting from an initial parameter set of SD/KP and SH/C, in the preprocessing part a grid is automatically generated and all boundary conditions are set.The CFD package TASCflow allows the pre-and postprocessing and the flow solver to be controlled by text-based input files.This important feature allowed a simple integration of the CFD package into the optimization software.
After the generation of an initial solution, a CFD calculation is conducted and the discharge coefficient C D of the seal is evaluated in the postprocessing part by determining the leakage mass flow rate and the pressure ratio.The coefficient C D serves as the input parameter for the optimization algorithm.The optimization algorithm decides whether the current design is optimal or the optimization procedure has to be continued by generating a new seal geometry.

OPTIMZATION ALGORITHM
Within the optimization environment, the optimization algorithm is the essential for the determination of the optimal seal geometry.In the present work, the simulated annealing method (Kirkpatrick et al., 1983) was chosen.The name of the method refers to the physical process of annealing.In contrast to classical approaches like gradient-based methods, simulated annealing is a heuristic method, which allows finding global optima of the evaluated functions.Starting from a current design with a design vector x = (SD/KP, SH/C) T , a randomly distributed new design in the neighborhood of x with x new = x + x is evaluated.If the new design is better than the current design, the new design is accepted as current design.Simulated annealing also accepts a poor new design with a certain probability.A poor design is accepted if a random number between 0 and 1 is lower then a calculated value of P with Therefore, there is a corresponding chance for the algorithm to leave a local extremum in favor of finding a better one.With regard to the basic idea of simulated annealing, the parameter T is termed as temperature.In the implemented algorithm, T is reduced every 30 evaluations (CFD calculations) by multiplying the actual value of T with 0.85.The rate of convergence and the probability of finding the global optimum are determined by this cooling schedule.As a starting value, T was chosen to be 0.3.Due to the continuous reduction of T , the probability of accepting a poor new solution successively decreases.In contrast to this approach, Fielding (2000) recommended a constant value for T .In a large number of numerical experiments he investigated the influence of T on the necessary number of function evaluations.He noticed that an optimally chosen fixed temperature tends to outperform implementations with a fast cooling schedule.
The random term x, which mathematically quantifies the expression "in the neighborhood of the current design", was specified by Here, the term rand is equally distributed random number in the range of 0 . . . 1.The indices min and max term the lower and upper acceptance limits of the design parameters.Term D allows one to narrow the search range.While at the beginning of the optimization process a value of D = 1 is preferred to cover a wide range of scanned geometries, in the further process smaller values of D are favorable for an accurate inspection of a smaller region next to the current solution.Thus, a value of D = 1 is chosen at the start of the optimization and later on D is reduced every 30 evaluations by multiplying the actual value of D with 0.8.
The outlined optimization algorithm is a rudimentary implementation of the simulated annealing optimization method.The CFD calculations were carried out on two processors in parallel.The complete optimization procedure required seven days.Modern variants of this strategy as presented by Ingber (1996) or Besnard et al. (1999) reach optimal solutions with a distinctly reduced computational effort.Furthermore, simple modifications like a temporary storage of computed seal efficiencies for a given geometry could avoid multiple evaluations of almost identical seal geometries.In this context, an interpolation of already predicted seal efficiencies to a new, similar geometry could reduce the computational effort (Torczon and Trosset, 1998).

OPTIMIZATION PROCESS
The boundary conditions for the CFD calculations in the optimization process were the same as those applied in the preliminary study, so the seal pressure ratio was approximately π = 1.08.Figure 6 shows the trend of the optimization parameters during the optimization process.At the beginning of the optimization the parameter D and T of the implemented simulated

FIGURE 6
Sequence of the optimization parameters.

FIGURE 7
Sequence of computed discharge coefficients.
annealing method were relatively large, so the trend of the design parameters SH/C and SD/KP reflects a random search process.Subsequently, the values of T were reduced and the probability of accepting a design deterioration declined.After about 200 evaluations the step distance leveled off at a value of about SD/KP = 0.8 while the step height still strongly fluctuated.This behavior is caused by the more distinct sensitivity of the discharge coefficient to the step position in comparison to the step height.After 600 evaluations the step height began to aspire to a value of about SH/C = 4.9.The progression of the discharge coefficient is shown in Figure 7.During the first 350 evaluations, the discharge coefficient showed a randomly fluctuating behavior.In the further process the fluctuations subsided and tended to a minimal discharge coefficient of C D = 0.404.This optimal seal had a step height of SH/C = 4.97 and a step position of SD/KP = 0.84, which is the border of the defined parameter domain.The worst observed sealing efficiency was a seal with a step position of SD/KP = 0.84 and a step height of SH/C = 1.38.The discharge coefficient for this seal was C D = 0.47.

DISCUSSION OF THE RESULTS
As a main result, Figure 8 shows a contour plot of the discharge coefficient C D dependent of the step position SD/KP and the step height SH/C.The plotted dots mark the computed seal geometries.Generally, it can be observed that the heuristic optimization algorithm covered almost the complete parameter domain.In the area of SD/KP < 0.65, the vertical contour lines show that the C D -value is nearly independent of the step height.In the area of SD/KP > 0.65, there is a distinct step height influence.Here the sealing efficiency is increasing from small to large step heights, which means the discharge coefficient decreases.
Besides the discharge coefficient, every CFD-calculation conducted during the optimization process yields a complete data set about the seal flow, which is used here for a more detailed discussion of the results.Four exemplary seal geometries (termed A to D) were chosen from the parameter domain.

FIGURE 8
Influence of the design parameters on the discharge coefficient.
Table 1 defines the geometries of the selected configurations and they are also marked in Figure 8.
Figure 9 shows the streamlines of the computed flow fields.For the configurations A and B, with a step near to the upstream fin, the flow passes the gap and hits the step.The jet is then strongly deflected and forms a stagnation point on the chamber ground below the step.Inside the labyrinth chamber the jet separates two counter-rotating vortices.With increasing step distance from A to C, the stagnation point moves along the chamber ground to and then up the downstream fin.If the stagnation point reaches the region near the fin tip respective to the next gap, part of the kinetic energy of the jet is directly transferred to the next chamber.This mechanism is known as the carry-over effect and particularly occurs in straight-through labyrinth seals, where it strongly impairs the sealing efficiency (Hodkinson, 1939;Vermes, 1961).In the flow field of configuration C, the beginning of carry over can be observed.This explains the maximum of the C D -value in this region (Figure 8).
An increasing step height results in a stronger deflection of the jet (see configurations A and B), which prevents the carry-over effect in configuration D. This configuration shows the lowest discharge coefficient of the discussed seals.Though there is no carry over for configurations A and B due to the step position,  A reason for this difference can be seen in the contact length between jet and wall after passing the gap.Along the jet high velocity gradients cause a raised turbulence level and high dissipation in the wall boundary layer as well as in the contact region to the chamber vortex.Regarding the velocity contour plot in Figure 10, it can be seen that the small distance of the step to the upstream fin restricts the jet length in configurations A and B. The large step distance combined with a sufficient step height to prevent carry over, led to the best sealing efficiency of configuration D.
Also from a construction point of view configuration D is advantageous.During operation due to thermal expansion, large axial displacements occur.Configuration D allows large axial displacements without contact of the rotor and stator of the seal.

FIGURE 10
Contour plot of the velocity in the first labyrinth chamber for configurations A to D.
To obtain an optimal step to fin position during operation, the axial displacements between downtime and base load conditions have to be taken into account when the seal is designed.

CONCLUSIONS
In a preliminary study, the leakage mass flow rate of a stepped labyrinth seal was investigated with several grid resolutions in order to assure grid independence of the results.The grid variation allowed the determining of a suitable grid resolution that results in an adequate discretization accuracy and at the same time a moderate computing effort.A comparison of the predicted discharge coefficient with experimental results from the literature showed a deviation of less than 2%.
A software environment for an automated shape optimization coupled with a commercial CFD flow solver has been applied to a stepped labyrinth seal.The grid resolution was chosen according to the results of the preliminary study.As optimization algorithm, the simulated annealing method was selected.The position and height of the labyrinth step were the variable design parameters for this study.The execution of the optimization yielded a seal geometry with a step distance to the upstream knife of 85% of the knife pitch.The optimal step height equaled five times the seal clearance.This configuration corresponded to the border of the defined parameter domain.In comparison to the baseline geometry of the preliminary study, the discharge coefficient of the optimized seal showed a leakage reduction of about 10%.
In the future, the chosen optimization algorithm of simulated annealing easily can be enhanced to three or more optimization parameters.The literature gives many suggestions to improve the efficiency of the basic simulated annealing method implemented in this study.It is expected that automated shape optimization as presented will become increasingly important in the design process of future turbomachines.

FIGURE 9
FIGURE 9Computed streamlines for configurations A to D.