Rough surface elastohydrodynamic lubrication and contact mechanics

Full text

(1):. L ICENTIATE T H E S I S. Rough Surface Elastohydrodynamic Lubrication and Contact Mechanics. Andreas Almqvist. Luleå University of Technology Department of Applied Physics and Mechanical Engineering, Division of Machine Elements :|: -|: - -- ⁄ -- .

(2) 2004:35. ROUGH S URFACE E LASTOHYDRODYNAMIC L UBRICATION AND C ONTACT M ECHANICS. A NDREAS A LMQVIST. Luleå University of Technology Department of Applied Physics and Mechanical Engineering, Division of Machine Elements. 2004 : 35 | ISSN : 1402 − 1757 |ISRN:LTU-LIC--04/35--SE.

(3) Cover figure:. A modeled surface topography pressed against a rigid plane, assuming linear elastic surface material. The theory describing the contact mechanics tool used to produce this result is given in Chapter 3. Title page figure:. Elementary surface features passing each other inside the EHD lubricated conjunction, see Fig. 6.5 for details.. ROUGH S URFACE E LASTOHYDRODYNAMIC L UBRICATION AND C ONTACT M ECHANICS c Andreas Almqvist (2004). This document is freely available at Copyright http://epubl.ltu.se/1402-1757/2004/35 or by contacting Andreas Almqvist, andreas.almqvist@ltu.se The document may be freely distributed in its original form including the current author’s name. None of the content may be changed or excluded without permissions from the author.. ISSN: ISRN:. 1402-1757 LTU-LIC--04/35--SE. This document was typeset in LATEX 2ε ..

(4) Abstract In the field of tribology, there are numerous theoretical models that may be described mathematically in the form of integro-differential systems of equations. Some of these systems of equations are sufficiently well posed to allow for numerical solutions to be carried out resulting in accurate predictions. This work has focused on the contact between rough surfaces with or without a separating lubricant film. The objective was to investigate how surface topography influences contact conditions. For this purpose two different numerical methods were developed and used. For the lubricated contact between rough surfaces the Reynolds equation were used as a basis. This equation is derived under the assumptions of thin fluid film and creeping flow. In highly loaded, lubricated, non-conformal contacts of surfaces after running-in, the load concentration no longer results in plastic deformations, however large elastic deformations will be apparent. It is the interaction between the hydrodynamic action of the lubricant and the elastic deformations of the surfaces that, in certain applications, enable the lubricant film to fully separate the surfaces. This is commonly referred to as full film elastohydrodynamic (EHD) lubrication. Typical machine elements that operate in the full film EHD lubrication (FL) regime include rolling element bearings, cams and gears. Unfortunately, a cost effective way of machining engineering surfaces seldom results in a surface topography that influence contact conditions in the same way as a surface after running-in. Such topographies may prevent the lubricant from fully separating the surfaces because of deteriorated hydrodynamic action. In this case the applied load is carried in part by the lubricant and in part by surface asperities and/or surface active lubricant additives. This could also be the case in lubricant starved contacts, which is a common situation in not only grease lubricated contacts but also in many liquid lubricated contacts, such as high speed operating rolling element bearings. The load sharing between the highly compressed lubricant and the surface and/or surface active lubricant additives is the reason why this lubrication regime is most commonly referred to as mixed EHD lubrication (ML). Machine elements that while running operate in the FL regime may experience a transition into the ML regime at stops or due to altered operating conditions. It is not possible to simulate direct contact between the surfaces using a numerical method based on Reynolds equation. A parameter study, of elementary surface features passing each other inside the EHD lubricated conjunction, was performed. The results obtained, even though no direct contact could be simulated, does indicate that a transition from the FL to the ML regime would occur for certain combinations of the varied parameters. At start-ups, the contact in a rolling element bearing could be both starved and drained from lubricant. In this case the hydrodynamic action becomes negligible in terms of load carrying capacity. The load is carried exclusively by surface asperities and/or surface aci.

(5) tive lubricant additives. This regime is referred to as boundary lubrication (BL). Operation conditions could also make both FL and ML impossible to achieve, for example, in the case in a low rpm operating rolling element bearing. The BL regime is in this work modeled as the unlubricated frictionless contact between rough surfaces, i.e., a dry contact approach. A variational principle was used in which the real area of contact and contact pressure distribution are those which minimize the total complementary energy. A linear elastic-perfectly plastic deformation model in which energy dissipation due to plastic deformation is accounted for was used. The dry contact method was applied to the contact between four different profiles and a plane. The variation in the real area of contact, the plasticity index and some surface roughness parameters due to applied load were investigated. The surface roughness parameters of the profiles differed significantly..

(6) Contents I. The Thesis. 1 Introduction 1.1 Elastohydrodynamic lubrication 1.2 Lubrication regimes . . . . . . . 1.3 Surface topography . . . . . . . 1.4 Objectives . . . . . . . . . . . . 1.5 Outline of this thesis . . . . . .. 1 . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. 3 3 4 5 5 7. 2 Full Film EHD Lubrication 2.1 Deterministic roughness models . . . . . . . . . 2.2 Governing equations . . . . . . . . . . . . . . . 2.3 The Block-Jacobi method . . . . . . . . . . . . . 2.3.1 Ordinary Jacobi . . . . . . . . . . . . . . 2.3.2 Coupled systems of equations . . . . . . 2.4 A brief overview of the multilevel technique . . . 2.4.1 Grid levels . . . . . . . . . . . . . . . . 2.4.2 Intergrid transfer operators . . . . . . . . 2.4.3 Two level solver for the Poisson equation 2.4.4 Two-dimensional functional operators . . 2.5 Dimensionless formulation . . . . . . . . . . . . 2.6 Discrete formulation . . . . . . . . . . . . . . . 2.7 Solution method . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. . . . . . . . . . . . . .. 9 9 11 15 15 16 17 17 18 19 20 21 23 24. 3 Dry Elasto-Plastic Contact 3.1 Statistical roughness models . . 3.2 Deterministic roughness models 3.3 Numerical solution techniques . 3.4 Governing equations . . . . . . 3.5 Spectral analysis . . . . . . . . 3.6 Dimensionless formulation . . . 3.7 Discrete formulation . . . . . . 3.8 Solution method . . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. 27 27 28 28 29 30 31 32 33. . . . . .. . . . . . . . .. . . . . .. . . . . . . . .. iii. . . . . .. . . . . . . . .. . . . . .. . . . . . . . .. . . . . .. . . . . . . . .. . . . . .. . . . . . . . .. . . . . .. . . . . . . . .. . . . . .. . . . . . . . .. . . . . . . . ..

(7) 4 Surface Characterization 4.1 Theoretical model . . . 4.2 Measured topographies 4.3 Parameter study . . . . 4.4 Conclusions . . . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. 35 35 35 36 39. 5 Reynolds vs. CFD 5.1 The CFD approach . . . . . . . . . 5.2 CFD - Governing equations . . . . . 5.3 The model problem . . . . . . . . . 5.3.1 Interpolation of solution data 5.3.2 Error estimation . . . . . . 5.4 The results of the comparison . . . . 5.5 Discussion and concluding remarks. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. . . . . . . .. 41 41 42 42 43 44 46 46. 6 Simulations of Rough FL 6.1 The different overtaking situations . . . . . . . . . . . . . . . . . . . . . 6.2 The Dent-Ridge overtaking . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 51 51 54 58. 7 Simulations of the dry contact 7.1 Varying the applied load . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 61 62 64. 8 Concluding Remarks. 65. 9 Future Work. 67. II Appended Papers. 69. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. A A.1 Introduction . . . . . . . . . . . . . . . . . . . . . . A.2 Governing equations . . . . . . . . . . . . . . . . . A.2.1 Boundary conditions and cavitation treatment A.3 Numerics . . . . . . . . . . . . . . . . . . . . . . . A.3.1 The numerics for the CFD approach . . . . . A.3.2 The numerics for the Reynolds approach . . A.3.3 Error estimation . . . . . . . . . . . . . . . A.3.4 Interpolation of solution data . . . . . . . . . A.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . A.6 Conclusions . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. 71 74 74 76 76 77 77 78 79 80 81 84. B.1 Introduction . . . B.2 Theory . . . . . . B.2.1 Equations B.2.2 Numerics. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. 87 90 91 91 92. B . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . .. . . . ..

(8) B.2.3 Error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C.1 C.2 C.3 C.4 C.5. Introduction . . . . . . . Theory . . . . . . . . . . Surface characterization . Results . . . . . . . . . . Conclusions . . . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. 93 94 98 103 106 107 108 109 112.

(9)

(10) Preface This licentiate thesis comprises the results from numerical simulations of both lubricated and unlubricated contacts, specifically on the influence of surface topography on contact conditions. The work has been carried out at the Department of Applied Physics and Mechanical Engineering, the Division of Machine Elements, Luleå University of Technology and the thesis is based on the three papers, A, B and C found in part II. Graduate studies are time consuming which means little time for other activities. This makes it hard work being the father of three wonderful children. My greatest gratitude is therefore given to my wife Ulrika who by being a such fantastic mother and taking care of our children has made these studies possible. For all those inspiring moments with my children I am also very grateful. My parents, grand parents and other close relatives and friends also deserve my gratitude for all their support. I would like to thank my supervisor Dr. Roland Larsson for introducing me to the subject of tribology and especially the subject of EHL and for the inspiring and encouraging meetings we have had and hopefully will continue to have. I would also like to thank all my colleagues at the Division of Machine Elements for their support, their co-operation and for introducing new concepts in tribology to me. I must also say that I am grateful to be a part of this particular division because of the friendly and relaxed atmosphere. Special thanks goes to my friend and former colleague Dr. Torbjörn Almqvist for his co-operation, enthusiasm and all the stimulating discussions we have had. Another special thanks goes to Mr. Fredrik Sahlin for his co-operation, and his support and help in scripting especially in Perl and LATEX 2ε . From the Department of Mathematics I would like to thank my associate supervisor Dr. Inge Söderqvist for his support in the field of Scientific Computing. I would also like to thank Mr. Reynold Näslund for stimulating and encouraging discussions and his help with mathematics in general. Finally, I would like to thank my sponsors; Fortum, Indexator, SKF Statoil, Volvo Car Corp.and the national research programme HiMeC and the national graduate school in scientific computning, NGSSC, both financed by the Swedish Foundation for Strategic Research (SSF).. vii.

(11)

(12) Nomenclature Pa −1. α. Pressure-viscosity coefficient. A. Functional space. µ¯. Dimensionless dynamic viscosity. ρ¯. Dimensionless density. ∆T. Dimensionless step size in time. s. ∆t. Step size in time. s. ∆X. Dimensionless step size in x. m. ∆x. Step size in spatial coordinate x. m. ∆z. Step size in spatial coordinate z. m. γ˙. Shear rate. E. Error. ε. Reynolds equation coefficient. s−1. εspatial Discretization error in space εtime. Discretization error in time. εx. Discretization error in space. εy. Discretization error in space. . Total potential complementary energy. I. Intergrid transfer operator. Λ. Wavelength. L. Functional operator. m. m s−1. u (x, y, z) Velocity field x (x, y, z) Spatial coordinates. m ix.

(13) µ. Dynamic viscosity. Pa s. µ0. Dynamic viscosity at ambient pressure. Pa s. ν. Poissont’s ratio. Ω. Integration domain 2D / 3D contact. ωx. Spatial frequency. 1/m. ωy. Spatial frequency. 1/m. O (n). Mathematical order of the number n. Φ. Solution variable. ψi. Topography of surface i. ψavg. Mean surface height. ρ. Density. kg m−3. ρ0. Density at ambient pressure. kg m −3. Σ. Dimensionless stress function. σ. Total stress tensor (A). Pa. τ. Lubricant shear stress. Pa. τ0. Eyring stress. Pa. τ1. Lubricant shear stress at surface 1. Pa. τm. Lubricant midplane shear stress. Pa. τxz. Lubricant shear stress. Pa. f. Fourier transformation of f. 1 ∇. Dimensionless reduced wavelength. ξ. Film thinning effect measure. A. Dent / Ridge Amplitude. ai. Matrix coefficient. aξ. Contact semi-width in the rolling direction. Ai j. Matrices. b. Hertzian half-width b =. bi. Matrix coefficient. bξ. Contact semi-width in the transversial direction. m / m2. m. m. (8wRx ) / (πE ). m. m. m.

(14) C1. Constant = 5.9 10 8. C2. Constant = 1.34. ci. Matrix coefficient. d. Elastic deformation. m. de. Elastic deformation. m. di. Matrix coefficient. dp. Plastic deformation. m. E. Modulus of elasticity. Pa. E. Effective mod. of elast. 2/E = (1 − ν21 )/E1 + (1 − ν22)/E2. Pa. F. Discrete right hand side (Reynolds equation). fi. Matrix coefficient. G. Discrete right hand side (Film th. equation). G. Shear modulus of elasticity of the lubricant. g. Determinant of the metric tensor. gi. Matrix coefficient. gs. The gap between the undeformed surfaces. H. Dimensionless film thickness. h. Film thickness. m. h0. Integration constant. m. Hs. Hardness of the softer material. Pa. H00. Dimensionless integration constant. hc. Central film thickness. m. hmin. Minimum film thickness. m. K. Elastic deformation integral kernel. L. Moes dimensionless speed parameter. li. Amplitude variation of dent/ridge. Lx. Length of the Fourier window in the x-dir.. m. Ly. Length of the Fourier window in the y-dir.. m. L1. Discrete operator (Reynolds equation). Pa. m.

(15) L2. Discrete operator (Film thickness equation). M. Mass flux per unit width. M. Moes dimensionless load parameter. mi. Wavelength variation of dent/ridge. P. Dimensionless pressure. p. Pressure. Pa. P0 ph. Constant in the viscosity expression Hertzian pressure p h = (wE ) / (2πRx ). Pa. pmax. Maximum pressure. Pa. Ra. Average roughness. µm. Rd. Radius of the lower surface (A). m. Ri. Radius of the surface i. m. Rk. Kurtosis. Rq. Root mean square (RMS) roughness. Ru. Radius of the upper surface (A). m. Rx. Reduced radius of curvature in x-dir. 1/R x = 1/R11 + 1/R22. m. Rz. Average maximum height. Rsk. Skewness. S. Non-Newtionian slip factor. s. Slide-to-roll ratio 2(u 1 − u2)/(u1 + u2 ). T. Dimensionless time. t. Time. u. Solution to the one-dimensional Poisson equation. u1. Surface velocity. m s −1. u2. Surface velocity. m s −1. us. Sum of velocities u s = u1 + u2. m s−1. W. Applied load, 2D / 3D contact. Pa m −1 / Pa. w. Applied load. X. Dimensionless spatial coordinate. kg m −1 s−1. µm. µm. s. Pa m −1 m.

(16) x. Spatial coordinate in the rolling direction. m. xc. Centre of dent / ridge x c (t) = xs − u t. m. xs. Initial placement of dent / ridge. m. z. Spatial coordinate across the film. m. z. Spatial coordinate across the film. m. z1. Spatial coordinate across the film. m. zvisc. Pressure-viscosity index.

(17)

(18) Part I. The Thesis. 1.

(19)

(20) Chapter 1. Introduction Machines consists of machine elements and their safe and efficient operation relies on carefully designed interfaces between the machine elements. The functional design of interfaces covers geometry, materials, lubrication and surface topography. Incorrect choice of these design parameters may lead to both lowered efficiency and shortened service life. Poor geometry could lead to large stress concentrations which in turn may lead to rapid fatigue. Large stress concentrations also implicitly imply a temperature rise because of energy dissipation due to plastic deformations. The choice of materials is also of great importance, for example, electrolytic corrosion may drastically reduce service life. Contact fatigue contact to low ductility would not only lower the service life but could lead to third body abrasion due to spalling which in turn could end up lowering the service life of other components. A lubricant serves several crucial objectives. When its main objective is to lower friction, the action of additives is of concern. If the interface is subjected to excessive wear, the lubricant’s ability to from a separating film is even more crucial. In this case the bulk properties of the lubricant have to be carefully chosen. All real surfaces are rough at some scale and the surface topography certainly influences the contact between machine elements. The design parameters mentioned are also mutually dependent. That is, they affect the way the others influence the operation of the system. For example, a change in geometry could require to another choice of materials. The new choice of material may change the objectives of the lubricant and force the operation into another lubrication regime. All these four design parameters are of great importance, however in this thesis work it is surface topography which is the main focus.. 1.1 Elastohydrodynamic lubrication Elastohydrodynamic (EHD) lubrication is the type of hydrodynamic lubrication in which the elastic deformations of the contacting surfaces cannot be neglected. This is often the case in non-conformal (concentrated) contacts. For example, the contact between the roller and the raceway in a typical roller bearing as shown in Fig. 1.1 will operate in the elastohydrodynamic lubrication regime. The actual contact zone for a rolling bearing is, in general, elliptic in shape. Depending on the design parameters previously mentioned and the actual running conditions, the shape of the ellipse will change as shown in Fig. 1.2. 3.

(21) CHAPTER 1. INTRODUCTION. 4. Figure 1.1: A typical rolling element bearing In any case, the contact region is small and the concentrated load implies high pressure which will lead to large elastic deformations and could also lead to plastic deformation. For a bearing in operation, it is the large elastic deformations that causes fatigue which in turn can lead to shortened service life due to for example spalling. When the contact is starved of lubricant, or when running conditions do not allow for a hydrodynamic action that fully separates the surfaces, the risk for plastic deformation increases. x. y. Roller width. Figure 1.2: Elliptic contacts (x/aξ )2 + (y/bξ)2 = 1 If bξ exceeds the minimum width of the raceway and the roller, then the contact will be truncated which could lead to increased stresses in the material; at least for the unlubricated contact. In the case of a contact with b ξ /aξ > 4 but with bξ still less than the minimum width of the raceway and the roller, the centerline in the rolling direction can be approximated to a line contact, Evans et al. [1].. 1.2 Lubrication regimes The lubrication regimes in EHD lubrication are commonly divided into: Boundary Lubrication (BL), Mixed Lubrication (ML) and Full Film Lubrication (FL). In the BL regime the lubricant’s hydrodynamic action is negligible and the load is carried directly by surface asperities or by surface active additives. In the ML regime the load is carried by the.

(22) 1.3. SURFACE TOPOGRAPHY. 5. lubricant’s hydrodynamical action, the surface active additives and/or directly by surface asperities. When the hydrodynamic action of the lubricant fully separates the surfaces and the load is carried totally by the lubricant film the contact enters the FL regime. In terms of traction the modified Stribeck curve shown in Fig. 1.3 gives a good visualization of the different regimes.. ML. FL. ξ. BL. ln (µu/p). Figure 1.3: A modified Stribeck curve for EHD lubrication. 1.3 Surface topography Depending on the lubrication regime, the surface topography will have a different influence on operation. In the BL regime the topography is preferably chosen so that friction is reduced but without increasing the rate of wear. In the ML regime the objective of the topography is to support the hydrodynamic action of the lubricant, to enable bonding of the surface active additives and to minimize friction in the contact spots without increasing wear. In the FL regime, traction may be reduced by carefully chosen topographies and, even though there is no direct contact, the topography must also prevent fatigue which lead to excessive wear in from of spalling.. 1.4 Objectives This work concerns the contact between rough surfaces; with or without a separating lubricant in between. The main objective was to investigate how the surface topography influences different contact conditions. The long term goal is to develop engineering tools, to predict different contact conditions, that takes into account the surface topography. For this purpose two different numerical methods were developed. A method of surface characterization was also developed, to be used to facilitate the interpretation of the data acquired when using deterministic roughness models. The FL regime was initially set in focus and a theoretical model based on Reynolds equation for line contacts was considered. The model was to account for transient effects,.

(23) CHAPTER 1. INTRODUCTION. 6 i.e., a time dependent model, and include: • Linear elastic surface materials • Surface roughness, measured or modeled • Non-Newtonian rheology. The Reynolds equation derived by Conry et al. [2] was chosen to handle the non-Newtonian rheology. A coupled solver was developed in which pressure and film thickness are solved simultaneously. The numerical method uses a multilevel technique to facilitate the solution process. A schematic picture of the theoretical model is given in Fig. 1.4.. Load R1 E1 ,ν1 µ0 ,α. hy. rap. g po. To. E2 ,ν2 R2. Figure 1.4: A schematic picture of the theoretical FL model used in this work The objective set for this thesis was to develop, verify and use this tool for systematic investigations of the local effects caused by elementary surface features passing each other inside the lubricated conjunction. To complement this work, a contact mechanics model was considered in order to be able to investigate the influence of surface topography for contacts operating in the BL regime. The theoretical model is schematically depicted in Fig. 1.5. The material model used was initially a linear elastic one. However, to be able to simulate running-in the material model was extended to allow for linear elastic-perfectly plastic materials. The specification of this model can be summarized as follows:.

(24) 1.5. OUTLINE OF THIS THESIS. 7. Load. H1. R1 E1 ,ν1. y. ph. gra. po To. E2 ,ν2. H2 R2. Figure 1.5: A schematic picture of the theoretical BL model used in this work. • 2D and 3D contacts • linear elastic-perfectly plastic surface material • Surface roughness, measured or modeled A model accounting for the energy dissipation due to plastic deformation was also considered, see Tian and Bhushan [3]. The numerical solution technique is, however, based on the model of Stanley and Kato [4] with some modification necessary to account for plastic deformations. The BL regime is here modeled as the unlubricated, contact mechanics problem, accounting for plastic deformations, and include the possibility of simulating the running-in process of real surfaces.. 1.5 Outline of this thesis The thesis is based on the three papers A, B and C found in part II. In Part I the content of these papers is extended with more a complete theoretical description and derivations and the inclusion of material not addressed in the papers. Part I is written in the form of a monograph..

(25) 8. CHAPTER 1. INTRODUCTION. Chapter 2, of this thesis introduces the topic of FL including the theory of the Reynolds approach used in this work. It includes a shortened derivation of the non-Newtonian Reynolds equation according to Conry et al. [2]. The full non-linear integro-differential system of equations including expressions for compressibility, viscosity, elastic deflection, film thickness and force balance are also described. The numerical ‘Block-Jacobi’ method used to solve the coupled system consisting of the Reynolds equation and the film thickness equation may also be found here. A multilevel technique used to facilitate the numerical solution process of the one dimensional Poisson problem and a possible modification to facilitate the numerical solution of coupled systems are then briefly described. The line contact problem is restated in dimensionless form and the discrete analog is derived. Finally, a description of the numerical solution process is given. The theory that constitutes the basis for the dry elasto-plastic contact method is found in Chapter 3 together with an introduction to the contact mechanics problem in tribology. Different models of the elastic contact are then discussed in connection to the present model. A description of how spectral analysis can be used in order to determine the elastic deformation integral is then given. The dimensionless and discrete formulation and the numerical solution process, for the dry elasto-plastic method, may be also found in this chapter. The interpretation of the data acquired when using a deterministic model, simulating either FL or the dry contact, may be facilitated by making use of the surface characterization technique developed in Chapter 4. This technique is based on truncated Fourier series and is applied to a sample topography. In Chapter 5, the Reynolds based approach described in Chapter 2 is compared with a Computational Fluid Dynamics (CFD) based approach. This comparison consist of simulations of one stationary and two different transient problems. The results are encouraging from several viewpoints: verification of the codes, the possibilities to further develop the CFD approach given by Almqvist and Larsson [5], and the justification of using a Reynolds approach under the running conditions chosen. In this comparison, the roughness has the form of a single ridge passing through the EHD lubricated conjunction. The topic of roughness in the FL regime more thoroughly addressed in Chapter 6. In a sliding contact, where the surface topography of both surfaces is, a continuously changing effective surface roughness occurs. For a line contact, the surface topography is superimposed by elementary surface features, i.e., dents and ridges. When passing each other inside the conjunction, these surface features cause local effects. In this chapter the influence of these local effects on the film thickness and the pressure is the subject of investigation. In Chapter 7 the method for modeling the dry elasto-plastic contact between rough surfaces is used to simulate the contact conditions of four different topographies. The simulations are restricted to two-dimensions and profiles taken from the measurements of four surfaces form the basis of the load parameter study performed. The content of the second part of this thesis, part II, is simply the three published papers that form the basis of the thesis. These are as published except for some additional corrections and proofreading..

(26) Chapter 2. Full Film EHD Lubrication Due to increasing demands on performance the lubricant film thickness in EHD lubricated contacts has decreased over the years, and will probably continue to do so. However, the surface machining processes are financially constrained which slows down the decrease in surface roughness and sometimes even causes an increase in commercially available machine elements. As a result, the ratio of film thickness to surface roughness will continue to decrease, which in turn will affect the performance of the EHD lubricated contact. Performance increases must also be matched with improved reliability hence a detailed understanding of EHD lubrication of rough surfaces deserves attention. In this chapter, the theoretical basis of the EHD lubrication model is described. First the adopted deterministic roughness model is described along with previous work done using such a model. A derivation of the one-dimensional non-Newtonian Reynolds-Eyring equation according to Conry et al. [2] is also given. The integro-differential problem for the isothermal EHD lubricated line contact consist of the Reynolds equation derived, an equation for the film thickness including the elasticity of the rollers, an equation for force balance and two semi-empirical equations for the fluid density and the fluid viscosity. These equations are restated in dimensionless form and then discretized. The coupled solution method, which allows for simultaneous solution of pressure and film thickness, referred to as the ‘Block-Jacobi’ method, is explained. A brief overview of the multilevel technique used to facilitate the solution method is also given.. 2.1 Deterministic roughness models There are many ways to model surface topographies. Mathematically, the models are either statistic or deterministic. Throughout this work the latter approach is used. In many theoretical studies of EHD lubrication, using deterministic models, one surface is considered smooth and the other as being rough. This is a suitable approximation when modeling a rolling contact or a contact in which the roughness of one of the surfaces is of minor importance. Because of the assumption of parabolic surfaces (Hertzian theory) it is possible to simply sum the roughness of the two contact surfaces and model the contact of the effective roughness and a perfectly smooth surface. This approach will be referred to as ‘one-sided roughness’. 9.

(27) CHAPTER 2. FULL FILM EHD LUBRICATION. 10. Many extensive, systematic studies have been carried out using the one-sided roughness assumption. Based on the knowledge that the surface and hence roughness deform inside the EHD lubricated conjunction, Lubrecht et al. [6] investigated how sinusoidal Fourier components deform as a function of wavelength and the contact operating conditions, including slide-to-roll ratio. They combined the results from the numerical simulations of the Newtonian line contact in a single formula,. where. Ad 1 = ˜ ˜2 Ai 1 + 0.125∇1 + 0.04∇ 1. (2.1). 3/4 1 = Λ M √1 , ∇ b L1/2 s. (2.2). that describes the relationship between the reduced amplitude A d /Ai and the initial wavelength Λ, the Hertzian half width b, Moes parameters M and L and the slide-to-roll ratio s. The relation described by Eq. (2.1) is visualized in Fig. 2.1. According to Eq. (2.2), 1. Ad /Ai. 0.8. 0.6. 0.4. 0.2. 0 0.1. 1. 10. 100. ˜1 ∇. 1 Figure 2.1: A d /Ai as a function of ∇ a shorter wavelength (Λ) or larger slide-to-roll ratio (s) will lead to a shorter reduced di 1 and thus a less reduced amplitude A d /Ai . In such cases the mensionless wavelength ∇ risk of film breakdown is high. The part of Eq. (2.2) describing the influence of: • applied load (w), • surface mean velocity ((u 1 + u2 )/2), • reduced radius of curvature of the contacting bodies (R x ), • the effective modulus of elasticity (E ), • bulk viscosity (µ 0 ), • the pressure-viscosity coefficient (α).

(28) 2.2. GOVERNING EQUATIONS. 11. concealed in M, L and b will also affect the reduction in roughness amplitude. If the 1 is expressed in the physical input parameters above: parameter ∇ 1 = ∇. pi1/2 Λ w1/4 3/4. 2 Rx. 1/2. E 1/4 µ0. 1/2. us. α1/2 s1/2. ,. (2.3). it can be seen that a decrease of w, or an increase in any of the parameters u s , Rx , E , µ0 , α will lead to a decreased reduction in roughness amplitude. In an EHD lubricated contact, that is subjected to sliding, the surface topographies of both the surfaces are of significance and a continuously changing effective surface roughness occurs. A number of investigations have been carried out on this subject, for example Evans et al., Tau et al., Chang, Venner and Morales-Espejel [1, 7, 8, 9]. The one-sided roughness approach is not valid in such situations and a ‘two-sided’ roughness treatment is needed. The topic of two-sided roughness is addressed in this work and results of an extensive parameter study may be found in Chapter 6.. 2.2 Governing equations A modified Reynolds equation, based on on the Eyring theory of non-Newtonian flow, is derived in one dimension. Johnson and Tevaarwerk [10] proposed the nonlinear constitutive equation for a lubricant under isothermal conditions given by Eq. (2.4). τ 1 dτ τ0 ˙γ = + sinh , (2.4) G dt µ τ0 where G is the shear modulus of elasticity of the fluid, τ is fluid shear stress, µ is the dynamic viscosity, and τ 0 is the Eyring shear stress. This is a Maxwell rheological model where the total shear strain rate is the sum of an elastic term and a nonlinear viscous term based on the Eyring’s theory of viscosity. The modified Reynolds equation is derived from on the Eyring equation, the nonlinear viscous portion of Eq. (2.4), and under the assumption of plain strain rate. Fig. 2.2 shows a fluid element in the thin lubricating film between two solids. (It should be noted that the velocity of the solids at the surface is approximated by the velocity u in the x-direction.) In this case the equation of equilibrium is the x-direction takes the form: ∂p ∂τxz = , (2.5) ∂z ∂x where p is fluid pressure. If the film thickness is denoted by h (x), then according to Fig. 2.2, 0 ≤ z ≤ h (x) and z = z − z1 . Assuming p = p (x), integration of Eq. (2.5) with respect to z yields: dp (2.6) τxz = τ1 + z , dx where τ1 is the shear stress at surface 1. Substituting Eq. (2.6) into the constitutive equation (˙γ = (τ0 /µ) sinh(τxz /τ0 )) yields: . τ1 + z dp du τ0 dx = sinh . (2.7) γ˙ = dz µ τ0.

(29) CHAPTER 2. FULL FILM EHD LUBRICATION. 12 z. u2. z. z. z1 u1 x. x. Figure 2.2: A description of the line contact region Assuming the velocity of the lower surface (z = 0) is u1 and that the viscosity does not vary across the film ,i.e., µ = µ (x), integration of Eq. (2.7) with respect to z gives the following expression for the velocity profile in the x-direction: . z τ1 + s dp τ0 dx sinh ds, (2.8) u x, z = u1 + τ0 0 µ which after evaluation of the integral expression becomes: . τ1 + z dp τ20 τ1 dx − cosh . u x, z = u1 + dp cosh τ0 τ0 µ. (2.9). dx. Introduction of the midplane shear stress: τm = τ1 +. h dp 2 dx. (2.10). and a dimensionless function, defined as: Σ=. h dp 2τ0 dx. makes it possible to rewrite Eq. (2.9) as: τ0 h 1 2z τm τm u x, z = u1 + −Σ 1− −Σ . cosh − cosh 2µ Σ τ0 h τ0. (2.11). (2.12). Application of the boundary condition u (x, h) = u 2 and utilizing hyperbolic relations gives: τm τm τ0 h u2 = u1 + + Σ − cosh −Σ cosh 2µΣ τ0 τ0 (2.13) τm τ0 h = u1 + sinh sinh Σ. µΣ τ0.

(30) 2.2. GOVERNING EQUATIONS. 13. Eq. (2.13) can be rearranged to allow for determination of the midplane shear stress as follows: τm µ (u2 − u1) Σ . (2.14) sinh = τ0 τ0 h sinh Σ The mass flux per unit width, M (x), is defined as: M (x) =. h 0. ρu x, z dz .. (2.15). After substitution of Eq. (2.12) into Eq. (2.15) and integration, the expression for mass flux per unit width becomes: τm τm τ0 h h + Σ − sinh −Σ − M (x) = ρu1 h + ρ sinh 2µΣ 2Σ τ0 τ0 (2.16)

(31) τm −Σ . h cosh τ0 Expanding and rearranging the terms in the brackets of Eq. (2.16) gives . τm τ0h2 sinh Σ M (x) = ρu1 h + ρ sinh + 2µ Σ τ0

(32) τm sinhΣ − Σ coshΣ cosh . Σ2 τ0. (2.17). Substitution of Eq. (2.14) into Eq. (2.17) yields: M (x) =. =. (u1 + u2 ) ρh + 2 ρτ0 h2 sinhΣ − Σ coshΣ τm cosh 2 2µ Σ τ0 (u1 + u2 ) ρh + 2. =. . . (2.18). τm ρτ0 sinhΣ − Σ coshΣ Σ cosh 2µ Σ3 τ0 (u1 + u2 ) ρh + 2

(33). τm dp ρh3 3 (sinhΣ − Σ coshΣ) cosh 12µ Σ3 τ0 dx h2. The dimensionless quantity inside the brackets of Eq. (2.18) is the non-Newtonian slip factor, which after using the hyperbolic relation cosh 2 (x) − sinh2 (x) = 1, and substitution of Eq. (2.14) is given by:. 2 µ (u2 − u1) Σ 3 (sinh Σ − Σ coshΣ) 1+ . (2.19) S (x) = Σ3 τ0 h sinh Σ.

(34) CHAPTER 2. FULL FILM EHD LUBRICATION. 14. The equation of continuity takes the following form: d dx. h(x) 0. ρu x, z dz = 0.. (2.20). or in terms of the mass flux per unit width: d M (x) = 0, dx. (2.21). thus requiring constant mass flux along the x-direction. Substitution of Eq. (2.18) into Eq. (2.21) using the expression of the non-Newtonian slip factor given by Eq. (2.19) finally yields the stationary one-dimensional ReynoldsEyring equation: d ρh3 dp (u1 + u2) d S (ρh) . (2.22) = dx 12µ dx 2 dx The representative lubricant stress τ 0 , characterizes the transition from Newtonian to non-Newtonian fluid behavior. When using the Eyring model an infinitely large τ 0 characterizes a Newtonian fluid, and by using L’Hospital rule it can be shown that S approaches unity. In this case Eq. (2.22) reduces to the conventional Reynolds equation. The value of S is always equal to or greater than one, which means that the effective viscosity µ/S is always less or equal to µ. It is possible to restate Eq. (2.22) in transient form, incorporating squeeze effects. The one-dimensional transient Reynolds-Eyring equation reads: ∂ ∂ ρh3 ∂p (u1 + u2) ∂ S (ρh) + (ρh). (2.23) = ∂x 12µ ∂x 2 ∂x ∂t Eq. (2.23) combined with an equation for the film thickness, an equation of force balance, and two semi-empiric equations for the density and the viscosity respectively, constitutes the basis of the EHD lubricated line contact problem studied in this work. The equation for film thickness is given by: h (x,t) = h0 (t) +. x2 + d (x,t) + ψ1 (x,t) + ψ2 (x,t) , 2Rx. (2.24). where h0 is an integration constant, R x is the reduced radius of curvature in the x-direction given by 1/R x = 1/Rx1 + 1/Rx2 , d is the elastic deformation of the contacting solids and ψ1 and ψ2 the topography of surfaces 1 and 2 respectively. The elastic deformation in one-dimension is given by: d (x,t) = −. 4 πE . ∞ −∞. ln |x − x|p x ,t dx ,. (2.25). where E is the effective modulus of elasticity given by 2/E = (1 − ν21 )/E1 + (1 − ν22 )/E2 . There are several different ways to model viscosity. In this thesis, the Roelands expression is used throughout and is given by: p zvisc αP0 −1 + 1 + . (2.26) µ (p) = µ0 exp zvisc P0.

(35) 2.3. THE BLOCK-JACOBI METHOD. 15. where α, P0 and zvisc are mutually dependent: zvisc (ln (µ0 ) + 9.67), α. P0 =. (2.27). The Dowson and Higginson equation is the semi-empiric expression used to model the density C1 + C2 p ρ (p) = ρ0 . (2.28) C1 + p At all times, the force balance condition (2.29). I. p (x,t) dx = w (t). (2.29). determines the integration constant h 0 (t). Also, the cavitation condition (2.30) is used to ensure that all negative pressures obtained during the solution process are removed. p ≥ 0.. (2.30). 2.3 The Block-Jacobi method In this section a generalized form of the Jacobi method, taking into account the coupling of the integro-differential system consisting of the Reynolds- and the film thickness- equation, is derived. To increase comprehension, the ordinary Jacobi method is revisited first.. 2.3.1 Ordinary Jacobi Let L (x) be general non-linear discrete operator L : R n → Rn . One way to solve the problem: ⎤ ⎡ ⎡ ⎤ 0 L1 (x) ⎥ ⎢ .. ⎥ ⎢ L (x) = 0 ⇔ ⎣ ... (2.31) ⎦ = ⎣ . ⎦, Ln (x). 0. is to apply the Newton-Raphson algorithm. In this case the determining system for the changes, δx, is ⎤⎡ ⎡ ∂L ⎤ ⎤ ⎡ 1 1 · · · ∂L δx1 L1 (x) ∂x1 ∂xn ⎥ ⎢ . ⎥ ⎥ ⎢ . . . . .. ⎥⎢ (2.32) Jδx = −L (x) ⇔ ⎢ ⎦ = − ⎣ .. ⎦. . . ⎦ ⎣ .. ⎣ .. ∂Ln ∂Ln δxn Ln (x) ··· ∂x1. ∂xn. Unfortunately, the need to compute the full Jacobian matrix J of the operator L and of solving the full linear system for highly resolved problems leads to an algorithm of high complexity. In the case of rough EHD lubrication, high resolution is needed and the Newton-Raphson algorithm is therefore not suitable. However, if the operator is either linear by nature or in a linearized form, i.e., L (x) ≡ Ax − b. (2.33).

(36) CHAPTER 2. FULL FILM EHD LUBRICATION. 16. where A is a n × n matrix and x and b are n × 1 matrices, then the ordinary Jacobi method could be applied. If ⎤ ⎡ a11 · · · a1n ⎥ ⎢ . . .. A = ⎣ ... (2.34) ⎦ . . an1. ···. ⎡. and. a11 ⎢ .. diag(A) = ⎣ . 0. ann ··· .. . ···. ⎤. 0 .. .. ⎥ ⎦,. (2.35). ann. then the changes δx are obtained by solving the following diagonal matrix system diag(A) δx = −L (x) .. (2.36). Since only the diagonal is used, the complexity of the Jacobi method is far less than that of the Newton-Raphson method.. 2.3.2 Coupled systems of equations Let us assume that a non-linear operator L is the discrete representation of a coupled system of two analytical equations, i.e.,.

(37) L1 (v1 , v2 ) , (2.37) L (v1 , v2 ) ≡ L2 (v1 , v2 ). where. L1 (v1 , v2 ) L2 (v1 , v2 ).

(38). =. f1 f2.

(39) .. (2.38). and (v1 , v2 ) represents the two dependent variables of the system. If this system is linearized at a point ( vˇ 1 , vˇ 2 ) it is possible to write it in general matrix form as:

(40)

(41).

(42) f A11 (ˇv1 , vˇ 2 ) A12 (ˇv1 , vˇ 2 ) v1 L (v1 , v2 ) ≡ − 1 . (2.39) A21 (ˇv1 , vˇ 2 ) A22 (ˇv1 , vˇ 2 ) v2 f2 If the ordinary Jacobi method is applied to the system L (v 1 , v2 ) = 0, the solution would be obtained by solving Eq. (2.36) with

(43). v1 , (2.40) x= v2. A=. A11 (ˇv1 , vˇ 2 ) A12 (ˇv1 , vˇ 2 ) A21 (ˇv1 , vˇ 2 ) A22 (ˇv1 , vˇ 2 ). and b=. f1 f2.

(44) (2.41).

(45) .. (2.42). In this case the coupling between the two equations would not be properly addressed and the solution method would most probably not converge..

(46) 2.4. A BRIEF OVERVIEW OF THE MULTILEVEL TECHNIQUE. 17. One way to include the coupling between the equations in the solution process is to improve the Jacobi method by introducing a slightly more expensive step δx = [δv 1 , δv2 ]T . Consider the system Eq. (2.39), where A l m are n × n matrices, and v l and fl are n × 1 matrices. A step that address the coupling may be solved from the system.

(47)

(48)

(49) J˜ 11 (ˇv1 , vˇ 2 ) J˜ 12 (ˇv1 , vˇ 2 ) L1 δv1 − , (2.43) δv2 L2 J˜ 21 (ˇv1 , vˇ 2 ) J˜ 22 (ˇv1 , vˇ 2 ) where J˜ l m = diag(Al m ) , l = 1, 2 m = 1, 2. Because of the block matrix structure, this method will be referred to as the ‘Block-Jacobi’ method. This makes a simple and not so expensive solution of Eq. (2.43) easy to implement. However, the most important feature of the the ‘Block-Jacobi’ method is that the coupling between the equations is considered.. 2.4 A brief overview of the multilevel technique Numerical linear system solvers like the Gauss-Seidel relaxation process and different types of Jacobi iterative processes are often not affordable due to their relatively slow rate of convergence. A multilevel technique applied to solve differential equations is based on the use of a conventional smoother, such as the Gauss-Seidel or Jacobi methods. The smoother itself is conventionally used to solve the problem with specified tolerance and resolution, i.e. using one discrete set of points referred to as a grid for the relaxation process. A decrease in the computational time for such a smoother can be achieved by either accepting a larger error or a lower resolution i.e. a coarser grid. Using multilevel methods it is actually possible to obtain a solution with both specified error and resolution without carrying out all the computational work on the fine grid level corresponding to the requested resolution. Instead the computations are restricted and performed on coarser grids and then prolonged back to the specified resolution. A brief overview of the multilevel technique used in the numerical solution is presented here. Multilevel techniques both for solving differential equations and to facilitate the solution procedure for different integral equations are described in detail by Venner and Lubrecht [11]. The brief description given here is based on a model problem consisting of the one-dimensional Poisson problem given in Eq. (2.46). A two level solution procedure is visualized for the one-dimensional functional operator defined for the Poisson problem. A description of how to modify this solution procedure for a two-dimensional functional operator in order to solve the coupled problem consisting of the Reynolds equation and the film thickness equation is then given.. 2.4.1 Grid levels A grid level is defined as a grid with number of nodes determined by the particular level. In this thesis work, the number of grid nodes n ∆x k on the grid level k with mesh size ∆x are obtained as follows: k n∆x (2.44) k = 2 + 1..

(50) CHAPTER 2. FULL FILM EHD LUBRICATION. 18 The mesh size ∆x may be chosen as:. ∆x = (b − a)/ n∆x k −1. (2.45). where the solution domain is given by the interval [a, b]. This means that if the grid level k corresponds to the mesh size ∆x then the grid level k − 1 corresponds to the mesh size 2∆x, i.e., a coarser grid.. 2.4.2 Intergrid transfer operators The intergrid transfer operators transfers the functions between different grid levels. The transfer to a coarser grid level is referred to as restriction and the transfer to a finer grid level is referred to as prolongation. These operators can be defined as follows: Definition 1 (Restriction operator) Let def .. ∆x v∆x i = v(xi ),. x∆x i ∈ [a = x1 , x1 + ∆x, . . . , x1 + i∆x, . . . , xn = b], ∆x = (b − a)/ n∆x k −1 , where n∆x k is the number of nodes and ∆x is the corresponding mesh size. Then the restriction operator I2∆x ∆x transfers discrete functions v∆x to the next coarser grid with mesh size 2∆x =. b−a n2∆x k−1 −1. cording to ∆x v2∆x = I2∆x ∆x v. Definition 2 (Prolongation operator) Let def .. v2∆x = v(x2∆x i i ), ∈ [a = x1 , x1 + 2∆x, . . ., x1 + i2∆x, . . ., xn = b] , x2∆x i 2∆x =. b−a n2∆x k−1 − 1. ,. where n2∆x k−1 is the number of nodes and 2∆x is the corresponding mesh size. Then the prolongation operator I∆x 2∆x transfers discrete functions v2∆x to the next finer grid with mesh size ∆x according to 2∆x v∆x = I∆x 2∆x v. ac-.

(51) 2.4. A BRIEF OVERVIEW OF THE MULTILEVEL TECHNIQUE. 19. 2.4.3 Two level solver for the Poisson equation The one-dimensional Poisson problem is given by Eq. (2.46). d2 v = f (x) , v (0) = 0, v (1) = 0. dx2. (2.46). It is possible to restate the problem in operator form, i.e., L (v). =. d2 v , dx2. =. f (x) , v (0) = 0, v (1) = 0.. def .. (2.47) L (v). Using second order finite differences, the problem defined in Eq. (2.47) is discretized into: Li (v). =. vi−1 − 2vi + vi+1 (∆x)2. , (2.48). Li (v). =. fi , i ∈ {2, . . . , n − 1}, v1 = 0, vn = 0.. The iterative solution process of this problem can be facilitated by a multilevel technique based on only two grid levels. For this purpose, the Gauss-Seidel relaxation process is assumed to be the smoother. Based on an initial guess on the finest grid level (v ∆x ) the solution can be obtained according to: • Use a few iterations of the Gauss-Seidel method to obtain an approximate solution v˜i ∆x and then compute the residuals r i∆x = Li∆x v˜∆x − fi∆x , where Li∆x (v) =. ∆x ∆x v∆x i−1 − 2vi + vi+1. (∆x)2. and f i∆x is equal to the discrete right hand side function f in Eq. (2.46). 2∆x = ∆x • Restrict the solution and the residuals to the coarse grid, v 2∆x = I2∆x I ∆x v˜i , and rI 2∆x ∆x I∆x ri , where I is the coarse grid index. • Calculate the coarse left hand side L i2∆x v2∆x of the Poisson equation to obtain the corresponding coarse right hand side, f I2∆x = LI2∆x v2∆x − rI2∆x .. • Use a few iterations the Gauss-Seidel method to obtain the approximate solution, v˜I 2∆x , to the system LI2∆x v2∆x = fI2∆x . • Calculate the coarse grid error, E I2∆x = v˜I 2∆x − VI , where VI = v˜∆x 2I . 2∆x and correct the solution, vˆ ∆x = v˜ ∆x + • Prolong the coarse grid error, E i∆x = I∆x i i 2∆x EI ∆x Ei . ∆x • Use a few Gauss-Seidel iterations to obtain the approximate solution, vˆi , on the fine grid..

(52) 20. CHAPTER 2. FULL FILM EHD LUBRICATION. Note the FAS right hand side, given by f I2∆x = LI2∆x v2∆x − rI2∆x on levels j where def . j < k and k is the finest level, is not equal to the discretized right hand side f I2∆x = f x2∆x i of Eq. (2.48). This two level multilevel technique can be extended to a more general multilevel technique, i.e., the multilevel V-cycle containing m levels. That is, the restriction process continues until the coarsest level is reached and then the prolongation process corrects and transfers the solution back to the finest level. (Venner and Lubrecht [11]).. 2.4.4 Two-dimensional functional operators The coupled integro-differential problem that consists of the Reynolds equation given by Eq. (2.23), and the film thickness equation given by Eq. (2.24) viz. ∂ ρh3 ∂p ∂ (u1 + u2) ∂ S (ρh) − (ρh) = 0, − ∂x 12µ ∂x 2 ∂x ∂t (2.49) 2 x h (x,t) − h0 (t) − d (x,t) = + ψ1 (x,t) + ψ2 (x,t) , 2Rx is two-dimensional in the functional space A 2 . For example, the two coordinates x and y can be thought of as one point (tuple) in R 2 , i.e, (x1 , x2 ) ∈ R2 . In the same way the pressure p and the film thickness h is a tuple in A 2 , i.e, (p, h) ∈ A 2 . This means that if Eq. (2.49) is restated in operator form, the operators map functions from A 2 to R, i.e., L : A 2 → R. In comparison, the one-dimensional functional operator previously defined for the onedimensional linear Poisson problem map functions from A to R, i.e., L : A → R. The operator system equivalent to Eq. (2.49) can be stated as L1 (v1 , v2 ). =. f1 (x,t) , (2.50). L2 (v1 , v2 ). =. f2 (x,t) ,. where (v1 , v2 ) represents (p, h), the operators L 1 and L2 are defined as the left hand x2 + ψ1 (x,t) + ψ2 (x,t). Note that the sides of Eq. (2.49) and f 1 (x,t) ≡ 0 and f 2 (x,t) = 2R x deformation, (d (x,t)), is a part of the operator, (L 2 ), because of the explicit dependence of pressure, here represented by v 1 . The multilevel technique applied to the solution of the one-dimensional Poisson problem can be modified to allow for two-dimensional functional operators. In this case, the smoother (section 2.3) is applied to the system Eq. (2.50) which after a few iterations produce the approximative solution (v˜ 1 i , v˜2 i ). The residuals rl i are then calculated according to: ∆x − fl∆x rl∆xi = Ll i v∆x 1 , v2 i , (l = 1, 2), to allow for the determination of the FAS right hand sides 2∆x 2∆x v2∆x − rl2∆x fl2∆x 1 , v2 I = Ll i I on the next coarser grid level. The coarse grid errors (E l I ) are then calculated by 2∆x El2∆x I = v˜l I − Vl I ,.

(53) 2.5. DIMENSIONLESS FORMULATION. 21. where Vl I = v˜∆x l 2I . These errors are prolonged in order to correct the solution (v˜ 1 , v˜2 ). Together with the description for the two-level technique for the one-dimensional Poisson problem, this description outlines the multilevel technique used to facilitate the numerical solution of Eq. (2.49).. 2.5 Dimensionless formulation The smooth Newtonian EHL line contact is dependent on the six input parameters w, E , us = (u1 + u2 ) , Rx α, η0 . For a smooth Reynolds-Eyring approach two extra input parameters are used; (u1 − u2 ) which could be introduced as a slide-to-roll ratio s = 2 (u 1 − u2 ) / (u1 + u2) and τ0 for the Eyring shear stress. For a rough Reynolds-Eyring approach the additional input parameters needed depend on how the surface topographies are modeled, e.g., for a single ridge/dent modeled by an amplitude and a wavelength only there are 8 + 2 = 10 input parameters. In order to restate the EHL line contact problem in dimensionless form the following set of dimensionless parameters where introduced, (Hamrock [12]): T = t/ (2b/us) ,. X = x/b,. H = h/ b2 /Rx , P = p/ph , ρ¯ = ρ/ρ0, where. b=. (2.51). µ¯ = µ/µ0,. 2w 8wRx , ph = = πE πb. wE 2πRx. (2.52). If this set of dimensionless parameters is used to restate the EHD lubricated line contact problem in dimensionless form, it becomes clear that the smooth incompressible problem is, mathematically, a two-parameter problem, see Moes [13], i.e., the solution is dependent on the parameters M and L, viz.. 2b3 ph 21/2 w π = 1/2 . (2.53) M= 1/2 1/2 4Rx µ0 us Rx E 1/2 µ us 0. . 2π E 3/4 µ0 us α = . (2.54) 1/4 M 21/4 Rx only. In this case, the Moes parameters given by Eq. (2.53) and Eq. (2.54) enables systematic parameter studies. That is, there is a one-to-one mapping between a given Moes pair (M, L) and parameters such as p max , hmin , hc , etc. It is also possible to use these parameters as determining factors for any two other independent parameters. Using the transformation given by Eq. (2.51) also prevents the introduction of truncation errors when solving the problem numerically. That is, because of the particular choice of transformation the values of the dimensionless parameters are restricted so as not to assume extremely large nor extremely small values. For example, the Hertzian pressure in a rolling element bearing is typically O 109 Pa whereas the film thickness is typically O 10−9 − 10−6 m and in terms of dimensionless quantities these parameters are O 100 . L = αph. 1/4 1/4.

(54) CHAPTER 2. FULL FILM EHD LUBRICATION. 22. In this work, the compressible solution is sought. This introduces an expression for the pressure-density relation which affect the solution. The density at ambient pressure, (ρ 0 ) does, however, not affect the solution The determining parameters are here represented by λ and α¯ which in terms of input parameters is given by: λ=. 3π2 µ0 us Rx E 6us R2x µ0 = ph b3 w2. and α¯ = αph =. αwE . 2πRx. (2.55). (2.56). Using the dimensionless parameters given by Eq. (2.51), Eq. (2.52) and the expressions Eq. (2.55) and Eq. (2.56) the resulting dimensionless Reynolds-Eyring equation (Eq. (2.23)) is transformed into: 3 ¯ ∂P ∂ ∂ ∂ ρH ¯ + ¯ S (ρH) (ρH) . (2.57) =λ ∂X µ¯ ∂X ∂X ∂T which is, explicitly dependent on the first determining parameter λ. Applying the same transformation, the dimensionless form of the film thickness equation (Eq. (2.24)) reads: H (X, T ) =. H0 (T ) +. X2 1 − 2 π. ∞ −∞. ln |X − X |P X , T dX (2.58). +. ¯ 2 (x,t) ¯ 1 (x,t) + ψ ψ. which, if ψi ≡ 0, i = 1 ∧ 2 is an equation totally independent of input parameters. Note that the dimensionless parameter H 0 (T ) has, in addition to the dimensionless form, a contribution from the dimensionless restated elastic deformation, i.e., 1 ∞ ln (b) P (X, T ) dX. H0 (T ) = h0 (T ) / b2 /Rx − π −∞. (2.59). The force-balance equation Eq. (2.29) transforms into:. I. P (X, T ) dX =. π . 2. (2.60). which is also an equation totally independent of input parameters when restated in dimensionless form. The Roelands expression for the viscosity depends on the determining parameter α¯ and the parameter p h when transformed into dimensionless form, viz. ¯ 0 αP Pph zvisc µ¯ (P) = exp −1 + 1 + . (2.61) ph zvisc P0 where. ph zvisc (2.62) (ln(µ0 ) + 9.67), α¯ or P0 is a determining parameter depending on which one of them P0 =. Note that, either zvisc that is specified..

(55) 2.6. DISCRETE FORMULATION. 23. The dimensionless form of the Dowson-Higginson expression for the density is given by ρ¯ (P) =. C1 + C2 Pph , C1 + Pph. (2.63). where the constants C1 and C2 is chosen so that ρ −→ 1.34 as p −→ ∞. Together with the cavitation condition P ≥ 0, (2.64) this completes the set of dimensionless equations used in a Reynolds-Eyring approach for the EHD lubricated line contact problem.. 2.6 Discrete formulation The set of equations yielding the EHD lubricated line contact problem given in section 2.5 are here discretized. The discrete form of the time-dependent Reynolds Eq. (2.57) used in this work can be written as: 1 k k k k k k k εi− 1 Pi−1 − εi− 1 + εi+ 1 Pi + εi+ 1 Pi+1 − 2 (∆X). 2. 2. 2. 2. k k λ k k ¯ ¯k ¯k 2∆X 3ρi Hi − 4ρi−1 Hi−1 + ρi−2 Hi−2 − λ ∆T. where εki± 1 2. (2.65). k−1 = 0, H ρ¯ ki Hik − ρ¯ k−1 i i. 3 εki + εki±1 ρ¯ ki Hik Sik k , εi = = . 2 µ¯ ki. (2.66). Which is of the second order in space (index i) and of the first order in time (index k). The following discrete form of the film thickness equation Eq. (2.58) is used: Hik − H0k − where Ki j. =. (Xi )2 1 + 2 π. ∑ Ki j Pkj + ψ¯ ki = 0,. (2.67). j. ∆X ∆X Xi − X j + ln Xi − X j + −1 − 2 2 ∆X ∆X ln Xi − X j − − 1 . Xi − X j − 2 2 . (2.68). The discretization technique used for the force-balance equation (2.60) is the trapetzodial rule. In the Roelands equation (2.61) and the Dowson-Higginson equation (2.63), used for determining the viscosity and density respectively, the dimensionless pressure P is simply replaced by its discrete representation Pi ..

(56) CHAPTER 2. FULL FILM EHD LUBRICATION. 24. 2.7 Solution method Equations (2.65) and (2.67) can be written in operator form (Eq. (2.38)) if the k − 1 term in Eq. (2.65) is moved to the right hand side of the equation. If the coefficients in Eq. (2.38) are linearized using the previously determined values of pˇ k and hˇ k obtained from any relaxation process, then it would be possible to rewrite the system in the form of Eq. (2.39), viz.

(57) k

(58) k−1

(59). A11 pˇ k , hˇ k A12 pˇ k , hˇ k p f1 = . (2.69) hk A22 pˇ k , hˇ k A21 pˇ k , hˇ k fk2 where the time dependence (k) is written out explicitly. As an example, for second order discretization in both the X- and the T - directions, the matrix A 11 is a tri-diagonal matrix with one sub and one super diagonal, A 12 is a tri-diagonal matrix with two sub diagonals. The matrix A21 depends on the discretization order of the integral expression for the elastic deformation and is a full matrix. The diagonal matrix A 22 is however, independent of the discretization order. If the Block-Jacobi method is applied to solve Eq. (2.69) it is straightforward to rewrite the system in indexes form, i.e., in the form of a 2 × 2 matrix system.

(60)

(61)

(62) ai bi pi f1 i = , (2.70) ci d i hi f2 i with coefficients: ai. =. bi. =. k k ε + ε 1 1 i+ 2 (∆X)2 i− 2 3 1 + −λ ρ¯ ki 2∆X ∆X 1. −. ci. =. ∆X ln. bi. =. 1. ∆X 2. (2.71). . and f1 i. =. − −. f2 i. =. . 1 2. (∆X) λ . k k εki− 1 Pi−1 + εki+ 1 Pi+1 2. . 2. λ k−1 k−1 k k − 4ρ¯ ki−1 Hi−1 − ρ¯ ki−2Hi−2 ρ¯ Hi 2∆X ∆T i. H0k +. Xi2 1 − 2 π. (2.72). ∑ Ki j Pkj − ψ¯1ki − ψ¯2ki . j

(63) =i. The solution of this system can be vectorized and thus very easily implemented in any higher level programming language. Unfortunately, the nature of the EHD lubricated line contact problem introduces some difficulties that have to be dealt with. The solution process consist of an inner relaxation loop for Eq. (2.70) and an outer loop to converge the.

(64) 2.7. SOLUTION METHOD. 25. force-balance criterion. When the solution process is to be facilitated by the multilevel technique described in section 2.4.4 the force-balance condition is relaxed only at the coarsest grid level. The convergence of the force-balance equation is very important in order to obtain sufficient accuracy of H 0 within each time step. This can be seen from the pressure solution of the coupled system Eq. (2.73), viz. Pi =. di f i − bi gi , a i d i − b i ci. (2.73). where the nominator contains the time dependence of the Reynolds equation. The nominator of Eq. (2.73) expresses numerically the time derivative, in dimensionless form, of the film thickness and thus the time derivative of H 0 (T ) . That is, a small ∆T will enhance the influence of an error in H 0 (T ) and the physical explanation of this behavior is the squeeze mechanism represented by ∂/∂T in Eq. (2.57). Moreover, this is not only important when using a coupled solution method; a serial/segregated solution method will be affected in the same way. At the outlet of the EHD lubricated conjunction the cavitation condition Eq. (2.64) is imposed. Thus, the solution of the Reynolds equation, which admits a negative pressure, is not valid here and the resulting non-Newtonian parameter S will drop to zero. Setting S = 1, for ∂P/∂X = 0 in the cavitated region, leads to an effective viscosity µ/S which equals the ambient viscosity µ 0 . This gives an almost smooth extension of S in the region of zero pressure. The importance of such an extension relates to the Block-Jacobi relaxation of Eq. (2.70), since Eq. (2.64) implies P = 0 in the cavitation zone and that the whole domain is considered at all iterations. In the solution process, this correction successively moves the starting point of the cavitation zone to a certain, converged, point..

(65)

(66) Chapter 3. Dry Elasto-Plastic Contact Machining processes aim to produce surface topographies to a given cost, that have a positive influence on contact conditions. For some applications product life may be of greatest importance in which case a machining process must be chosen so that wear is minimized. In other applications, low friction may be of primary importance and one of the tasks of the surface topography might consequently become to minimize the real area of contact. However, this may significantly reduce service life since a small contact area causes a very high pressure and thus an increased risk of wear. This chapter includes theory for both 2D and 3D dry elasto-plastic contacts. A contact mechanics problem can be governed by the minimum potential complementary energy theory which is an approach based on a deterministic roughness model. That is, pressures and displacements can be obtained by solving a minimum value of an integral equation, i.e., a variational problem, see Kalker [14], for rough surface contact mechanics. Combined with an equation of force-balance, the equation for potential complementary energy including constraints provides the basis for the solution of the dry elasto-plastic contact. It should be mentioned that it is possible to use pressure and gap constraints together with the Boussinesq equation, i.e., Eq. (3.5) to solve the elastic problem, Lubrecht and Ioannides [15]. For both the models described above it can also be argued that, if the pressure constraint includes an limiting value p ≤ H s , where Hs is the hardness of the softer surface material, and that the load is such that only a very small fraction of the resulting pressure profile is truncated by this limiting value, then this is one way of modeling the dry contact assuming elastic perfectly plastic surface material. This way of modeling elastic perfectly plastic surfaces does not, however, allow for computation of the resulting plastically deformed surface since it does not include in the energy dissipation due to plastic deformation. Using the variational approach, however, it is possible to account for energy dissipation and thus allow computation of the resulting plastically deformed surface. This model is meant to be used for prediction of parameters determining different contact conditions, such as the real area of contact, and not for determining the sub-surface stress field in detail.. 3.1 Statistical roughness models Modeling topographies of real surfaces is a difficult task because of their random structures. Various approaches to modeling the dry contact between surfaces have been re27.

(67) 28. CHAPTER 3. DRY ELASTO-PLASTIC CONTACT. ported in the literature. Greenwood and Williamson [16] (GW), considered topographies consisting of hemispherically shaped asperities of uniform radius. They assumed that the asperities had a Gaussian distribution, in height, about a mean plane. Greenwood and Tripp [17] extended this approach to handle the contact between two rough surfaces. Further improvements to the GW approach were added by Whitehouse and Archard [18], Nayak [19], Onions and Archard [20], Bush et al. [21, 22] and Whitehouse and Phillips [23, 24]. The input needed for these types of models is deduced from surface measurements. More precisely, all of these models are based on asperity curvature. This parameter depends on the measurement resolution (Poon and Bhushan [25]) and a mean value of the asperity curvature is often a poor approximation since it can vary considerably in between individual asperities. It is, therefore, a difficult task to find the correct input value for the model. Moreover, use of only a few parameters to describe a surface generates a one-tomany mapping, i.e., the same set of parameters can be deduced for surfaces obtained by completely different machining processes.. 3.2 Deterministic roughness models In line with improvements in affordable computing power, deterministic approaches to the dry contact problem have been developed. These are becoming more and more realistic as computational speed increases. However, simplified models of the material and/or the topography are still needed. This is necessary to minimize the computing time and to be able to study different effects independently. For example, the Hertzian contact assumes linear elastic, frictionless materials and parabolic surface profiles for which an analytical solution is available. That is, the contact between a single asperity, approximated by a parabolic profile, and a rigid plane may be studied analytically. Westergaard [26] showed that it is possible to solve analytically the contact between a single frequency sinusoidal surface (mathematically: sin (2πnx)) in contact with a rigid plane. On a micro scale these analytical solutions represent quite realistic topographies. However, the coupling to asperities at other scales is neglected and a more complete representation is therefore needed. In this work the total surface content, i.e., all data of the measured or modeled surface topography is used to determine different contact conditions such as the real area of contact. It should be mentioned that a deterministic roughness model benefits from surface characterization (Chapter 4) which can be used tp deduce the data and possibly facilitates the interpretation of the results obtained.. 3.3 Numerical solution techniques In addition to the different ways of modeling the mating materials and topography, there are numerous deterministic numerical techniques that can be applied to solve contact problems between rough surfaces. Lubrecht and Ioannides [15] applied multilevel techniques to facilitate the numerical solution of the elastic contact problem. Ren and Lee [27] applied a moving grid method to reduce storage of the influence matrix when the conventional matrix inversion approach is used to solve the contact problem of linear elastic bodies with rough surfaces. Björklund and Andersson [28] extended the conventional matrix inversion approach by incorporating friction induced deformations. This was done by using the assumption of linear elastic material. Ju and Farris [29] introduced an FFT-based method to.

(68) 3.4. GOVERNING EQUATIONS. 29. solve the elastic contact problem. Stanley and Kato [4] combined this FFT-based method with a variational principle to solve both the 2D and the 3D contact problem of rough surfaces. It should be noted that the contact between realistic topographies, under relatively small loads, leads to plastic deformations. See for example, Tian and Bhushan [3] who based their theoretical model on a variational principle for both linear elastic and linear elastic perfectly plastic materials. In contrast to the model given by Stanley and Kato, the present model incorporates linear elastic-perfectly plastic materials and the energy dissipation due to plastic deformation is accounted for. In this way, not only the in-contact topography and the corresponding pressure distribution but also the unloaded plastically deformed topography is obtained.. 3.4 Governing equations Assuming frictionless linear elastic contact, the variational problem including constraints to be solved ([14]) is given by (3.1). A general formulation is introduced here that admits both 2D and 3D problems. . 1 p de dΩ − p ψ dΩ , min () = min p≥0 p≥0 2 Ω Ω (3.1). S. p dS = W,. where p is the pressure, d e is the normal, elastic deformation, ψ is the gap between the undeformed surfaces and W is the applied load. The mathematical model given by Eq. (3.1) has one equivalent in the system p = 0,. h > 0,. p > 0,. h = 0,. S. (3.2). p dS = W,. where h = de + ψ + const, see Lubrecht and Ioannides [15]. Restated, this model may be expressed as: h = de + ψ + const, S. p dS = W,. (3.3). p ≥ 0. However, some modification of Eq. (3.1) is necessary to solve the contact problem assuming linear elastic - perfectly plastic surface materials, model depicted in Fig. 3.1..

(69) CHAPTER 3. DRY ELASTO-PLASTIC CONTACT. 30. σ. σs. Elastic. Perfectly plastic. ε. Figure 3.1: Material model for linear elastic - perfectly plastic surfaces To account for the energy dissipation due to plastic deformations, the variational statement given by Eq. (3.1) has to be restated (see Tian and Bhushan [3]) and the variational approach reads: . 1 min () = min p de dΩ − p (ψ − d p) dΩ . 0≤p≤Hs 0≤p≤Hs 2 Ω Ω (3.4). p dS S. = W,. where Hs is the hardness of the softer material and d p is the normal, plastic deformation. For two elastic half spaces, assuming plain strain, the total normal surface displacement de (x) for a given pressure distribution can be obtained from Eq. (3.5). de (x) =. ∞. −∞. K (x − s) p (s) ds + const.. (3.5). where bold symbols indicate vector valued parameters, e.g., x = x in 2D and x = (x, y) in 3D and the integral kernel K is given by K (x − s) = − and K (x − s, y − t) =. 4 ln|x − s| , πE . 2 1 πE (x − s)2 + (y − t)2. (3.6). (3.7). for 2D and 3D contacts respectively.. 3.5 Spectral analysis In the case of a sinusoidal pressure p (x) of one frequency only it is possible to find the exact, closed form solution for the normal surface displacement d e (x), see for example.

No results found