On the Two-Scale Modelling of Elastohydrodynamic Lubrication in Tilted-Pad Bearings

: A two-scale method for modelling the Elastohydrodynamic Lubrication (EHL) of tilted-pad bearings is derived and a range of solutions are presented. The method is developed from previous publications and is based on the Heterogeneous Multiscale Methods (HMM). It facilitates, by means of homogenization, incorporating the effects of surface topography in the analysis of tilted-pad bearings. New to this article is the investigation of three-dimensional bearings, including the effects of both ideal and real surface topographies, micro-cavitation, and the metamodeling procedure used in coupling the problem scales. Solutions for smooth bearing surfaces, and under pure hydrodynamic operating conditions, obtained with the present two-scale EHL model, demonstrate equivalence to those obtained from well-established homogenization methods. Solutions obtained for elastohydrodynamic operating conditions, show a dependency of the solution to the pad thickness and load capacity of the bearing. More precisely, the response for the real surface topography was found to be stiffer in comparison to the ideal. Micro-scale results demonstrate periodicity of the flow and surface topography and this is consistent with the requirements of the HMM. The means of selecting micro-scale simulations based on intermediate macro-scale solutions, in the metamodeling approach, was developed for larger dimensionality and subsequent calibration. An analysis of the present metamodeling approach indicates improved performance in comparison to previous studies.


Introduction
Tilted-pad bearings are found in lubricated rotating devices as they provide separation under heavy loads while keeping the power losses at a minimum [1]. It is the hydrodynamic action that generates pressure in the lubricant and gives the lifting force, and the bearing is said to operate under hydrodynamic lubrication (HL) conditions, if the pressure is not large enough to deform the bounding surfaces. Titled-pad bearings typically have a facing comprised of a compliant material, e.g., Babbitt and polytetrafluoroethylene (PTFE). This means that they operate under elastohydrodynamic lubrication (EHL) conditions where the pressure that is generated in the lubricant deforms the facing and sometimes also the backing, and significant fluid-structure interaction (FSI) occurs. The modelling approach including elastohydrodynamically generated FSI, is also well-established in the literature under the assumption of perfectly smooth surfaces [2][3][4]. For this purpose, it is conventional to use the Reynolds equation to model the hydrodynamic pressure generation, and to use infinitesimal strain theory to describe the deformation caused by the lubricant pressure. However, in lubricated contacts, surface micro-scale features, e.g., surface roughness, cannot be neglected as it contributes to the tribology at the interface. This also becomes more and more pronounced as the thickness of the film, separating the contacting bodies, becomes smaller and smaller. For example, Etsion [5] indicated that surface texturing can lead to a reduction of friction in EHL contacts and Greenwood and Johnson [6] showed that surface waviness generates ripples in the resulting pressure distributions. Other authors have developed deterministic models which aim to incorporate roughness into the simulations of EHL by modelling the features together with the bearing domain, see an example Reference [7][8][9][10][11][12]. These methods encounter significant computational limitations, due to requirements on grid density associated with the separation of scales between roughness and the contact region. Therefore, no general approach is yet accepted and authors have developed other approaches to modelling surface features in EHL, see e.g., the recent review paper by Vakis et al. [13].
One such approach is homogenization, in which information pertaining to a set of individual micro-scale surface features are characterized over a range of variables and subsequently coupled into a model for the macro-scale flow of the lubricant. Periodicity in the lubricant flow and surface features ensures that the homogenization process produces results representative of the mechanics at both scales. Stochastic theories such as those reported by Tzeng and Saibel [14] and by Christensen [15] are early ways of dealing with the effects of roughness by means of this type of two-scale approach in hydrodynamic lubrication. In his work Reference [16], Elrod was very early to present a homogenized Reynolds type of equation, obtained by multiple-scale analysis. Patir and Cheng [17,18] are precursors of an approach where the terms of the Reynolds equation are multiplied by flow factors to incorporate the effect of the surface features in the bearing fluid flow. The development of homogenization methods for HL has been investigated for a range of applications including compressible flows [19], mixed/boundary lubrication [20,21], tilted-pad bearings [22], porous media [23], inertial effects [24,25], cavitation [26] and soft contacts [27,28]. The rigorous mathematical derivation of these homogenization approaches is significantly different to that of Patir and Cheng [17,18] but exhibit similar responses under certain conditions [29]. Homogenization methods are preferable because they are capable of accurately describing the micro-scale effects in a macro-scale model without the associated computational requirements of deterministic modelling [30]. The requirement to include surface deformation, cavitation and other tribological phenomena has led to the derivation of more complex homogenization-based methods; see e.g., Reference [31][32][33].
One set of homogenization methods applicable to EHL are the Heterogeneous Multiscale Methods (HMM) [34,35], which describe a set of general techniques that facilitate a problem to be solved over multiple scales. The approach is applicable where there is at least an order of magnitude difference between these scales, as is the case when considering the bearing contact region and surface features. Gao and Hewson [36] were the first to derive the HMM in application to EHL and solved the influence of surface topography in sliding line contacts. In Reference [37], de Boer et al., subsequently generalized this two-scale method and adapted the approach for lubricant rheology, non-Newtonian behavior and inertial fluid flow. Part of the two-scale method is the use of metamodeling as a means of coupling the problem scales; this was investigated to allow surface features to be optimized for a minimum coefficient of friction [38]. Further studies based on the twoscale method have investigated line contact geometries [39] and cavitation at the scale of topography [40]. The HMM has also been applied to lubricated contacts by Pérez-Ràfols et al. [41] to investigate stochastic roughness in pressure-driven flows and metal-to-metal seals.
Hydrodynamic cavitation occurs in lubricated contacts when the pressure drops below the ambient value, causing air trapped in the lubricant to vaporize and form bubbles [1]. Conventionally, this does not occur in tilted-pad bearings where the fully-convergent shape of the film thickness prohibits such pressure being obtained. However, where there are divergent features, such as between the roller and the raceways in rolling element bearings or at the scale of surface topography, cavitation will be onset and must therefore be considered in such flows. Modelling of cavitation has been dominated by the well-established Elrod model [42] but the implementation and accuracy of such Heaviside methods are limited and authors have sought alternatives. More recently, the concept of complementarity has been applied to facilitate a more representative transition between the liquid and vapor phases in the contact, see e.g., Reference [43,44].
In connection to this, Pérez-Ràfols et al. [45] presented a very elegant transformation between incompressible and arbitrarily compressible pressure driven flow. de Boer and Dowson [46] recently solved the cavitation problem based on a moving mesh type approach, although film reformation was not implemented. Gao et al. [40] used hyperbolic tangents to reduce the fluid viscosity and density to capture the change in phases over a surface feature and explore micro-cavitation. This concept was more robustly explored by Söderfjäll et al. [47], who solved flows with textured surface features in oil control rings and introduced a model capable of describing the transition in phases accurately. The cavitation model of Söderfjäll et al. [47] is implemented in this article as a means of capturing the effects of micro-scale cavitation.
This article develops the two-scale method to EHL based on the HMM [36][37][38][39][40] to explore threedimensional tilted-pad bearings. This includes modelling cavitation at the scale of the surface features and investigating the difference between ideal and real surface topographies on the solutions generated. The two-scale method is compared to the well-established homogenization methods [22] under HL conditions showing identical results. Subsequently the choice of the EHL model is analyzed in which the effects of material thickness and load capacity are explored. Two different surface topographies are investigated which facilitate ideal sinusoidal surfaces and real experimentally measured surfaces are directly compared. This showed that real surface topography has a stiffer response than in the ideal case. Simulations at the scale of topography indicated how periodicity was maintained and the effect of cavitation in the resulting distributions. Additionally, an analysis of the metamodeling approach used in coupling the problem scales is undertaken showing how the solution evolves and is calibrated throughout the solution procedure.

Materials and Methods
In this section the theoretical formulation for the two-scale method for modelling surface topography in tilted-pad bearings is outlined. Descriptions of both the macro and micro-scale models are given along with details for the metamodeling approach used to couple the scales and the corresponding solution procedure.

Macro-Scale Model
The macro-scale model considers the EHL of a three-dimensional tilted-pad bearing in which pressure generated in the lubricant causes deformation of the contacting surfaces, this interaction is fully-coupled to reach a specified load-carrying-capacity. The macro-scale geometry is depicted in Figure 1, showing an axially-symmetric (periodic) section of the bearing and the collar. The inner radius is 1 and the outer radius is 2 , the bearing facing has a sector angle of and the sector angle of the periodic section . The collar rotates with the angular velocity about the -axis and lubricant, assumed to perfectly stick to the collar's surface, is dragged into the divergent gap between the collar and the inclined bearing facing surface. At the macro-scale a fully separating fluid film is considered, thus assuming that there is no asperity contact nor any starvation at the interface. Additionally, the surface of the rotating collar is modelled as perfectly smooth and the bearing facing surface includes topography. Isothermal conditions are also applied. Both the bearing backing and the collar has a thickness of , and modeled, as they were made of stainless steel, with material data that is listed in Table 1. The bearing pad facing is comprised of PTFE.
where ∇= ( , , ) is the gradient operator for the , , Cartesian coordinate directions, is the identity matrix, = ( , , ) is the vector of deformation components, is the strain tensor and is the stress tensor. The shear modulus = 2(1 + ) ⁄ and Lamé's first parameter are related to the solid material Young's modulus and Poisson ratio as indicated, these must be specified for the different materials in the model. Two different cases are considered: (i) One which models the pad and collar simultaneously; and (ii) another which only considers the pad. The latter case is preferred since the computational requirements to model the simple geometry outweigh the benefits of including the collar in the simulation under certain conditions, as demonstrated by an example later in the article.
For Case (i), the boundary conditions are specified as follows. The lower surface of the collar backing is fixed such that deformation is zero = 0. The upper surface of the pad is loaded with load-per-unit-area * such that • = − * , where is the surface normal vector. The boundaries of the collar representing axial-symmetry of the domain are coupled with a periodic condition such that their deformations are identical . All remaining boundaries are free to deform. It is of note that the load-per-unit-are * relates to the lubricating pressure applied to the pad surface, see Section 2.1.4. A visualization of these boundary conditions is presented in Figure 2a. In Case (ii), the lower surface of the pad is specified using a fixed constraint = 0 and the remaining boundary conditions are the same as Case (i), see Figure 2b. This assumes that in Case (i) the magnitude of deformation of this surface is negligible in comparison to that of the upper surface of the pad, thus allowing the collar backing to be removed from the model.
where and are the coordinate direction mass fluxes in the lubricating region, is the film thickness, 0 is the macro-scale fluid density, 0 is the macro-scale fluid viscosity, and are the components of rotational speed, and 1, , 2, , 1, and 2, are flow factors.
Incompressible and isoviscous conditions are considered such that the lubricant density and viscosity are constant in the macro-scale model. The rotational speed components are specified by defining a cylindrical coordinate system at the origin with the vertical component aligned to that of the Cartesian system: = √ 2 + 2 ; = tan −1 ( ⁄ ). Given an angular velocity of about the vertical axis of rotation yields = cos( ) and = − sin( ). Boundary conditions for the fluid flow problem are required on the outer perimeter of the upper surface of the pad, at each location the pressure is specified as that of the atmosphere = 0 (zero gauge).

Film Thickness
The macro-scale film thickness is given by the sum of the contact separation 0 , the pad geometry , and the normal mechanical deformation of the pad surface = −( • ) . See Equation (7), where in this paper is described by Equation (8), with 1 being the pad differential film thickness. Subsequently, Equation (9) provides a relation between and which is the film thickness used to couple scales of the problem, where the local deformation of the pad is . This local deformation is modelled at the micro-scale, see Section 2.2.

Load-Per-Unit-Area
The load-per-unit-area * and the pressure differ because of the multiscale approach, see Section 2.2. At the macro-scale, * is used to determine deformation and load capacity, whereas is the solution to the macro-scale fluid flow problem Equations (4)- (6). An additional relation is made between these two variables by introducing a further flow factor , see Equation (10). * = In the case where all flow factors are identical to unity, e.g., ≡ 1, the problem reduces to the perfectly smooth problem. However, when the flow factors deviate from unity, the solution incorporates the effect of topography in the macro-scale model. The flow factors are functions of 6 macro-scale variables, which couple the macro with the micro-scale, namely the pressure ; the pressure gradients and ; the film thickness ; and the surface velocities and . In other words, = ( , , , , , ) and the reader is referred to Section 2.3 for the details on how the flow factors are determined.

Load Carrying Capacity
The load carrying capacity of the contact is also modelled at the macro-scale. The requirement is that the bearing should hold a specified load (given as an input) as described by Equation (11), where is the surface area of the pad. To achieve this force-balance requirement, the contact separation 0 is considered as a dependent variable, which takes a unique value given a specified load large enough not to collapse the film completely. A monotonic increase in leads to a monotonic decrease in 0 , and vice versa, and this relation is the basis for a quasi-static solution procedure to achieve the required load. The procedure is outlined in Section 2.4.

Micro-Scale Model
The micro-scale model considers the EHL of a domain independent of the tilted-pad bearing, in which a homogenized pressure , pressure gradients and , film thickness and surface velocities and defined at the macro-scale determine the micro-scale parameters. The model is based on the HMM [32,33] and is valid so long as two assumptions are upheld: (i) the size of the micro-scale domain is an order of magnitude, or more less than that of the macro-scale domain; and (ii) periodicity in geometry and flow conditions are maintained in the micro-scale.
The geometry is specified by considering two length scales and that describe the bounds of the micro-scale domain = {0 ≤ ≤ , 0 ≤ ≤ } , where and are the micro-scale Cartesian coordinate directions. The value of each is chosen based on a wavelength relative to the size of the macro-scale model , = ( , ), where , = 2 − 1 quantify the largest length scale of the macro-scale domain in each of the coordinate directions. The wavelength must be selected such that there is at least an order of magnitude difference between the macro and micro length scales, e.g., < 10 −1 . This assumption limits the size of surface topographies which can be investigated using the two-scale approach. However, the assumption of < 10 −1 is consistent with the scaling used in deriving the Reynolds equation. Conversely, as the wavelength increases the assumption of near parallel surfaces in contact becomes further from that represented by the two-scale mechanics, the model breaks down. Given periodicity in the geometry and fluid flow as the wavelength tends to zero the solution to the micro-scale problem becomes that described exactly by the Reynolds equation.
It is possible to include phenomena in the micro-scale model which are not described by the macro-scale mechanics alone, in this formulation we are concerned with surface topography and the resulting effects of cavitation.

Fluid Mechanics
Micro-scale fluid flow is modelled by considering the conservation of mass and momentum in the lubricant according to the Reynolds equation. Mathematically, the fluid transport at the microscale domain is described by Equations (12)- (14), where and are the micro-scale mass flux (per unit depth) in each coordinate direction, is the micro-scale pressure, ℎ is the micro-scale film thickness, is the micro-scale fluid density, and is the micro-scale fluid viscosity. The solution to Equations (12)- (14) are achieved by specifying periodic boundary conditions for the pressure at the extents of the micro-scale domain ( , ) = (0, ) + ∆ , ( , ) = ( , 0) + ∆ , where ∆ and ∆ are the micro-scale pressure differences in each of the coordinate directions. These conditions are imposed such that the distribution of pressure on each set of boundaries are the same, but their magnitudes differ by a constant value, as defined by the macroscale pressure gradients ∆ = , ∆ = .
In combination with the periodic boundary conditions, a pressure point constraint must be included to solve the problem. This is defined by two requirements: (i) the average pressure in the micro-scale domain should represent the macro-scale pressure; and (ii) there should be no pressure defined on the boundary which is negative, due to the requirements of cavitation. To satisfy the first requirement, a constraint pressure 0 is specified according to Equation (15), which is ideally imposed at the origin of the micro-scale domain. If this condition is upheld and it results in a negative pressure on any of the boundaries, a solution will then not be physically representative when cavitation is considered, see Section 2.2.2. Therefore, a minimum required pressure is calculated = min( 0 , 0 + ∆ , 0 + ∆ , 0 + ∆ + ∆ ) and the pressure point constraint is given by Equation (16).

Cavitation
Cavitation is considered in the micro-scale model because, unlike the macro-scale model where the shape of the contact is fully convergent, the micro-scale geometry will have converging-diverging geometrical features relating to surface topography. Therefore, under low pressure conditions, cavitation will occur in the diverging regions, and if left unmanaged, will produce non-physical negative pressures. It is a conventional to approach the cavitation problem by rapidly reducing the fluid density and viscosity as the pressure drops below a zero gauge to represent a transition from liquid to vapor phases. In this article, Söderfjäll's et al. [47] cavitation model is adopted as described by Equation (17), where is the degree of saturation of the fluid and is a threshold pressure. The micro-scale fluid density and viscosity are subsequently expressed by Equations (18) and (19), respectively, where is threshold constant and takes a value approaching zero. The values of and used in this study are the same as Söderfjäll's et al. [47] and are given in Table 2.

Film Thickness
The micro-scale film thickness is described by Equation (20), where ℎ is a function defining the surface topography and is the micro-scale deformation. This deformation is equal to that of a spring with stiffness-per-unit-area under a load-per-unit-area of * , see Equation (21).
The function describing topography must be periodic on the micro-scale boundaries ℎ ( , ) = ℎ (0, ), ℎ ( , ) = ℎ ( , 0) but can take any value within the domain. Additionally, ℎ must have an average value over the domain of zero such that the average initial film thickness (without deformation) is equal the macro-scale film thickness . The load-per-unit-area is subsequently determined from the average micro-scale pressure, as given by Equation (22), * ( , , , , , ) = 1 ∬ and is a function of the 6 variables used to couple the problem scales. The pressure and load-perunit-area * differ because, at the micro-scale pressure, it is not uniformly distributed over the domain due to the presence of topography and cavitation; as such, there are different contributions made in coupling the fluid flow and solid deformation in the two-scale model. In the model derived here, the micro-scale deformation is a constant value over the domain, to include deformation as a function of the micro-scale domain, an equivalent solid thickness can be derived from the material and stiffness properties of the problem, as previously investigated by the authors in Reference [38,40]. This effect has been neglected due to the large increase in computational costs, where only a 0.5% change in the micro-scale surface topography would be expected under the conditions imposed [38,40].

Homogenization
To couple the micro-scale data to the macro-scale information, it is homogenized and mapped across the scales. The homogenized variables become functions of the 6 parameters used in defining the micro-scale model, resulting in Equation (21) for * , and Equations (23)-(24),  (4), (5) and (10). Therefore, the flow factors describe fluctuations in the mass flux and the load-per-unit-area at the macro-scale due to the inclusion of surface topography and cavitation in the micro-scale model.

Flow Factors
A metamodeling approach is used to calculate the flow factors that couple to the macro and micro-scale models. The Moving Least Squares (MLS) method is employed to derive a multidimensional metamodel of the macro-scale variables = ( , , , , , ) to determine each of the flow factors . Underpinning the calculation of flow factors is the requirement for a range of predetermined micro-scale experiments to be known such that the metamodels can be assessed. This forms the metamodel building stage in which a Design of Experiments (DOE) is used to determine which micro-scale experiments to collect data for. As part of the metamodel building, a Cross Validation (CV) procedure is used to tune the metamodel fit to the known data. Since the macro-scale solution procedure is an iterative process (see Section 2.4), additional micro-scale data is required and the DOE evolves with solution time according to the instantaneous difference in macro-scale parameters. The procedure outlined develops the approach de Boer et al. [39] implemented for investigating highly loaded line contacts and aims to improve both accuracy and efficiency for higher dimensionality when using metamodels as a means of coupling the problem scales.

Moving Least Squares
In each of Equation (4), (5) and (10), the left-hand-side terms are written as a series summation of flow factors multiplied by terms that are only functions of the variables in . Subsequently, by using these equations as the basis functions, a regression analysis allows to be assessed as coefficients, given a known set of data from the DOE.
MLS is a form of conventional least squares regression in which the coefficients do not remain constant but instead become functions of the space in which the metamodel is assessed = ( ). A weighting parameter γ is introduced to change the influence of known data points in the regression analysis. A Gaussian decay type function is conventionally used such that the weight diminishes with increasing Euclidian distance from the assessment location as described by Equation (25), where is the closeness of fit and ɍ 2 is the square of the normalised Euclidian distance between the 'th known experiment and assessment location , and is the number of data points in the DOE. The normalised Euclidian distance is given by Equation (26), where is the normalised 'th component dimension of , and , is the normalised 'th component dimensionof for the 'th known experiment. The normalised variables of are obtained from , , = (( , , ) − min( , )) (max( , ) − min( , )) ⁄ .
An example of the MLS metamodel building is presented in Equation (27) for the variable , where additional subscripts denote the position in the DOE. The metamodel is assessed numerically as a system of overdetermined equations to give the flow factors 1, and 2, as a function of the space described by , a matrix operation is required for each assessment location and there is no analytical form of the metamodel. This formulation is derived in a similar manner for the remaining metamodels, required as part of the two-scale method, and * .

Closeness of Fit
The closeness of fit controls the rate of decay of the MLS weighting functions and therefore provides a means of tuning the response of the metamodel to minimise the error between the predicted and known values. In the case where = 0 the weights become equal to unity and the MLS metamodel reduces to a conventional least squares regression, under these conditions the resulting coefficients are constant. In the opposing case where → ∞ the influence of all data points on the response becomes zero and the metamodel does not provide a measureable quantity.
A value of must be chosen for each of the metamodels required as part of the two-scale method, this is achieved by conducting CV studies on the known DOE data. In this process, a range of values are inspected and the one that produces the lowest absolute error between the known data points and that predicted by the metamodel is selected as the true value. After the true value is known the metamodel can be assessed to produce the flow factors as a function of the design variables.
The CV procedure implemented in this article is known as the Leave-One-Out method. In this method, a range of the closeness of fit parameters are selected. Subsequently, and for each value of , all data points in the DOE are removed in turn and the metamodel built from the remaining DOE values. Then, the absolute error between the metamodel predication and the removed DOE value is assessed. By repeating and averaging all the DOE points, the mean error is given, and by varying , the best value can be selected, that is, the one which minimises the mean error. This final value of is the one for which the metamodel is assessed to determine the flow factors.

Design of Experiments
The flow factors are calculated at the macro-scale from MLS metamodels of the homogenized micro-scale data and are updated as the solution procedure progresses, see Section 2.4. For each interim macro-scale solution, a check is performed on the DOE to determine whether additional micro-scale experiments are required. This check adds DOE data points based on whether the new solution is far from the previously known information, such that the flow factors will be accurately described in the current region of interest.
At any instance within the macro-scale solution procedure, the homogenized variables describe a curve within a 6-dimensional space, as defined by = ( , , , , , ). The curve is given from Equation (28), where is the arc length of the curve . The curve is discretized by dividing it into 1000 evenlyspaced points along . Then, based on the normalised Euclidian distance between the curve and the current DOE points, additional points are added to the DOE if they are more than 1% from it. In the case where no data points exist (first iteration of the macro-scale solution procedure), then the entire curve is added. For each of the added points, the corresponding micro-scale simulations are assessed and the MLS metamodels are updated. In the subsequent CV procedure, the range of the closeness of fit is chosen based on the current known value and increased/decreased by 50%, or in the case where no value (or a zero value) of has been determined, 0 ≤ ≤ 3000 is implemented.

Macro-Scale Solution Procedure
Macro-scale solutions were calculated using the software Comsol Multiphysics [48] in which the domain and governing equations were discretized using the Finite Element (FE) method. In the first instance, solutions without topography were solved to provide initial estimates for the cases where topography was included. To do this, Equations (1)-(3) were solved in the solid phase and Equations (4)- (6) were solved in the fluid phase. The interaction was full-coupled through the boundary load applied to the pad surface and the flow factors were equal to unity = 1. Initially, an arbitrary value of 0 = 5 μm was specified (this value was chosen to allow a comparison to the hydrodynamic case without deformation) and was subsequently changed to balance the required load capacity, according to Equation (29), where is the solution time and is a positive constant. At each increment of 0 in Equation (29), the initial solution to pressure became the solution of the previous time step, for the first time step, a zero initial solution was assumed. The total solution time 0 is chosen arbitrarily and the magnitude of is changed to achieve numerical convergence over 0 ≤ ≤ 0 . In the case without topography, the flow factors are each specified as = 1 over the entire macro-scale domain. Referring to Section 2.1, Case (i), inclusive of the pad and collar together, were discretized using 65,000 2nd order tetrahedral elements in total. A sensitivity study of this mesh demonstrated a less than 1% improvement in predicting the maximum lubricating pressure when a greater number of elements were used. Similarly, for Case (ii) in which the pad was modelled alone, a mesh of 7800 2nd order tetrahedral elements were shown to give the same level of accuracy as Case (i). To allow the results produced to be compared to the literature, see Almqvist [22], a hydrodynamic case without deformation was considered. Case (iii) therefore modelled the pad surface alone in two-dimensions and there was no coupling to the solid, in this model the domain was discretized using 1200 2nd order triangular elements. Again, with no significant change in the maximum pressure calculated when using a more refined mesh. All computations were performed on an Intel Xeon CPU (Intel Corporation, Santa Clara, CA, USA) with 4 cores operating at 3.3 GHz.

Micro-Scale Effects
Micro-scale effects, as described by flow factors, were included in the macro-scale solution by using the results without topography as an initial solution. Based on the 6-dimensional curve describing this solution = ( , , , , , ), the DOE and corresponding micro-scale solutions were calculated and the metamodels built, see Section 2.5. Subsequently, the flow factors were calculated as functions of the initial macro-scale solution and kept constant throughout the remainder of the procedure. The macro-scale solver was then computed for the current flow factors and the new solution generated was checked for convergence in the maximum normalized pressure predicted. The criteria implemented measures the difference in maximum pressures measured between two iterations and divided by the maximum pressure of the previous iteration. If there was no convergence, the DOE check is conducted and additional micro-scale simulations conducted were required. When the size of the DOE was increased, the metamodel tuning parameters were updated by undertaking a further CV study. After the metamodel was updated, the initial solution and flow factors were updated based on the current solution and the process was repeated until convergence was reached. The calculation of flow factors and metamodel building was undertaken using Matlab [49]. All solver tolerances were set to 10 −3 to ensure the accuracy of the returned solutions.

Macro-Scale Parameters
The geometry and operating conditions used to describe the macro-scale calculations undertaken in this article are given in Table 1. These have been chosen to provide meaningful comparisons to the work of Almqvist [22] who studied the homogenized influence of surface topography in hydrodynamic tilted-pad bearings.

Micro-Scale Solution Procedure
Micro-scale simulations were calculated by inheriting the parameterized variables from the macro-scale , , , , and such that for any combination of them, a solution to , and * can be obtained. Comsol Multiphyics [48] was used to solve Equations (12) which is consistent with the constraint pressure specified in Equation (16).
For the micro-scale domain, a mesh of 6000 2nd order triangular elements were used for all surface topographies and the cases were considered. This was found to produce grid independent results with a change of less than 1% in the maximum pressure observed when using a more refined domain. The same CPU, as described in Section 2.4, was used to calculate the micro-scale solutions; the solver tolerance was specified as 10 −3 for convergence.

Surface Topography
Two definitions of surface topography are investigated in the remainder of this article: (i) an ideal topography; and (ii) a real topography. Each represents a periodically repeating surface feature over the macro-scale domain, changes in topography from one location to the next can be achieved at the macro-scale so long as the variation is continuous. However, this is beyond the scope of the current paper. The opposing boundary values of the surface topography function are identical and the length scales must be chosen so that they are consistent with the required scale separation of the HMM used in deriving the model.
The function used to describe the ideal surface topography is given by Equation (31), where is the amplitude and ℎ is consistent with that implemented by Almqvist [22]. In the case of real topography, experimental data for a PTFE pad was measured using a vertical scanning interferometer. This provided height data for the real surface as a function of the Cartesian coordinate directions. Subsequently, the experimental data was processed by implementing a 2dimensional Fourier transformation in Matlab [49], the resulting power spectrum was truncated to keep the top 10% of frequencies and the inverse transformed to create a smooth continuous function representing the real data. To provide a periodic function, the real surface data was mirrored about both the x and y-axes. Lastly, the mean value of the function over x and y was determined and taken away from the mirrored function to give ℎ for the experimental data, this facilitates that the microscale domain represents the macro-scale film thickness without deformation.

Micro-Scale Parameters
To define surface topography and the micro-scale model, additional parameters were introduced and values for these are given in Table 2. These have been selected to generate meaningful solutions, but also to directly compare with the results of Almqvist [22]. Note that the value of ε is selected such that the required separation of scales is maintained. The sizes of the micro-scale domain and are subsequently derived, for which the surface topography function ℎ is scaled accordingly.

Results
In this section results are presented for the two-scale model to demonstrate the following:

Verification of the Two-Scale Model under Hydrodynamic Conditions
A comparison of the results generated by the two-scale method under hydrodynamic conditions and the homogenization approach of Almqvist [22] was undertaken. For this purpose, simulations were calculated under the conditions described by Case (iii) in Section 2.4 without the load balance considered and 0 = 5 μm. At the micro-scale, the ideal topography was used and cavitation was not modelled, this allowed a direct comparison between the two sets of results. Figure 3a,b presents the macro-scale pressure distributions from the two-scale method resulting from this analysis, where two different scale separations are investigated. Both values of satisfy the requirements of the HMM in defining the two scales and are both presented for the same conditions by Almqvist [22]. The pressure distributions for the two scale separations are shown to be very close together in Figure 3a,b, with an 87 Pa difference being the maximum value between any two elements on the pad surface in the corresponding simulations. This result was expected since changing the size of the micro-scale domain under hydrodynamic conditions with no cavitation leads to a direct scaling of the homogenized variables. Additionally, a visual comparison between these results and those presented by Almqvist [22] shows that they are similar. This verifies the two-scale model against the well-established homogenization methods for lubricated contacts under hydrodynamic conditions and indicates the scope of the two-scale method for modelling EHL including the effects of surface topography.

Verification of the Choice of Macro-Scale EHL Model
To establish the choice of macro-scale EHL models, a study into the effect of modelling the pad and backing together versus modelling the pad alone was undertaken. In this analysis, the thickness of the pad was varied between 0.1 mm and 10 mm with the backing thickness kept constant at = 10 mm. The contact separation was kept constant at 0 = 5 μm and the load balance was not considered. A range of calculations were then conducted according to Case (i), and as described in Section 2.4, surface topography was not included in the model and as such all flow factors in the model were set to unity = 1. Average deformations were then calculated according to Equation (32), where is the average deformation of the backing, is the average deformation of the pad and is the surface deformation of the backing. The ratio of these deformations is plotted against the ratio of thicknesses and is presented in Figure 4a. As the pad decreases in thickness, the proportion of backing deformation increases in comparison to that of the pad itself, implying that as the pad becomes thinner, the requirement to model the collar becomes more important. Additionally, for a pad of thickness = 2 mm, and a backing thickness = 10 mm, the ratio of deformations is approximately 0.01, which can be considered small enough to neglect and assume that the backing deformation is negligible. This verifies the choice to model the pad under EHL conditions without including the backing, as described by Case (ii) in Section 2.4. This assumption is carried forward throughout the remainder of the results presented. A visualization of the deformations generated in pad and backing models, under these conditions, is presented in Figure 4b. The result has been mirrored in the axial line of symmetry to show the full bearing and the volume was deformed according to a scale factor of 500 with the deformation field exaggerated and visually observed in the solution.
A further consideration to make is that the choice of macro-scale EHL model is load dependent, with larger loads causing larger deformations of both the pad and collar, thus moving further from the assumption that the backing is fully constrained. Two further simulations were undertaken at the macro-scale using the pad and backing model (for the thicknesses described previously), neglecting surface topography but including the load balance. Two loads were chosen, = 10 N and = 30 N, to show the dependence of increasing the load on the resulting film thickness and deformation distributions. The results are presented in Figure 4a,b for these two loads, respectively, showing deformation of the pad surface. A measurement is introduced in Equation (33) to quantify the differences in the average trailing edge film thickness , where the integration is performed on the boundary [ , 2 ⁄ ], as defined in the polar coordinate system.  The low load result, depicted in Figure 5a, corresponds more or less exactly to a direct scaling of the pressure distribution on the pad, whereas the high load result, depicted in Figure 5b, corresponds more to a Boussinesq type relation to the pressure (see Figure 6a). Further, in the low load case, the average trailing edge film thickness and contact separation are = 17.56 μm and 0 = 17.34 μm, respectively, whereas in the high load case, they are = 6.01 μm and 0 = 4.49 μm. In the low load case, the trailing edge and contact separation are very close together, indicating that the trailing edge does not become out of shape. In the high load case, there is a significant difference between the two measurements that implies that the trailing edge is deforming in such a way that it becomes out of shape under the load. This Boussinesq type relation is exhibited in the high load case only, we have therefore chosen to implement a pad thickness of = 2 mm and load of = 30 N in the remainder of this article so as to demonstrate the conventional solution to tilted-pad bearings additionally, including the effects of surface topography.

Solutions of Macro-Scale EHL Model Inlcuding Surface Topography
Using the two-scale method for EHL outlined results that were generated at the macro-scale, incorporating the effects of both the ideal and real topographies at the micro-scale. Figure 6a,b present the pressure and film thickness distributions for the case without topography. Subsequently, Figure  6c,d present the differences between the pressure and film thickness distributions without topography, and includes the ideal topography, and Figure 6e,f present the differences between the pressure and film thickness distributions without topography and includes the real topography. For all simulations considered in this section, the separation of scales was selected as = 2 −4 and the load = 30 N. It is observed that there are significant differences in both the pressure and deformation distributions at the macro-scale when considering different topographies at the micro-scale. These differences are shown to be greater than the order of magnitude of the calculated values themselves. Figure 6c,e show that the ideal topography generates a larger difference in pressure than in the real topography when compared to the smooth surface case. This corresponds to a larger deviation of film thickness from the smooth surface case, as illustrated in Figure 6d,f. It is of note that the average trailing edge film thickness for the ideal and real topographies were found to be = 9.12 μm and = 8.57 μm, respectively, compared to = 6.90 μm for the case without topography. This demonstrates a stiffer response of the real topography, compared to the ideal, which is compounded by the pad becoming more out of shape (larger deformation of the pad surface). Overall these two observations lead to the hypothesis that a real surface topography leads to a significant difference in hydrodynamic performance, compared to an ideal one, and that neither of these solutions are represented well by smooth surface assumptions.
Convergence at the macro-scale was measured by the maximum normalized pressure between iterations of the solution procedure. The history of this criterion is given in Figure 7a,b for the two topographies investigated. It is seen that the real topography converges in less iterations than the ideal, and implies that deviations from the case without topography are smaller in this case. The convergence histories each demonstrate that the solutions generated are not trivial and require the iterative process outlined. As the initial solutions are updated, each iteration and the metamodels are concurrently modified and the results demonstrate the convergence of both the metamodel's accuracy and macro-scale prediction.

Example Solutions of the Micro-EHL Model
To  The distributions presented show that for the ideal topography there is an exact symmetry of the micro-scale pressure over the domain for both the conditions investigated. This generates a symmetry of the degree of saturation due to the direct coupling between the two variables. For the ideal topography, it is also shown that between the two sets of conditions presented in Figures 8 and  9, there is a rotation of the pressure and degree of saturation solutions corresponding to the difference in velocity applied to each set of simulations. The film thickness distributions presented for the ideal topography are identical between the two conditions, which corresponds to the same average pressure over the domain in each of the results produced.
For the real topography, there is a smaller variation of pressure compared to the ideal topography for the same conditions imposed. Additionally, the plane of symmetry identified in the ideal case is also produced for the real case. This smaller variation in pressure is linked to the film thickness distributions, which have much larger scales in the ideal case than for the real case, this means that the pressure must vary significantly more than the ideal topographical features than the real topographical features. Further, in the real case, it is shown that the degree of saturation is distributed in a way that corresponds to the feature patterns described by the film thickness and this complexity is linked to how pressure is distributed over the surface features. Again, for the real case, there is no significant difference between the film thicknesses under each of the conditions because the average pressures are similar. However, unlike the ideal case, there is no direct rotation of the solutions with the direction of velocity, this is linked to the complexity of the flow as it passes over the real surface features, as compared to the flow over ideal surface features.
Overall, the example micro-scale solutions demonstrate the difference between models inclusive of surface topography and cavitation, as compared to the smooth surface assumptions used to define the macro-scale model. The homogenization of parameters at the micro-scale links directly to the values implemented in the macro-scale model, mapping the differences illustrated in this section to that described in the bearing mechanics.

Metamodeling
To determine the effectiveness of the metamodel in predicting the micro-scale data at the macroscale, a check was performed at various locations on the pad for both the ideal and real topographies. These points were selected arbitrarily as the center of the pad and the locations of (macro-scale) maximum pressure and minimum film thickness in the case without topography. This check compared the percentage difference between the metamodel predictions at the macro-scale by the metamodels (as presented in Section 3.3) to the exact values obtained at the micro-scale. The results of this study are given in Table 3 for the ideal topography and Table 4 for the real topography. Both results show that the metamodel predictions are with 1% of the actual values for all the cases and topographies considered. This is an improvement in performance as compared to previous works [39], despite an increase in dimensionality of the problem, which is attributed to the curvilinear discretization and DOE checks conducted as part of this study. A further study was conducted to analyze the performance of the metamodeling as part of the macro-scale solution procedure. Figure 10a,b show the variation of the size of the DOE with the number of macro-scale iterations for the ideal and real topographies, respectively. Correspondingly, Figure 11a,b show the change in metamodel closeness of fit parameters with the number of macroscale iterations for the different topographies. As the macro-scale solution procedure progresses and the solvers reach convergence in the pressure field, there is also convergence in both the size of the DOE and the closeness of fit parameters. Throughout the first few iterations, the size of the DOE growth is reasonably linear, with almost the full 1000 points added each instance; as the pressure field converges, the solutions on each iteration are closer together and therefore less data points are added. In these cases, the total number of DOE points was significantly different with 11,930 needed for the ideal topography case and 4462 for the real topography case, this is the result of many more iterations being required in the ideal case. The closeness of fit parameters are shown to change significantly over the first few iterations, but as the pressure field converges, so do each of the values, corresponding to the number of DOE points reaching a steady value. When no points are added to the DOE, the closeness of fit parameters remains unchanged.

Discussion
In this article, a two-scale method for modelling surface topography in the EHL of tilted-pad bearings is derived. The approach is based on the HMM and facilitates the homogenized effects of micro-EHL to be coupled into the macro-EHL mechanics using flow factors generated from metamodeling techniques. This work furthers the methods previously investigated by the authors in order to include three-dimensional bearings and a direct comparison of ideal and real surface topographies.
A study was conducted comparing the results generated by the two-scale method with those presented by Almqvist [22] under hydrodynamic conditions. Results were given for two different values of the separation of scale parameters considering the ideal topography. Macro-scale pressure distributions compared directly to those of Almqvist [22] and were identical for the different separation of scales examined. This verified the two-scale method for modelling surface topography under hydrodynamic conditions against the well-established approach in the literature and the scope for the method to analyze EHL, inclusive of the homogenized effects of surface topography.
The choice of macro-scale EHL model was verified by comparing results generated by simulations including the pad and collar together against simulations of the pad alone. This was undertaken to determine the significance of implementing a zero-deformation boundary condition on the pad backing. The ratio of backing to the pad average deformation was shown to increase with a ratio of backing to pad thickness such that the thinner pads generated larger backing deformations. The backing thickness was chosen as = 10 mm and the pad thickness = 2 mm, where the ratio of pad to backing average deformation was 0.01. This influence of the backing deformation was assumed to be small enough to be neglected in the remainder of EHL studies conducted and a zerodeformation boundary condition for the backing was used. For larger loads and thinner pads, this assumption cannot be made and the full pad and collar model must be simulated to correctly describe the deformation in the bearing.
Macro-scale EHL simulations were carried out that compared the results generated without topography to those inclusive of the effects of the homogenized ideal and real topographies by the two-scale method. Simulations including the ideal topography generated larger differences in pressure and film thickness than the simulations including the real topography. Additionally, the average trailing edge film thickness deviated further in the ideal case than the real case from that without topography. Overall this demonstrated a stiffer response of those simulations inclusive of the real topography in comparison to those inclusive of the ideal topography under the conditions. It is hypothesized that the real topography has a stiffer response and this generates a better performing hydrodynamic wedge in comparison to the ideal topography, which leads to an improved load carrying capacity for the bearing.
Micro-scale EHL simulations were undertaken to demonstrate the range of solutions generated at this scale when including surface topography, the examples selected were representative of conditions at the macro-scale and illustrated interesting solutions. Periodicity in the distributions of pressure, film thickness, and the degree of saturation was demonstrated and satisfied the assumptions of the HMM underlying the two-scale method. Symmetry in the solutions under both ideal and real conditions were also observed owing to the periodicity required and the conditions imposed. The responses of the real topography showed more variation compared to the ideal topography and was shown to relate to the difference in complexity of the surface features. However, in the ideal case, the difference in pressure was much larger than the real case, owing to the larger scales of the surface features.
An analysis of the metamodeling techniques used to couple the scales of the problem was also undertaken. A check of the macro-scale EHL solutions generated by the metamodels was compared to the exact values obtained at the micro-scale. Four different locations on the pad surface were chosen, and for both the ideal and real surface topographies, there was no difference of more than 1% between them. It was also demonstrated that as the macro-scale solution procedure progressed, there was almost a linear increase in the size of the DOE used, meaning that the intermediate solutions differ significantly from the previous iterations. As the size of the DOE is increased and more is known about the design space then the intermediate solutions become closer together and less points are added. Eventually the sizes converge, which also corresponds to convergence in the macro-scale distributions, and the closeness of fit parameters used in tuning the metamodel response. Overall the metamodeling procedure was shown to be more accurate and efficient than implementations in previous works [39], despite the increase in dimensionality.

Conclusions
The two-scale method for tilted-pad bearings derived in this paper facilitates the homogenized effects of surface topography at the micro-scale to be modeled at the macro-scale. The method is based on the HMM and furthers past studies [36][37][38][39][40] to include a three-dimensional bearing and a comparison of the ideal and real surface topographies. Results showed that the two-scale method compared well to the established literature in modeling hydrodynamic lubrication [22] and the pad can be modeled without the backing under EHL conditions when the ratio of pad to backing thickness and load capacity are suitably chosen. Macro-scale EHL results demonstrated a stiffer response of the real topography than the ideal topography when compared to the case without. It was hypothesized that this mechanism could explain why such surfaces lead to a better hydrodynamic performance than in the without or ideal topography case. Micro-scale EHL results illustrated the significant differences in flow features for the two topographies considered and that the distributions generated at this scale are consistent with the assumptions of the HMM. Metamodeling solutions showed how the size of the DOE increased with the progression of the macro-scale solution procedure, and that, as convergence was reached in the resulting distributions, there was also convergence in the size of the DOE and the metamodel tuning parameters.
The method derived will be furthered to analyze additional flow phenomena conventionally examined in lubricated contacts-such as compressibility, piezoviscosity, and non-Newtonian behavior of the lubricant, as well as including derivatives across the film and inertial effects. All of which have previously been included for line contact problems [38,40]. Surface topography will be parameterized and the effects on homogenized variables will be analyzed across the scales; this will allow macro-scale variations in topography and the performance of different surface features to be parametrically compared using the approach.