Numerical Analyses of Static Characteristics of Liquid Annular Seals Based on 2D LBM-LES Model

Static characteristics and leakage flow rates of liquid annular seals have great influences on the hydraulic efficiency of turbomachinery. In this paper, a two-dimensional (2D) mathematical model for predicting the leakage flow rates and static characteristics of liquid seal is established, based on the lattice Boltzmann method (LBM) combined with the D2G9 velocity model for incompressible fluid and large eddy simulation (LES) turbulence model, in which the transformation equation of reference pressure is developed with the Bernoulli equation. Moreover, the proposedmodel is validated by comparing with the experimental results, calculation results based on the finite volume method (FVM), and the results based on the empirical method of three seals under different operating conditions. )e comparisons show that the maximum deviation in leakage prediction of the calculating model based on 2D LBM is 4%, and this calculating model will effectively improve the leakage prediction accuracy of the seals compared with the FVM and theoretical method.


Introduction
In turbomachinery, there exist several sets of liquid annular seals, including neck-ring seals, interstage seals, and balancepiston seals. ese seals prevent leakage flow while increasing the volumetric efficiency of the machinery. At present, there are mainly three methods to investigate the static and dynamic characteristics of these seals, including empirical formulas, bulk-flow method, and computational fluid dynamics (CFD) method. However, these three methods are all based on the hypothesis of continuous media, which means that the methods are all developed from a point of macroscopic view. With the development of the lattice Boltzmann method (LBM), many researchers started to use this method to solve different fluid flow problems of compressible and incompressible fluids. Compared with the widely used finite volume method (FVM), LBM has the advantages of simple algorithm, accurate calculation, and strong adaptability of boundary conditions, which also has provided an effective method for calculating fluid flow from a point of mesoscopic view. And in recent years, many researchers have analyzed the static and dynamic characteristics of sliding bearings, which are structurally similar to liquid annular seals based on LBM. Kucinschi and Afjeh [1] analyzed the internal flow characteristics and the fluid-film lubrication details within model sliding bearings based on 2D LBM. Kim et al. [2] have used LBM to study the flow conditions in a nano-sized air bearing, and Ramirez et al. [3] used this method to analyze the switching flow characteristics between different channels in a micron-sized air bearing. e simulation results showed that LBM has high accuracy in simulating the fluid flow of the sliding bearing with mini clearance.
Although the structures of sliding bearings are similar to those of liquid annular seals, there are major differences between them. e internal flow within sliding bearings is laminar flow with small Reynolds number, while the internal flow of seals in the centrifugal pumps is highly turbulent with large Reynolds number. Recently, Posa and Lippolis [4], Li et al. [5], and Si et al. [6] analyzed the internal flow characteristics of centrifugal pumps considering liquid annular seals channels based on different turbulent models. Andres et al. [7], Saber and Abdou [8,9], and Nagai et al. [10] proposed the theoretical calculation method for static and dynamic characteristics of different kinds of labyrinth seals and investigated the effects of operating conditions on these characteristics. As the maximum Reynolds number that fit the single-relaxation-time LBM (SRT-LBM) is 10 4 and the maximum Reynolds number that fit the multiple relaxationtime LBM (MRT-LBM) [11] is 4 × 10 4 , analyses on the fluid flow with large Reynolds number based on LBM are often performed with the turbulent model, such as LBM-k-ε model established by Succi et al. [12], mixed length algebraic turbulent model established by Teixeira [13], and LBM-LES turbulent model established by Hou et al. [14]. However, the LBM-k-ε turbulent model shows poor performances on describing flow separation, reattachment, and recovery effects [15], and mixed length algebraic turbulent model requires a large number of empirical coefficients, which make the two models not so widely used.
In contrast, the LBM-LES turbulence model combines the advantages of LBM and LES, in which the turbulent viscosity coefficient is calculated based on the Smargorinsky eddy viscosity model of LES. Premnath et al. [16] investigated the turbulent kinetic energy for the separation and attachment phenomenon of complex turbulent flow within a back-step model based on the LBM-LES turbulence model. Schneider et al. [17] used the LBM-LES turbulence model to solve the inlet flow characteristics of a centrifugal pump firstly. Chang et al. [18] analyzed the heat transfer between square cavity flow and back-step flow with the LBM-LES turbulent model. Comparisons between the numerical and experimental results showed that the method could describe the turbulent heat transfer flow well. Hamane et al. [19] simulated the flow within a 2D cylindrical channel with a large Reynolds number based on the LBM-LES turbulent model. Eitel et al. [20] introduced hierarchically refined lattice technology to the LBM-LES turbulent model and contrasted simulations of the cylinder flow model with different Reynolds number. e refinement of lattice reduced the number of lattice and significantly improved the computational efficiency. Nadim et al. [21] analyzed the flow of airfoil using LBM which combined the Smagorinsky SGS model, while parallel calculations were also used in the paper.
With the development of LBM, static characteristics and leakage analyses of liquid annular seals from a mesoscopic point of view become possible. In this paper, a two-dimensional (2D) mathematical model for predicting the leakage flow rates and static characteristics of liquid seal is established, based on the lattice Boltzmann method (LBM) and large eddy simulation (LES) turbulence model. Because the internal fluid of the seals is incompressible, the D2G9 velocity model for incompressible fluid is applied. Nonequilibrium extrapolation scheme is used for the boundary conditions of the wall.
Moreover, on comparisons of the experimental results with the results of LBM, FVM and empirical formulas are conducted and analyzed to verify the effectiveness and accuracy of LBM in calculating the results of static characteristics and leakages of seals.

Simulation Method
As the internal flow of the seal is mainly affected by the pressure-induced force, the lattice Boltzmann equation (LBE) can be expressed as follows: In equation (1), M is the rate of change from initial state to final state of the particle distribution function f, which is called collision operator. e basic thought of the Boltzmann equation is not to determine the motion state of each molecule but to find the probability of each molecule in a certain state and to obtain the macroscopic parameters of the system through statistical methods [22].
e distribution function f acts to describe the number of molecules in a certain position in a certain range of speed at a particular time. e macroscopic values of velocity and pressure can be obtained by f, which is related to discrete velocity model, equilibrium distribution function, and its evolution equation [23].

Governing Equations.
e governing equation of the LBM calculating model can be written as equation (2) under the SRT-LBM model. It has second-order accuracy at low Mach number. In equation (2), c is collision frequency, τ is relaxation time, and the relationship between c and τ is c � 1/τ. τ is mainly related to the lattice kinematic viscosity ] 0 as shown in equation (3).
For the LBM calculating model, the distribution function evolution process is the most important component and the evolution process consists of two steps: the collision step that is described as equation (4) and the streaming step, which is described as equation (5).  [25] studied the compressible fluid flow which is around NACA0012 airfoil by the D1Q4 velocity model. Yang analyzed the 3D viscous compressible flows [26] by using the lattice Boltzmann flux solver. Meanwhile, for incompressible fluid flow calculation, Qian et al. [27] established the D2Q9 velocity model for weakly compressible. In this model, fluid flow is driven by density gradient, and the relationship between pressure P and ρ is shown as equation (6). On the basis of the D2Q9 velocity model, Guo et al. [28] developed the D2G9 velocity model for incompressible fluid flow. e fluid flow is driven by pressure gradient in D2G9. e equilibrium distribution function satisfies Σf eq α � const. e macroscopic equations of the D2G9 velocity model can be derived by Chapman-Enskog expansion [29], and the equations can be written as equations (7) and (8) [28].
zu zt e form of equations (7) and (8) is equivalent to the incompressible Navier-Stokes equations.
In this paper, the D2G9 velocity model is used to solve the fluid flow within the liquid annular seals. e velocity directions of the model are shown in Figure 1. e speed configurations of the D2G9 model are shown as equation (9).where e a is the discrete velocity, c is the lattice speed and it is defined as c � Δx/Δt, Δx is the lattice step, and Δt is the time increment.

Equilibrium Distribution Function.
e equilibrium distribution function of the D2G9 model is shown as equation (10). e equilibrium distribution itself is calculated from the macroscopic values which themselves are low order moments of f [28].

Turbulence Model.
According to the lattice Reynolds number equation shown as shown in equation (14), the lattice kinematic viscosity ] 0 will decrease with the increase in Reynolds number. Since the flow within seals is highly turbulent with a Reynolds number of more than 10 4 , the relaxation factor τ is close to 0.5, which will make the computation difficult to converge. To solve this problem, the LBM calculating model combined with the 2D LES Smagorinsky model is established specially for the static characteristics and leakage of liquid annular seals. As shown in equation (15), the total lattice viscosity ] t consists of two parts: the lattice kinematic viscosity ] 0 and the Smargorinsky eddy viscosity ] t,LES . e value of ] 0 can be described in equation (14), and the value of ] t,LES is described in equation (16), in which Cs is the Smagorinsky constant, Δ is the mean lattice spacing defined as Δ � (ΔxΔy) 1/2 , and s ij is the strain rate tensor that could be obtained by equation (17). In this paper, the value of Smagorinsky constant Cs is fixed at 0.1 [18]. Besides, near the walls, ] t,LES is damped, which means ] t,LES will decrease and approach to zero as the wall is encountered.

Boundary Conditions.
Boundary conditions have great influences on accuracy, stability, and computational efficiency of LBM simulation. In this paper, velocity inlet, open outlet, and nonequilibrium extrapolation scheme are used as the initial boundary conditions.

Velocity Boundary Condition.
e modified velocity boundaries based on bounce back part of the nonequilibrium extrapolation scheme are shown as equation (18). Extrapolation is used to calculate the outlet functions defined in equation (19) [30].

Wall Boundary Condition.
Nonequilibrium extrapolation scheme [30] is chosen as the wall boundary condition. Pressure is used as the dynamic variable for eliminating compressibility errors in the simulation. As shown in Figure 2, Node A, Node O, and Node C are located at the boundaries, and Node E, Node B, and Node D are located in the flow field. e distribution function of the boundary Node O can be calculated by the following equation: 2.4. Code Program. Calculation code for leakage flow rate of liquid annular seals is accomplished by using C++ on the basis of 2D LBM combined with the governing equations, turbulence model, and boundary conditions described above. e overall program block diagram is shown in Figure 3.

Simulation Parameters
In order to verify the method proposed in this paper, the predicted leakage flow rates of three different liquid annular seals under different operating conditions based on LBM are compared with FVM results, experimental results, and theoretical results of empirical formulas. Figure 4 illustrates the basic geometry and coordinate system used in this paper. e geometric parameters of the three seals are listed in Table 1, and the operating conditions are illustrated in Table 2.

Parameters of LBM.
e geometric parameters of model seals should be transformed to reference parameters with lattice unit [34], and the reference parameters must be the same value in the same model. e reference length H r , reference velocity u r , and reference density ρ r are defined as follows: In equation (21), where c s ' and c s are the sound speed in physical unit and lattice unit, respectively, in this study, the value of u r is 2678.44 m/s. Andρ′ and ρ are the density in physical unit and lattice unit, respectively; the value of the ρ r is 992 kg/m 3 . e reference length H r is defined by grid resolution; for example, if the lattice number of the model gets 9047 and the physical length H′ is 13.13 mm, the reference length H r should be 1.451 × 10 −6 m. Besides, the Set up parameter including lattice velocity U 0 , reynolds number Re, smagorinsky constant Cs, and lattice space ∆ Calculate the equilibrium distribution function of the pressure variable, as shown in Equation (8) Calculate turbulent viscosity vt, LES in large eddy simulation, Equation (14) Calculate boundary conditions, Equation (17) and (18) Calculate macro velocity and macro pressure, Equation (10) and (11) Collision, Equation (4) Streaming, Equation (5) No Yes

Error analysis
End and output data file Set up the initial macroscopic parameters including the grid dimensions Lx and Ly of the model, grid density ρ, inlet pressure P in , and outlet pressure P out  Shock and Vibration 5 transformation equation of lattice pressure P and reference pressure P r is shown as follows.
In equation (22), P c is the stagnation point pressure in the lattice unit, and it must take the same value for the stimulating of same model seal. e value of P c generally takes less than 1 but more than 0.09 for accelerating convergence.
Moreover, as the lattice Reynolds number Re lattice shown in equation (14), it is equal to the macroscopic Reynolds number Re, and it is possible to choose proper values for U lattice and ] 0 under the condition that Re lattice equals to Re [35]. Recently, ] 0 remains consistent in different stimulation cases for the same model. erefore, the value of ] 0 in this work remains the same for the same model seal when stimulations for different operation conditions are conducted. e simulation parameters of the model seals are listed in Table 3.
A grid independence study for each of the three models using seven different sizes of lattices has been carried out, as shown in Figure Figure 5 is calculated based on equation (23), where Q n is the numerical result and Q e is the experimental result:

Parameters of FVM.
e calculation results based on FVM which is realized by FLUENT software are compared with the results based on the proposed method in this study. Transient 2D model, 2D coupled solver, and implicit difference scheme are applied. e LES Smagorinsky model is introduced to describe the fluid stress term in the control equations and is discretized with the first-order upwind style. e convergence condition is that the residuals of all the equations are less than 1 × 10 −5 .

Parameters of Empirical Formulas.
As the empirical formulas shown in equation (24) are the most commonly used method for calculating the leakage flow rates of the seals in engineering, in this work, the theoretical results of empirical formulas also act as an important part of the comparisons. In the empirical formulas, C is the seal clearance, µ is the flow coefficient, L is the seal length, S is the open area of seal clearance, and Q is the leakage of seal. As shown in the above equations, µ is related to the input rounded coefficient, and the friction coefficient number is determined by the Reynolds number and the surface roughness. e coefficient λ generally takes the value between 0.04 and 0.06, and the coefficient η takes the value between 0.5 and 0.9.

Results and Discussion
e comparisons of predicted leakage flow rates of three different liquid annular seals under different operating conditions based on LBM, FVM results, empirical formulas, and experiences are conducted and discussed in this section. Table 4, the simulation results of Model 1 based on LBM, FVM, and empirical formulas are compared with the experimental results conducted by Childs [28]. Since the clearance of the annular seal is small, enlarged velocity contour diagrams at the seal exit  Figure 6. From Figure 6, it can be found that the fluid velocity increases as the pressure difference increases. In the seal clearance, the fastest fluid velocity is at the center of the clearance. e comparisons show that the empirical formulas are the most time efficient method but with the lowest prediction accuracy of leakage flow rates. e maximum error of the empirical formulas is nearly 15%. LBM and FVM have a good accuracy with a maximum error of less than 3.86% and 6.42%, respectively. For this model, the accuracy of LBM is higher than that of FVM. In contrast with the complex modeling and meshing process of FVM, for different model seals using LBM, only a few initial parameters such as geometric parameters, inlet velocity, and pressure difference need to be changed, which has made LBM more time efficient and more suitable for engineering application.

Results of Model 2.
Experimental results are conducted by Darden et al [32], and the calculation results which are calculated by the three methods for leakage flow rates of Model 2 are shown in Figure 7 and Table 5. It is observed that the leakage rates based on empirical formulas were under-predicted with an error of more than 10% compared to the experimental results. e prediction error of the leakage rate based on FVM is 2.85% under the pressure difference of 9.29 MPa but with error of more than 5% under the other two operating conditions. Among the three methods, LBM has obvious advantages in leakage prediction accuracy, especially under high pressure differences. e prediction error of LBM-LES decreases from 3.48% to 2.32% as the pressure difference increases, while the error of FVM shows an opposite variation.     [33], and the calculation results based on three methods for predicting leakage flow rates of Model 3 are shown in Figure 8 and

Conclusions
In this paper, a 2D LBM calculating model which consists of 2D LES turbulence model, D2G9 velocity model, and nonequilibrium extrapolation scheme is established for predicting of leakage flow rates and static characteristics of liquid annular seals. e leakage rates of three different plain annular seals under three different operating conditions are calculated, and the results of them are compared with the experimental results to prove the validity and accuracy of the proposed model. Comparisons show that all of the prediction errors of the 2D LBM calculating model are less than 5%, some of them even less than 2%, which indicates that the calculating model satisfies the requirements of engineering application better compared with the other commonly used methods. e accuracy of the calculating model for predicting the leakage flow rates of seals is mainly related to the clearance and pressure differences of the model seals. e accuracy decreases with the increasing clearance of the model seals which increases with the increase in pressure difference. e transformation equation of reference pressure Pr is developed based on the Bernoulli equation for incompressible fluids within annular seals. Reference velocity, reference density, and stagnation point pressure Pc in lattice unit are introduced in this transformation equation, and the optimal value range of Pc for leakage rates predictions of annular seals is proposed.

Nomenclature c s ′:
Physical velocity of sound ρ′: Physical density H′: Physical length P′: Physical pressure u: Physical velocity C: Clearance of seals L: Length of seals e: Eccentricity of the rotor D: Diameter of seals P in : Supply pressure of seal P out : Exit pressure of seal Ω: Rotor whirl velocity ω: Rotor angular velocity ] 0 : Lattice kinematic viscosity