DESIGN OPTIMISATION OF LABYRINTH SEALS USING LES

Labyrinth seals are extensively used in gas turbines to control leakage between components. In this research, the effects of geometry on the sealing performance are investigated. To obtain the best sealing performance, an investigation is undertaken into the possibility of optimising labyrinth seal planforms using a genetic algorithm (GA). Large Eddy Simulation (LES) is used for its ability to accurately capture the complex unsteady behaviour of this type of flow. Three hundred LES calculations are carried out. By making use of a large number of processors, an optimum geometry can be achieved within design cycle timescales. The optimised design shows a leakage reduction of 27.6% compared to the baseline geometry. An immersed boundary method (IBM) is used with LES to generate complex geometries on a background Cartesian grid. The GA is inherently parallel, and this enables the exploitation of the reliability and accuracy benefit of LES as demonstrated. Mathematics Subject Classification. 35L05, 35L70. Received July 24, 2020. Accepted December 8, 2020.


Introduction
High speed turbomachinery produces power from high temperature and pressure flows. A crucial factor in increasing the efficiency of these machines is reducing leakage. Consequently, labyrinth seals are used in gas turbines to separate regions of different pressure to control such leakage. They operate by dissipating as much of the kinetic energy of the fluid as possible. They have been explored extensively experimentally, and the effect of a variety of different geometries has been investigated, such as clearance [2,26,42,45], tooth configuration [2,6,25,37], and the groove geometry [18,23,24]. Although these important geometrical parameters have been identified, no single parameter has been found to enable an optimum seal design to be obtained. Combining positive trends can be counter productive. Aside from experimental methods and empirical relationships, CFD (Computational Fluid Dynamics) methods, based on the Reynolds-Averaged Navier-Stokes (RANS) equations has long been used to aid improvements. Detailed measurements within seals are generally lacking, making flow-based design impossible.
For labyrinth seals, there is a general lack of publicly available detailed experimental data to make comparisons with. For example, turbulent stress data would be beneficial. The physical size and complexity of seals makes data acquisition challenging. To increase our understanding of these flows, we therefore pursue computational methods. A range of RANS and LES have previously shown RANS to predict incorrect leakage trends and has significant variation in error (up to 30%) with turbulence model [32]. Poor linear eddy-viscosity RANS performance can be attributed to extreme flow acceleration, separation, streamline curvature and turbulence anisotropy. Such flow phenomena are naturally resolved using LES. Given the complex relationship between design parameters, optimisation methods used in conjunction with LES simulations allows multiple advanced designs to be examined with confidence in predictions. Even though LES is computationally expensive compared to RANS, it provides detailed insight into physical processes due to greatly reduced turbulence modelling impact. Here, we perform near DNS resolution LES. Each simulation is however of relatively small scale, enabling distribution across PCs with a short turn-around time and low cost relative to experiments.

Optimisation
Rhode et al. [22] attempted to optimise labyrinth seals through experiments, with the seal width, step width, step height, and clearance size set as the variables. An optimal seal was obtained which resulted in a 60% improvement in leakage compared to the baseline design. Experimental methods, however, have two key drawbacks -they require substantial user involvement, and they may miss large portions of the design space. There is another extreme optimisation method which involves exhaustively searching the entire design space. Although this is, conceptually, a straightforward process to automate, and no part of the design space is missed, for more than trivial problems it is likely to be impossibly expensive.
Between these two extremes lie a range of self-guiding optimisation methods, which aim to sensibly and automatically explore the design space while avoiding the extreme cost of an exhaustive search. Coupled with computational physics approaches, these can be very effective tools for largely automated design optimisation. With the increase in computational power that is available, these self guiding optimisation methods have been used widely, and have made automated optimisation realisable. It is essential to select the appropriate optimisation algorithm as it will directly affect the optimisation performance. A detailed study of optimisation methods for a variety of applications is given in Jin et al. [10].
The hill climbing algorithm is probably the simplest of the optimisation methods. Hill climbing is an iterative algorithm which attempts to find a better solution by comparing the current objective function value with the surrounding nodes, and moves the design in the direction which produces an improvement. By repeatedly doing this until no more improvement is possible, an optimised design is found. Hill climbing is an efficient way of finding a solution, however, in more complicated design spaces with many peaks and troughs, it is easy for it to get trapped in a local optimum, with the solution dependent on the location of the initial design. Schramm et al. [27] optimised labyrinth seals to minimise the leakage in stepped geometries using a simulated annealing (SA) optimisation algorithm. Two variables were taken into consideration in the optimisation process: the step position and the height of the step. The best sealing performance was achieved when combining a large step height and a large step distance to the upstream knife. However, when the range of the parameters is high, SA can become less effective. The SA algorithm is able to accept solutions which are worse than the current, with the expectation that a more optimum solution exists beyond the current one [12]. This decreases its efficiency at finding a local optimum, but gives some protection against becoming trapped in a locally, rather than a globally, optimal design. Hill climbing and related methods are sequential and hence not suitable for parallelisation, which makes them prohibitively expensive in real time for LES. Obayashi and Takanashi [19] employed a GA to design turbine blades more efficiently, enabling convergence towards the best possible solution through a process of evolution. The use of GAs is becoming increasingly common. For example, Pierret [20] explored the ability of GAs to optimise a compressor blade design, which was found to perform well. Watson and Tucker [38] optimised the performance of film-cooling slots for cutback trailing edges using a GA. Six hundred LES simulations were conducted. The optimised design showed a considerable improvement over the previous experimental geometries. physical trends. The flow characteristics of labyrinth seals lend themselves well to solution using LES. Tyacke et al. [33] noted that LES has the advantage of providing deeper insight into physical processes than RANS when investigating these types seals. Thus, LES is used in this study to conduct a large number of simulations. LES is typically more than two orders of magnitude more expensive than RANS [33], and complex geometries can push grid requirements and mesh construction times higher. IBM is used here to reduce mesh construction effort. In the current case, a single larger mesh is used, resulting in over-resolved fluid regions. It is hence an ideal approach, as hundreds of labyrinth seals will be tested in the optimisation process.
The design optimisation was performed using a GA. Significant benefits of using a GA are that a wide design space can be explored, reducing the possibility of designs becoming trapped in local minima and providing a large data set for lower order modelling. Also, the GA allows almost perfect parallelism. Thus, in this way, the advantages of LES, IBM and GA are coupled together, and are capable of providing highly accurate and efficient design optimisation.

Governing equations
The incompressible governing equations for LES are written in the following form: x, y and z correspond to the streamwise, cross-stream and spanwise directions respectively (here, equivalent to axial, radial and tangential components). For LES, the spatial filtering operation separates the large resolved filtered scales and small scales which are not resolved. A subgrid scale (SGS) model is used to model the overall dissipative effect of the small scales. The SGS model also aims to sensibly model the interaction between the resolved and unresolved scales. β represents a pressure gradient source term used to drive the axially periodic flow.

Turbulence modelling
In LES, the energy transfer could involve nonlinear interactions between the resolved and unresolved scales. The complex flow in labyrinth seals involves flow separation, high streamline curvature, and numerous other unsteady phenomena. Therefore, a nonlinear SGS model is employed to model turbulence anisotropy and energy backscatter to make the prediction more accurate. However, on their own, the nonlinear models do not generally dissipate enough energy from the flow, which can contribute towards solution instability. To remedy this, a mixed nonlinear SGS model [32] is used in this study which combines the linear Yoshizawa model [43] and the nonlinear Lagrangian averaged Navier-Stokes (LANS)-α model [4]. In their review paper, Meneveau and Katz [16], conclude a priori tests in a variety of flow conditions show that mixed models are, on the whole, superior to pure eddy viscosity models. Such nonlinear models are also simpler to implement and more computationally efficient than models such as dynamic Smagorinsky [7]. The mixed-LANS-α model has been used successfully for a variety of applications including turbulent channel flow [11], plane jets [14], turbine cooling ducts [35] and electronics [36] involving similar complex, low Reynolds number, turbulent flows.

Yoshizawa model
The SGS stress τ ij can generally be expressed as: where L and N L represent the linear and nonlinear terms, respectively. Here, the linear contribution to the SGS stress is provided by the Yoshizawa [43] model as L = 2k SGS δ ij /3 − 2µ t S ij . Here k SGS is the subgrid kinetic energy, δ ij is the Kronecker delta, S ij = 0.5(∂u i /∂x j + ∂u j /∂x i ) and µ t = ρC µ ∆k 1/2 SGS , where C µ = 0.07. The nonlinear term is discussed in Section 2.2.2.
This gives a linear and purely dissipative SGS model that introduces history and non-local effects by use of an extra transport equation for the SGS kinetic energy k SGS .
where P k SGS = 2µ t S ij S ij and ε = C ε k 1.5 SGS /∆ are production and dissipation terms, respectively. The LES filter length ∆ = l ε,LES = l µ,LES relates the smallest resolved scales with the largest unresolved scales. The constants C ε and C µ are 1.05 and 1.07, respectively.

LANS-α model
As presented by Liu et al. [14] the nonlinear terms, N L, for the LANS-α model can be expressed as: where α = ∆ and the repeated term denotes an implicit summation. This nonlinear term is derived from the convection of SGS vorticity by the resolved scales.

Solver details
The governing equations are solved using the NEAT code as presented by Tucker [31]. It has been shown to be well suited to LES of this flow when compared to a variety of solvers [32]. This is a structured, finite volume, parallel CFD solver. Spatially second order central differences are used for the calculation of fluxes, and the Crank-Nicolson scheme is used to integrate the solution in time. The SIMPLE algorithm is used to couple the velocity and pressure fields. As periodic boundary conditions are employed in x and z directions, a cyclic TDMA algorithm is used in these two directions. TDMA is used in the y direction. To parallelise the code, OpenMP is implemented to speed up the simulation as provided by Tyacke [34].

Genetic algorithm
A GA is used to optimise designs automatically. GAs were inspired by Darwin's evolutionary theory, and were first proposed by Professor Holland from Michigan University in 1975 [9]. Through the biological process of evolution, it enables gradual convergence towards the best possible solution for an optimisation problem by employing natural selection. GAs can search with robustness, have a global optimisation ability, and are not limited by the need to calculate derivatives and associated problems. With a proper setting of constraints, they are able to rapidly explore a large design space whilst also quickly rejecting suboptimal parameter sets. As each solution can be evaluated separately, it lends itself well to parallel architectures. The large data sets generated can also be used to inform design practices, lower order modelling and to provide a deeper understanding of factors impacting complex flow physics.
For the cases laid out here, a custom procedure was developed to link the main CFD code (including the IBM) with the GA. Figure 1 shows a flow chart of this process to give an overview of the coupling between the codes.
The implementation of GA can be considered as consisting of four basic stages.  1. Prior setup. Before the optimisation begins, the design space is encoded as a stream of data, referred to in the literature as a genome. To determine the extent of the design space, a systematic study of the key dependent variables was conducted. Each member of the initial set of designs is stochastically generated by giving a random value to each gene within its genome. In this case the initial generation consisted of 100 members. Typical members of this initial population are shown in Figure 2. Subsequent generations consisted of 30 members. 2. Evaluation of the fitness parameters.
The fitness parameter describes how well a design performs. Designs which have better fitness parameters are more likely to pass on their genes to the next generation. In this case, the fitness parameters were evaluated by performing an LES on each design. Each design was then ranked by its fitness. In terms of the design space, the individuals with the best performance tend to be clustered in regions where the objective function is more likely to achieve high values. 3. Selection, crossover, and mutation.
The key part of the GA is the manner in which the genomes of the next generation are determined from those of the previous population. The first stage in this process is selection, in which the best performing members of the population are chosen to be parents based on their objective functions. Once the parents have been selected, their genes are recombined in a process known as crossover to produce offspring. For each parent's gene, the genotype had a 50% chance to pass on to the child's genome.
To allow further exploration of the design space, mutations of the child's genome could also happen. A complete mutation, where a gene is completely lost and replaced with a random value occurs in 0.1% of cases. A swap, in which the genotypes of two neighbouring genes are swapped in the genome, occurs in 0.5% of cases. 4. Termination condition and output.
In this study, the GA is run until N gen is reached, at which point the optimised design is output.
Unlike the hill climbing or simulated annealing algorithms, each generation of individuals produced by the GA can be evaluated simultaneously. In this situation the whole process can be conducted in a perfectly parallel fashion. The evaluation of each generation of thirty individuals took, in total, about 60,000 core hours. However, in this work only 8 cores were used for each individual seal simulation. By making use of the inherent parallelism of the heuristic, the real time for each generation was reduced to around ten days.

Immersed boundary method
The IBM combines with the GA to improve design efficiency. IBM has the benefit of using simple rectilinear grids to model different complex geometries, which can reduce meshing effort. The IBM employs the ghost cell method of Ferziger and Peric [30]. Validation of the IBM is provided in Dai [3] using flow in an inclined channel, lid-driven inclined cavity flow and flow over a circular cylinder. Figure 3a shows the basic IBM method. Cells near geometric boundaries are taken as ghost cells (triangular symbols) and internal cells require no treatment (squares). The flow is reconstructed at the geometric boundary, O by quadratic reconstruction of known flow variables (crosses) outside the boundaries. Figure 3b shows linear reconstruction for illustrative purposes. Quadratic reconstruction uses the 5 nearest points following [30]. For example, all velocity components at walls are zero, hence internal ghost cell values, G, can be interpolated from external values X1, X2, to enforce no-slip and impermeability conditions. Interpolation coefficients can become large where X1, X2 are near O, hence new ghost cells are chosen to alleviate the issue.
Ghost cells are updated for every sweep of the solver algorithm to allow the flow properties of the ghost cell to continually impose the position of the boundary. The ghost cells are never more than one sub-iteration behind the main fluid simulation. Thus an accurate surface is preserved.
The pressure can be updated according to the solution of the local flow field without any interference from the IBM. Although acceptable results can be obtained, areas of high pressure can be found within the surface. This is because ghost cell values are forced directly to suit the local flow, which results in continuity violations. Here, we take the pressure as a distance-weighted combination of nearest nodes. This eliminates large pressure spots within the surface and hence models surfaces with greater accuracy. Further implementation details can be found in Preece [21].

Computational setup
This optimisation study starts with a simple geometry, as shown in Figure 4a, studied by Gamal and Vance [6]. This baseline seal has been studied using a variety of LES models in Tyacke et al. [32]. Validation data used was for a 4-tooth seal. For the LES, a quarter-length periodic domain and associated mean pressure gradient was used. The periodic domain also alleviates numerical stiffness by transferring variables either side of the seal blade, requiring the assumption of incompressible flow. Key turbulence is generated near walls and convected within cavities and is hence all low-Mach. Typical compressible solvers introduce significant numerical smoothing near walls, often scaling with 1/(M a 2 ) and would hence lead to poor LES resolution. Reduction of such smoothing using low-Mach pre-conditioning has led to instability in previous studies by the authors, and the current setup has proven to be superior [32].
For the cases tested, the dimensionless mass flow coefficient, f , is calculated as given in Gamal and Vance [6] (originally presented in [44]) containing the validation data: whereṁ is the mass flow rate through the seal, R is the gas constant, and A i is the clearance area. The seal clearance δ = 0.1 mm, the seal blade height H = 12.7 mm and the seal blades are H/4 wide giving a total domain length of 1.25H. Neither the lower shaft nor the upper seal are rotating. For seal A, [6] Re H = ρu bulk H/µ = 690. Periodic boundary conditions were employed in both the streamwise and spanwise directions, and the noslip condition is applied at walls. For the periodic streamwise direction, the mean pressure gradient, is fixed, matching [6] (β = 1.387 × 10 3 kPa/m). This maintains a constant mass flow rate with an equivalent pressure ratio P R = P in /P out to measured conditions (P in,exp = 207, 000 Pa, P out,exp = 101, 342 Pa, T in,exp = 294.95 K, m exp = 0.004617 kg/s).
Several other studies consider incompressible flow [1,8,17,28], and we take this starting point whilst validating against available measured data. Results for the baseline seal match the measurements of Gamal and Vance [6] taken at the lowest pressure ratio PR of 2.04 equivalent, to within approximately 1% as shown previously by Tyacke and Tucker [32] using the same setup. LES has therefore been shown to be a reliable method to study these types of flows. Once validated the optimisation simulations were made at a constant pressure ratio of 2.04.

Mesh, resolution and runtime
The baseline mesh is 359 × 399 × 42 nodes in the x, y, and z directions respectively. These correspond to an average grid spacing of ∆x + = 12, ∆y + ≈ 1, ∆z + = 25 in the seal tooth region. In the cavity zone, where the shear is less intense, the average grid spacing is ∆x + = 5, ∆y + = 15, ∆z + = 12. At walls, ∆ + normal ≈ 1. These spacings are within LES-Quasi-DNS resolution requirements, and the total mesh size is around 6 million cells. Figure 5 shows the structured mesh of the computational domain in this study with every 4th line plotted. In the x-and y-directions, more cells are applied around the seal blade and in the near wall region to capture the turbulent shear layer and small scales near the wall. The CFL number, which is defined as u∆t/∆x, is set to CFL ≈1. For the seals with an inclination angle, the mesh has been refined in the x direction and contains 512 × 399 × 42 nodes in the x, y, and z directions, and the total mesh size is around 8.6 million cells.
A key time scale is the time taken for the development of the correct mass flow rate from the initial conditions. This can be inferred from Figure 6. Here t * = L gap /U max represents one through flow time for the gap. The flow is developed by 15t * , and then statistics are averaged for at least another 15t * . Around 140 hours using 8 CPU cores are required for each 15t * .

Geometry of the labyrinth seal
Some key cases, shown in Figure 4, have been run to investigate the effect of various geometrical parameters. These geometries contain different inclination angles (θ), numbers of grooves (n), groove widths (l), groove geometries, and corresponding positions. The width of the groove was determined based on the tooth width (L) and number of grooves (n). As given in Table 1, for the baseline seal, the clearance size δ is set as δ/H ≈ 0.008 and the groove height h fixed at L/10. This is of a scale similar to the mean flow recirculation at the clearance gap entrance and ensures for inclined cases, the grooves do not pass through the tooth sides creating floating  geometry. For Seal B, the grooves were located at the start and end of the tooth and the groove width was set as (1/12)L. For Seal E the grooves were located at two positions along the tooth, the corresponding l was set as (1/6)L; for Seal F the teeth were located fully along the groove with l set as (1/10)L.

Mass flow results
To help better define the initial GA population, the mass flow results for the nine example seals are analysed through varying the imposed pressure ratios, to investigate the effect of various parameters on the seal performance. The mass flow coefficient, f , is selected as the objective function. Figure 7 compares the pressure Table 1. Geometrical parameters of the labyrinth seal. ratio against the mass flow coefficient for these nine example seals, at pressure ratios of 1.5, 1.75, and 2.04. For all cases, the mass flow increases almost linearly as the pressure ratio increases supporting validity of results. The calculated results also show a strong dependence of the objective function on the seal configuration. As can be seen from Figure 7, the mass flow coefficient of the seals with grooves is much smaller than that of the baseline, Seal A, suggesting that the leakage performance of the grooved seals is better. It is also noted that the more grooves, the better the leakage performance that can be achieved. The leakage mass flow coefficient of Seal F reduced by 16% compared to that of the baseline Seal A for a given pressure ratio. While seals D to F have better leakage performance, seals B and C fail to show a positive effect, with Seal C performing particularly poorly.
In order to highlight the difference in vortex patterns that different seal layouts can produce, Figure 8 presents the instantaneous vorticity magnitude contours for each of the example seals. The flow fields shown in Figure 8 for seals with grooves (D to F) are similar. For these geometries, the flow through the gap enters the grooves, which causes higher loss and a large flow blockage area. This complicated flow structure caused by the grooves results in a larger pressure drop. For a given pressure ratio, then, there is less leakage for the seals with grooves. For Seal F, with grooves, the flow passing through the clearance separates into three recirculating regions. In Seal E, two vortices rather than three were formed. As a result, the vorticity for Seal F is more intense. Consequently, it shows the best performance in reducing the leakage.
Fasheh [5] observe that an inclination angle would lead to lower leakage rates than a straight tooth profile. This finding is supported by the results of this study, as the straight seal is again found to have greater leakage than the upstream inclined seals. The leakage performance between the three inclined seals is quite similar, though as the inclination angle increases, there is slightly improved performance (Fig. 7). Seal I (with the greatest inclination angle) reduced the leakage by 3.2% relative to the baseline straight seal (Seal A) at a pressure ratio of PR = 2.04. As shown in Figure 8, Seal I produces a more distinct deflection of the leakage flow. After the fluid has flowed through the first tooth, deviation of fluid from the mean axial flow direction is greater. This leads to increased recirculation strength and kinetic energy dissipation in the main cavity.

Effects of rotational speed
Waschka et al. [41] show that at low speeds the impact of rotation on sealing performance could be considered negligible. This is further explored here, invoking the assumption that the rotor radius 20δ, so that the Cartesian coordinate system can be maintained. Figure 9 compares the mass flow results based on different rotational Reynolds numbers. Here, the rotational Reynolds number is defined as Re ω = ρω R R r 2 /µ, where ω R represents the rotational speed, and R r = 0.0509 represents the seal tip radius.
The rotor surface is set to be the rotating wall boundary. Simulations were conducted with the rotational Reynolds number in the range of 0-7000 with a constant sealing clearance. From the results, as shown in Figure 9, the leakage flow rate remains approximately constant. Thus to keep the number of simulations -and mesh requirements -manageable, the effect of rotational speed is not considered in this optimisation exercise. However, Lu et al. [15] show that the rotational speed and its physical impact is aligned with that of clearance. Hence, variations in clearance are, in a sense, a surrogate for rotational speed.

Effect of the inclination angle direction
The inclination angle of the tooth profile is allowed to be either positive or negative. Two example cases (Seal G, Seal J) were selected to investigate the effect of the inclination angle direction. They were generated from the baseline seal, Seal A, with Seal G and J being given a +5 • inclination angle and a −5 • inclination angle, respectively.
A comparison of the mass flow coefficient of these two geometries at different pressure ratios is shown in Figure 10. For reference, the baseline Seal A is also included. It is clear that seals with positive inclination angles have better sealing performance than the baseline seal. Moreover, for Seal G (+5 • ), the increase over the baseline effect was approximately double that for Seal J (−5 • ). Figures 11 and 12 present the vorticity magnitude contours of the two examples seals.
From the above test results and previous studies in the literature, five geometrical parameters are determined. They are: the inclination angle (θ); clearance size (δ); tooth width (L); number of grooves (n); and groove width (l). The parameters and their allowable ranges are given in Table 2. A schematic indicating the corresponding model and geometrical parameters is given in Figure 13. These geometrical parameters of the seal geometry are then used to generate hundreds of individuals as the input design variables.

Optimisation results
An initial generation is randomly generated. By iterating through the genetic optimisation process, an optimised design is obtained. By analysing the evolutionary process, it is possible to see how the optimisation advanced, as well as the optimal values found within each generation.

The progress of the evolution
The average performance of each generation is presented in Figure 14. As the population in the first generation are randomly generated, the corresponding performance is poor. Thus there is considerable scope for improvement in the subsequent generations. After the first generation, an average reduction of 8.5% of the mass flow rate was achieved when compared with the initial geometry configuration. The rate of the improvement begins to slow, as further improvements have a stronger dependence on the random mutations, and higher peaks are increasingly rare. For the second generation the average reduction was up to 4.0%. Figure 15 shows the development of the relationships between the mass flow coefficient and the blockage ratio (δ/L), where the symbol "+" represents the first generation, "•" shows the third generation and "♦" the fifth. The performance gains over these three generations can be seen from this figure, as well as the tendency the optimiser had to reduce the blockage ratio. These effects are especially clear between the first and the third generations.

Empirical relations
As well as producing the optimal design, the calculations carried out on the way can also give useful information. Here, the relationships between the sealing performance and various design parameters are presented using the populations from the final three generations. By this time, the gene pool has become largely saturated with "good" genes. The correlation lines on each scatter plot are determined using a Theil-Sen median slope estimator. This estimator can be computed efficiently, and is insensitive to statistical outliers. As defined by Theil [29], in each correlation line the slope is chosen as the median slope among all lines through pairs of points. Once the slope m has been determined, the y-intercept b is given by the median y-intercept of each of these gradients.
As expected, the correlations between the overall leakage performance and each geometrical parameter are fairly weak (as shown by the gradient and variance from Thiel-Sen predictor summarised in Tab. 3). This is because, when conducting the optimisation process using GA, there is no requirement to alter only one parameter at one time. In fact, the capability of the GA to select combinations of weakly dependent parameters is a particularly useful attribute. Figures 16-19 plot the mass flow coefficient against clearance size, inclination angle, tooth width and groove width, respectively. The Theil-Sen median slope, m, and variance from Thiel-Sen predictor, σ 2 m , for these parameters are shown in Table 3. The slope is plotted using thick dashed lines in Figures  16-19.
Of all the figures of merit explored, the clearance size (δ) is found to be the critical parameter with the strongest correlation (σ 2 m = 0.00128), as shown in Figure 16. This is not unexpected, as with larger clearance more fluid is less strongly influenced by the tooth. As a result, a larger portion of the kinetic energy is carried over to the next cavity without being dissipated and accordingly the mass flow coefficient increases. Figure 16 also shows the standard deviation of the points around Thiel-Sen median line (fine dashed lines). Unfortunately, the scope to change δ is limited, often due to running clearances. Hence other parameters become more influential. The second and third most important parameters are the tooth width (L) and the groove width (l), with σ 2 m = 0.00187 and 0.00228 respectively, as shown in Figures 17 and 18. It is interesting to note that as the generation increases, there is a tendency to move down the Thiel-Sen gradient. This is particularly apparent in Figure 17, where different generations can clearly be seen. The inclination angle (θ) shows relatively weak correlation with σ 2 m = 0.00344, as presented in Figure 19

Optimised design
The most successful seal tooth configuration from the final generation is obtained. Table 4 gives the geometrical parameters for the optimised configuration. Also shown are the geometrical parameters of the initial baseline configuration. The geometry of the optimised design can be seen from Figure 20.
The mass flow coefficient of this optimised labyrinth seal is found to be significantly better than that of the baseline seal. A significant reduction of 27.6% in the mass flow coefficient was observed. By simply qualitatively       comparing the configuration of the best performing labyrinth seal and that of the baseline seal, it is clear why the latter has a significantly higher mass flow coefficient. The best performing geometry consists of three grooves which are quite evenly spaced in the streamwise direction. These generate substantially (50%) higher turbulence relative to the baseline seal. This significantly increases loss and hence the effort required to force the fluid between the cavities. The high speed wall jet from the baseline Seal A is more likely to impinge directly onto the tooth of the next downstream cavity.

Flow fields
In order to better understand the capability of the optimised design to control flow leakage, the flow field is briefly analysed. The instantaneous vorticity magnitude is compared for the optimised seal and the baseline seal in Figure 20. As can be seen, at the exit of the seal a wall jet is generated. Fluid from the cavity is highly Figure 22. Comparison of Q-criterion iso-surface between the baseline Seal A and optimised design.
accelerated into the seal gap. On exit, it is decelerated rapidly, forming a turbulent shear layer next to the lower wall. Also, it can be seen that the flow penetration into the cavity is moderate and the flow decelerates in the corner regions of the cavities. The shear layer sizes of Seal A and the optimised seal are presented using the arrows in the averaged vorticity magnitude contours in Figure 21. The size of the shear layer (expressed as a half width) for the optimised seal is approximately 8.75% thicker than that for Seal A. This increased thickening is due to more mixing and hence loss, which gives the improved sealing.
To identify coherent structures more clearly, contours of the Q-criterion are plotted in Figure 22. The flow spreads after passing through the seal gap. After passing through the clearance, the main flow path comes to the next cavity. The flow then begins to accelerate again, with vortices becoming stretched near the upper part of the seal clearance, and large recirculations are generated with larger coherent structures. Small vortices are also generated in the boundary layer. Hence a mixed flow is formed together with the larger vortical structures, which further improves the sealing performance.

Conclusions
The present study investigated the effect of geometry on the leakage mass flow coefficient for around three hundred different labyrinth seals. The clearance size, inclination angle, number of grooves, tooth width, and groove width were specified as design variables. RANS calculations have previously been shown to be inadequate, whereas the flow structures can be reliably produced using LES.
Multiple geometries were generated using the GA, controlled by the selected geometrical parameters and manufacturing constraints. Geometries are generated in sets, representing a population for each evolutionary generation. The range of geometries is represented using an IBM on a background Cartesian grid. This avoids the need to laboriously mesh by hand, or to trust in automatic meshing software and manually deal with the inevitable failures. As a further benefit, the computational cost of each of the individual simulations could be reduced by the use of an efficient structured solver. IBM gives a much easier way to resolve the wide range of complex geometries with a simple mesh.
The best performing seal is that with the lowest mass flow coefficient, f , for the applied pressure gradient. In comparison to the baseline geometry, the design optimised by the GA showed a leakage reduction of 27.6%. The optimised design consisted of a seal inclined towards the flow direction with grooves at its tip. Along with this, among the five variable design parameters, the effect of the clearance size, δ, had the strongest correlation (σ 2 m = 0.00128) on the mass flow coefficient. The parameter with the second strongest correlation was the tooth width, L, with σ 2 m = 0.00187. The use of LES for large scale design exercises has only recently come within the reach of modern computational systems. With careful problem selection, such as the labyrinth seals shown here, the use of a GA based optimising heuristic can achieve significant design improvements in relatively short real-time scales even on a modest local cluster or group of PCs.