A numerical model for design and optimization of surface textures for tilting pad thrust bearings

A B S T R A C T A numerical model based on the Reynolds equation to study textured tilting pad thrust bearings considering mass-conserving cavitation and thermal effects is presented. A non-uniform and adaptive ﬁ nite volume method is utilized and two methods are compared and selected regarding their ef ﬁ ciency in handling discontinuities; speci ﬁ cally placing additional nodes closely around discontinuities and directly incorporating discontinuities in the discrete system. Multithreading is applied to improve the computational performance and three root-ﬁ nding methods to evaluate the bearing equilibrium are compared; namely Newton-Raphson method, Broyden's method with Sherman-Morrison formula and a continuation approach with fourth-order Runge-Kutta method. Results from the equivalent untextured bearing are utilized to accelerate the computation of the textured bearing and results are validated by comparison with CFD data.


Introduction
Surface texturing is becoming a promising method for enhancing the performance of hydrodynamic bearings in terms of increasing the oil film thickness and reducing the frictional loss for a safer and more efficient bearing operation. However, successful industrial applications of textured bearings are still limited. One of the main challenges is the dependency of optimum texturing parameters on the type of contact and the operating conditions [1]. A poor texture selection may even lead to a deterioration of the bearing performance. This makes the design of optimized texture patterns a challenging task, which generally requires the utilization of advanced computational models due to the large number of parameters involved. Hence, a successful application of surface texturing relies to a great degree on fast and robust mathematical models that allow an accurate evaluation of the impact of surface textures on the performance of bearings under a wide range of conditions.
The key task in the theoretical analysis of hydrodynamic bearings is the solution of the Reynolds equation to obtain the pressure field, which after integration yields the bearing's main performance parameters, such as load carrying capacity, friction and power loss. While solving the Reynolds equation is quite straightforward for conventional bearings, a number of issues are encountered when simulating textured bearings. For example, texturing can result in the development of multiple cavitation zones and consequently a mass-conserving mathematical treatment of cavitation becomes necessary [2,3]. Also, the fine meshes generally required to capture the complex geometry of textured bearings result in significantly increased computation times. Furthermore, textures introduce numerous discontinuities in the film thickness distribution, which if untreated, can lead to considerable discretization errors. One of the most popular discretization methods in the field of hydrodynamic lubrication is the finite volume method (FVM) due to its simplicity and massconserving properties. Unlike methods based on the weak solution of the Reynolds equation, e.g. finite element methods, the FVM is based on boundary flux approximations, i.e. derivatives at film discontinuities directly depend on the mesh size. Consequently, discontinuities should be treated in order to avoid large discretization errors or high computation times caused by fine meshes. Two ways to deal with discontinuities in finite difference based approaches are available: A local mesh refinement [4,5] and a direct incorporation of discontinuities in the discrete system as proposed by Arghir et al. in 2002 [6]. However, these methods have not been evaluated previously regarding their capability of decreasing discretization errors or reducing computation times. Despite the benefits, discontinuities are rarely directly handled in finite difference based numerical approaches, resulting in unnecessarily fine meshes and high computation times.
Another key step in the analysis of hydrodynamic bearings is the evaluation of the bearing equilibrium, i.e. the specific film geometry that balances the applied load. This generally requires solving the Reynolds equation multiple times for different film geometries. Due to the increased complexity in solving the Reynolds equation for textured bearings, effective methods for finding the bearing equilibrium are crucial. While numerous root-finding methods to evaluate the bearing equilibrium are available, the majority of numerical studies are based on the Newton-Raphson method due to its simplicity and quadratic convergence. However, this method requires the determination of the Jacobian matrix at each iteration and an initial film thickness guess sufficiently close to the actual solution in order to converge. Other methods, such as Broyden's method or continuation methods may provide enhanced stability and computational performance when applied instead of the Newton-Raphson method or in combination with the Newton-Raphson method.
The aim of this work was the development of a fast and robust numerical model to analyse the influence of surface texturing on the performance of tilting pad thrust bearings. To allow for parametric studies and the optimization of texture designs, the model is optimized in terms of computational speed and robustness. The model is based on a finite volume discretization of the Reynolds equation while considering massconserving cavitation and thermal effects. Two methods of handling discontinuities (local mesh refinement and the direct incorporation in the discrete system) and three different root-finding methods (Newton-Raphson method, Broyden's method and a continuation method) are compared and selected based on computation speed and numerical stability. Computation times are decreased by utilizing results from the equivalent untextured bearing and results are validated through comparison with data from commercial CFD published in literature.

Bearing geometry and film thickness
A point-pivoted tilting pad thrust bearing and details of its pads are shown in Fig. 1.
Nomenclature a coefficient for the discrete system (m.s or kg/s) a coefficient for viscosity temperature relationship A control volume face dimension (m or rad) B interface Bernoulli coefficient (Pa) B P Bernoulli coefficient (kg/s) b coefficient for the discrete system (m.s) c p lubricant specific heat (J/kg/K) D damping parameter d f discontinuity coefficient e p ; e e ; e t tolerance value for pressure, equilibrium and temperature solver F nonlinear system for equilibrium solver f interpolation factor G homotopy function h local film thickness (m) h p film thickness at pivot (m) h texture texture depth (m) ii total number of nodes in radial direction J Jacobian matrix J P jump coefficient (kg/s) jj total number of nodes in circumferential direction k con convection parameter m; n coefficients for viscosity temperature relationship n pad number of pads n r number of textures in radial direction n θ number of textures in circumferential direction p local pressure (Pa) p cav cavitation pressure (Pa) Q volumetric flow rate (m 3 /s) q mass flow rate (kg/s) r radial coordinate (m) r i inner pad radius (m) r o outer pad radius (m) r p radial coordinate of pivot (m) T temperature ( C) T f friction torque (Nm) T K temperature ( K) u average fluid velocity (m/s) w 0 applied specific load (MPa) x solution vector for equilibrium solver x; y Cartesian coordinates (m) α relative texture extend in circumferential direction α r ; α θ pitch and roll angle (rad) β relative texture extend in radial direction Γ diffusion coefficient (m.s) δr radial distance from centre of pressure to pivot (m) δW difference in load carrying capacity and applied load (N) δθ circumferential distance from centre of pressure to pivot (rad) ε p ; ε e ; ε t fractional residuals for pressure, equilibrium and temperature solver η lubricant dynamic viscosity (Pa.s) Θ fractional film content θ circumferential coordinate (rad) θ p circumferential coordinate of pivot (rad) θ pad pad angle ( ) λ homotopy parameter ν 40 ; ν 100 lubricant kinematic viscosity at 40 C and 100 C (cSt) ν cSt lubricant kinematic viscosity (cSt) ξ pressure drop coefficient Π frictional power loss (W) ρ lubricant density (kg/m 3 ) ρ r texture density in radial direction ρ θ texture density in circumferential direction ω rotational speed (1/s) ω p ; ω Θ relaxation parameter for pressure and fractional film content D þ pressurized regions D 0 cavitated regions F computational domain Subscripts and superscripts À; þ value just before and after discontinuity eff effective quantity i; j control volume indices ir quantity at pad inner radius k iteration number max maximum min minimum opt optimum out quantity at pad outlet sup supplied quantity W; E; S; N; P west, east, south, north, central nodal value w; e; s; n west, east, south, north boundary value