A Two-Scale Approach for Lubricated Soft-Contact Modeling: An Application to Lip-Seal Geometry

. We consider the case of soft contacts in mixed lubrication conditions. We develop a novel, two scales contact algorithm in which the ﬂuid-and asperity-asperity interactions are modeled within a deterministic or statistic scheme depending on the length scale at which those interactions are observed. In particular, the e ﬀ ects of large-scale roughness are deterministically calculated, whereas those of small-scale roughness are included by solving the corresponding homogenized problem. The contact scheme is then applied to the modeling of dynamic seals. The main advantage of the approach is the tunable compromise between the high-computing demanding characteristics of deterministic calculations and the much lower computing requirements of the homogenized solutions.


Introduction
Compliant contacts, most commonly known as soft-contacts, are very common in nature (e.g., cartilage lubrication, eye-eyelid contact) and technology (e.g., tires, rubber sealings, adhesives).It has long been stated that the friction and fluid leakage characteristics of wet soft-contacts are strongly related, among the other factors, to the local interactions occurring at the contact interface [1][2][3][4][5].In the case of randomly rough surfaces, the basic understanding of the role played by the asperity-asperity and fluid-asperity interactions, occurring over a wide range of roughness length-scales, has been largely investigated and debated in the very recent scientific literature [6][7][8][9][10][11][12].Given the (usual) fractal nature of random roughness, a number of interesting phenomena have been highlighted, as for example, the viscous-hydroplaning [6], the viscous flattening [9][10][11][12][13][14][15], the fluid-induced roughness anisotropic deformation [10,11], the local [10,11] and global [8,16] fluid entrapment, and many others.The way to deal with random roughness contact mechanics, despite being nontrivial and suffering of a certain description fragmentation, is however well described in the current scientific literature.
On the other side, nowadays bio-inspired research [17,18], together with the widely-spreading practice of surface engineering [19], is showing the many (mainly unexplored) opportunities offered by the physical-chemical ordered modification of surfaces in order to tailor targeted macroscopic contact characteristics, such as adhesion and friction.Bio-inspired adhesive research [20] is probably the best state of the art example of such research trend.However, investigating the combined effect of, let us say, quantized roughness and fluid action has not equally attracted the scientific community attention, apart from few experimental investigations [12,21,22] and basic theoretical investigations [23,24].This may be justified by the complexity of the numerical formulation of the problem, which is expected to not to present an analytical treatment.As a result, and to the best of authors' knowledge, the combined effect of lubricant action and single-scale (contact splitter) roughness has not been practically investigated by tribologists' community.
In this work we give our personal contribution to this research field and, in particular, we discuss a novel numerical scheme of soft mixed lubrication, which can be adopted to perform such investigations.In our model the roughness is split into two contributions.A threshold scale  ζ th = L 0 /λ th > 1 (where L 0 is the representative macroscopic size of the contact and λ th is the threshold roughness lengthscale) is identified.For length scales λ = L 0 /ζ > λ th , that is ζ < ζ th , the system is investigated by using a deterministic approach, whereas for

Advances in Tribology
the problem is treated by homogenizing the equations.This is of utmost help in performing numerical simulations, as for example in the case of microstructured or bio-inspired surfaces where the surface geometry roughness is characterized by a single scale texture (e.g., a pillars array) combined with random roughness at sub-micrometer scales.The main advantage of the proposed approach is, therefore, the strong reduction of numerical complexity.The paper is organized as follows.In Section 2 we describe the mixed lubrication model scheme (we refer to appendix for the details of numerics), whereas in Section 3 we report an example of application of the proposed model to the case of a lip sealing geometry operating in steady-sliding contact.For a detailed description on dynamic sealings modeling, the reader is referred to [25][26][27][28][29].

Problem Formulation
Here we consider a generic rough compliant solid in steady sliding contact with a rigid smooth counter surface, as shown in Figure 1(a).Due to the coupled action of asperity-fluid and asperity-asperity interactions, the compliant surface gains the actual (or deformed) configuration (Figure 1(b)).In soft-contacts, the shape difference between the deformed and initial configurations cannot be usually neglected, due to the occurrence of large deformations.Nevertheless, the small displacement assumption is usually adopted in the literature of soft-elastohydrodynamics [6,30,31] (soft-EHL), supported by a relevant experimental validation for ballon-flat contact geometries.However, if one analyses the contact at much shorter length scales, that is, at the asperities length scales, then one finds out that the asperities may be completely flattened because of the high contact and fluid pressures.Hence, neglecting the influence of large deformation would lead to strong inaccuracy in describing the evolution of the system at the micro-scales.For this reason, in this work, the small displacement assumption has been relaxed.Moreover, in some cases, for example, for dynamic sealings modeling, the precise calculation of the large tangential displacement of the rubber at the interface of the sliding contact is a must in order to capture key phenomena as the well known reverse pumping effect [26].
Consider now the schematic drawing of Figure 1(a).We assume that the whole texture belongs to the class of Reynolds roughness, that is, with |∇h| 2 1 (where h is the surface height distribution).In such a case, the thin film lubrication formulation is sufficiently accurate in describing the fluid dynamics at the sliding interface at all roughness lengths.The (ensemble) average local separation, u(y 1 , y 2 ) is calculated as: where x i and y i (the subscript i = 1 or 2 or 3, see Figure 1(b)) are the three coordinates of the lip surface points at the initial (x i ) and at the deformed (y i ) state respectively.The quantity u r represents the macroscopic contact shape and h r the undeformed low-pass filtered (i.e., for λ > λ th ) surface roughness.w i is the generic component of the displacement vector describing the average local surface deformation of the compliant body.In the classical EHL approach the displacement w i is calculated within linear elasticity framework [31], by adopting a boundary element approach which requires O(n 2 ) operation where n is the number of discretization points.In our case, however, we adopt a more general nonlinear rheology and the displacement components are calculated by employing a classical finite-element solver, which requires O(n 3 ) operations.Observe that the normal (locally averaged) stress σ 3 (y) is (see also [6]) where y = (y 1 , y 2 , y 3 ), σ f 3 (y) is the locally averaged fluid pressure and σ s3 (y) the locally average solid contact pressure, the latter coming from the asperity-asperity interactions occurring at length-scales λ < λ th .We observe that [6]: where A c (y)/A 0 is the local normalized area of solid contact, where τ f is the average fluid shear stress (see later for more details).We assume that the solid friction shear stress τ s is constant (and directed along the sliding direction), as it happens in the case for a rubber-inert substrate contact [6].
Clearly, the model can be easily extended to include different boundary friction conditions.The relation between the average solid contact pressure σ s3 and the local average interfacial separation u is obtained from Persson's theory of contact mechanics [32,33].In particular, given the local elastic properties of the material and the power spectral density (PSD, see, e.g., Figure 2) x of the surface (h is the surface roughness field, with average value h(x) = 0, q 0 = 2π/λ 0 , where λ 0 is the low frequency cut-off wavelength), the theory allows to calculate the local areal fraction of solidsolid contact α c = A c /A 0 as a function of the locally averaged solid-solid contact pressure σ s3 as where the quantity [∇h(x)] 2 1/2 is the root mean square gradient of the rough surface and is related to the PSD C(q) of the rough surface through the formula [∇h(x)] 2 = d 2 q q 2 C(q ).The relation between the local interfacial separation u and the solid-solid contact pressure is instead determined, in the adhesionless case, by requiring that the change of elastic energy per unit nominal contact area dU el at the interface equals the work done by the contact pressure to deform the elastic body: In ( 6) the elastic energy per unit contact area U el is calculated as a function of the contact pressure σ s3 (see [32,33] for more details).The final formula linking the local separation u and the solid-solid contact pressure is relatively complicated, but at large separation it simplifies and becomes [33] where the quantities β 1 and β 2 can be easily calculated once the PSD of the rough surface is known [33], and We consider now the case where a Newtonian fluid is sandwiched at the interface between the solids.We assume constant viscosity η, constant density ρ, and isothermal conditions.The adoption of a Reynolds roughness, and the assumption of a representative average interfacial separation value u/L 0 1, suggests that the fluid velocity varies slowly with the coordinates y 1 and y 2 compared to the variation in the orthogonal direction y 3 .In such a case, combining the equilibrium with the continuity equation, the homogenized problem formulation reads [9,34] where, as before, u and σ f 3 are the locally averaged (over a length scale given by the longest wavelength surface roughness component, i.e., λ th ) interfacial surface and fluid pressure, respectively.The fluid flow conductivity (tensor) σ eff and σ eff can be related, respectively, to the pressure flow factor (tensor) [7,9,35,36] φ p = 12ηu −3 σ eff and to the shear flow factor (tensor) φ s = u −1 σ eff , which we assume to be a function of the average interfacial separation u only.Here φ p has been calculated on the basis of the Bruggeman' effective medium theory and on the Persson's contact mechanics [34].Within this approach, the effect of solid contacts percolation, as well as of local fluid trapping, on the fluid flow can be taken approximately into account, as recently discussed in [16].For simplicity, we have assumed φ p = φ p I, and φ s = I, where I is the identity tensor.Moreover, by considering the steady sliding condition, the homogenized lubricant equation simplifies into: where U = 2U m is the sliding velocity.However, in order to include the effect of micro-cavitations occurring at large-scale roughness lengths, a JFO cavitation model (constant mixture pressure in the cavitation zone) is adopted.The Reynolds equation is then reformulated under a mass conservative equation valid throughout the cavitating/not cavitating domain:

Lip section Spring
Liquid side

Shaft
Air side y 1 y 2 ω Figure 3: A schematic of a typical lip seal construction.A garter spring enables the radial compression of the seal lip on the shaft surface, ideally preventing fluid leakage from the high pressure side (the fluid reservoir) to the low pressure side.In the running-in stage, the seal lip is macroscopically reshaped as a consequence of the initial lip wear process, which depends on the seal-shaft designed interference and on the solids chemical and physical affinity with the actual lubricant composition.After this bulk material removal (and depending on the viscoelastic rubber properties), a successful seal reaches a steady-state configuration with a well defined roughness pattern, strongly dependent on the actual shaft micro-geometry, and with an established thin lubricant film at the interface.
where we have adopted the cavitation index F (F = 0 in the cavitation zones, and F = 1 otherwise) and the dummy variable φ defined as where ρ 0 is the lubricant density, and G is the representative solid shear elastic modulus.Note that (10) must be solved on the deformed configuration, which is an unknown of the problem: therefore, reformulating the fluid equation in the initial configuration may result numerically convenient.Thus (10) has been rephrased by using the mapping rule The Jacobian J = ∂(x 1 , x 2 )/∂(y 1 , y 2 ) of the transformation can be calculated as the inverse J = g −1 of the Jacobian g = [g i j ] = ∂(y 1 , y 2 )/∂(x 1 , x 2 ) of the transformation Observe that the determinant of the Jacobian tensors must be necessarily larger than zero, that is, det g = (det J) −1 > 0. The above transformation enables the generation of an adaptive mesh grid which follows the changing of the surface shape during the deformation process, thus without loosing spatial resolution.The detailed numerics derivation is reported in the Appendix.
The boundary conditions to be applied in the resolution of (10) depend, clearly, on the particular soft-contact problem under investigation.The dimensionless formulation of (10) is where U = 6μUλ 0 /(2Gh 2 rms ) and α = 2h rms /λ 0 = q 0 h rms /π (note that u has been made dimensionless with h rms , whereas other lengths with λ 0 ).λ 0 corresponds in this case to the largest deterministic roughness wavelength.

Results
In this work we use the proposed model to investigate the lubricated contact of the lip seal schematically shown in Figure 3.Moreover, being λ 0 much smaller than the shaft diameter, it is possible to reduce the computational domain to only a small angular fraction of the lip surface.Therefore, ( 10) is solved with the constant pressure boundary conditions at the low pressure side p = p env and at the high pressure side p = p oil , and with the periodicity conditions on the circumferential y 1 -direction, with spatial period λ 0 .The two-scale mixed lubrication model has been numerically solved, as described in the Appendix.We have opted for a fractal self-affine isotropic geometry.For any self-affine fractal surface the statistical properties are invariant under the transformation In such a case it can be shown that for isotropic surface the PSD is where q 0 = 2π/λ, q = |q|, H is the Hurst exponent of the randomly rough profile, which is related to the fractal dimension D f = 3 − H.To apply our method we need to numerically generate the surface h r (x) in the low frequencies spectrum, that is, for ζ > ζ th , see, for example, Figure 2. To this end we have utilized the spectral method described in [37] where the roughness is described by a periodic surface in the form of Fourier series where q kh = (kq 0 , hq 0 ), x = (x, y).Since h r (x) is real we must have a −h,−k = a hk .Moreover for randomly rough surfaces the following relation must be satisfied a hk a lm = 0 with l / = ±h, m / = ± k, where the symbol • • • is the ensemble average operator.The PSD of Surface Equation ( 17) is from which it follows: For isotropic surfaces we have C(q) = C(q) which simply gives C(q hk ) = C(q 0 √ h 2 + k 2 ) and assuming self-affine fractal surface (see ( 16)) one obtains Hence the quantities |a hk | 2 can be determined once |a 11 | 2 and the Hurst exponent of the fractal surface are known.However to completely characterize the rough profile we still need the probability distribution of the quantities a hk .We first observe that the condition a hk a lm = 0 with l / = ±h, m / = ±k is satisfied if the phases ϕ hk of the complex quantities a hk are random numbers uniformly distributed between 0 and 2π.We also recall the condition a −h,−k = a hk also implies that |a −h,−k | = |a h,k | and that the quantities ϕ −h−k = −ϕ hk .So what we need now is only the probability distribution of |a h,k |.Of course there are several choices and the simplest one is to assume that the probability density function of It can be shown that this choice guarantees also that the random profile h(x) has a Gaussian random distribution.

Sealing at Nearly Zero Sliding Velocity.
Here we report on the model application to the case where the sliding velocity between the two mating surfaces is vanishing ( U → 0), that is, for the case of static seals.The self-affine roughness of the lip surface presents a long-distance cut-off frequency q 0 = 2π/200 μm −1 , fractal dimension D f = 2.4, h 2 1/2 = 10 μm, and small scale cut-off frequency q 1 = 1000q 0 .The threshold frequency adopted in the calculation is q th = 2π/λ th = 5q 0 .In the following, the seal is assumed linear elastic.The calculation is performed at a constant rigid penetration of the lip.In Figure 4 we show the locally averaged interfacial separation field u (note that x 1 is the circumferential direction, whereas x 2 is the axial direction).Interestingly the corrugation which can be observed in the figure is just a consequence of the deterministically included roughness (ζ < ζ th ), at larger frequency the surface appears smooth since the high frequency (ζ < ζ th ) content of the roughness has been included through Persson's statistic model, that is, by means of homogenization.The corresponding average solid contact pressure is shown in Figure 5(a).Observe first that, due to the presence of roughness, the contact is split into many contact patches in the whole lip apparent contact area and, correspondingly, normal stresses are concentrated into contact spots.Observe also that the average solid contact pressure is smooth at the contact borders, differently, from what instead expected in the case of smooth elastic contact (e.g., in the case of the Hertzian contact).This is due to the homogenized roughness contribution, which distributes, on a wider local contact area, the pressure acting on the single deterministic asperity.The average fluid x 2 (μm) (with x 1 = 0) pressure is shown in Figure 5(b).Observe that the fluid pressure presents steep variations (almost step-like) in those locations where the separation between the surfaces takes the smallest values.This is in perfect agreement with the critical junction theory of leakage in seals [8], which predicts that the hydraulic resistivities are only concentrated in few points where the fluid flow encounters strong restrictions.At these restrictions the pressure must present high gradients as clearly shown in Figures 5(b) and 6.In particular, in Figure 6 we report the average interfacial separation (black curve), the average solid contact pressure (green curve), the average fluid pressure (red curve), the total average normal stress (blue curve, slightly distinguished from the average solid contact pressure), and the solid contact pressure in the case of perfectly smooth contact (black dashed curve, useful for comparison), for x 1 = 0 (i.e., along the axial direction).Note that the largest pressure gradients occur in the very proximity of the largest values of σ s3 , that is, where the local average separation takes its minimum value.It is also interesting to observe the normalized solid contact area A c /A 0 for this static contact case, see Figure 7(a).Note that, due to the high local squeezing pressure, the asperities lowfrequency features (i.e., the roughness asperities described by the deterministical model) have been squeezed so much to coalesce in larger contact patches.As expected [8], the presence of asperity-asperity contacts increases the hydraulic resistivity at the interface and, in particular, largest values of fluid velocity occur at the local minimum in the average interfacial separations.
We have also carried out calculations assuming the seal obeys a Mooney-Rivlin model.The results are then compared with the elastic case.In Figure 8 we show the normalized area of real contact for the elastic (Figure 8(a)) and hyperelastic (Figure 8(b)) bulk rheology.Observe that for the linear case, the area of contact are slightly larger, as expected, then for the other case.This however, does not  strongly affect the leakage flow for the current geometry, as shown in Figure 9.However, interestingly, the nonlinear rheology is characterized by a larger value of surface area change, evaluated as ∂w 1 /∂x 1 + ∂w 2 /∂x 2 (see Figure 10), where we show in red colors that on part of the interaction domain there is a shrinking, that is, ∂w 1 /∂x 1 + ∂w 2 /∂x 2 < 0.
Observe that in plane surface displacements are responsible for the frequency-shift of roughness PSD, for example, in the case of Figure 10(b) we expect a relevant modification of the homogenized power spectral density shape, which can therefore affect the local contact mechanics and flow factors.
In the literature, surface displacements are not usually taken into account; however the present study shows that such a surface effect can be relevant for contact modeling.

Sealing at Non-Zero Sliding Velocity.
In this section we show the case of non-zero sliding velocity.The macroscopic lip geometry is the same as before, as well as the roughness PSD.However, this time the part of roughness that is described deterministically is characterized by only one length scale, that is, just one frequency is included.The remaining roughness has been therefore included within the homogenized approach.In Figure 11 we show the cavitation Interestingly, cavitation originates as expected in the low pressure side of the region (i.e., on the left side of the domain, whereas sliding velocity is directed from top to the bottom) and at the trailing edge of asperities.By increasing the sliding velocity, the cavitation extends from the low to the high pressure side.The cavitated areas may even coalesce at the largest sliding velocity, with the formation of cavitation  fingers shown in Figure 11(e).The adoption of our twoscale approach allows then to capture the complex solid contact and fluid-dynamics characteristics of real contact geometry, which can help engineers to have useful insights into the mixed-EHL lubrication condition occurring at the interface of soft-contacts, especially in terms of friction and leakage, with much less computation effort.Figure 12 shows the axially component of fluid flow at the interface.Red areas corresponds to the counter-gradient flow (i.e., the flow directed in the opposite direction compared to the Figure 13: Algorithm scheme.In the two-scales approach, the solution of the contact problem is split into two, coupled problems, respectively, the deterministic and the homogenized problem.Those are solved iteratively towards convergence. surface stress is determined by solving the homogenized part, which consists of the Persson's contact mechanics and of the homogenized fluid problem (previously discussed).
The deterministic part follows, where the macroscopic deformation problem is solved as a function of the previously calculated average surface stress field.To do so, in this work we have used the Ansys finite element code; in particular, the lip geometry, similar to that described in [27], has been meshed with tetrahedral structural elements of type 92.The macroscopic gap relation is then finally updated, determining the loop restart.The solver iterates until certain convergence criteria are satisfied.In our case, convergence is checked on the average interfacial separation field u(x 1 , x 2 ).Under-relaxation, with relaxation factors in the range [0.01, 0.1], is usually adopted to numerically damp the interfacial separation solution.

Nomenclature η:
Fluid viscosity ν: Poisson's ratio α: Reduced root-mean-square roughness α c : Local normalized area of solid contact u: Locally averaged interfacial separation β 1 : Parameter for the asymptotic interfacial separation law β 2 : Parameter for the asymptotic interfacial separation law λ th : Threshold roughness wavelength ρ 0 : Full film density σ eff : Shear flow conductivity σ 3 : Normal (locally averaged) stress σ eff : Pressure flow conductivity σ f 3 : Normal (locally averaged) fluid stress σ s3 : Normal (locally averaged) solid contact stress τ f i : Tangential (locally averaged) fluid stress i-component τ i : Tangential (locally averaged) stress i-component τ si : Tangential (locally averaged) solid contact stress i-component U: Reduced sliding velocity ζ: Roughness magnification or wavenumber ζ th = L 0 /λ th : Threshold roughness wavenumber A 0 : Representative area of interaction A c : Area of solid contact in A 0 C(q): Total roughness power spectral density C sto (q): Power spectral density of the statistically-calculated roughness E: Y o u n g ' sm o d u l u s E * : Reduced elastic modulus F: Cavitation index G: Shear elastic modulus g i j : Mapping rule Jacobian h r : Deterministic surface roughness height function J i j : Inverse of the mapping rule Jacobian L 0 : Representative macroscopic size of the contact U: Shaft sliding velocity U m : M e a nv e l o c i t y u r : Macroscopic contact shape function U el : Locally stored elastic energy w i : Generic component of the average surface displacement vector x i : Initial state reference y i : Actual or deformed state reference.

Figure 1 :
Figure 1: Contact scheme: (a) reference geometry; (b) deformed geometry obtained by solving the fluid/solid contact problem.

4 ]Figure 2 :
Figure2: A self-affine power spectral density of a generic surface roughness.Points correspond to that spectral content included deterministically in the calculations, whereas the continuous line represents the homogenized surface roughness.

Figure 4 :
Figure 4: Dimensionless average interfacial separation u/h rms as a function of contact position, reported in the undeformed reference.

fluid pressure σ f 3 Figure 5 :
Figure 5: Average fluid and solid contact pressure for the static linear elastic seal.

Figure 6 :
Figure 6: All solutions in arbitrary scaling.

Figure 7 (
b) shows a sample of the leakage flux lines calculations (red lines) superposed to the locally averaged fluid velocity intensity (in gray scale).
Sample of fluid flux lines

Figure 8 :
Figure 8: Normalized area of real contact for different bulk rheologies.Dashed contours represent the bearing area paths.