Similar articles

Abstract: To improve lubrication effect and seal performance, complicated geometrical hydrodynamic grooves or patterns are often processed on end faces of liquid lubricated mechanical seals. These structures can lead to difficulties in precisely estimating the seal performance. In this study, an efficient adaptive finite element method (FEM) algorithm with mass conservation was presented, in which a streamline upwind/Petrov-Galerkin (SUPG) weighted residual FEM and a fast iteration algorithm were applied to solve the lubrication equations (Reynolds equation). A mesh adaptation technique was utilized to refine the computation domain based on a residual posterior error estimator. Validation, applicability, and efficiency were verified by comparison among different algorithms and by case studies on seals’ faces with different groove structures. The study investigated the influence of the order of shape function and the mesh number on the leakage balance. Mesh refinement occurred mainly in cavitation zones when cavitation happened, otherwise it occurred in regions with a high pressure gradient. Numerical experiments verified that the proposed algorithm is a fast, effective, and accurate method to simulate lubrication problems in the engineering field apart from end face seals.

液体润滑端面密封质量守恒自适应有限元算法

研究目的：为液体润滑端面密封空化问题的数值模拟提供一种高效精确的计算方法。
创新要点：1.采用流线迎风SUPG有限元法求解具有对流扩散特征的控制方程；2. 针对未知量的互补关系提出了一种高效的数值迭代技术；3.提供了针对润滑液膜空化问题的网格自适应技术。
重要结论：1.相比文献算法，本文算法具有计算速度快，计算精度高的优点；2.空化区域和大压力梯度区域的计算网格得到加密。

关键词：空化；有限元法；端面密封；流线迎风伽辽金法SUPG法；质量守恒

1.  Introduction

 Liquid lubricated end face seals are mainly used to prevent a high pressure medium from escaping out of the machine. It includes a couple of stators and rotors respectively sitting on the house and axial of the rotary machinery. Between end faces of the couples, there exists liquid lubrication film to avoid the direct contact and severe wear. The end faces are often etched or machined with complicated geometrical hydrodynamic grooves to improve the lubrication state, decrease temperature rise and minimize leakage. These groove structures (Fig. 1) include spiral groove, circular groove, rectangle groove, and even surface textures. According to mechanical seals, the flow rate or leakage rate is a more important performance parameter than journal bearings or thrust bearings. However, it is difficult to precisely evaluate the leakage rate of mechanical seals because of the existence of complex hydrodynamic grooves and the cavitation effect resulting from Poiseuille flow and Couette flow of the lubrication film.




 According to the classical lubrication theory, the leakage rate is proportional to the film pressure gradient distribution along the normal direction of the leak boundary. Therefore, it is important to derive the exact film pressure distribution and pressure gradient from the lubrication governing equation—Reynolds equation. There are three major factors leading to difficulties in achieving a solution, i.e., (1) cavitation effect with mass conservation which implies prior uncertainties on the division of cavitation and full film zones, (2) convention dominated problems from elliptic-hyperbolic type equations which leads to unstable solutions under the high Peclet number (Zienkiewicz and Taylor, 2000), and (3) complicated hydrodynamic groove structures where the film thickness is discontinuous and high pressure differences between faces which cause more complex pressure distribution than bearings.

 Nowadays, plenty of mathematical models and numerical algorithms have been proposed since Jakobsson and Floberg (1957) and Olsson (1965) presented the JFO cavitation theory. It is now the only boundary condition preserving mass conservation of fluid films. The finite difference method (FDM) (Payvar and Salant, 1992), finite volume method (FVM), finite element method (FEM), boundary element method (BEM) (Yu and Keith, 1995), and space-time conservation element and solution element method (CE/SE) (Cioc and Keith, 2003; Cioc et al., 2003) were used to solve the linear or nonlinear equation group deduced from the Reynolds equation. It is widely accepted that Elrod and Adams (1975) and Elrod (1981) initialized the successful numerical algorithm. A universal equation was proposed to describe the lubrication film by assuming that the film was compressible in non-cavitation zone and liquid/gas mixture exists in the cavitation zone, which was called the α-g model. The FDM was chosen to solve the governing equation by central and upwind difference styles (Vijayaraghavan and Keith, 1990a; 1990c). To improve the computational efficiency, Brewe (1986) introduced a multi-grid technology into the Elrod algorithm and analyzed the performance of the bearings. Pavyer and Salant (1992) presented an FVM algorithm where convection items of the equation were handled by a back difference method and an alternating direction iterative (ADI) technique was applied to speed the computational convergence. Because the above α-g model gives an unrealistic fluid film fraction (larger than 1) in the full film region for the numerical reasons from the choice of proper bulk modulus, a p-θ model is frequently used in various studies, where the JFO cavitation boundary is converted into a complementarity condition. According to the model, Ausas et al. (2007; 2009) presented an effective relaxation scheme to update the solution of satisfying complementarity conditions. Due to the capability of FEM to create a complex boundary, it was first introduced into the prediction of cavitation regions with the JFO theory by Kumar and Booker (1991; 1994) and Boedo and Booker (1995), who used Murty (1988)’s algorithm to distinguish the film rupture and reformation boundaries. Bonneau et al. (1995) and Optasanu and Bonneau (2000) used this technique to analyze dynamically loaded bearings. Bayada et al. (1990; 1998; 2006) and Durany et al. (1997; 2006) treated the cavitation as a free boundary problem and gave a series of detailed verifications. They used characteristic-Galerkin FEM and a duality iteration method to stabilize the solution. Giacopini et al. (2010) applied a linear complementary problem (LCP) method to deal with the complementarity condition and validated the algorithm on the textured bearings and squeeze film dampers in a 1D model. Hajjam and Bonneau (2004; 2007) and Fatu et al. (2005) proposed a modified p-θ FEM algorithm where a high order shape function in a non-cavitation zone and a low order shape function in a cavitation zone were used to simulate the lip seal and dynamically loaded journal bearings. Schweizer (2009) rewrote the Reynolds equation and JFO cavitation boundary to an arbitrary Lagrangian-Eulerian (ALE) problem and the FEM formulation was presented to solve some numerical examples. Due to the prior uncertainty and irregularity of film rupture and reformation boundaries, an adaptive mesh-refinement technique is frequently used to improve the computational accuracy. Vijayaraghavan and Keith (1990b) applied a grid transformation and adaption technique in the cavitated journal bearings. Case studies verified a better accuracy and lower time-consumption. Moreover, Nilsson and Hansbo (2007) applied an adaptive FEM for hydrodynamic lubrication. A posteriori error estimator and adaptive algorithms were discussed in detail with a pity of a non-mass-conserving formulation used.

 In this paper, an effective mass-conserving finite element algorithm is provided for predicting the seal performance of liquid lubricated end faces etched with complex hydrodynamic grooves. A stream-upwind/Petrov-Galerkin (SUPG) weighted residual method was implemented to the compressible Reynolds equations and an FEM formulation was deduced. A fast and stable block iteration algorithm in (Shi and Salant, 1999; Shi and Paranjpe, 2002) was applied to deal with the JFO boundary condition. Film pressure and density distribution were solved simultaneously conforming to the complementarity condition. Then a classical residual posterior error was estimated, according to which the computation domain meshes were refined and the FEM formation was solved iteratively. Case studies were conducted to validate the effectiveness and efficiency of the algorithm for some conventional liquid face seals. Mass conservation characteristics and seal performance were compared by different orders of FEM shape functions and mesh adaptation techniques.

2.  Mass-conserving finite element model and mesh adaptation algorithm

2.1.  Governing equations

 According to the lubrication film between two end faces, by applying the conventional hypothesis of lubrication theory, where the flow is assumed to be laminar and isothermal and the inertial force is neglected, the film pressure p can be solved from the compressible Reynolds equation in the orthogonal coordinates form (1) , where h is the seal clearance or film thickness between two faces, and μ represents the dynamic viscosity of lubricants. U=(U_x , U_y ) describes the velocity of the rotor face in the x-y coordinates. When cavitation happened, fluid film can be subdivided into two regions, full film zone and cavitation zone. The common boundaries of the two zones are called film rupture boundary and reformation boundary according to whether fluid starts to cavitate or reforms. In the full film zone, the lubricants are viewed as incompressibility, so the density ρ is constant (equal to the density of fluid film ρ _L) and the pressure p is larger than the cavitation pressure p _c (in this study, we suppose p _c=0). In the cavitation zone, a liquid/gas mixture exists, so ρ<ρ _L and p=p _c. Then Eq. (1) can be rewritten as follows, which is an elliptic equation in the full film zone and a hyperbolic equation in the cavitation zone: (2) with a boundary condition on film rupture and reformation boundaries in the form of (3) where U_n is the velocity along the normal direction n.

 Boundary condition Eq. (3) preserves the mass conservation across the film rupture and reformation boundaries. If the film reformation boundary employs the same formula as the film rupture boundary, then the boundary is called a “Reynolds cavitation boundary”. For computational convenience, a new variable θ=ρ/ρ _L, named as the film fraction, is introduced, then Eq. (2) can be written in the following form in both the full film zone and cavitation zone (4) , with complementary condition (CC) in the form (5)

 Eq. (4) with Eq. (5) is called the p-θ cavitation model, which is an elliptic-hyperbolic equation similar to the transonic potential flow (Vijayaraghavan et al., 1991). This equation is also a convection-diffusion-reaction equation. The p-θ model can be converted to an LCP and solved by Murty (1988)’s algorithm or some suitable pivoting algorithm like in (Giacopini et al., 2010).

 Once the film pressure p is solved, the seal performance parameter such as the open force F _o and leakage rate Q can be easily calculated by (6) , where Ω is the computational domain, that is the fluid film region, and l is the leakage boundary of the end faces. The subscriptions i and o represent the inner boundary and outer boundary of the seal faces.

2.2.  Finite element variational formulation

 To compute Eq. (4) with Eq. (5) by the FEM, triangular mesh is used to discrete the whole domain. The groove boundary conforms to the mesh edge, by which the different film thickness can be strictly endowed to the groove zones and non-groove zones. A Galerkin method is chosen to handle the elliptic part of Eq. (4) and the SUPG stabilized method is performed to manage the hyperbolic part by adding a numerical diffusion term (artificial viscosity) along the streamline direction (Evans, 1997; Zienkiewicz and Taylor, 2000). So the following integral weak form of Eq. (4) can be obtained by the way of integration by parts (Hajjam and Bonneau, 2004; 2007), in addition of Eq. (5) on the transition boundaries between the full film and cavitation zone (7) , where w is a differentiable function defined in domain Ω, and τ ^SUPG is the stabilization parameter which has the following form (Evans, 1997; Zienkiewicz and Taylor, 2000), , where h _T is the characteristic length of each element in the direction of velocity U (U_x , U_y ). That means τ ^SUPG varies on each element.

 Let w=N, p=p _i N _i, and θ=θ _i N _i, where N is linear or quadratic depending on the order of the shape function. The weak variational form of Eq. (7) is written as (8) , where K is the stiffness matrix which can be determined by (9)

 To solve Eq. (8) with Eq. (5), it is convenient to introduce another universal variable Φ here. The form is (10) , where F is the zone indicator, which is 1 in the full film zone and 0 in the cavitation zone. That is to say, Φ is equal to pressure p in the full film zone and to film fraction θ in the cavitation zone. Substituting Eq. (10) into Eq. (8), we can obtain: (11) , where (12) , where I is an identity matrix, and C is a matrix where all entities are zeros except the diagonal entities filled with zone indicator F_i .

 An iteration procedure is performed to obtain the pressure and density distribution. Before iteration, matrixes K ^p and K ^θ are computed with the Gauss integration method from Eq. (9) by assembling all the element stiffness matrixes. Then the lubrication region is initially assumed to be in full film zone, that means F=1. Matrix A is calculated according to Eq. (12). Eq. (11) is solved by the generalized minimal residual method (GMRES) iteration technique and the derived variable Φ is checked whether Eq. (5) is satisfied. If F_j =1 and Φ_j <0, then let F_j =0; if F_j =0 and Φ_j >0, then let F_j =1. If Φ conforms to F, the convergent solution is gained; if not, the iteration procedure should be continued. The following case studies shows that 6–12 iterations should be conducted before convergent solutions are obtained. Note that different from Hajjam (2007)’s algorithm, the global stiffness matrix K of the present algorithm is computed only once during computation, so the algorithm is very fast and effective.

2.3.  Adaptive mesh-refinement technique

 Because the cavitation zone and non-cavitation zone are not distinguished until a stable solution is computed, the initial little uniform FEM meshes results in large computational errors. A remedy is the adaptive mesh-refinement technique, where the meshes are refined according to the residual posteriori error of the unknown variable Φ. The following notation similar to that in (Hajjam and Bonneau, 2004; 2007) is defined to avoid the jump of Φ across the film rupture and reformation boundaries: (13)

 Then the local residual a posteriori error η can be evaluated for every element by (14) , where , where T and L respectively are the set of elements and edges, δ_T is the longest radius of T, δ_l is the length of element edge l, and [∂D/∂n] is the jump of function ∂D/∂n across the edge. According to Eq. (14), both the posteriori error of pressure in the full film zone and film fraction in the cavitation zone can be simultaneously solved. Meanwhile, the film rupture and reformation boundary position error can be estimated although unknown variable Φis discontinuous.

 Once η is computed, we can use it to adapt and update the meshes. The new mesh size δ_{k +1} is given by (15) , where δ is used to control the triangular meshes sizes. is the mean of η in all element meshes. c ₁ and c ₂ are used to control the size of refinement or coarseness region. Here, c ₁=3 and c ₂=1/3 are chosen to avoid a generation of too thin or too dense meshes. From Eq. (15) we can deduce that the larger the local residual a posteriori error η is,i.e., the smaller the mesh size δ is, i.e., the meshes where the local error is large are refined.

3.  Validations and case studies

3.1.  Validations

 To validate the present algorithm, a comparison was conducted according to a classic journal bearing in (Elrod, 1981) and a seal face with one dimple in (Qiu and Khonsari, 2009), which are shown in Fig. 2 and Fig. 3, respectively. We can see that the present algorithm gives similar pressure distributions for the journal bearing and face seal, indicating that the present algorithm is valid.







 To further verify the efficiency of the algorithm, a comparison was made among different models and algorithms for a seal face etched with only one dimple (its radius is 1.5 mm). These algorithms are the present algorithm (JFO-FEM), finite difference method with JFO cavitation theory (JFO-FDM) from (Ausas et al., 2007; 2009) and FEM with Reynolds cavitation boundary condition (REYNOLDS-FEM). According to REYNOLDS-FEM, Uzawa’s algorithm (Ito and Kunisch, 2003) was applied to assure pressure satisfying the Reynolds boundary condition. Table 1 lists the value of the face geometrical parameters and operating parameters. Fig. 4 shows the pressure and film fraction distribution obtained from the above algorithms, respectively. We can see that JFO-FEM gives a similar pressure and density distribution to those from JFO-FDM. The cavitation zone (Figs. 4d and 4e) and maximum film pressure (Figs. 4a and 4b) show an excellent agreement between the two algorithms. Fig. 4c shows the pressure distribution from REYNOLDS-FEM. The maximum film pressure is larger than that in Figs. 4a and 4b from the JFO models. Furthermore, a larger cavitation zone in Fig. 4c is also found. All these results illustrate that the Reynolds cavitation boundary condition cannot provide accurate seal performance predictions.







 Fig. 4f gives a comparison among the above three algorithms on a cost time for computation t and the seal performance parameter, including open force F _o and flow rates, Q _i and Q _o. Flow rate and open force show an excellent agreement between JFO-FEM and JFO-FDM. Due to JFO model’s mass conservation, flow rate (leakage rate) Q _i is approximately equal to Q _o. This demonstrates the accuracy of the present model and algorithm. Although more meshes were used for the current algorithm (30 044 triangular meshes) than JFO-FDM (251×91 structured meshes), it performed in a very high efficiency which is almost 180 times higher than the later. The REYNOLDS-FEM algorithm also had a very high computation efficiency (the computation time is 6.92 s according to 30 044 triangular meshes) as JFO-FEM does, and it overestimates the open force F _o by 13.8% and yields an unequal leakage rate between Q _i and Q _o with a bias of 11.6%. All numerical results in Fig. 4 show that the current algorithm has a high efficiency and precision.

 To examine the applicability of the current algorithm, in the following sections, some case studies were conducted on the seal performance prediction of the conventional mechanical seals such as laser textured surface (LST), spiral groove surface (SGS), laserface surface (LFS), etc. These end face structures are shown in Fig. 1. The seal face geometrical parameters and operating conditions listed in Table 1 are used to evaluate the film pressure distribution and seal performance.

3.2.  Case study 1: LST mechanical seal

 Etsion and Michael (1994), Etsion and Burstein (1996), and Etsion et al. (1999) invented LST mechanical seals (LST-MS), of which the surface is etched with regularly distributed micro-dimples. The dimples act as numerous micro-bearings between seal faces by which the seal performance and resistance to wear are greatly improved. In (Etsion and Michael, 1994)’s mathematical model, the Reynolds cavitation boundary condition is always applied, which overestimates the open force and is a non-mass -conservative as mentioned above. Here, our JFO-FEM algorithm was used to recalculate the seal performance and pressure distribution.

 Fig. 5 shows the results from solving for LST-MS. Fig. 5a is the final FEM meshes after a few adaptations. We can see that the mesh refinement zones are clearly recognized in the cavitation region, which means larger local errors occur here. From Eq. (2) we know that governing equation describing the non-cavitation region is elliptic and that describing the cavitation region is hyperbolic. Meanwhile, the same order shape function of FEM elements are used in the lubrication region, so computation accuracy in the cavitation region is lower than in the full film region. Fig. 5b illustrates the film pressure distribution. It can be seen that the local maximum pressure exists at the convergent step of every dimple, while the cavitation occurs at the divergent step. Moreover, the global maximum pressure appears at the high pressure side. Fig. 5c shows the density distribution of LST-MS, where θ=1 means the zone belongs to the full film zone and θ<1 determines the cavitation zones. The cavitation zone area of every dimple increases along the radial direction from the outer side to inner side. This indicates that low pressure can easily lead to occurrence of a cavitation. To further improve the dynamic effect of micro-dimples, a partially distributed LST-MS (PLST-MS) is introduced by Etsion and Michael (1994). Figs. 5d–5f show the computation results of PLST-MS with 4×4 micro-dimples. Mesh refinement zones are also located in the cavitation region just as in LST-MS. The apparent ‘accumulation effect’ can be seen in Fig. 5e. The local maximum pressure of every dimple increases from the inner side to outer side and along the film flow circumference direction. Local minimum density (Fig. 5f) shows the same rules as the maximum pressure does. Comparisons were made on the seal performance with different mesh numbers, orders of FEM shape function N, and mesh adaptation. Simulation results of LST and PLST are listed in Table 2. We can see that the mesh number, order of N and mesh adaptation have little influence on the open force F _o, while the leakage rate is greatly influenced by these factors. The higher order shape function and larger mesh number give a more accurate evaluation of the balance between Q _i and Q _o, defined as 2|(Q _o−Q _i)/(Q _o−Q _i)|. JFO-FEM with a mesh adaptation technique with the first order of N gives about the same accuracy of the leakage rate as JFO-FEM without an adaptation technique with the second order of N. These results demonstrate that the proposed mesh adaptation JFO-FEM algorithm has high applicability and accuracy for analyzing seal performance and optimizing micro-dimples’ distribution of LST-MS.

3.3.  Case study 2: spiral groove end face seals (SGS)

 In this section, another case study was conducted according to spiral groove end face seals (SGS) where spiral grooves were etched to improve their hydrodynamic effect. Fig. 6 illustrates the computation region meshes with and without mesh adaptation and the final pressure distribution. We can see that the mesh refinement mainly happens at the steps and sharp corners of the spiral grooves where the high pressure gradient occurs. The hydrodynamic effect of spiral grooves leads to a high pressure distribution at the sharp corner and a low one at the entrance of the grooves. Under the current geometric structure and operating conditions used, no cavitation zone is predicted. Table 3 lists the influence of the mesh number, order of shape function N, and mesh adaptation on the computation accuracy of the SGS’ seal performance. Same as in LST-MS, a large mesh number, high order of N, and mesh adaptation technique give rise to high accuracy on the balance between Q _i and Q _o, while having little impact on the open force F _o. Compared with LST, SGS shows a more excellent hydrodynamic and severe leakage rate. This case study demonstrates that the current JFO-FEM algorithm and mesh adaptation technique are applicable and effective in precisely estimating seal performance of end face seals with no cavitation occurrence.

3.4.  Case study 3: laserface end face seals (LFSs)

 To further verify the high efficiency of the current algorithm, a more complex hydrodynamic combination grooves end face seals called laserface seals (LFSs) were investigated, including rectangular grooves and crescent-shaped grooves spaced alternately with each other (Fig. 1). The rectangular grooves, called inlet grooves, introduce the fluid into the seal face clearance. Crescent-shaped grooves, called return grooves, force the fluid to the outer periphery of the seal. The fluid film thickness is described by where a and b are parameters determining the shape of the return groove, l_y denotes the position of the groove and c represents the return groove deepness.

 Fig. 7 shows film pressure and density distribution of the periodic face structure obtained from JFO-FEM with mesh adaptation and JFO-FDM. Due to the hydrodynamic effect of the grooves, two pressure peaks are shown respectively in the inlet groove and return groove. At the same time, two cavitation zones are shown. This type of pressure profile shows the existence of high pressure gradients on the outer side of the seal face, which reduces the accuracy of the leakage rate Q _o. Thus, fine meshes or mesh adaptation techniques are recommended to be utilized here. Similar pressure and density distribution are shown in Figs. 7a and 7b except for existence of local spikes around the return groove cavitation zone in Fig. 7b and extra biased predicting cavitation zone in Fig. 7d, which are the results from JFO-FDM. These errors are primarily induced by the inconsistencies between the structured mesh of JFO-FDM and the complicated groove boundaries. To remove these mesh errors, unstructured meshes should be generated, whose edge is strictly consistent with the groove boundaries just as those used in the current JFO-FEM algorithm.




 Table 4 shows the influence of different algorithms, mesh numbers and order of N on the predicted seal performance. JFO-FEM with mesh adaptation techniques provides a most accurate leakage rate prediction. Numerical experiments show that JFO-FDM provides a high open force and high time consumption for computation. This case study demonstrates that the current algorithm is greatly free in dealing with end face seals with complicated groove patterns.

4.  Conclusions

 In this paper, an effective adaptive FEM algorithm with mass conservation is presented to analyze the seal performance of end face liquid seals with complicated hydrodynamic groove patterns. The classical p-θ cavitation model was used to describe the lubrication problem between end faces. The SUPG stabilization technique was brought into the FEM procedure to deal with the elliptic-hyperbolic equation. An effective numerical iterative algorithm was applied to handle the complementary relationship between pressure and density. The mesh adaptation technique was introduced into the algorithm to precisely predict the seal performance based on the residual posterior error estimation. Numerical result comparisons were conducted among the present algorithm, JFO-FDM, and REYNOLDS-FEM to analyze seal performance of end face seals with different groove patterns. Numerical experiments demonstrate that the adaptive technology gives accurate results with low order shape function. The SUPG technology stabilizes the numerical results and provides the smooth pressure and film fraction distribution. The refinement of mesh mainly occurs in the cavitation zone and large pressure gradient zones. The proposed algorithm is especially capable to simulate lubrication problems with complicated surface patterns.

---

* Project supported by the National Natural Science Foundation of China (Nos. 51005209 and 51375449)