An artiﬁcial viscosity to model the eﬀect of roughness

– This paper presents a new way to consider the eﬀect of roughness in ﬂuid ﬁlm lubrication. An artiﬁcial viscosity function is introduced to model the eﬀect of roughness. The parameters of the function are derived either from the homogenized solution or the stochastic solution. With this method it is possible to determine the average velocity proﬁle, which cannot be done with the usual homogenization or stochastic methods


Introduction
The objective of fluid film lubrication is to separate two surfaces in relative motion by a thin film of fluid to a Corresponding author: noel.brunetiere@univ-poitiers.fr avoid contact and wear and to reduce friction.Because of the thickness of the film, perturbations induced by the surface roughness became an issue a long time ago.Tzeng and Saibel [1] were among the first to propose a solution for the pressure distribution in an infinite bearing with random rough surfaces.In the same period of time, Christensen [2] presented an averaged Reynolds equation obtained by a stochastic model.These two approaches are, however, limited to striated rough surfaces, while the topography of real surfaces often varies in both directions.About ten years later, Patir and Cheng [3,4] proposed a new approach based on numerically derived flow factors which overcomes this limitation.The flow factors are included in an averaged Reynolds equation to consider the roughness induced perturbations on the fluid flow.This method has the advantage of not being limited to striated surfaces, and has thus been used in many applications.Many improvements can be found in the literature.Elrod [5] and then Tripp [6] performed a theoretical calculation of the flow factors by using an asymptotic expansion.Later, Harp and Salant [7] included the effect of inter-asperities micro-cavitation in this model.
An alternate approach based on homogenization was proposed by Bayada and Faure [8].The main idea consists in using two scales, the smaller being used to describe a periodical micro roughness.With this method, it is possible to consider all types of orientations of surface roughness, which is more difficult with the flow factors method.This approach has also been improved to include cavitation [9], as well as real measured surface roughness [10,11].
The improvement of computer performance now allows performing deterministic simulations where a real surface topography is included in the model without any averaging procedure [12,13].However this method is still limited to small surfaces compared to averaged roughness methods because of the computation burden.
In some applications, it is of interest to obtain not only the fluid pressure.For example, the resolution of the energy equation in the film requires also knowing the velocity distribution [14].However, stochastic and homogenized approaches are not relevant for determining the velocity profile.In the present paper, an artificial viscosity function is introduced to model the effect of roughness.This artificial viscosity can be compared to the turbulent viscosity used in fluid mechanics to model turbulent random eddies [15].The parameters of the viscosity function are derived either from the homogenized solution or the stochastic solution.One of the main advantages is the possibility of calculating the velocity profile based on usual through film viscosity integral functions [16].Another possible field of application could be to consider roughness effects in the bulk flow type equations used for inertial flows in thin films [17].This procedure can be used for 2D isotropic roughness (as the viscosity is isotropic) or with the infinitely long device assumption (1D Reynolds equation).

Averaged Reynolds equation incorporating roughness
The studied problem is presented in Figure 1.A smooth and flat surface is moving with a velocity U in the x direction.It is separated from a stationary rough surface by a distance h.The gap between the two surfaces is filled by an incompressible fluid of viscosity μ 0 .The local film thickness h is a function of the nominal film thickness h 0 describing the macroscopic geometry and of a roughness function f representing the local asperities height variation: where is a dimensionless parameter controlling the height of the roughness.This roughness is supposed to be isotropic.Depending on the type of roughness (periodical pattern or random distribution), the Reynolds equation governing the pressure distribution can be derived by a homogenization or by a stochastic approach.For an isotropic 2D roughness, the two approaches lead to the following Reynolds equation: (2) where p is the averaged or homogenized pressure.φ p and φ s are correcting factors used to model the effect of roughness.The averaged or homogenized shear stress on the smooth moving surface along the x direction can be written φ τ p and φ τ s are correcting factors used to model the effect of roughness.
For an infinitely long device, the following 1D Reynolds equation is valid: The shear stress equation remains unchanged.

1D Periodic roughness
In the case of a 1D periodic roughness profile, the homogenization leads to an analytical expression for the correcting factors.The film thickness can be expressed by: where x 1 is the local scale of the roughness.f is a periodic function of x 1 of wave length 1.Let us introduce the following notation: with this notation, it can be demonstrated that the correcting factors can be expressed as [8]: The shear factor φ s and the pressure friction factor φ τ p are equal.

2D Random roughness
The homogenization technique is able to deal with random roughness but the expression for the flow factors cannot be computed analytically in terms of f .Only numerical computation is available [8].In this case, it is better to use the stochastic method.If the roughness function f is isotropic with an exponential auto-correlation function and a standard deviation equal to one, it can be shown with the Tripp method [6] that As with the homogenization approach, it is found that φ s = φ τ p .In this particular case, these two factors are equal to the pressure factor φ p .Note that many engineered surfaces are found to have an exponential autocorrelation function [18].

Reynolds equation with variable viscosity
In this section, it is assumed that there is no roughness, f = 0, but that the viscosity of the fluid varies through the film thickness: where z is the through-film coordinate.This is a wellknown problem developed in the case of ThermoHydro-Dynamic (THD) lubrication in bearings [14,16].Let us introduce the following notation for the through film viscosity integrals: Using this notation, it can be shown that the Reynolds equation for smooth surfaces is (see Appendix A): The shear stress on the moving wall can be expressed as The velocity profile across the film thickness can be easily expressed:

Equivalence between averaged and variable viscosity Reynolds equations
The Reynolds and the shear stress equations with a variable viscosity are similar to those obtained with rough surfaces by homogenization or the stochastic approach.By identification, it is found that they are identical if These relations suggest that it should be possible to model the roughness effect by using an artificial viscosity varying through the film thickness.If there exists a viscosity function g such that the corresponding integrals I 0 , I 1 and I 2 satisfy Equations ( 18) to (20), then the pressure and shear stress will be identical.The description of the average flow due to roughness will be the same as the description of the flow with a variable viscosity along a smooth surface.Equations ( 18) to ( 20) are equivalent to Let us assume that Because of the three equations, the function must contain three parameters.The choice of G is not unique.However for each function satisfying Equations ( 7) to (9) or Equations (10) to (12), the pressure and shear stress fields are the same.The objective of this paper is not to find the best function or the function which provides the most physical solution.The objective is to show the possibility of using a variable viscosity to model the effect of roughness.In the following, a polynomial function which can be analytically integrated is chosen.In order to enforce that the viscosity is not modified close to the smooth wall, we impose the following constraints: Function G is thus defined by: It can be shown that: The parameters A, B are C are solutions of the system The determinant of the matrix is different from zero and the system always has a solution.The matrix can thus be inverted:

Numerical application to periodic roughness for an infinitely long device
In this section, it is assumed that the roughness is a "crenel" function: (33)

Example 1: Viscosity and velocity profiles
By solving the system (32), it is possible to find the values of A, B and C and then the function G (the solution is given in the Appendix B).The profile of this function is presented in Figure 2.This function is equal to 1 close to the smooth wall (bottom), increases slightly, and then decreases sharply when the distance to the rough wall is reduced.As is logical, the magnitude of the variations increases with the roughness height .G can become negative for greater than about 0.37.In this case, the viscosity function 1 + g, which is the inverse of G, would be singular.This situation must be avoided.This is a limitation to the artificial viscosity method.
The corresponding viscosity roughness functions 1 + g = μ μ0 are presented in Figure 3.The artificial viscosity is increased close to the rough wall when the roughness height increases.Moreover, the maximum effective viscosity value is shifted toward the center of the channel when the roughness height is increased.
The advantage of the viscosity approach compared to homogenization is that it is possible to compute the averaged velocity profile of the fluid in the lubricating film.Some examples of velocity profiles corresponding to a pure pressure flow and a pure shear flow are presented  in Figure 4 and compared to smooth wall solutions.Generally speaking, the fluid velocity is reduced because of the higher friction (and viscosity) induced by the roughness.For the pressure flow, the maximum of the speed is shifted toward the smooth wall (bottom) where the friction is lower.For the shear flow, an almost linear speed distribution is observed close to the smooth wall but with a steeper slope than for the smooth surfaces solution.

Example 2: 1D slider bearing
Although not previously mentioned, the present analysis is valid for a non-uniform gap h 0 (x).In this case, , A, B and C as well as G and g vary along the x direction.It is thus possible to use the present approach to study a 1D slider bearing.In this case, the nominal film thickness is defined by where L is the length of the contact and h o the outlet nominal film thickness.In this case the inlet film thickness is 3ho 2 .The roughness height is constant along the contact.
In Figure 5, the pressure profiles obtained with different methods are presented.For the deterministic approach, 10 waves were used for the contact.This exact approach leads to pressure oscillations along the contact because of the height variations.The homogenization method gives an excellent approximation of the average pressure variations.As expected, the present viscosity model gives exactly the same pressure profile as the homogenized pressure profile.If roughness is not considered in the simulation, the pressure is significantly underestimated.
The shear stress profiles on the moving wall obtained with the different approaches are presented in Figure 6.
Here also, the viscosity model gives exactly the same results as the homogenization, which is a good approximation to the deterministic solution.
The load generated in the contact and the friction force are presented as a function of the roughness height in Figures 7 and 8.The viscosity model and the deterministic model lead to results in close agreement, showing that the roughness tends to increase the load as well as the friction.The slight difference on the load could be due to the too small number of waves (10) used for the deterministic simulation.

Numerical applications to random roughness in 2D devices
In this section, the roughness function f is a Gaussian height distribution having an exponential autocorrelation function.The height of the roughness is = σ h0 where σ is the standard deviation of the roughness distribution.

Example 1: viscosity and velocity profiles
By solving the system (32), it is possible to find the values of A, B and C and then the function G (the solution is given in the Appendix).The profile of this function is presented in Figure 9.This function is equal to 1 close to the smooth wall (bottom), increases slightly, and then decreases sharply when the distance to the rough wall is reduced.As is logical, the magnitude of the variations increases with the roughness height .As with the crenel roughness profile, G can become negative for greater than about 0.33.In this case, the viscosity, which is the inverse of G, would be singular.This is a limitation to the viscosity method.However this situation corresponds  approximately to the appearance of contact between the two surfaces, for which the viscosity concept makes no sense.
The corresponding effective viscosity functions are presented in Figure 10.The artificial viscosity increases close to the rough wall when the roughness height increases.The viscosity distribution is different from the one obtained with a crenel roughness profile.For example, the location of the maximum of the effective viscosity is almost constant, whereas it changes with the roughness height for a striated surface (see Fig. 3).Some examples of velocity profiles corresponding to a pure pressure (Poiseuille) flow and a pure shear (Couette) flow are presented in Figure 11 and compared to smooth wall solutions.As with 1D roughness, the fluid velocity is reduced because of the higher friction (and  viscosity) induced by the roughness.For the pressure flow, the maximum of the speed is shifted toward the smooth wall (bottom) where the friction is lower.For the shear flow, an almost linear speed distribution is observed close to the smooth wall but with a steeper slope than that of the smooth surfaces solution.Because of the different viscosity distribution, the velocity profiles are different from those obtained with a crenel roughness.The viscosity in the vicinity of the rough wall is decreased leading to a higher speed than with smooth surfaces.Indeed, with a smooth surface, all the fluid has a zero speed at z = h 0 whereas, with rough surfaces, a part of the fluid has the possibility of flowing between the asperities ending at z ≈ h 0 + 3σ.

Example 2: 2D slider bearing
The different approaches are now used to study a 2D slider bearing.In this case, the nominal film thickness is  defined by where L is the length of the contact and h o the outlet nominal film thickness.The inlet film thickness is 3ho 2 .Note that the width of the contact is equal to length of the contact L. The standard deviation σ of the roughness is constant in the contact.For this comparison, five different rough surfaces were used.These surfaces are Gaussian and isotropic.The correlation lengths of the surfaces vary from 0.022L to 0.053L.They were numerically generated by the Hu and Tonder algorithm [19].
In Figure 12, the pressure profile obtained in the center line of the bearing with the deterministic (averaged over 5 surfaces) and viscosity methods are presented.The viscosity method gives a good approximation of the pressure increase due to roughness.A better correlation would certainly be obtained by using more than 5 surfaces.Note that the solution given by the viscosity approach is exactly the same as the one given by the stochastic approach, which is thus not presented here.The shear stress profiles on the moving wall obtained with the different approaches are presented in Figure 13.The viscosity model provides a higher shear stress than the smooth surfaces solution, which is in agreement with the shear stress averaged over the 5 deterministic solutions.
The load generated in the contact and the friction force are presented as a function of the roughness height in Figures 14 and 15.Generally speaking, the load and the friction increase with the roughness height.However, the evolution can be different from one surface to the other.The viscosity model is a good approximation of the results obtained with the five different rough surfaces.

Conclusion
In this paper a new method to treat the roughness effect in fluid film lubrication has been presented.The idea is to introduce an artificial viscosity function.The parameters of this function can be set so as to obtain similar results as those calculated with the homogenization method or the flow factors method.The great advantage is that it is possible to derive the velocity profile, which is not possible with the usual approaches.The method has been satisfactory compared to deterministic solutions.The approach is, however, limited to isotropic roughness.An anisotropic viscosity function could be introduced to deal with more complex roughness patterns.Moreover, a first extension to the case of two rough surfaces is proposed in the Appendix D of the paper.

A The Reynolds equation with a variable viscosity function
The Reynolds equation is derived for the configuration given in Figure 1.The simplifed Navier-Stokes equation for a thin viscous fluid film is given by [16]: Since p does not vary with z, the two first equations can be integrated twice to find the velocity components V x and V y : The continuity equation integrated through the film thickness is: Considering that V x (h) = V y (h) = 0, the equation can be modified in this way: By replacing the velocity components V x and V y by their expressions, the following Reynolds equation is finally obtained:

D Extension to two rough surfaces
In many situations, the two surfaces involved in the contact exhibit some level of roughness.For two rough surfaces, the situation is more complex, as several scales exist in both space and time.However the homogenized situation can also be associated to a Reynolds equation with known coefficients written between two smooth surfaces.They are obtained by first making a space homogenization (with time acting as a fixed parameter) and then by making some local time averaging [20].In the case of random isotropic roughness, simple expressions for the flow factors can be found.Let us assume that 1 is the relative roughness height of the bottom moving surface and 2 is the relative roughness height of the top stationary surface: An inverse viscosity function G can be computed for each surface using the exact solution presented in the Appendix.When expressing the polynomial solution, it is necessary to invert the z direction for the bottom surface: If the top surface is smooth, 2 = 0, the viscosity is It can be demonstrated that the correcting factors solution with 2 = 0 and 1 = 0 is recovered.When the two surfaces are rough, it is possible to propose several viscosity function combinations.The first one  Due to the complex expression for the inverse of the viscosity, it is not possible to obtain a simple analytical expression for the integrals I n .Another way consists in using the product of the contributions: (D.8)This gives a very simple expression for the inverse of the effective viscosity.These two solutions will be compared to the exact flow factors expression (Eqs.(D.1) to (D.3)).For this, 1 and 2 are varied in the range from 0 to 0. 1 is lower than 0.33 ( 2 < 0.1).This value corresponds to the appearance of contact between asperities.Above this value, some significant deviations are observed for φ s and φ τ s .This is not surprising since the artificial viscosity function for two rough walls was arbitrarily and not theoretically derived.
An interesting point is that the two ways to combine viscosity functions (product and sum) lead to almost the same results and both of them can be used.As can be seen in Equation (D.9), this means that the product of g 1 by g 2 is small compared to unity.Indeed these functions are   [1 + g 1 (z)] [1 + g 2 (z)] = 1 + g 1 (z) + g 2 (z) + g 1 (z) g 2 (z) ≈ 1 + g 1 (z) + g 2 (z) ( D .9 )

Fig. 2 .
Fig. 2. G function distribution for a crenel roughness profile and different values of the roughness height .

Fig. 3 .
Fig. 3. Roughness function for a crenel roughness profile and different values of the roughness height .

Fig. 6 .
Fig. 6.Shear stress distribution on the moving wall of a 1D slider bearing with a crenel rough wall, (L) = 0.3.

Fig. 7 .
Fig. 7. Generated load in a 1D slider bearing with a crenel rough wall as a function of the roughness height (L).

Fig. 8 .
Fig. 8. Friction force in a 1D slider bearing with a crenel rough wall as a function of the roughness height (L).

Fig. 9 .
Fig. 9. G function distribution for a random roughness and different values of the roughness height .

Fig. 10 .
Fig. 10.Effective viscosity function for a random roughness and different values of the roughness height .

Fig. 12 .
Fig. 12. Pressure distribution in the center line of a 2D slider bearing having a random rough rough wall, (L) = 0.2.

Fig. 13 .
Fig. 13.Shear stress distribution on the moving wall of a 2D slider bearing having a random rough rough wall, (L) = 0.2.

Fig. 14 .
Fig. 14.Generated load in a 2D slider bearing with a random rough wall as a function of of the roughness height (L).

514-page 7 NFig. 15 .
Fig. 15.Friction force in a 2D slider bearing with a random rough wall as a function of of the roughness height (L).

Fig. D. 1 .
Fig. D.1.Comparison of the resulting pressure flow factors for two rough surfaces with the analytical expression.

Fig. D. 2 .
Fig. D.2.Comparison of the resulting shear flow factors for two rough surfaces with the analytical expression.

3 . 2 2 + 2
The effective viscosity function is numerically integrated and the resulting flow factors are compared to the analytical expressions in Figures D.1 to D.3.A good agreement between the viscosity method and the exact correcting factors is obtained when the combined roughness =

Fig. D. 3 .
Fig. D.3.Comparison of the resulting shear friction factors for two rough surfaces with the analytical expression.

Fig. D. 4 .
Fig. D.4.Artificial viscosity of the upper wall (μ1 with 1 = 0.25), of the lower wall (μ2 with 2 = 0.2), of the two walls by product and sum methods.
close to zero in the half channel where the other function is different from zero (see Fig. D.4).