1. Introduction
The task of sealing technology is to prevent the undesired leakage of liquids or gases from or the ingress of dirt into machines and technical systems or to reduce such leakage to a minimum. According to the relative movement between the components to be sealed, seals are divided into two groups: “static seals” and “dynamic seals”. Static seals include, for example, the frequently used O-ring, while dynamic seals include, for example, mechanical seals, piston rings, labyrinth seals, and radial shaft seals. While static sealing is mastered to a large extent, the sealing of dynamic contact surfaces continues to pose challenges. For sealing dynamic shaft interfaces without overpressure in machine housings (e.g., gearboxes), radial shaft seals are frequently used. Therefore, this type of seal is used millions of times in technical systems and is one of the seals dominating the market in terms of quantity and sales [
1].
In order to comply with stricter legal CO
2 emission limits, electrified powertrain concepts are increasingly being realized. The electric engines used in these systems are operated at very high speeds (20,000–50,000 min
−1 [
2]) to increase their efficiency and economy. However, these high speeds introduce new challenges in ensuring proper sealing at shaft interfaces.
With increasing circumferential speed, radial shaft seals cause significant friction losses, which on the one hand lead to a non-negligible reduction in efficiency. On the other hand, the friction losses result in high local heat generation. Due to the small contact area between the sealing lip and the shaft (contact width 100–300 µm) [
1,
3], the elastomer material is subjected to high thermal loads, which can lead to material changes, increased wear and consequently to significant leakage. During operation, the relative movement of the rough surfaces generates a hydrodynamic pressure in the sealing contact, which leads to a (partial) separation of the surfaces. The seals currently used develop a lubricant film thickness of an average 1 micrometer. The friction is dominated by the high shear of the oil [
1]. Therefore, friction reduction offers a chance to extend the performance limits with regard to the current speed limitations and increase efficiency within existing applications.
Approaches for friction reduction are based on the manipulation of the tribological system “seal-oil-shaft”. One of the approaches aims at increasing the hydrodynamic pressure build-up in the sealing contact to increase the lubricant film thickness. In this case, a sufficient backflow of the fluid to be sealed must be provided to prevent the inherent increase in leakage. Since the lubricant is usually defined by the application, the increase in hydrodynamic pressure build-up must be induced by manipulating the shaft or the seal surface. For this purpose, microscopic deterministic or stochastic structures are applied to the surfaces.
The beneficial effect of a microstructured sliding surface on the macroscopic frictional behavior of tribologically stressed surfaces has long been experimentally proven in many applications. For axial plain bearings [
4] and piston rings [
5,
6] as well as for mechanical seals [
7,
8] and metal to metal face seals [
9], larger lubricant film thicknesses and thus reduced friction are achieved by a specific microstructuring of the sliding surfaces. Few studies have been conducted to examine the influence of microstructuring on elastomer lip seals. In [
10,
11], a structure is developed and investigated to provide a back-pumping capability for a PTFE lip seal. For this purpose, a deepening structure was applied to the side oriented opposite to the shaft. This structure affects the stiffness of the seal and deforms during assembly in such a way that an inhomogeneous contact pressure results, which deflects the fluid towards the fluid side and thus creates a pumping effect.
A methodology for friction reduction of radial shaft seals is considered in this paper, which was presented by the authors in [
12,
13] and is based on microstructuring and surface treatment of the elastomer sealing sliding surface.
Figure 1 illustrates this methodology using the Stribeck curves of a conventional (grey line) and a friction-optimized (black line) radial shaft seal.
Deterministic microstructures are applied to the sealing sliding surface (a) to increase the hydrodynamic pressure build-up and thus the film thickness. As a result, a load-bearing lubricant film is generated at lower sliding speeds. This leads to a shift in the friction regime from solid to mixed friction or from mixed to fluid friction, which leads to an overall friction reduction. However, since even slight wear leads to an impairment of the structural function due to the microscopic size of the structure, it is essential to protect the elastomer surface against wear. Therefore, a plasmapolymeric coating is applied to the structured sealing sliding surface (b). This coating additionally reduces friction in a solid friction regime.
The design and optimization of a low-friction radial shaft seal poses a major challenge when it comes to minimizing friction while ensuring no-leakage. In the scope of this paper, a simulation-based method is employed, since experimental optimization is linked to high costs and time. Within this method, the focus is put on the fluid mechanical optimization of the microstructure. It is shown that a computation of the fluid flow should be based on the complete Navier–Stokes equation and that cavitation phenomena have to be considered. Furthermore, flow simulations under varying structural geometry are performed and the results are used for structural optimization.
2. Simulation Setup
The optimization of the microstructure is carried out by investigating its influence on the local fluid flow using Computational Fluid Dynamics (CFD) simulations. For this purpose, a simulation model of a microstructure is derived from the radial shaft seal system, which is described in
Section 2.1. The underlying simulation method is then presented in
Section 2.2. A brief description of the mathematical principles required for the calculation of the fluid flow is given in
Section 2.3. Finally, the derived simulation model is shown in detail in
Section 2.4. In order to draw quantitative conclusions on the basis of this model, an estimation of the numerical error is given in
Section 2.5.
2.1. Radial Shaft Seal
The design of radial shaft seals is specified in ISO 6194 [
14], among others.
Figure 2 shows a schematic assembly situation of a type A radial shaft seal. While the housing is at rest, the shaft performs a rotary motion with constant direction of rotation. The primary function is the dynamic sealing between the rotating shaft and the static sealing lip. The secondary function is the static sealing between the housing and the outer diameter of the seal.
Static sealing is provided by fitting the seal into the housing, which is supported by a metal case. The outer diameter of the seal can be elastomer-coated (type A) or non-elastomer-coated (type B). The dynamic sealing is provided by the sealing edge, which is formed by two conical surfaces. During assembly, a radial force is caused by the difference in diameter between the shaft d and the seal, whereby the diameter of the seal is smaller than the diameter of the shaft. This difference ranges between 0.01–0.06 · d [
15]. A garter spring is located on the sealing membrane above the sealing edge, which contributes to the radial force. The result is a circumference-related radial force of about 0.1–0.15 N/mm. The local contact pressure distribution is determined on the one hand by the distance between the center planes of the garter spring and the sealing edge and on the other hand by the front side lip angle α and the back side lip angle β. The distance between the center planes lies in the range 0.4–0.7 mm. The front side lip angle α is within a range of 40–60° and is therefore always greater than the back side lip angle β, which is within a range of 20–35° [
1,
15]. This results in an asymmetrical pressure distribution shown in
Figure 2, which causes the backflow effect of a radial shaft seal [
1]. The mean surface pressure p
m is given on the basis of the circumference-related radial force and the width of the sealing edge, which is in the range of 100–300 µm, by a typical value of p
m = 10,000 mBar (p
m = 1 N/mm
2) [
3].
2.2. Methodology
The design and optimization of the low-friction radial shaft seal is based on a two-part simulation method, which is shown schematically in
Figure 3.
In order to consider the macroscopic system behavior, the sealing system (shaft, seal, and housing) is modelled as a flexible multibody system (FMBS) taking into account micro-elastohydrodynamic lubrication (EHL). In this case, the local lubricant film thickness h is determined considering the kinematics of the system and the deformation of the elastomer material. The material deformation is caused by the local pressure p in the sealing gap, which can be divided into a fluid pressure p
f and a solid contact pressure p
c. For operating conditions at low circumferential speed or low lubricant viscosity, only a slight fluid pressure results, causing the seal to operate in mixed friction regime [
1,
16]. Therefore, the model takes into account both fluid pressure p
f (considering surface roughness according to Patir–Cheng [
17]) and solid contact pressure p
c (microcontact model according to Greenwood–Williamson [
18]). The detailed design and structural optimization is based on separate Computational Fluid Dynamic (CFD) simulations, on which this article focuses. The optimized structure is in turn integrated into the FMBS model.
2.3. Governing equations
To analyze the flow behavior of a fluid, e.g., in a microstructure, a mathematical description of the general behavior of a fluid is required. This description is based on the conservation laws for mass, momentum and energy. Assuming an isothermal flow, the conservation law for energy can be neglected. The two remaining conservation laws for the mass (Equation (1)) and the momentum (Equation (2)) form the system of equations for calculating the flow variables [
19].
In these equations, t denotes the time, ρ the density, U the velocity, p the pressure, and τ the stress tensor.
In the case of fluid flow in microstructures, the pressure can drop below the saturation pressure of the lubricant in regions with an increasing film thickness. In such a case, the lubricant liquid changes into the gas phase (vapor cavitation). To account for this two-phase flow (liquid, vapor), the transient, compressible solver cavitatingFoam of the open source software package openFOAM [
20] is used. The underlying cavitation model tries to maintain a constant pressure in cavitating regions (isobaric cavitation model). As pressure drops below saturation pressure, liquid is converted into vapor. If pressure rises again, vapor is converted into liquid [
21]. The cavitatingFoam solver is based on the homogeneous equilibrium model (HEM). The HEM assumes that the liquid–vapour mixture is in a mechanical and thermodynamic equilibrium. Thus, the mixture can be described as a single fluid whose properties are determined based on the vapor fraction γ. The vapor fraction can be calculated from the density of the mixture ρ, the densities of the liquid ρ
l, and the vapor ρ
v at saturation pressure p
sat (Equation (3)). If only liquid is present, the value is γ = 0, if only vapor is present, the value is γ = 1.
In addition, it is assumed that the densities of the liquid and the vapor depend solely on pressure, such that both state variables can be linked by means of a barotropic equation of state (Equation (4)), where ψ denotes the compressibility of the mixture and ρ
l,0 denotes the density of the liquid at zero pressure.
Different models are available for calculating the compressibility of the mixture ψ. In this paper, the Chung model [
22] is used due to its more stable behavior compared to the linear model [
21]. The viscosity of the mixture μ is obtained by linear interpolation from the viscosity of the liquid μL and the vapor μ
v using the vapor fraction γ. The resulting system of equations for calculating the flow variables is solved numerically using the Pressure-Implicit with Splitting of Operators (PISO) algorithm [
23].
2.4. Model
The calculation of the fluid flow within a microstructure is based on the conservation laws shown above. This requires a model of the region to be investigated and the definition of boundary conditions. The aim of the calculations is to investigate the effects of a structural variation on the resulting film thickness and thus on friction. In order to keep the calculation times low, a reduction of the flow region as well as a reduction of the geometric complexity are additionally required.
The structure considered in this paper is based on previously published experimental results [
13]. In these experiments, not yet optimized, deepening paraboloidal structures were applied to elastomer samples by laser micromachining. In tribometer tests, a friction reduction up to 50% has been achieved compared to an unstructured and uncoated reference. The geometry derived from the experimental investigation is shown in
Figure 4.
Assuming a uniform flow behavior in the circumferential direction of the sealing lip, a simulation model of a single microstructure with cyclic boundary conditions (unit cell) is sufficient. The geometry of the model is defined by the diameter D, the structural depth t, the length of the unit cell L, and the lubricant film thickness h
0. In order to identify similar geometries and thus reduce the amount of required simulations, the geometry can also be characterized by the dimensionless parameters structural density S = A
Structure/A
Cell, relative structural depth α = t/h
0, and the aspect ratio λ = D/t. Additional design parameters of the radial shaft seal such as the front side lip angle α, the back side lip angle β or the contact width b between the sealing lip and the shaft are not needed to characterize the geometry of the microstructure. However, these parameters are indirectly incorporated into the structural optimization in
Section 3.3, as they influence the mean hydrodynamic pressure p
m.
A two-dimensional fluid flow along the main flow direction in circumferential direction is assumed. The backflow effect of the radial shaft seal, causing an axial flow, is thus not considered in the simulations. Its influence on the hydrodynamic pressure build-up is negligible due to a low flow velocity as compared to the main flow. In addition, the fluid flow is assumed to be incompressible and isothermal. Solid contact between the sealing lip and the shaft is not considered. The fluid properties are based on PAO base oil used in [
13] with a kinematic viscosity ν
40 = 55.5 mm
2/s and a density ρ
15 = 0.834 g/cm
3.
The boundary condition for the velocity is set to a no-slip condition for the surfaces of the moving shaft (sliding velocity U) and the stationary seal. For density and pressure, a zero-gradient condition applies to the shaft and the seal. A cyclic boundary condition is defined at the input and output for all flow variables. Both the seal and the shaft are modelled as solid bodies.
Since no closed solution of the system of equations comprising the conservation equations is available, a numerical approximate solution must be found. Therefore, the conservation equations need to be discretized, i.e., they have to be transferred from a continuous to a discontinuous description. In this work, the Finite Volume Method (FVM) is applied [
24]. In this method, the flow region is divided into a limited number of finite volumes (cells) and the flow variables are calculated at a limited number of vertices (usually in the cell centers). The discretization is based on transferring the differential conservation equations, given in (1) and (2), into an integral formulation and approximating the resulting single terms [
24,
25]. The approximation of surface integrals requires the interpolation of values from cell centers to face centers. The Laplacian and the Gradient terms are discretized using a linear interpolation scheme that is second-order accurate. To ensure numerical stability, the velocity gradient is additionally limited such that the interpolated values at the face centers are limited by the values at the adjacent cell centers. The Divergence terms are discretized using an upwind scheme that sets the value of the face centers according to the value of upstream cell centers [
26]. For spatial discretization, a parameterized, block-structured mesh with hexahedral cells only is used. An estimation of the discretization error caused by the spatial discretization is given in the next Section.
2.5. Error Estimation
The numerical fluid flow calculation based on the conservation laws represents an approximation to the real flow. The deviation between the real flow and the numerical approximation is referred to as the total error of the model ε
tot. In order to draw correct conclusions from the approximation, it is necessary to be aware of the total error, which can be described as the sum of a model error ε
mod and a numerical error ε
num [
27]. The model error ε
mod is defined as the difference between the real flow and the exact solution of the conservation laws. Since in most cases neither the real flow can be measured exactly nor the analytical solution of the conservation laws is known, the model error cannot be specified exactly. The numerical error ε
num is defined as the difference between the exact solution of the conservation laws and the iteratively determined solution of the discretized equations. It is composed of the iteration error ε
it and the discretization error ε
h. Assuming a converged solution for each time step, the iteration error can be neglected and the numerical error can be estimated by the discretization error (ε
it << ε
h).
The discretization error ε
h can be quantified by examining the convergence behavior of the flow variables with increasing mesh resolution. In this paper, the Grid-Convergence-Index (GCI) method developed by Roache [
28] is used. This method is also recommended for reporting uncertainty due to discretization by the fluids engineering division of the American Society of Mechanical Engineers (ASME) [
29]. The resolution of the different meshes considered is increased by a uniform increase in the number of cells. Two exemplary meshes with a low and a medium resolution are shown in
Figure 5.
The simulations are performed without taking cavitation into account using the steady-state, incompressible solver simpleFoam available in openFOAM, which uses the SIMPLE algorithm [
30]. Operating conditions critical to the discretization error are encountered at maximum Reynolds number Re = U h
0 / ν. Therefore, operating conditions with maximum speed, lubricant film thickness and temperature (minimum viscosity) are selected. The simulations are based on a use case with a shaft diameter of d = 35 mm. The speed is set to n = 20,000 min
−1 (see Introduction), which equals a circumferential speed of U
max = 40.8 m/s. As maximum temperature T
max = 100 °C is set, the lubricant kinematic viscosity is ν
min = 5.8 mm
2/s. A constant film thickness of h
0,max = 3 µm is applied and as deviating boundary conditions a value of 0 bar is set for the pressure at the inlet and outlet. Additional mesh studies were carried out for deviating operating conditions with lower Reynolds numbers. Only a small influence of the operating conditions is found, whereas the discretization error is consistently smaller (see
Table A1).
The convergence behavior of the averaged hydrodynamic pressure p
m on the seal surface is used to quantify the discretization error. A graphical overview of the simulation results is shown in
Figure 6 together with the corresponding computing times t. The variable h denotes the discretization parameter, which can be interpreted as a mean cell edge length, such that h decreases with increasing mesh resolution, i.e., with an increasing number of cells.
With an increasing mesh resolution, the pressure p
m decreases and tends towards a finite value for h approaching 0 (number of cells → ∞). Using the GCI method, the results shown in
Table A2 are obtained. The discretization error, expressed as GCI, decreases as expected with increasing mesh resolution and the computing time increases disproportionately. Since the numerical effort increases when cavitation is taken into account, a compromise between a solution that is as exact as possible and a reasonable computing time must be found. Thus, a discretization error of 2.49% for the computational mesh #7 with a number of 18522 cells is considered sufficiently accurate.
3. Results
The simulation model introduced in
Section 2 is used to analyze the fluid flow within the microstructure. In
Section 3.1 the influence of high circumferential speeds on the local fluid flow is examined. Subsequently, in
Section 3.2, the necessity of taking cavitation into account is discussed. Considering these information, a variation of the geometric parameters of the microstructure is carried out in
Section 3.3 and the influence on the resulting lubricant film thickness is considered.
3.1. Influence of Inertial Effects
The fluid mechanical analysis of microstructured surfaces is commonly based on solving the Reynolds equation, which can be derived from the conservation laws shown under certain assumptions [
31]. One of the essential assumptions is the neglect of inertial forces. An investigation of their influence for microstructured surfaces is given in [
32]. In this paper, the authors show that with an increasing influence of inertial effects, an additional hydrodynamic pressure build-up results, which leads to a gap-opening force. They conclude that a fluid mechanical analysis based on solving the Reynolds equation is unsuitable for Re > 1. In [
33], this statement is extended by additionally including the aspect ratio of the microstructure λ. For a rectangular, two-dimensional microstructure, they derive a graphic from which the applicability of Reynolds equation can be assessed. In addition, the consideration of all physical effects influencing the fluid flow in the microstructure is of decisive importance in order to be able to predict the suitability of a microstructure, since it generally can also have a detrimental influence [
34,
35].
Considering the most critical operating conditions with regard to Reynolds number (see
Section 2.5) with a sliding speed U
max = 40.8 m/s, a kinematic viscosity ν
min = 5.8 mm
2/s and a lubricant film thickness h
0,max = 3 µm, a range of Re = 0…21.1 results. For these extreme values, the results of a simulation using the complete Navier–Stokes equation are shown in
Figure 7.
Figure 7a shows the resulting velocity field along with corresponding streamlines.
Figure 7b shows the pressure field along with contour lines of constant pressure. In this simulation, no cavitation was considered.
Analyzing the local velocity field within the microstructure reveals the development of a vortex. The center of the vortex shifts downstream with increasing sliding speed, which effects the hydrodynamic pressure field. Based on the isobars shown, a variation of fluid pressure across the film thickness can be observed. This variation already occurs at low speeds and becomes more severe with increasing sliding speeds. Additional simulations were carried out under a variation of Reynolds number to examine the influence of inertial forces (see
Table 1).
The results of the corresponding simulations are summarized in
Figure 8 and
Figure 9, which show the pressure distribution acting on the seal surface along the relative length of the unit cell x
rel (see
Figure 7).
In
Figure 8, the pressure curves for simulations 1–3 calculated with the complete Navier–Stokes equation (black lines) are compared with those based on the analytical solution of the Reynolds equation (grey lines). The calculation of the pressure distribution based on the solution of the Reynolds equation assumes a constant pressure across the film thickness. This assumption leads to a deviation between the calculated pressure distributions already for a Reynolds number Re = 1. The deviation increases with increasing Reynolds number. As shown in [
36], even at moderate Reynolds numbers Re = 5, significant deviations between the pressure distributions of more than 20% occur when comparing the pressure acting on the shaft surface. At higher Reynolds numbers Re = 21, the influence of inertial forces increases, leading to more significant differences of more than 60%. The increasing deviation with rising Reynolds number is more clearly visible in
Figure 9.
In the figure above, it can also be seen that the location of point p = 1000 mBar shifts upstream (towards negative xrel) and the pressure extrema increase with increasing Reynolds numbers. This leads to an asymmetrical pressure distribution, causing a force with a gap-opening effect. By solving the Reynolds equation, this force is not detectable, since it is caused by inertial effects neglected in that equation. By taking the inertia-induced pressure build-up into account, a fully load-bearing lubricant film is already predicted at lower sliding speeds. Thus, the beneficial effect of microstructuring for friction reduction is already apparent in applications with lower speeds. In addition, the pressure build-up enforces the positive effect with increasing sliding speeds. For these reasons, the further analysis of the fluid flow as well as the structural optimization is carried out by simulations based on the complete Navier–Stokes equation.
3.2. Cavitation
The investigation of the influence of inertial effects has shown that, depending on the operating conditions, the fluid pressure can drop below saturation pressure with increasing film thickness at the beginning of the microstructure. In the real fluid, a change from the liquid to the gas phase occurs under such conditions (cavitation). For structural optimization, this phenomena must be taken into account, hence the further simulations are carried out using cavitatingFoam solver described in
Section 2.3.
Exemplary simulation results are shown in
Figure 10 and
Figure 11 for the operating conditions underlying
Figure 7. The fluid parameters used can be taken from
Table A3. As deviating boundary conditions, a value of 1 bar is set for the pressure at the inlet and outlet. In
Figure 10a, the vapor fraction within the microstructure is shown. From that Figure, it can be seen that at the beginning of the microstructure, liquid is converted into vapor. Downstream, a phase change from vapor to liquid occurs with increasing fluid pressure, thus the cavitation zone is limited to the beginning of the microstructure. The size of the cavitation zone depends on the operating conditions as well as on the film thickness and the geometry of the microstructure. The cavitation zone directly influences the developing velocity field. The vortex shown in
Figure 10b develops only in the area filled with liquid.
Figure 11 illustrates both the pressure field in the microstructure (a) and the pressure distribution acting on the seal surface (b), both of which are also influenced by the cavitation zone.
Within the cavitation zone, the pressure is limited to the defined saturation vapour pressure of psat = 50 mBar. Furthermore, a variation of the pressure across the film thickness remains visible. As a result of the cavitation zone, the point p = 1000 mBar shifts downstream. However, the asymmetry in the pressure distribution remains due to the saturation pressure, which represents a lower limit for the pressure.
3.3. Structural Optimization
In order to achieve maximum friction reduction for the radial shaft seal, an optimum geometric shape of the microstructure has to be found. The optimum shape is defined by a maximum hydrodynamic pressure build-up to increase the lubricant film thickness, thus reducing friction. To find an optimum structure, simulations are carried out by varying the structural depth t and the structural density S for a constant diameter of the microstructure D =30 µm, according to the experimental investigations shown in [
13]. An isothermal fluid flow is assumed considering cavitation. The operating conditions and fluid parameters are summarized in
Table 2.
For the determination of the lubricant film thickness resulting from the hydrodynamic pressure build-up induced by the microstructure, a sequence of stationary simulations is run to reduce the computing time. Starting from two specified film thicknesses h
0,1 and h
0,2, the film thickness h
0,i is determined in an iterative process according to Brent’s method [
37] for which an equilibrium between a defined external load and the load-bearing capacity induced by the microstructure is reached. The external load is expressed as the mean pressure p
m acting on the seal surface. An illustrative result obtained from the evaluation of these simulations is shown in
Figure 12.
The figure shows that both the structural density S and the aspect ratio λ of the microstructure have a considerable influence on the resulting film thickness h
0. A reduction of the structural density up to a value of S ≈ 10% leads to an increase in the film thickness. For lower values, h
0 decreases again. However, the increase in film thickness becomes smaller as the value of S decreases. For a constant structural density, an increasing aspect ratio leads to an increase of the film thickness. From a ratio dependent on the structural density, however, this behaviour reverses and the film thickness decreases again with a further increase in the aspect ratio. For the operating conditions mentioned, a microstructure with a structural density in the range S = 10–30% and an aspect ratio in the range λ = 15–30 (D = 30 µm → t = 1–2 µm) leads to a maximum film thickness of h
0 ≈ 1.4 µm. The friction in the simulations underlying the optimization shown above results from fluid friction only, since solid contact between the seal and the shaft has not been modelled so far. An estimation of the friction regime using the specific lubricant film thickness h
rel indicates that for the surface roughnesses measured in [
13], pure fluid friction is present for h
0,crit > 0.8 µm (h
rel = 4). At these low film thicknesses, fluid friction is dominated by fluid shear outside the microstructure. The friction decreases disproportionately with increasing film thickness (F
F ~ 1/h
0), so that the structure leading to a maximum film thickness also leads to minimum friction.