Ecole Sigle

Doctorale

Renée El

Melhem Bât

Pascal Blaise

R Gourdon

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

M Philippe
email: philippe.normand@univ-lyon1.fr

Mme Emmanuelle
email: emmanuelle.canet@univ-lyon1.fr

M Luca

Jean-Yves Buffière
email: jean-yves.buffiere@insa-lyon.fr

Stéphanie Cauvin

M Jocelyn Bonjour

Véronique Guichard

J Y Toussaint

M Christian Montes

Keywords: Elastohydrodynamic lubrication, Numerical simulation, Piezoviscous elastic regime, Amplitude Reduction Theory, Friction variation i Lubrification élastohydrodynamique, Simulation numérique, Régime élasto piézovisqueux, Théorie de réduction d'amplitude, Variation de frottement measured surface roughness . . . . . . . .

The friction of interfacial surfaces greatly influences the performance of mechanical elements. Friction has been investigated experimentally in most studies. In this work, the friction is predicted by means of numerical simulation under an elastohydrodynamic lubrication (EHL) rough contact condition.

The classical Multigrid technique performs well in limiting computing time and memory requirements. However, the coarse grid choice has an important influence on code robustness and code efficiency to solve the rough problem. In the first part of this work, a coarse grid construction method proposed by Alcouffe et al. is implemented in the current time-independent EHL Multi-Grid code. Then this modified solver is extended to transient cases to solve the rough contact problem.

The friction curve is usually depicted as a function of "Λ ratio", the ratio of oil film thickness to root-mean-square of the surface roughness. However this parameter is less suitable to plot friction variations under high pressure conditions (piezoviscous elastic regime). In the second part of this work, the friction coefficient is computed using the modified EHL code for many operating conditions as well as surface waviness parameters. Simulation results show that there is no single friction curve when the old parameter "Λ ratio" used. Based on the Amplitude Reduction Theory, a new scaling parameter depends on operating condition and waviness parameters is found, which can give a unified friction curve for high pressure situation.

For more complex rough surfaces, a power spectral density (PSD) based method is proposed to predict friction variations in the third part of this work. The artificial surface roughness is employed to test the rapid prediction method firstly. Good agreement is found between the full numerical simulation and this rapid prediction. Then the rapid prediction method is applied to analyze the friction variation of measured surface roughness. A comparison is also made between predictions and experiments.

Both the new scaling parameter and the friction increase predicted by the PSD method show good engineering accuracy for practical use.

Résumé

Le frottement à l'interface des surfaces influence les performances des éléments mécaniques. Le frottement a été étudié expérimentalement dans la plupart des études. Dans ce travail, le frottement est prédit à l'aide d'une simulation numérique dans des conditions de contact rugueux avec une lubrification élastohydrodynamique (EHL).

La technique classique Multigrille fonctionne bien pour limiter le temps de calcul et les besoins en mémoire. Cependant, le choix de la grille grossière a une influence importante sur la robustesse du code et son efficacité pour résoudre le problème brut. Dans la première partie de ce travail, une méthode de construction de grille grossière proposée par Alcouffe et al. est implémenté dans le code EHL Multigrille indépendamment du temps. Ensuite ce solveur modifié est étendu aux cas transitoires pour résoudre le problème de contact avec rugosité.

La courbe de frottement est généralement représentée en fonction du « Λ ratio », le rapport entre l'épaisseur du film d'huile et la valeur moyenne quadratique de la rugosité de la surface. Cependant, ce paramètre est moins approprié pour tracer les variations de frottement dans des conditions de haute pression (régime élasto piézo-visqueux). Dans la deuxième partie de ce travail, le coefficient de frottement est calculé à l'aide du code EHL modifié pour de nombreuses conditions de fonctionnement ainsi que pour les paramètres d'ondulation de surface. Les résultats de la simulation montrent qu'il n'y a pas de courbe de frottement unique lorsque l'ancien paramètre « Λ ratio » est utilisé. En se basant sur la théorie de la réduction d'amplitude, un nouveau paramètre de dimensionnement qui dépend des conditions de fonctionnement et des paramètres d'ondulation est trouvé, ce qui peut donner une courbe de frottement unique pour les situations de haute pression.

Pour les surfaces rugueuses plus complexes, une méthode basée sur la densité spectrale de puissance (PSD) est proposée pour prédire les variations de frottement dans la troisième partie de ce travail. La rugosité artificielle de la surface est utilisée pour tester d'abord la méthode de prédiction rapide. Un bon accord est trouvé entre la simulation numérique complète et cette prédiction rapide. La méthode de prédiction rapide est ensuite appliquée pour analyser la variation de frottement de la rugosité de surface mesurée. Une comparaison est également faite entre les prédictions et les expériences.

List of Figures

= 0.5 × (u 1 + u 2 ) m/s U dimensionless speed parameter U = (η 0 u r )/(E R x ) ∆U slide-to-roll ratio ∆U = δu/u r = (u 1 -u 2 )/u r U rat slip parameter U rat = u 1 /u r v 1 , v 2 Poisson ratios of the two contact bodies w normal load N W 2 2d dimensionless load parameter W 2 = w/(E R 2 x )
x coordinate in the rolling direction m

X dimensionless coordinate X = x/

Background

Energy is the most important of all resources, which is needed to support economic and social progress and build a better quality of life. According to the report of IEEJ Outlook 2019 (Institute of Energy Economics, Japan), the world primary energy consumption will continue growing from 2018 to 2050. Most of this growth comes from non-OECD (non-Organization for Economic Cooperation and Development) countries, where demand is driven by strong economic growth (shown in Figure 1.1). Fossil energy like oil, coal and natural gas are still the largest energy source of the world, and the reserves of fossil energy are limited. Meanwhile, the over-consumption of fossil fuel leads to the over-release of carbon dioxide (CO 2 ), methane, oxynitride (NO x ) and particulates into the air, which disturb the natural balance of the atmosphere. The rapid rise in the carbon dioxide contributes to the serious global warming problem. NASA (National Aeronautics and Space Administration) reported that the global surface temperature has been persistently increasing since the late 19th century (shown in Figure 1.2). An investigation from the U.S. Energy Information Administration (EIA) shows that the industrial and transportation sectors have consumed the most energy (shown in Figure 1.3) as well as produced the most CO 2 emissions (show as Table 1.1).

However, a substantial amount of energy is not put to useful purposes. Researches [START_REF] Holmberg | Global energy consumption due to friction in passenger cars[END_REF][START_REF] Holmberg | Global energy consumption due to friction in paper machines[END_REF][START_REF] Holmberg | Global energy consumption due to friction in trucks and buses[END_REF][START_REF] Holmberg | Global energy consumption due to friction and wear in the mining industry[END_REF] show that a considerable amount of energy in industrial and transportation is consumed to overcome friction. For instance, energy consumed to overcome friction over the total energy consumption in heavy-duty vehicles is 33%, in paper machines is 32%, in passenger cars is 33% and in mineral mining industry is 40%. Recently, the increasing environment awareness requires efforts to improve energy efficiency and reduce CO 2 production. Correct lubrication between engineering provides sufficient separation of the roughness present on the contact surfaces, which contributes to reducing friction losses. Better understanding and control of friction in mechanical components has the potential to offset large energy savings and CO 2 emission reduction [START_REF] Lee | Tribological opportunities for enhancing america's energy efficiency[END_REF]. Studies [START_REF] Holmberg | Global impact of friction on energy consumption, economy and environment[END_REF][START_REF] Holmberg | Influence of tribology on global energy consumption, costs and emissions[END_REF][START_REF] Wong | Overview of automotive engine friction and reduction trends-effects of surface, material, and lubricant-additive technologies[END_REF] estimated that with the implementing advanced tribological technologies, such as using new contact surface, materials, lubricants, energy losses due to friction and wear could potentially be reduced by 40% in the long term (15 years) and by 18% in the short term (8 years )and CO 2 emissions globally can also reduced by 1,460 MtCO 2 (million tonnes CO 2 ) in the short term and by 3,140 MtCO 2 in the long term. only for meeting the increasing requirements of energy efficiency and CO 2 emissions reduction but also for providing a theoretical tool in element design and optimization.

Literature review

Elastohydrodynamic lubrication (EHL) is the type of lubrication for frictional pairs having elastic contact under very high pressure and forming lubricant film in non-conformal contacts, such as rolling bearings, gears, human synovial joints and so on [START_REF] Lugt | A review of elasto-hydrodynamic lubrication theory[END_REF]. Lubricant film and surface roughness play an important role for improving reliability and effectiveness of mechanical parts as well as reducing friction losses. The majority of the published work on the influence of surface roughness on friction has been experimental, the minority of theoretical work has been done on friction prediction. This section represents the literature review on the methods to solve the rough EHL contact problem as well as friction in rough EHL contacts.

Methods to solve the rough contact problem

Typically, the EHL model consists of five equations [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF], in which the Reynolds equation Equation 1.1 is a partial differential equation and the film thickness equation Equation 1.2 contains an integral term, both equations are required to be solved simultaneously, making these equations very complex. There are many approaches that can be used to solve this EHL model: the inverse method [START_REF] Dowson | A numerical solution to the elasto-hydrodynamic problem[END_REF], the Newton-Raphson method [START_REF] Wolff | The application of newton-raphson method to thermal elastohydrodynamic lubrication of line contacts[END_REF], the homotopy method [START_REF] Nurgat | Solving ehl problems using iterative, multigrid, and homotopy methods[END_REF], the finite element method [START_REF] Oh | Numerical solution of the point contact problem using the finite element method[END_REF], the Multigrid method [START_REF] Lubrecht | Multigrid, an alternative method for calculating film thickness and pressure profiles in elastohydrodynamically lubricated line contacts[END_REF][START_REF] Lubrecht | Multigrid, an alternative method of solution for two-dimensional elastohydrodynamically lubricated point contact calculations[END_REF] and the Navier-Stokes approach [START_REF] Almqvist | The navier-stokes approach for thermal ehl line contact solutions[END_REF]. So far, the Multigrid algorithm has been considered as one of the most efficient methods and applied frequently to EHL problems.

∂ ∂x ( ρh 3 12η ∂p ∂x ) + ∂ ∂y ( ρh 3 12η ∂p ∂y ) -u r ∂(ρh) ∂x = 0 (1.1) h(x, y) = h 0 + x 2 2 + y 2 2 -r r (x, y) + 2 π 2 +∞ -∞ +∞ -∞ P (x , y ) (x -x ) 2 + (y -y ) 2 d x d y (1.2)
Where p represents the pressure, h denotes the film thickness, h 0 is the mutual approach and u r = (u 2 + u 2 )/2 is the mean velocity of two contact surfaces. ρ and η are the density and viscosity of the lubricant, respectively. The x axis is aligned with the direction of the mean velocity ū, and the y axis is perpendicular to the x direction.

In engineering, no surface is perfectly smooth, the order of magnitude of the surface roughness is often the same as or greater than that of the film thickness predicted by smooth contact conditions [START_REF]History of EHL Development[END_REF]. Therefore, the surface roughness should be considered. Generally, there are two methods to treat the rough lubrication problem numerically. CHAPTER 1. INTRODUCTION One approach is called the "stochastic" method. Early work was conducted on the hydrodynamic lubrication (HL) problem. Theoretical analysis of the implementation of the stochastic theory on rough HL contact problem was described by Christensen [START_REF] Christensen | The hydrodynamic lubrication of rough bearing surfaces of finite width[END_REF][START_REF] Christensen | Some aspects of the functional influence of surface roughness in lubrication[END_REF]. Then Patir and Cheng [START_REF] Patir | An average flow model for determining effects of threedimensional roughness on partial hydrodynamic lubrication[END_REF][START_REF] Patir | Application of average flow model to lubrication between rough sliding surfaces[END_REF] proposed an average flow model to determine the effects of surface roughness on rough-lubricated contacts. In this model the Reynolds equation is simplified as an average Reynolds equation using a independent flow factor (shown in Figure 1.4). Since the pioneering studies on the "average flow model", a number of authors have extended and generalized this work. Tripp [START_REF] Tripp | Surface roughness effects in hydrodynamic lubrication: the flow factor method[END_REF] re-computed the flow factors using a perturbation expansion of the pressure in a nominal parallel film. When small roughness amplitude is considered, the flow factors calculated in Reference [START_REF] Tripp | Surface roughness effects in hydrodynamic lubrication: the flow factor method[END_REF] agree well with that of Patir and Cheng [START_REF] Patir | An average flow model for determining effects of threedimensional roughness on partial hydrodynamic lubrication[END_REF]. Hu and Zheng [START_REF] Hu | Some aspects of determining the flow factors[END_REF] studied the influence of boundary conditions, grid systems and statistics of rough surfaces on the flow factors. Lunde and Tonder [START_REF] Lunde | Pressure and shear flow in a rough hydrodynamic bearing, flow factor calculation[END_REF] calculated the flow factors for an isotropic rough bearing and found that the boundary conditions of the selected bearing part can not affect the flow passing through. Subsequently, Zhu and Cheng [START_REF] Zhu | Effect of surface roughness on the point contact ehl[END_REF] extended the flow factors method in the point EHL contact problem. Some authors [START_REF] Harp | An average flow model of rough surface lubrication with inter-asperity cavitation[END_REF][START_REF] Bayada | An average flow model of the reynolds roughness including a mass-flow preserving cavitation model[END_REF] applied the flow factors to deal with the cavitation problem. Letalleur et al. [START_REF] Letalleur | Average flow model of rough surface lubrication: flow factors for sinusoidal surfaces[END_REF] studied the flow factors for two rough cases: smooth-rough stationary case and rough-rough unstationary case. Sahlin et al. [START_REF] Sahlin | Rough surface flow factors in full film lubrication based on a homogenization technique[END_REF] developed a novel method using a homogenization technique to compute the flow factors. However, in stochastic model, roughness asperities are mainly treated as rigid. Another way is to incorporate the surface roughness term r r (x, y, t ) in the film thickness equation (shown as Equation 1.2) and to solve the system of equations directly. Due to the limitation of computation of speed and storage space, preliminary research [START_REF] Goglia | The effects of surface irregularities on the elastohydrodynamic lubrication of sliding line contacts. part i-single irregularities[END_REF][START_REF] Goglia | The effects of surface irregularities on the elastohydrodynamic lubrication of sliding line contacts. part ii-wavy surfaces[END_REF][START_REF] Houpert | Elastohydrodynamic lubrication calculations used as a tool to study scuffing[END_REF][START_REF] Venner | Surface roughness effects in an ehl line contact[END_REF][START_REF] Hooke | The behaviour of low-amplitude surface roughness under line contacts[END_REF] studied the steady state line rough contact problem, where the surface roughness is time-independent and one-dimensional model was considered. Later on, the stationary two dimensional rough contact problem [START_REF] Lubrecht | The numerical solution of the elastohydrodynamically lubricated lineand point contact problem using multigrid techniqes[END_REF][START_REF] Lubrecht | The influence of longitudinal and transverse roughness on the elastohydrodynamic lubrication of circular contacts[END_REF][START_REF] Kweh | Micro-elastohydrodynamic lubrication of an elliptical contact with transverse and three-dimensional sinusoidal roughness[END_REF][START_REF] Kweh | Simulation of elastohydrodynamic contacts between rough surfaces[END_REF] were carried out. With the increasing development of computational technique, transient cases were studied by many authors [START_REF] Ai | A transient ehl analysis for line contacts with measured surface roughness using multigrid technique[END_REF][START_REF] Chang | On the pressure rippling and roughness deformation in elastohydrodynamic lubrication of rough surfaces[END_REF][START_REF] Chang | A study of elastohydrodynamic lubrication of rough surfaces[END_REF][START_REF] Greenwood | The behaviour of transverse roughness in sliding elastohydrodynamically lubricated contacts[END_REF][START_REF] Greenwood | The behaviour of transverse roughness in ehl contacts[END_REF][START_REF] Lubrecht | Aspects of two-sided surface waviness in an ehl line contact[END_REF][START_REF] Morales | Elastohydrodynamic lubrication of smooth and rough surfaces[END_REF][START_REF] Osborn | Time dependent line ehd lubrication using the multigrid/multilevel technique[END_REF][START_REF] Venner | Numerical simulation of the overrolling of a surface feature in an ehl line contact[END_REF][START_REF] Venner | Multilevel solution of the EHL line and point contact problems[END_REF][START_REF] Venner | Numerical simulation of a transverse ridge in a circular ehl contact under rolling/sliding[END_REF][START_REF] Venner | Numerical simulation of waviness in a circular ehl contact, under rolling/sliding[END_REF]. Based on the previous studies on transient rough contact problem, Venner and Lubrecht et al. [START_REF] Venner | Amplitude reduction of waviness in transient ehl line contacts[END_REF][START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF][START_REF] Lubrecht | Waviness amplitude reduction in ehl line contacts under rolling-sliding[END_REF][START_REF] Venner | Amplitude reduction of small-amplitude waviness in transient elastohydrodynamically lubricated line contacts[END_REF][START_REF] Lubrecht | Elastohydrodynamic lubrication of rough surfaces[END_REF][START_REF] Jacod | Amplitude reduction of waviness in elastohydrodynamic lubrication using an eyring fluid model[END_REF][START_REF] Venner | An engineering tool for the quantitative prediction of general roughness deformation in ehl contacts based on harmonic waviness attenuation[END_REF][START_REF] Chapkov | Roughness amplitude reduction under non-newtonian ehd lubrication conditions[END_REF][START_REF] Wang | Amplitude reduction in ehl line contacts under rolling sliding conditions[END_REF] published a series of papers on the "Amplitude Reduction Theory" describing the relation between the surface roughness deformation and operating conditions. They found that under very high pressure situations (piezoviscous elastic regime), the surface roughness will deform, and this deforma-CHAPTER 1. INTRODUCTION tion depends on operating conditions. A master curve (shown in Figure 1.5) describes this relation quantitatively. Then Šperka et al. [START_REF] Šperka | Experimental study of real roughness attenuation in concentrated contacts[END_REF] verified the "Amplitude Reduction Theory" by measuring the deformed surface roughness on an optical test rig, the comparison of the measured results and the predicted results can be found in Figure 1.6. Recently, some extension work on the rough contact problem were addressed [START_REF] Holmes | Transient elastohydrodynamic point contact analysis using a new coupled differential deflection method part 1: theory and validation[END_REF][START_REF] Evans | Deterministic mixed lubrication modelling using roughness measurements in gear applications[END_REF][START_REF] Zhang | Influence of lubrication starvation and surface waviness on the oil film stiffness of elastohydrodynamic lubrication line contact[END_REF][START_REF] Gu | A comparative study of tribological performance of helical gear pair with various types of tooth surface finishing[END_REF]. One of the most important studies about the influence of surface roughness on friction will be represented in the next subsection. The surface roughness term r r (x, y), which is often considered to be the same order of magnitude as the oil film thickness, incorporated in Equation 1.2 makes the coefficient (ρh 3 )/(12η) in the Reynolds equation jump orders of magnitude, which leads to a significant variation in EHL equations. From a mathematical point of view, the coefficient (ρh 3 )/(12η) is continuous, while in the numerical simulation, the coefficient causes a strong discontinuity in discrete CHAPTER 1. INTRODUCTION Reynolds equation. Work by Alcouffe et al [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] proposed an efficient way to overcome the above discontinuity through constructing the coarse grid in a Multi-grid code. The present work employs the Multigrid techniques [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF] to solve the rough contact problem. In addition, the idea provided by Alcouffe et al [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] is also applied to improve the code robustness and code efficiency.

Friction in rough EHL contact problem

Most heavily loaded machine elements are working under elastohydrodynamically lubricated conditions. Understanding the frictional behavior in such contacts plays an important role for reducing friction, preventing wear as well as improving service life. To reveal the relation between fluid properties and friction, research has been conducted as follows: Crook [START_REF] Crook | The lubrication of rollers iv. measurements of friction and effective viscosity[END_REF] used a disc machine for measuring frictional traction. It was shown that the rolling friction (the traction due to rolling) is independent of load and simply proportional to the film thickness in the elasto-hydrodynamic regime. Johnson and Cameron [START_REF] Johnson | Fourth paper: shear behaviour of elastohydrodynamic oil films at high rolling contact pressures[END_REF] measured traction in a rolling contact disc machine and results showed that the traction first increases and then decreases when the sliding speed increases. Johnson and Roberts [START_REF] Johnson | Observations of viscoelastic behaviour of an elastohydrodynamic lubricant film[END_REF] observed the visco-elastic behavior of film thickness through measuring shear forces on a rolling-contact test rig. Evans and Johnson [START_REF] Evans | Regimes of traction in elastohydrodynamic lubrication[END_REF] constructed traction maps depending upon pressure, temperature as well as shear rate for different fluids, where different areas represent different traction behavior. Zhang et al. [START_REF] Zhang | Ehl analysis of rib-roller end contact in tapered roller bearings[END_REF] studied the elliptical contact between rib face and roller end in tapered roller bearings by means of a full numerical simulation. They found that the elastic deformation has a non-negligible influence on the friction coefficient. Yu and Medley [START_REF] Yu | Influence of lubricant additives on friction in a disc machine[END_REF] studied the influence of lubricant additives on friction via a side-slip disc machine. They concluded that the limiting shear stress, which is a useful parameter for predicting friction, is affected by the lubricant additives. Jacod et al. [START_REF] Jacod | A generalized traction curve for ehl contacts[END_REF] predicted the coefficient of friction over a wide range of operating conditions and obtained a single generalized friction curve based on a full numerical simulation for a non-Newtonian EHL contact model. Vicente et al. [START_REF] Vicente | The frictional properties of newtonian fluids in rolling-sliding soft-ehl contact[END_REF] explored friction in rolling-sliding, soft-EHL contacts numerically and experimentally. Numerical calculations of the Couette friction are in good agreement with measured results. Very recently, Liu et al [START_REF] Liu | Starved lubrication of a spur gear pair[END_REF] calculated the friction coefficient in a gear contact interface numerically, based on a thermal starved EHL model. They found that the maximum friction coefficient appears at the engaging-in point where a considerable slide-to-roll ratio exists. Björling et al. [START_REF] Björling | The effect of ageing on elastohydrodynamic friction in heavy-duty diesel engine oils[END_REF] measured the friction under EHL conditions on a ball-on-disc test rig for aged and fresh oils. Results showed that there is no difference in friction. In addition, Zhang [START_REF] Zhang | Effect of base oil structure on elastohydrodynamic friction[END_REF] measured the EHL friction for a wide range of base fluids and compared the friction values for five different operating conditions. The study underlined the importance of molecular structure of the base fluid in determining the EHL friction.

Studies [START_REF] Emmens | The influence of surface roughness on friction[END_REF][START_REF] Wu | A friction model of partial-ehl contacts and its application to power loss in spur gears[END_REF][START_REF] Lee | Modeling of the friction caused by lubrication and surface roughness in sheet metal forming[END_REF] showed that surface roughness has a significant impact on the friction behavior of lubricated surfaces. A useful tool to investigate the frictional behavior between rough surfaces is the classical Stribeck curve, showing the friction coefficient is a function of a ratio of the averaged oil film thickness to the combined surface roughness. The original research about the Stribeck curve dates back to the 19th century. In 1879, Thurston gave precise values of the friction coefficient and he was probably the first person to report that the friction coefficient passed through a minimum as the load increased [START_REF] Thurston | Friction and Lubrication: Determinations of the Laws and Coëfficients of Friction by New Methods and with New Apparatus[END_REF][START_REF] Dowson | History of tribology[END_REF]. Twenty years later, Stribeck [START_REF] Stribeck | Die wesentlichen eigenschaften der gleit-und rollenlager[END_REF][START_REF] Stribeck | Die wesentlichen Eigenschaften der Gleit-und Rollenlager: Untersuchung einer Tandem-Verbundmaschine von 1000 PS[END_REF] systematically published results of a carefully conducted and wide-ranging series of experiments on journal bearings, which are frequently referred to as 'the Stribeck curve' (shown in Figure 1.7). G ümbel [START_REF] Gümbel | Das problem der lagerreibung[END_REF] organised Stribeck's experimental results in a single curve by plotting the friction against the parameter ηω/ p, where η is the lubricant viscosity, ω is the angular velocity of the shaft and p is the load per unit length. At the same time, Hersey [START_REF] Hersey | The laws of lubrication of horizontal journal bearings[END_REF] conducted CHAPTER 1. INTRODUCTION experiments on journal bearings and plotted the friction coefficient against the load, speed, temperature, viscosity and rate of oil supply. He showed that hydrodynamic friction should be a function of ηn/p in which n is the rotational speed and p is the pressure. Many years later, Wilson and Barnard [START_REF] Wilson | The mechanism of lubrication[END_REF] replotted the Stribeck curve by introducing a new variable i.e. zn/p, where the lower-case z stands for the lubricant viscosity. Subsequently, McKee [START_REF] Mckee | The effect of running-in on journal bearing performance[END_REF] provided a similar dimensionless group Z N /P . Vogelpohl et al. [START_REF] Vogelpohl | Die stribeck-kurve als kennzeichen des allgemeinen reibungsverhaltens geschmierter gleitflächen[END_REF] incorporated the boundary and fluid friction coefficient and showed a transition from the hydrodynamic lubrication regime to the mixed lubrication regime. All of the work mentioned above is performed under low pressure conditions, in the isoviscous rigid regime [START_REF] Johnson | Regimes of elastohydrodynamic lubrication[END_REF]. The situation for non-conforming contacts, such as those occurring in rolling element bearings, gears and cams, is somewhat different [START_REF] Spikes | Mixed lubrication-an overview[END_REF]. Shotter [START_REF] Shotter | Experiments with a disc machine to determine the possible influence of surface finish on gear tooth performance[END_REF] experimentally showed that the friction increases with the surface roughness. Tallian and his co-workers [START_REF] Tallian | Lubricant films in rolling contact of rough surfaces[END_REF][START_REF] Tallian | Paper 14: Partial elastohydrodynamic lubrication in rolling contact[END_REF] proposed a ratio |ξ 0 | between the elastohydrodynamic film thickness and the composite root mean square roughness to represent the mixed elastodydrodynamic regime ( 1 < |ξ 0 | < 4 ). Poon [START_REF] Poon | Third paper: frictional behaviour of lubricated rollingcontact elements[END_REF] was concerned with the transition from the boundary to the mixed regime with a dimensionless parameter 1 ξ 2 and the transition from mixed to full EHL region with 2 < ξ 2.4 by using electrical-conductivity measurements. Bair and Winer [START_REF] Bair | Regimes of traction in concentrated contact lubrication[END_REF] plotted the reduced traction coefficient as a function of a lambda ratio by performing sliding-rolling experiments. They found that when the lambda ratios is less than 2 the contact moves into the mixed regime. In general, the Stribeck curve can be divided into three regimes [START_REF] Stachowiak | Engineering tribology[END_REF]: λ > 3 represents the full-film regime, 1 λ 3 is the mixed EHL regime and λ < 1 indicates the boundary regime. However, study [START_REF] Cann | The lambda ratio-a critical reexamination[END_REF] shows that this lambda ratio is not a suitable parameter to determine lubrication states especially when some aspects such as non-Newtonian, thermal and transient effects are considered. Transition locations from mixed to boundary lubrication regime or from full-film to mixed lubrication regime are still ambiguous. Therefore, an appropriate grouping including the speed, film thickness and roughness is required. Schipper [START_REF] Schipper | On the transitions in the lubrication of concentrated contacts[END_REF] suggested a so-called Lubrication number L , which takes viscosity, speed and pressure into consideration, to detect the variation of the friction coefficient. Recently, Gelinck [START_REF] Gelinck | Calculation of stribeck curves for line contacts[END_REF] extended Johnson's model [START_REF] Johnson | A simple theory of asperity contact in elastohydro-dynamic lubrication[END_REF] to calculate the coefficient of friction for the whole mixed EHL regime. Lu and Khonsari [START_REF] Lu | The stribeck curve: experimental results and theoretical prediction[END_REF] CHAPTER 1. INTRODUCTION examined the behavior of the Stribeck curve theoretically and experimentally on a journal bearing and found a good agreement. Wang et al. [START_REF] Wang | Simulations and measurements of sliding friction between rough surfaces in point contacts: From ehl to boundary lubrication[END_REF] presented a numerical approach developed on the basis of deterministic solutions of mixed lubrication to evaluate sliding friction. Meanwhile, they measured the sliding friction on a commercial test rig. Both results were plotted against sliding velocities and also showed good agreement. Kalin [START_REF] Kalin | The stribeck curve and lubrication design for nonfully wetted surfaces[END_REF] investigated changes of the Stribeck curve when one or two surfaces in the contact are non-fully wetted. Afterwards, Kalin [START_REF] Kalin | Non-conventional inverse-stribeck-curve behaviour and other characteristics of dlc coatings in all lubrication regimes[END_REF] tested the variations of the friction coefficient with diamond-like carbon coatings (DLC). Zhang [START_REF] Zhang | Friction prediction of rolling-sliding contact in mixed ehl[END_REF] developed a numerical approach assuming the asperity interaction friction is proportional to the contact area to predict the mixed EHL friction coefficient. Bonaventure [START_REF] Bonaventure | Transition between mixed lubrication and elastohydrodynamic lubrication with randomly rough surfaces[END_REF] and his co-authors conducted rolling-sliding experiments with random surface roughness, they found that the onset of ML occurs at a higher entrainment product η 0 u e (in which η 0 is inlet viscosity and u e is entrainment speed) and a relevant roughness scalar parameter was obtained to predict the onset position.

Most of the work on Stribeck curve was done by experiments. Current study employs the Amplitude Reduction Theory [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF] to study the frictional behavior in piezoviscous elastic regime [START_REF] Dowson | The piezo-viscous fluid, rigid solid regime of lubrication[END_REF] by means of numerical simulation.

Research aims and Outlines

Research aims

Long term successful operation of mechanical devices greatly depends on correct lubrication of the mechanical elements to provide sufficient separation of the roughness present on the contact surfaces. However, lubrication provides another important role, reducing friction between rough contact surfaces.

The objective of the present research project is to develop an efficient and robust Multi-Grid-based algorithm to study the frictional behavior between rough contact surfaces. Current MultiGrid codes show the required efficiency, but are not sufficiently robust to treat the rough surface problem in a general way. Difficulties may lie in the following aspects:

(i) The efficient construction of the coarse grid of EHL Multi-Grid model to guarantee the code robustness and code efficiency of impact, rough surface EHL contact problems.

(ii) Tests of the increased robustness of the new EHL Multi-Grid solver.

(iii) Implementation applied to test the code robustness and code efficiency of rough surface EHL contact problems.

(iv) The extension of the developed Multi-Grid lubrication code to transient contact problems.

(v) The computation of the friction coefficient of rough contact surfaces.

(vi) The unification of friction curves that differ according to operating conditions.

(vii) The extension of the lambda ratio parameter predicting the transition from mixed to fullfilm regimes.

CHAPTER 1. INTRODUCTION

Outlines

According to the research aims listed in the previous sub-section. The layout of this thesis is as follows:

Chapter 1: This chapter first emphasizes the important role of friction played in energy consumption and environmental issues. Subsequently, a literature review on the methods to solve rough contact problems and friction in rough contact problems is illuminated. The objective and structure of the present thesis are given in the last section.

Chapter 2: This chapter represents the numerical model and algorithm for solving the transient rough EHL contact problem. The governing equations for transient EHL model are introduced first. Then the method proposed by Alcouffe et al [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] is employed to construct transfer operators as well as coarse grid operator. Finally the Multi-Grid method [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF] is implemented.

Chapter 3: In this chapter, lubricant rheological models are illustrated in the first place. The relative friction coefficient, an indicator for the full-film-mixed lubrication regime transition, is proposed in methodology section. Then the relative friction coefficient is calculated numerically for isotropic as well as anisotropic harmonic surface roughness respectively. Finally, a single friction curve is obtained using a new "lambda ratio" parameter.

Chapter 4: A rapid analytical prediction method using the power spectral density [START_REF] Jacobs | Quantitative characterization of surface topography using spectral analysis[END_REF] is proposed to study a more complex surface topography in this chapter, firstly. Then an artificial surface roughness is employed to test this rapid prediction method. Finally, the prediction method is applied to predict friction for measured rough surfaces.

Introduction

Multi-grid methods have been used successfully to treat Elastohydrodynamic lubrication (EHL) problems in the past [START_REF] Venner | Advanced multilevel solution of the ehl line contact problem[END_REF][START_REF] Venner | Higher-order multilevel solvers for the ehl line and point contact problem[END_REF][START_REF] Venner | Multigrid techniques: a fast and efficient method for the numerical simulation of elastohydrodynamically lubricated point contact problems[END_REF]. However, when taking the surface roughness into account, film thickness and viscosity jump violently, both of them are strongly discontinuous parameters in discrete equations and will influence code robustness and code efficiency. The paper by R. Alcouffe [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] proposed an efficient way to solve this problem through constructing the coarse grid in a Multi-grid code. In this chapter, the Multigrid method is applied to solve the transient EHL model, and the algorithm outlined in Reference [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] is also implemented.

Transient EHL model 2.2.1 Governing equations

The lubrication of rough surfaces in EHL contacts is inherently a highly transient process. Study [START_REF] Ai | A transient ehl analysis for line contacts with measured surface roughness using multigrid technique[END_REF] shows that the surface roughness induced by the transient effect has a remarkable influence on the pressure and film thickness profiles. For the time-dependent problem [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF], the Reynolds equation is given as:

∂ ∂x ( ρh 3 12η ∂p ∂x ) + ∂ ∂y ( ρh 3 12η ∂p ∂y ) poiseuille -u r ∂(ρh) ∂x couette - ∂(ρh) ∂t transient = 0 (2.1)
with p = 0 on the boundaries and the cavitation condition p 0 everywhere. Where p is the pressure, h is the film thickness whose expression is shown as Equation 2.2 and u r = (u 1 +u 2 )/2 is the mean velocity (u 1 and u 2 are the velocities of two contact surfaces respectively). The CHAPTER 2. NUMERICAL MODEL direction of the x axis is as same as that of the mean velocity u r , the y axis is perpendicular to x axis and t is time.

The equation used to describe the gap between the two contact bodies is the film thickness equation:

h(x, y, t ) = h 0 (t ) + x 2 2R x + y 2 2R y -r r (x, y, t ) + 2 πE +∞ -∞ +∞ -∞ p(x , y , t ) (x -x ) 2 + (y -y ) 2 d x d y elastic deformation (2.2)
in which r r (x, y, t ) stands for surface roughness. R x and R y represent the reduced radius of curvature in x and y direction respectively. h 0 denotes the rigid body approach. E is called the reduced elastic modulus and its expression can be found below. The elastic deformation term is calculated with the approach named multilevel multi-integration [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF][START_REF] Lubrecht | A fast solution of the dry contact problem and the associated sub-surface stress field, using multilevel techniques[END_REF].

2

E = 1 -v 2 1 E 1 + 1 -v 2 2 E 2
and E 1 and E 2 are the elastic moduli of the two contact bodies. v 1 and u 2 are the Poisson ratios.

In order to have a load balance. The integral of the pressure distribution should be equal to the applied load w.

+∞ -∞ +∞ -∞ p(x , y , t )d x d y = w(t ) (2.3)
In the Reynolds equation (2.1), ρ is the density and η is the viscosity of the lubricant. Both of them are functions of pressure. A simply density pressure relation is given by Dowson and Higginson [START_REF] Dowson | Elasto-hydrodynamic lubrication: the fundamentals of roller and gear lubrication[END_REF]:

ρ(p) = ρ 0 5.9 × 10 8 + 1.34p 5.9 × 10 8 + p (2.4)
where ρ 0 is the atmospheric density. The simplest viscosity pressure relation is proposed by Barus [115]:

η(p) = η 0 exp(αp) (2.5)
in which η 0 is the atmospheric viscosity and α is the pressure viscosity coefficient. However, this exponential Barus equation usually predicts a higher viscosity value when the pressure is very large. A more realistic relation is derived by Roelands [START_REF] Roelands | Correlational aspects of the viscosity-temperature-pressure relationship of lubricating oils[END_REF]:

η(p) = η 0 exp[(ln(η 0 ) + 9.67)(-1 + (1 + p p 0 ) z )] (2.6)
where η 0 is the atmospheric viscosity and z is the pressure viscosity index, typically z = 0.6 and p 0 = 1.98 × 10 8 Pa.

Dimensionless equations and parameters

To simplify the equation system and generalize the EHL model, the equations described above are made dimensionless using dimensionless variables based on the Hertzian dry contact solution [START_REF] Johnson | Contact Mechanics[END_REF]. For the dry point contact case, the pressure distribution profile required for contact deformation reads:

p(x, y) = p h 1 -(x/a h ) 2 -(y/a h ) 2 if x 2 + y 2 ≤ a 2 h 0 otherwise (2.7) CHAPTER 2.

NUMERICAL MODEL

with a h the radius of the contact area:

a h = 3 3wR x 2E (2.8)
and p h is referred to as the Hertzian pressure:

p h = 3w 2πa 2 h .
(2.9)

Then the dimensionless variables are introduced to simplify the EHL model:

X = x/a h Y = y/a h P = p/p h H = hR x /a 2 h η = η/η 0 ρ = ρ/ρ 0 T = u r t /a h ᾱ = αp h (2.10)
Substituting the dimensionless variables in Equation 2.1 yields:

∂ ∂X ( ρH 3 η λ ∂P ∂X ) + ∂ ∂Y ( ρH 3 η λ ∂P ∂Y ) - ∂( ρH ) ∂X - ∂( ρH ) ∂T = 0 (2.11) with X ∈ [X a , X b ] and Y ∈ [-Y a , Y a ]. Where λ = (12u r η 0 R 2 x )/(a 3 p h ). And the boundary condi- tions are P (X a , Y a ) = P (X a , Y ) = P (X b , Y ) = P (X , Y a ) = P (X , -Y a ) = 0. The cavication condition is P (X , Y , T ) ≥ 0.
The dimensionless film thickness equation becomes:

H (X , Y , T ) = H 0 (T ) + X 2 2 + Y 2 2 -RR(X , Y , T ) + 2 π 2 +∞ -∞ +∞ -∞ P (X , Y , T ) (X -X ) 2 + (Y -Y ) 2 d X d Y (2.
12) where H 0 (T ) is determined by the dimensionless force balance equation:

+∞ -∞ +∞ -∞ P (X , Y , T )d X d Y = 2π 3 (2.13)
The dimensionless density equation for a compressible lubricant reads: ρ(P ) = 5.9 × 10 8 + 1.34p h P 5.9 × 10 8 + p h P (2.14)

The dimensionless forms of viscosity equations are:

Barus: η = exp( ᾱP ) (2.15)
and Roelands:

η = exp((ln(η 0 ) + 9.67)(-1 + (1 + p h p 0 P ) z )). (2.16)
Beside the dimensionless variables mentioned in Equation 2.10, two dimensionless numbers are often used to reduce the number of parameters, they are referred as Moes dimensionless parameters [START_REF] Moes | Discussion on a contribution by K. Jakobsen and H. Christensen[END_REF][START_REF] Moes | Design charts for optimum bearing configurations: 1-the full journal bearing[END_REF]. For point contact they are defined as [START_REF] Moes | Design charts for optimum bearing configurations: 1-the full journal bearing[END_REF]:

M = w E R 2 x ( 2η 0 u r E R x ) -3/4 (2.17) CHAPTER 2. NUMERICAL MODEL and L = αE ( 2η 0 u r E R x ) 1/4 . (2.18)
For convenience, the Moes parameters can be used to re-write the parameters ᾱ and λ:

λ = ( 128π 3 3M 4 ) 1/3 ᾱ = L π ( 3M 2 ) 1/3 (2.19)
Hamrock and Dowson [START_REF] Hamrock | Isothermal elastohydrodynamic lubrication of point contacts: Part 1-theoretical formulation[END_REF] introduced three parameters to simplify the study of film thickness. For point contact they are written as:

W = w E R 2 x U = η 0 u r E R x G = αE (2.20)

The finite difference scheme

The second-order self-adjoint elliptic partial differential equation considered by Alcouffe [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] is

-∇ • (D(x, y, t )∇U (x, y, t )) + σ(x, y, t )U (x, y, t ) = f (x, y, t ) (x, y) ∈ Ω (2.
- ∂ ∂X (D ∂P ∂X ) - ∂ ∂Y (D ∂P ∂Y ) = ∂( ρH ) ∂X + ∂( ρH ) ∂T (X , Y ) ∈ Ω (2.22)
In the present work, the calculation domain

Ω is a rectangle [X a , X b ] × [-Y a , Y a ].
This domain is covered with a uniform grids with a system of straight lines parallel to the coordinate axes. The mesh size in the two directions is hx = (X b -X a )/N x and h y = 2×Y a /N y , in which N x and N y are the number of mesh points in both directions. Figure 2.1: Mesh point (x i , y j ) and it's related mesh region r i , j CHAPTER 2. NUMERICAL MODEL To derive the difference scheme, a mesh region r i , j (shown as Figure 2.1) defined by the lines x = x ih x /2 , x = x i + h x /2 , y = y jh y /2 and y = y j + h y /2 for each mesh point (x i , y j ) is selected. In terms of each mesh point (x i , y j ), P (x i , y j , t k ) := P i , j ,k is unknown, now integrating Equation 2.22 over the corresponding mesh region r i , j :

-

r i , j [ ∂ ∂X (D ∂P ∂X ) + ∂ ∂Y (D ∂P ∂Y )]d xd y = r i , j [ ∂( ρH ) ∂X + ∂( ρH ) ∂T ]d xd y (2.23)
According to Green's Theorem [START_REF] Riley | Mathematical methods for physics and engineering[END_REF], Equation 2.23 can be expressed as:

- c i , j [(D ∂P ∂X )d y -(D ∂P ∂Y )d x] = r i , j [ ∂( ρH ) ∂X + ∂( ρH ) ∂T ]d xd y (2.24)
where c i , j is the boundary of r i , j and the integration path along this boundary is anticlockwise. Supposing f (x i , y j , t k ) := f i , j ,k , the double integrals of the right hand side of Equation 2.24 can be simply approximated by means of

r i , j f (x, y, t )d xd y . = f i , j ,k • a i , j . (2.25)
where a i , j = hx• h y is the area of the rectangle region r i , j shown in Figure 2.2. Referring again to Figure 2.2, the line integal of Equation 2.24 over the four boundaries of r i , j is approximated by means of central differences as:

- c i , j [(D ∂P ∂X )d y -(D ∂P ∂Y )d x] . = (h y)[D i +1/2, j ,k ( P i , j ,k -P i +1, j ,k hx ) + D i -1/2, j ,k ( P i , j ,k -P i -1, j ,k hx )] + (hx)[D i , j +1/2,k ( P i , j ,k -P i , j +1,k h y ) + D i , j -1/2,k ( P i , j ,k -P i , j -1,k h y )] (2.
A i , j ,k (P i , j +1,k -P i , j ,k )+A i , j -1,k (P i , j -1,k -P i , j ,k )+ B i , j ,k (P i +1, j ,k -P i , j ,k ) + B i , j -1,k (P i -1, j ,k -P i , j ,k ) = F i , j ,k (2.27)
where

A i , j ,k = -( 1 2 )( hx h y )(D i , j ,k + D i , j +1,k ) B i , j ,k = -( 1 2 )( h y hx )(D i , j ,k + D i +1, j ,k ) F i , j ,k = (hx • h y) f i , j ,k
In terms of f i , j ,k , the same discrete schemes used in Reference [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF] is adopted. At this point, the right hand side of Equation 2.27 can be taken as:

F i , j ,k . = h y(1.5 ρi,j,k H i , j ,k -2 ρi-1,j,k H i -1, j ,k + 0.5 ρi-2,j,k H i -2, j ,k ) + hx • h y ht (1.5 ρi,j,k H i , j ,k -2 ρi,j,k-1 H i , j ,k-1 + 0.5 ρi,j,k-2 H i , j ,k-2 ) (2.28)
where ht is the mesh size in time domain. A more detailed derivation of the above difference scheme can be found in Reference [START_REF] Mcguire | Matrix structural analysis[END_REF].

Transfer operators

Intergrid transfers are used for connecting the fine grid with the coarse grid. After a number of relaxations the error on the fine grid is smooth enough to be approximate on the coarse grid. Hence a restriction operator I H h is needed to transfer the approximated solution P h and the residual r h . When the low frequency errors have been eliminated on the coarse grid, it is necessary to define a new error υ h (υ h = P h -P h ) on the fine grid to correct the fine grid approximate solution P h . The classical bi-linear interpolation works quite well for most load cases. However when D jumps by orders of magnitude, Alcouffe [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] proposed a more efficient interpolation operator and this type of operator allows D∇P to be continuous over the whole calculation domain and gives a more reasonable physical representation on the coarse grid [START_REF] Noutary | A robust piston ring lubrication solver: Influence of liner groove shape, depth and density[END_REF].

Interpolation

Having defined the coefficients A and B in Equation 2.27, it is time to define the interpolation operator. In matrix form, the interpolation can be represented as:

υ h = I h H υ H (2.29)
where υ h and υ H are the fine grid and coarse grid error vectors respectively. I h H is the interpolation operator and the superscripts h and H stand for the fine grid and the coarse grid respectively. The new coarse grid construction method proposed by Alcouffe et al. [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] is used here, the interpolation process will be illustrated as follows: The first step is to interpolate the fine grid points (black points shown in Figure 2.3 (b)) coinciding with the coarse grid points (green points shown in Figure 2.3 (a)):

υ h i F, j F,k = υ H iC , jC ,k (2.30)
where υ h is the error on the fine grid, and υ H is the error on the coarse grid. Subscripts (i F, j F, k) and (iC , jC , k) are applied for illustrate the mesh points on the fine grid and on the coarse grid at the k th time step respectively. The second step is to obtain the middle points, represented as blue dots in Figure 2.3 (c), on the fine grid. Along horizontal lines, the expression for middle points is:

υ h i F +1, j F,k = (B h i F, j F,k υ H iC , jC ,k + B h i F +1, j F,k υ H iC +1, jC ,k ) (B h i F, j F,k + B h i F +1, j F,k ) (2.31)
A similar expression can be derived for vertical lines, which reads:

υ h i F, j F +1,k = (A h i F, j F,k υ H iC , jC ,k + A h i F, j F +1,k υ H iC +1, jC ,k ) (A h i F, j F,k + A h i F, j F +1,k ) (2.32)
Finally, the central point represented as a red point in Figure 2.3 (d) on the fine grid, which is obtained as:

υ h i F +1, j F +1,k = (A h i F +1, j F +1,k υ h i F +1, j F +2,k + A h i F +1, j F,k υ h i F +1, j F,k + B h i F, j F +1,k υ h i F, j F +1,k + B h i F +1, j F +1,k υ h i F +1, j F +1,k )/ (A h i F +1, j F +1,k + A h i F +1, j F,k + B h i F, j F +1,k + B h i F +1, j F +1,k ) (2.33)
The above pointwise description (from Equation 2.30 to Equation 2.33) can be replaced by the matrix expression Equation 2.29, in which the matrix is large and complex. A simply way to describe this matrix is by using a stencil notation. As was shown in Figure 2.4, in the interpolation process, the stencil provides weighting factors for dividing the coarse grid value in point (iC , jC , k) to the coinciding fine grid point (i F, j F, k) as well as its 8 adjacent points. Observing those pointwise expressions, the contribution of the coarse grid point to the 9 corresponding fine grids can be written as a stencil I h H in Equation 2. [START_REF] Venner | Surface roughness effects in an ehl line contact[END_REF].

I h H =         NW h i F, j F,k N h i F, j F,k N E h i F, j F,k W h i F, j F,k C h i F, j F,k E h i F, j F,k SW h i F, j F,k S h i F, j F,k SE h i F, j F,k         (2.34)
where CHAPTER 2. NUMERICAL MODEL

C h i F, j F,k = 1, N h i F, j F,k = A h i F, j F,k A h i F, j F,k + A h i F, j F +1,k , E h i F, j F,k = B h i F, j F,k B h i F, j F,k + B h i F +1, j F,k , S h i F, j F,k = A h i F, j F -1,k A h i F, j F -2,k + A h i F, j F -1,k , W h i F, j F,k = B h i F -1, j F,k B h i F -1, j F,k + B h i F -2, j F,k
N E h i F, j F,k = [B h i F, j F +1,k A h i F, j F,k /(A h i F, j F,k + A h i F, j F +1,k )] A h i F +1, j F +1,k + A h i F +1, j F,k + B h i F +1, j F +1,k + B h i F, j F +1,k + [B h i F, j F,k A h i F +1, j F,k /(B h i F, j F,k + B h i F +1, j F,k )] A h i F +1, j F +1,k + A h i F +1, j F,k + B h i F +1, j F +1,k + B h i F, j F +1,k , NW h i F, j F,k = [B h i F -1, j F +1,k A h i F, j F,k /(A h i F, j F,k + A h i F, j F +1,k )] A h i F -1, j F +1,k + A h i F -1, j F,k + B h i F -2, j F +1,k + B h i F -1, j F +1,k + [B h i F -1, j F,k A h i F -1, j F,k /(B h i F -1, j F,k + B h i F -2, j F,k )] A h i F -1, j F +1,k + A h i F -1, j F,k + B h i F -2, j F +1,k + B h i F -1, j F +1,k , SE h i F, j F,k = [B h i F, j F -1,k A h i F, j F -1,k /(A h i F, j F -1,k + A h i F, j F -2,k )] A h i F +1, j F -1,k + A h i F +1, j F -2,k + B h i F, j F -1,k + B h i F +1, j F -1,k + [B h i F, j F,k A h i F +1, j F -1,k /(B h i F, j F,k + B h i F +1, j F,k )] A h i F +1, j F -1,k + A h i F +1, j F -2,k + B h i F, j F -1,k + B h i F +1, j F -1,k , SW h i F, j F,k = [B h i F -1, j F -1,k A h i F, j F -1,k /(A h i F, j F -1,k + A h i F, j F -2,k )] A h i F -1, j F -1,k + A h i F -1, j F -2,k + B h i F -2, j F -1,k + B h i F -1, j F -1,k + [B h i F -1, j F,k A h i F -1, j F -1,k /(B h i F -1, j F,k + B h i F -2, j F,k )] A h i F -1, j F -1,k + A h i F -1, j F -2,k + B h i F -2, j F -1,k + B h i F -1, j F -1,k
.

Injection

In general, the restriction operator matrix is the transposed matrix of the interpolation operator [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF]:

I H h = (I h H ) T (2.35)
In order to derive the restriction operator, a basis function e l i , j whose value is 1 at the point (i , j ) on the l th grid is employed.

W H = (I H h W h ) (iC , jC ) = [(I h H ) T W h ] (iC , jC ) = i F j F 〈I h H e H iC , jC , e h i F, j F 〉 〈e H iC , jC , e H iC , jC 〉 W h i F, j F (2.36)
in which W H and W h are vectors of unknowns on the coarse and fine grid respectively. In terms of the orthogonal basis e l i , j , the dot products are 〈e H iC , jC , e H iC , jC 〉 = HxHy and 〈e h i F, j F , e h i F, j F 〉 = hxhy.

Supposing the vectors of unknowns on coarse grid and on fine grid are marked as

W H = [u H iC , jC ,k ] and W h = [u h i F -1, j F +1,k , u h i F -1, j F,k , u h i F -1, j F -1,k , u h i F, j F +1,k , u h i F, j F,k , u h i F, j F -1,k , u h i F +1, j F +1,k , u h i F +1, j F,k , u h i F +1, j F -1,k ].
According to Equation 2.36, one can obtain:

u H iC , jC ,k = hx • h y H x • H y (C h i F, j F,k u h i F, j F,k + NW h i F, j F,k u h i F -1, j F +1,k + SW h i F, j F,k u h i F -1, j F -1,k + SE h i F, j F,k u h i F +1, j F -1,k + N E h i F, j F,k u h i F +1, j F +1,k + N h i F, j F,k u h i F, j F +1,k + S h i F, j F,k u h i F, j F -1,k + W h i F, j F,k u h i F -1, j F,k + E h i F, j F,k u h i F +1, j F,k ).
(2.37) CHAPTER 2. NUMERICAL MODEL Hence, the stencil of the restriction operator I H h is:

I H h = hx • h y H x • H y         nw h i F, j F,k n h i F, j F,k ne h i F, j F,k w h i F, j F,k c h i F, j F,k e h i F, j F,k sw h i F, j F,k s h i F, j F,k se h i F, j F,k         (2.38)
where

c h i F, j F,k = 1,
n h i F, j F,k = A h i F, j F,k A h i F, j F,k + A h i F, j F +1,k , e h i F, j F,k = B h i F, j F,k B h i F, j F,k + B h i F +1, j F,k
,

s h i F, j F,k = A h i F, j F -1,k A h i F, j F -2,k + A h i F, j F -1,k , w h i F, j F,k = B h i F -1, j F,k B h i F -1, j F,k + B h i F -2, j F,k , ne h i F, j F,k = [B h i F, j F +1,k A h i F, j F,k /(A h i F, j F,k + A h i F, j F +1,k )] A h i F +1, j F +1,k + A h i F +1, j F,k + B h i F +1, j F +1,k + B h i F, j F +1,k + [B h i F, j F,k A h i F +1, j F,k /(B h i F, j F,k + B h i F +1, j F,k )] A h i F +1, j F +1,k + A h i F +1, j F,k + B h i F +1, j F +1,k + B h i F, j F +1,k , nw h i F, j F,k = [B h i F -1, j F +1,k A h i F, j F,k /(A h i F, j F,k + A h i F, j F +1,k )] A h i F -1, j F +1,k + A h i F -1, j F,k + B h i F -2, j F +1,k + B h i F -1, j F +1,k + [B h i F -1, j F,k A h i F -1, j F,k /(B h i F -1, j F,k + B h i F -2, j F,k )] A h i F -1, j F +1,k + A h i F -1, j F,k + B h i F -2, j F +1,k + B h i F -1, j F +1,k , se h i F, j F,k = [B h i F, j F -1,k A h i F, j F -1,k /(A h i F, j F -1,k + A h i F, j F -2,k )] A h i F +1, j F -1,k + A h i F +1, j F -2,k + B h i F, j F -1,k + B h i F +1, j F -1,k + [B h i F, j F,k A h i F +1, j F -1,k /(B h i F, j F,k + B h i F +1, j F,k )] A h i F +1, j F -1,k + A h i F +1, j F -2,k + B h i F, j F -1,k + B h i F +1, j F -1,k , sw h i F, j F,k = [B h i F -1, j F -1,k A h i F, j F -1,k /(A h i F, j F -1,k + A h i F, j F -2,k )] A h i F -1, j F -1,k + A h i F -1, j F -2,k + B h i F -2, j F -1,k + B h i F -1, j F -1,k + [B h i F -1, j F,k A h i F -1, j F -1,k /(B h i F -1, j F,k + B h i F -2, j F,k )] A h i F -1, j F -1,k + A h i F -1, j F -2,k + B h i F -2, j F -1,k + B h i F -1, j F -1,k . CHAPTER 2. NUMERICAL MODEL
Subsequently, the solution P H and the residual r H on the coarse grid are restricted as:

P H iC , jC ,k = I H h P h i F, j F,k = hx • h y H x • H y [c h i F, j F,k P h i F, j F,k + n h i F, j F,k P h i F, j F +1,k + s h i F, j F,k P h i F, j F -1,k + w h i F, j F,k P h i F -1, j F,k + e h i F, j F,k P h i F +1, j F,k + nw h i F, j F,k P h i F -1, j F +1,k + sw h i F, j F,k P h i F -1, j F -1,k + se h i F, j F,k P h i F +1, j F -1,k + ne h i F, j F,k P h i F +1, j F +1,k ].
(2.39) Special attention is needed when restricting the residual r h . The factor (hx • h y)/(H x • H y) in Equation 2.38 is used to restrict functions mathematically, which is not used to restrict flows physically. For Reynolds equation, the restriction operator I H h must be replaced by J H h when restricting the right hand side of the discrete Reynolds equation [START_REF] Noutary | A robust Reynolds solver for textured surfaces in the piston ring cylinder liner contact[END_REF] .

J H h = H x • H y hx • h y I H h =         nw h i F, j F,k n h i F, j F,k ne h i F, j F,k w h i F, j F,k c h i F, j F,k e h i F, j F,k sw h i F, j F,k s h i F, j F,k se h i F, j F,k         (2.40) yielding: r H iC , jC ,k = J H h r h i F, j F,k . (2.41)

Coarse grid operator

After defining the transfer operators, the coarse grid operator will be described here. The coarse grid operator is formed by a restriction equation:

J H h (L h υ h ) = J H h r h (2.42)
In which υ h = I h H υ H and equation 2.42 can be represented as:

J H h L h I h H (υ H ) = J H h r h (2.43)
Thus, the coarse grid operator can be defined as:

L H = J H h L h I h H (2.44)
Assuming the stencil of L l at point (i , j ) on level l and at k t h time step is:

L l =        L l nw (i , j , k) L l n (i , j , k) L l ne (i , j , k) L l w (i , j , k) L l c (i , j , k) L l e (i , j , k) L l sw (i , j , k) L l s (i , j , k) L l se (i , j , k)        (2.45)
The construction of the coarse grid operator can be seen in Appendix A.

The coarse grid operator defined here involves 9 points, while as equation 2.27 shows an operator referring to 5 points on the finest grid, reads:

L h (i F, j F, k) =         0 A h i F, j F,k 0 B h i F -1, j F,k -A h i F, j F,k -A h i F, j F -1,k -B h i F -1, j F,k -B h i F, j F,k B h i F, j F,k 0 A h i F, j F -1,k 0         (2.46) CHAPTER 2. NUMERICAL MODEL

Relaxation

Employing the grid operator L l defined before, the discrete Reynolds equation for point (i , j , k) of the grid l at the k t h time step is shown as follows:

L l c P l i , j ,k + L l n P l i , j +1,k + L l w P l i -1, j ,k + L l s P l i , j -1,k + L l e P l i +1, j ,k + L l nw P l i -1, j +1,k + L l sw P l i -1, j -1,k + L l se P l i +1, j -1,k + L l ne P l i +1, j +1,k = (h l • h l )( ρH ) l x + (h l • h l )( ρH ) l t (2.47)
with P l i , j ,k = 0 for points on the boundary and P l i , j ,k > 0 for the cavitation condition. In Equation 2.47, h l stands for the mesh size on the l t h grid. In terms of the couette and the transient term, their discrete forms are:

( ρH ) l x =      ρl i , j ,k H l i , j ,k -ρl i -1, j ,k H l i -1, j ,k h l if i = 1 1.5 ρl i , j ,k H l i , j ,k -2.0 ρl i -1, j ,k H l i -1, j ,k +0.5 ρl i -2, j ,k H l i -2, j ,k h l if i 2 (2.48) and ( ρH ) l t =            0 if nt = 0 ρl i , j ,k H l i , j ,k -ρl i -1, j ,k H l i -1, j ,k ht if nt = 1 1.5 ρl i , j ,k H l i , j ,k -2.0 ρl i -1, j ,k H l i -1, j ,k +0.5 ρl i -2, j ,k H l i -2, j ,k ht if nt 2 (2.49)
with

H l i , j ,k = H 0 (k) + X 2 i ,k 2 + Y 2 j ,k 2
+ i j K l l i ,i , j , j ,k P l i , j ,k (2.50)
For the discrete Equation 2.47, Reference [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF] shows that an iterative approach of combining the Gauss-Seidel line relaxation and the Jacobi distribution line relaxation is pretty stable and efficient. Hence, in this work the same method is applied. Figure 2.5 shows the relaxation process represented in Reference [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF], where ξ l is the local coefficient and its definition can been seen as Equation 2.51.

ξ l i , j ,k = ρ(P l i , j ,k )(H l i , j ,k ) 3 η(P l i , j ,k ) λ (2.51) with ξ l i ±1/2, j ,k = (ξ i , j ,k + ξ i ±1, j ,k )/2 and ξ l i , j ±1/2,k = (ξ i , j ,k + ξ i , j ±1,k )/2.
Whether the Gauss-Seidel or the Jacobi distribution line relaxation is used, for each grid point (i , j , k) a new approximation P l i , j ,k to P l i , j ,k is computed by:

P l i , j ,k = P l i , j ,k + ωδ l i , j ,k (2.52)
with

ωδ l i , j ,k =    ω g s δ l i , j ,k
Guass-Seidel

ω j a [δ l i , j ,k - (δ l i +1, j ,k +δ l i -1, j ,k +δ l i , j +1,k +δ l i , j -1,k ) 4
] Jacobi distribution where ω g s and ω j a are the relaxation factors for the Guass-Seidel line relaxation and the Jacobi distribution line relaxation respectively. The changes δ l i , j ,k for the line relaxation in the X direction can be obtained from:

A j δ l j = r l j (2.53) CHAPTER 2. NUMERICAL MODEL in which δ l j is a vector of changes δ l i , j ,k and r l j is the residual vector r l i , j ,k , both are n x element vectors. A j is a matrix of coefficients A j i ,m , whose derivation can be found in Appendix B.

Implementation of the Multi-Grid method

So far, pre-preparations for implementing the Multi-Grid method to solve the transient EHL equation system have already been prepared. Final step is to organize all steps together. Figure 2.6 shows a simple two-level "V" cycle which is used to illustrate the implementation of the Multi-Grid method. The corresponding steps are as follows:

• Several relaxations on grid h to obtain an approximate solution P h by:

L h 〈 P h 〉 = F h , (2.54)
compute residual on grid h below:

r h = F h -L h 〈 P h 〉. (2.55)
• Compute Alcouffe's coefficients A, B and D, prepare the transfer operators I H h and I h H , construct the coarse grid operator L H . CHAPTER 2. NUMERICAL MODEL Figure 2.6: Implementation of the Multi-Grid method with a two level "V" cycle.

• Coarsen pressure P h and residual r h to grid H :

P H = I H h 〈 P h 〉 r H = J H h 〈r h 〉 (2.56)
with the right hand side term F H expressed as:

F H = L H 〈 P h 〉 + r H . ( 2.57)
• Relax coarse grid solution P H on grid H through:

L H 〈 P H 〉 = F H . (2.58)
if grid H is the coarsest grid, the mutual approach H 0 should be updated.

• Interpolate and correct P h from grid H to h using:

P h = P h + I h H 〈 P H -I H h P h 〉. (2.59)
• Relaxations on grid h.

CHAPTER 2. NUMERICAL MODEL

For simplicity, the current code uses "V" cycles to solve the EHL equation system. In terms of the transient EHL model, the transient term ∂ ρH /∂X should be arranged in calculation cycles.

Here use the equation below to denote the transient Reynolds equation with the first order backward discretization of transient term:

L l 〈P l 〉 -(h l • h l )( ρH ) l x = (h l • h l )( ρl i , j ,k H l i , j ,k -ρl i -1, j ,k H l i -1, j ,k ht ).
(2.60)

In order to construct the transient "V" cycle, Equation 2.60 will be rewritten as:

L l 〈P l 〉 -(h l • h l )( ρH ) l x -(h l • h l )( ρl i , j ,k H l i , j ,k ht ) = (h l • h l ) (-ρl i -1, j ,k H l i -1, j ,k ) ht . (2.61)
Observing Equation 2.61, one finds that the right hand side term of the Reynolds equation is

(h l •h l ) (-ρl i -1, j ,k H l i -1, j ,k ) ht
. As shown in Figure 2.7, the right hand side term of the Reynolds equation for the (nt ) t h time step is the value of (

h l • h l ) (-ρl i -1, j ,k H l i -1, j ,k ) ht
for the (nt -1) t h time step.

Conclusion

In this chapter, the coarse grid construction method proposed by Alcouffe et al. [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] is implemented in the current EHL Multi-Grid code. Instead of the bi-linear interpolation stencil used in the class code [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF], the new interpolation stencil applying Alcouffe's method is derived. Then the new restriction stencil is obtained through the Equation 2.35. A nine point coarse grid operator is also constructed by the Galerkin method in Appendix A. Finally, the implementation of the Multi-Grid method is introduced briefly.

Chapter 3

Friction influence of harmonic surface waviness

Introduction

As mentioned in Chapter 1, the friction of interfacial surfaces greatly influences the performance of mechanical elements. Studies show that the surface roughness amplitude, wavelength and anisotropy affect friction. Thus this chapter investigates the effect of surface waviness on the friction evolution, in which the surface waviness is defined by Equation (3.1). In this chapter, a relative friction coefficient is proposed to indicate the transition from the full-film to the mixed lubrication regime. The definition of the relative friction coefficient is illustrated in Section 3.3. Subsequently, the influence of operating conditions and surface waviness topography on friction are studied. Finally, a single friction curve is obtained depending on a new "Lambda ratio" parameter including operating conditions as well as surface waviness parameters. The waviness is defined as:

RR(X , Y , T ) = Ai × 10 -10[max(0, X - X λx /a h ) 2 ] cos(2π (X -X ) λ x /a h )cos(2π Y λ y /a h ) (3.1)
where X = X st +U rat × T with U rat = u 1 /u r , Ai is the initial amplitude of the surface waviness, λ x and λ y are the wavelengths in x and y direction respectively. At the same time, the parameter r (r = λ x /λ y ) is used to determine the surface waviness anisotropy. For convenience, the parameter λ is defined as λ = min(λ x , λ y ). The exponential term is used to avoid discontinuous derivatives when the waviness moves into the calculation domain.

Lubricant rheological models

The lubricant rheology greatly determines the friction in a tribological contact [START_REF] Ponjavic | Local rheology of lubricants in the elastohydrodynamic regime[END_REF]. Once the rheology of the lubricant is given, the friction force can be obtained by integration of the CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS shear stress over the contact area [START_REF] Olver | Prediction of traction in elastohydrodynamic lubrication[END_REF]. Figure 3.1 shows a typical shear stress-shear rate relationship for an EHL contact. In the low shear strain rate range, the lubricant is Newtonian and the shear stress increases linearly with the shear strain rate. This linear model is valid for predicting hydrodynamic friction [START_REF] Bair | The high pressure high shear stress rheology of liquid lubricants[END_REF]. However, in the EHL regime of lubrication, the Newtonian assumption is insufficient, especially when large slide-to-roll ratios are employed. Studies [START_REF] Bair | Shear thinning correction for rolling/sliding elastohydrodynamic film thickness[END_REF][START_REF] Kumar | Full ehl simulations using the actual ree-eyring model for shear-thinning lubricants[END_REF][START_REF] Bair | The rheological assumptions of classical ehl: What went wrong?[END_REF] show that the lubricant in the inlet zone shows shear-thinning behavior (i.e. the fluid viscosity decreases with increasing shear stress rate) which influences the film thickness. Some frequently-used non-Newtonian fluid models are listed as follows:

γ = 1 G e d τ d t + τ 0 η sinh( τ τ 0 ) (3.2)
in which G e is the elastic shear modulus, τ is the stress, τ 0 is referred to as the Eyring stress and γ is the shear strain rate. The strain rate shown as Equation 3.2 consists of two components where the first elastic term is always neglected [START_REF] Huge | History, origins and prediction of elastohydrodynamic friction[END_REF]. Thus Equation 3.2 reduces to:

γ = τ 0 η sinh( τ τ 0 ). (3.3)
Observing Equation 3.3, when the Eyring stress τ 0 approaches infinity, the limit of Equation 3.3 becomes the Newtonian constitutive equation, i.e.: lim

τ 0 →∞ τ 0 η sinh( τ τ 0 ) = τ η . ( 3.4)
•

Bair-Winer model

For larger strain rates the lubricant exhibits a limiting shear stress τ L which is the threshold of the shear stress. Bair and Winer [START_REF] Bair | A rheological model for elastohydrodynamic contacts based on primary laboratory data[END_REF] modified the "Maxwell" model applying this limiting shear stress, yields:

γ = 1 G ∞ d τ d t + τ L η ln(1 - τ τ L ) (3.5) CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS
where G ∞ is the limiting elastic shear modulus.

• Carreau-Yasuda model (power-law model)

Carreau and Yasuda [START_REF] Yasuda | Shear flow properties of concentrated solutions of linear and star branched polystyrenes[END_REF] offered a more general model whose equation is:

η = η 0 [1 + (λ t γ) 2 ] (n-1)/2 (3.6)
where λ t is the time constant for the fluid and n -1 is the power-law slope.

• Circular fluid model

In order to overcome the difficulty lies in incorporating the Bair-Winer model into the Reynolds equation, Lee and Hamrock [START_REF] Lee | A circular non-newtonian fluid model: Part i-used in elastohydrodynamic lubrication[END_REF] suggested an appropriate lubricant rheological circular model, it reads:

γ = τ η [1 -( τ τ L ) 2 ] -1/2 (3.7)
• Actual Ree-Eyring model Bair [136] proposed an actual Ree-Eyring model for shear-shinning lubricants [START_REF] Kumar | On the role of lubricant rheology and piezo-viscous properties in line and point contact ehl[END_REF]:

τ = N i =1 x i τ i sinh -1 (λ i γ) (3.8)
in which λ i = η/τ i is a characteristic time of the fluid and x i is a weighting factor ( N i =1 = 1).

Methodology

Relative friction coefficient

Even though the previous section described several non-Newtonian rheological model, in this section, we will use a Newtonian model to describe friction variations resulting from pressure variations. Two arguments can be used.

• When the shear stress variations are small, a linearised model can be used. See the following example sinh.

• A linear model allows an FFT based sum over all wavelengths which a non-linear model does not.

Example:

γ = τ 0 η sinh τ τ 0 = τ 0 2η (e τ/τ 0 -e -τ/τ 0 ) = τ 0 2η (1 + τ τ 0 + τ 2 2τ 2 0 + ... -1 + τ τ 0 - τ 2 2τ 2 0 + ...) ≈ τ η CHAPTER 3.

FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

For Newtonian lubricant, the above equation can be rewritten as:

τ = η d u d z (3.9)
where u is the velocity in rolling direction, which is defined as [START_REF] Rahnejat | Tribology and dynamics of engine and powertrain: fundamentals, applications and future trends[END_REF][START_REF] Wen | Principles of tribology[END_REF]:

u = 1 2η ∂p ∂x (z 2 -zh) + z h (u 2 -u 1 ) + u 1 .
(3.10) Substituting Equation 3.10 into Equation 3.9, yields:

τ = 1 2 ∂p ∂x (2z -h) + η u 2 -u 1 h (3.11)
where two parts contribute to the shear stress, the first term is a parabolic part due to the Poiseuille flow and the second term is a linear part due to the Couette flow. Thus the viscous shear force on the lower surface (z = 0) is obtained by integrating Equation 3.11 [START_REF] Bonaventure | Influence of random surface roughness on friction in elastohydrodynamic, mixed and boundary lubrication[END_REF]:

f = - ( ∂p ∂x • h 2 )d xd y + (η • u 2 -u 1 h )d xd y. (3.12)
In the full-film EHL regime, the friction force is dominated by the viscous shear force, the dimensionless friction force is:

F friction = - ( ∂P (X , Y , T ) ∂X • H (X , Y , T ) 2 )d X d Y + (η • SRR H (X , Y , T ) )d X d Y (3.13)
in which SRR = (u 2u 1 )/ ū is the slide-to-roll ratio.

Experimental results [START_REF] Evans | Measurement and mapping of the rheological properties of elastohydrodynamic lubricants[END_REF] also found that the Poiseuille force is lower than the Couette force. Hence the dimensionless friction force is simplified as:

F (X , Y , T ) = η[P (X , Y , T )] • SRR H (X , Y , T ) d X d Y . (3.14)
Different lubricant oils give different curves of friction coefficient. For the sake of simplicity, a relative friction coefficient is proposed to detect friction variations, which is defined as:

µ r µ s (T ) = F r w (T )/ F s w = F r F s (T ) (3.15)
where subscripts r and s are used to distinguish the rough and smooth case.

Numerical solution

The numerical simulation uses the domain -2.5 X 1.5 and -2.0 Y 2.0 with 513 × 513 equal-spaced points. The time step is selected equal to the spatial mesh size on the finest grid, i.e. with ∆T = hx = h y = 0.0078125. Meanwhile, the calculation starts with X st = -2.5 and the surface topography moves into the high pressure zone with the velocity of the rough surface u 1 .

The monitoring time should be long enough so that 'steady oscillations' of the results occur. The present work considers small-amplitude roughness and a small slip parameter is selected i.e. U r at = 1.01(SRR = 0.02). This small slip assumption and and small amplitude allow us to use the Amplitude Reduction Theory [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF] for pure rolling, as shown in [START_REF] Šperka | Experimental study of real roughness attenuation in concentrated contacts[END_REF].

CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

The numerical code used in the present work is validated by comparing the current results with results from Reference [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF]. Figure 3.2 shows the relative deformed amplitude (Ad /Ai ) as a function of f (r )∇ 2 . This figure shows that results from the current model with SRR = 0 (blue squares) fall onto the master curve on Reference [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF]. The choice of the mesh size influences computing time as well as precision. A large mesh size leads to large discretization errors and a small mesh size causes long computation times. This work applies an intermediate mesh size (513 × 513 points). Table 3.1 suggests that it yields a friction precision better than 1%.

Time-dependent solution

Defined by Equation 3.1, Figure 3.3 shows the top view of the surface waviness for three cases: the isotropic surface waviness of r = 1 and λ/a h = 0.5 (a), the longitudinal surface waviness of r = 2 and λ/a h = 0.5 (b) and the transverse surface waviness of r = 0.5 and λ/a h = 0.5 (c). For the transient case, in the high pressure zone the Poiseuille term vanishes, the Reynolds equation reduces to a transport equation:

- ∂( ρH ) ∂X - ∂( ρH ) ∂T ≈ 0. (3.16)
In this case, the final film thickness and pressure profiles depend on the profiles of steady state case as well as on the inlet disturbances [START_REF] Morales | Elastohydrodynamic lubrication of smooth and rough surfaces[END_REF]. show that the propagation speed of the pressure increase is the same as that of the film thickness variation. This is due to the small slide-to-roll ratio, the velocity of the rough surface is nearly the same as the entrainment velocity ū. In addition, the wavelength of the pressure increase and the film thickness variation are the same as the wavelength of corresponding initial surface waviness.

Figure 3.7 shows the dimensionless central film thickness H c rough (0, 0, T ) as a function of time T for the above three rough cases. The initial amplitude of the surface waviness is Ai = 0.01215, while from those figures, one can observe that the amplitude of the central film thickness is smaller. This means that the surface roughness is deformed under current operating conditions. For those cases in Figure 3.7, the deformed amplitudes are: Ad = 0.0057 for r = 1, Ad = 0.00375 for r = 2 and Ad = 0.0057 for r = 0.5. Reference [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF] studied this deformation over a large range of operating conditions, and a single equation is obtained as Equation 3.17.

Ad Ai = 1 1 + 0.15 f (r )∇ 2 + 0.015( f (r )∇ 2 ) 2 (3.17)
where

f (r ) = e 1-1/r , if r > 1 1, otherwise r = λ x /λ y and ∇ 2 = (λ/a h ) M /L with λ = min(λ x , λ y ).

Effect of operating conditions

The relative friction coefficient can be plotted as a function of the classical parameter "lambda ratio" i.e. H c/Ai (Ai is varying and H c is fixed) for a specific operating condition: Figure 3.9. This figure shows that as H c/Ai increases, the relative coefficient decreases monotonically.

CHAPTER 3. FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS

Early work on conformal condition shows that with the "Λ ratio" increasing, the friction coefficient decreases to a minimum value (the transition position from the full-film to the mixed lubrication regime) first then increases. The following work will detect the influence of the operating conditions M , L and surface wavelength on the relative friction coefficient for the isotropic case. The relative friction coefficient is affected not only by the operating conditions M and L, but also by the wavelength of the surface waviness. Figure 3.12 describes the effect of the wavelength λ/a h on the relative friction coefficient for M = 1000 and L = 10. As λ/a h increases, the relative friction coefficient decreases, each single curve has a small left-shift (shown in Figure Figure 3.13 presents the relative friction coefficient as a function of the classical "lambda ratio" H c/Ai for many different operating conditions. This figure shows that for each operating condition, a very smooth curve is obtained, however, one does not obtain a single curve as for the low pressure case using the parameter H c/Ai . According to the Amplitude Reduction Theory [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF][START_REF] Venner | An engineering tool for the quantitative prediction of general roughness deformation in ehl contacts based on harmonic waviness attenuation[END_REF], for very high pressures, the surface roughness will deform. Instead of using this simple parameter Ai or a measured surface roughness parameter σ, it is better to use the deformed parameter Ad . where θ 2 = L -1.1 M 0.33 (λ/a h ) 0.67 (H c/Ai ). The physical justification of this scaling parameter can be understood from a simplified analysis given in Appendix C.

Effect of surface anisotropy

The previous section demonstrated the influence of the operating conditions and the wavelength on the relative friction coefficient for the isotropic case. Figure 3.8 shows that the surface anisotropy affects the relative friction coefficient under the same operating conditions. This section describes the influence of the surface anisotropy on the relative friction coefficient in more detail. in order to distinguish the surface roughness anisotropy function in Equation 3.17, in this work, the anisotropy function is marked as ff(r ). The transverse case is studied for M = 1000, L = 10 and λ x /a h = 0.5. Figure 3.17 shows the relative friction coefficient (µ r /µ s ) as a function of H c/Ai for different r (0 ≤ r ≤ 1) values. Results for the transverse case are more complicated. It can be observed from Figure 3.17 that as r increases (r varies from 0 to 0.33), the relative friction coefficient decreases gradually (shown as blue lines) while for r varying from 0.4 to 1, the relative friction coefficient increases (shown as black lines).

Longitudinal and transverse wavy cases

Once again, a function can be found to scale all results and form a single curve as shown in Figure 3.

Purely longitudinal wavy case

Study [START_REF] Venner | Numerical simulation of waviness in a circular ehl contact, under rolling/sliding[END_REF] showed that the transient purely longitudinal case remains a stationary problem. This phenomenon can also be seen in the current study. The influence of the waviness anisotropy on friction can be found in Table 3.2. It is reported that with r increasing from 1 to 32 the relative friction coefficient (µ r /µ s ) decreases, while the case r = ∞ predicts the highest relative friction coefficient value.

θ * 2 = M 0.4 L -1.2 [-1.1(λ y /a h ) 2 + 2.1(λ y /a h ) + 0.2](H c/Ai

Conclusion

In this chapter, the relative friction coefficient is predicted by means of numerical simulation for a large range of transient EHL operating condition and surface waviness anisotropy (r = λ x /λ y ). A relation is derived for the relative friction coefficient for isotropic and anisotropic harmonic surfaces. For the isotropic wavy case, all relative friction coefficient values can be unified onto a single curve applying a dimensionless parameter θ 2 . For anisotropic wavy cases, a roughness anisotropy function ff (r ) is used to scale all results together, when the same operating conditions are considered. Finally, results of all cases can be combined into a single equation:

µ r µ s = 1 + 0.56[ff (r )θ 2 ] -2 + 0.23[ff (r )θ 2 ] -4 (3.23)
where ff (r ) = (2.7r 2 + 0.1)/(r 3 + 1.7r 2 + 0.1), if r < 1 1.23 -0.23(1/r ) 1.4 , otherwise and θ 2 = M 0.33 L -1.1 (λ/a h ) 0.67 (H c/Ai ) with λ = min(λ x , λ y ).

However, from a lubrication point of view, the purely longitudinal wavy case (r = ∞) has a very different frictional behavior. A separate scaling parameter θ * 2 is obtained for the purely longitudinal wavy case.

Early work on conformal contacts showed that the onset of mixed lubrication regime was roughly around a "lambda ratio" equal to 3, in which the "lambda ratio" is defined as the ratio of oil film thickness to the combined surface roughness. The current work reveals that besides the classical "lambda ratio" parameter, the operating conditions as well as the surface topography play an important role on the onset of the mixed lubrication regime.

Chapter 4 Friction of complex rough surfaces

Introduction

The previous chapter has already predicted the friction increase for simple harmonic surface waviness, while the range of wavelengths is limited. In reality, the real surface roughness is more complicated. It contains very different wavelengths and amplitudes. On the other hand, the full numerical simulation is time-consuming for real rough surfaces. In this chapter a rapid prediction method based on the roughness power spectral density (PSD) is provided to predict the friction increase due to the roughness. Section 4.2 will describe this method in the first place. Section 4.3 tests this method for an artificial fractal surface roughness. Finally, section 4.4 predicts the transition of the full-film to the mixed lubrication regime by applying this method for real measured surface roughness.

Power spectral density friction method 4.2.1 PSD friction model

The power spectral density (PSD) is a mathematical tool that decomposes a rough surface into harmonic components of different frequencies [START_REF] Jacobs | Quantitative characterization of surface topography using spectral analysis[END_REF], which enables the pressure increase to be calculated analytically for each frequency component. Subsequently, the friction variations for the whole rough surface can be obtained. At last, the relative friction coefficient is obtained. The calculation process is as follows:

A rough surface topography r r x,y can be expressed in the frequency domain by means of the Fourier transform: r r q x ,q y = ( 4

N x N y )
x,y (r r x,y )e -i (q x x+q y y) (

where r r x,y is the discrete form of the surface roughness r r (x, y), q x and q y are the wavenumbers in the x and y direction respectively. In general, Equation 4.1 is computed by the fast Fourier transform (FFT) algorithm.

Combine Equation 4.1 and Equation 3.17, the deformed surface roughness r r d q x ,q y in the frequency domain is:

r r d q x ,q y = ( A d A i ) q x ,q y • |r r q x ,q y |. (4.2)
According to the relation between the pressure and the elastic deformation of the waviness given in Appendix D, the pressure increase in the frequency domain follows the expression below δp q x ,q y = πE 2k(r )λ (r r q x ,q yr r d q x ,q y ) (

where λ is defined as λ = min(2π/q x , 2π/q y ). With the inverse discrete Fourier transform, the pressure increase in the space domain is obtained:

δp x,y = ( 4
N x N y )
q x ,q y δp q x ,q y • e i (q x x+q y y) . (

According to the friction force Equation 3.14, the ratio of the shear stress τ r /τ s can be derived as:

τ r (x, y) τ s (x, y) = η r (x, y) η s (x, y) • h s (x, y) h r (x, y) = η r (x, y) η s (x, y) • h s (x, y) h s (x, y) -a d (x, y) (4.5)
It is easy to obtain the shear stress distribution τ s (x, y) for the smooth surface case, where the pressure distribution for the smooth surface case can be replaced by a semi-elliptical pressure distribution:

p s (x, y) = p h 1 -(x/a h ) 2 -(y/a h ) 2 , if x 2 + y 2 ≤ a 2 h 0, otherwise. (4.6)
Afterwards, the pressure distribution for roughness cases is computed by p s + δp. The shear stress distribution τ r (x, y) for a rough surface case is obtained as:

τ r (x, y) = η r (x, y) η s (x, y) • h s (x, y) h s (x, y) -a d (x, y) • τ s (x, y). ( 4.7)
Finally, the shear forces for both the smooth case and the rough case are computed by integrating the shear stress τ s (x, y) and τ r (x, y), respectively. The relative friction coefficient is then calculated according to Equation 3.15. A detailed description of the prediction process of the relative friction coefficient is shown in

Model validation

To validate the model described in subsection 4.2.1, the relative friction coefficient evaluated from a full numerical simulation is compared with that predicted by PSD under the same operating conditions. The same calculation domain, mesh size and slide-to-roll ratio as in subsection 3.3.2 are used here. An artificial fractal rough surface is chosen to validate the model mentioned in the previous subsection and the operating condition parameters are listed in Table 4.1. 4.2 presents the friction ratio as a function of the number of mesh points for the two methods. The results predicted by the two schemes are basically identical. The ratio of the friction coefficients predicted by the PSD method changes slightly (< 0.08%) with decreasing mesh size. However, in the full numerical simulation, a large mesh size leads to a relative large error. This is because some high frequency components of the rough surface can not be correctly represented on such a large mesh size. In this chapter, the precision of the numerical results simulated by a 513 × 513 points is considered acceptable.

The artificial surface roughness

In this section, an artificial fractal surface roughness is selected to test the rapid friction prediction method.

Surface roughness power spectrum

Many parameters like the root-mean square (R q ), standard deviation (σ), skewness (Sk), kurtosis (K ) et al. are employed to describe a measured surface roughness [START_REF] Bhushan | Modern tribology handbook, two volume set[END_REF]. The power spectral density (PSD) is perhaps the most used method for the surface description. Reference [START_REF] Stover | Optical scattering: measurement and analysis[END_REF] gives the calculation of the PSD for a surface roughness:

C 2D q x ,q y = (L x L y ) -1 |r r q x ,q y | 2 (4.8)
where L x and L y are the length of the measured surface roughness profiles in the x and y directions. A 2D surface roughness is shown in Figure 4.2 (a) and it can be represented as a 2D PSD C 2D in frequency space. When the surface is isotropic, its power spectral density is radially symmetric (reported in Many real rough surfaces show self-affine characteristic. A surface is self-similar when it has the same statistical properties as a magnified version of itself [START_REF] Barabási | Fractal concepts in surface growth[END_REF].The power spectral density of a self-affine surface often follows the following equation [START_REF] Jacobs | Quantitative characterization of surface topography using spectral analysis[END_REF][START_REF] Persson | On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion[END_REF] and shown in Figure 4.3:

C i so (q) =        C 0 , if q L ≤ q < q r C 0 q -2-2H , if q r ≤ q < q s 0, otherwise (4.9)
in which q = q 2 x + q 2 y , q L is the long-wavelength cut-off wave vector, q r is the long-wavelength roll-off wave vector and q s is the short-wavelength cut-off wave vector. H is the Hurst exponent, which is related to the fractal dimension D f of the surface [START_REF] Barabási | Fractal concepts in surface growth[END_REF] :

H = 3 -D f .
Once the power spectral density of a rough surface is obtained, other parameters like the RMS roughness h rms , RMS slope h rms and RMS curvature h rms can be derived [START_REF] Nayak | Random process model of rough surfaces[END_REF].

Friction increase prediction of a rough surface

The friction increase is determined by the non-linear viscosity variations from pressure variations caused by roughness deformation. In this subsection, these variables are plotted by employing the artificial surface roughness with q r = 50000 m -1 shown in Figure 4 Subsequently, the pressure increase can be computed by Equation 4.3 and Equation 4.4. Figure 4.7 gives the pressure increase distribution in frequency and space domains, respectively. Then the pressure distribution as well as film thickness distribution for the rough case can be obtained as p s + δp and h sr r d respectively. Shear stress for both cases are predicted by η(u 2u 1 )/h, which are shown in Figure 4.8. Those values are unrealistically high, and more complex rheological models are required for an absolute shear stress prediction, but here we are only interested in the relative shear. The shear stress for rough case is higher than that of smooth case, this is because pressure variations caused by roughness deformations make a big contribution to the viscosity. Finally, the friction force is calculated by Equation 3.14 and the relative friction coefficient is given by Equation 3.15. For this rough surface, friction forces for CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES the smooth case and the rough case are F s = 1.0606×10 6 N and F r = 1.3624×10 6 N, respectively. The relative friction coefficient is µ r /µ s = 1.2846.

Comparison between the EHL simulation and the PSD prediction

The chosen artificial surface topography is generated by means of fractals without a roll-off region. The friction prediction only occurs in the high pressure zone (X 2 + Y 2 ≤ 1.0), hence in the full numerical EHL simulation this generated surface roughness is located in the high pressure zone (shown in Figure 4.9 (b)). In order to make results keep periodical, the "patch" also needs to be periodical shown in Figure 4.9 (c). As shown in the previous subsection, one can find that the deformed amplitude is a fundamental intermediate variable. It is necessary to compare this parameter in the first place. The resulting deformed micro-geometry, of which the original surface topography is shown in Figure 4.9 (a) for a full numerical EHL simulation and a PSD prediction, are presented in Figure 4.10. In terms of numerical simulation results, the deformed micro-geometry r r d is obtained by h s -h r and removing data outside the high-pressure zone (X 2 +Y 2 ≤ 1). Once again, for this specific surface, the operating conditions are the ones given in Table 4.1 where M = 1000 and CHAPTER 4. FRICTION OF COMPLEX ROUGH SURFACES L = 10. Both the deformed surface topographies r r d are shown in the same high pressure region. The maximum Hertzian pressure reaches 1.66 GPa and the maximum surface roughness height deformed significantly from 1 × 10 -7 m to 3.5 × 10 -8 m. In addition, it is shown that the height distribution of the deformed surface roughness from the EHL simulation and the PSD prediction are very similar. Magnitudes of the central line r r d (x, 0) and p(x, 0) for both methods are also similar. However, the full numerical simulation shows a more "smooth" results in terms of the deformed surface roughness, this is because small wavelength components can not be represented correctly in full numerical simualtion. On the other hand, the Amplitude Reduction Theory is valid for small roughness amplitude values, this is the reason for some components shown in subfigure(d) whose maginitude is a little higher. In addition, the pressure distribution predicted by the PSD method is the combination of a semi-elliptical pressure distribution and the pressure deflection i.e. p s + δp x,y . Twenty artificial random rough surfaces are generated (shown in Figure 4.12) with the same input parameters i.e. the standard deviation σ = 5 × 10 -8 m, lengths of final topography L x = L y = 8.29×10 -4 m, roll-off wave number q r = 0 m -1 and Hurst exponent=0.8. These generated roughness are used to compare results from the full numerical simulation as well as the PSD prediction.

The averaged relative friction coefficients for these artificial random rough surfaces are given in Table 4.3, showing that the two different prediction methods give close results. The average deviation is around 8%. It seems that the averaged relative friction values simulated by the full numerical simulation are higher than those predicted by PSD. This is beacuse in PSD prediction pressure spike effects are not taken into consideration. Here we need to notice that when the number of mesh points is selected as 513 × 513, the calculation time of the full numerical simulation is almost 3 days, this is because more than 2000 time steps are needed. Meanwhile the calculation time of the PSD prediction is only 15 minutes.

Measured surface roughness

The surface roughness used in this section is the surface roughness measured from discs employed in friction experiments in Reference [START_REF] Bonaventure | Influence of random surface roughness on friction in elastohydrodynamic, mixed and boundary lubrication[END_REF].

Friction prediction under a specific operating condition

In this subsection, the measured surface roughness dART is collected within an area of 1 cm 2 and having a sampling interval of hx = h y = 3.653 µm [START_REF] Bonaventure | Influence of random surface roughness on friction in elastohydrodynamic, mixed and boundary lubrication[END_REF]. This surface roughness is measured from a finished disc made of AISI 52100 steel with a 60 mm diameter and an 8 mm thickness. The surface roughness is then corrected by removing large scales (through a high pass-filtering operation), this corrected surface roughness and the original measured surface roughness are shown in Figure 4.13. The root-mean-square of this surface roughness is R q = 0.038µm.

Since the area of high pressure happens in an area of Hertz contact zone, small square windows with the length of 2a h are used to extract effective prediction areas (shown in Figure 4.14). For a measured surface roughness, more effective prediction areas are obtained. For each extracted sub-surface, the friction prediction process described in subsection 4.2.1 is applied to predict its relative friction coefficient under a specific operating condition. Measured operating conditions and lubricant parameters are listed in Table 4 For this measured surface roughness dART under this specific operating condition, the relative friction coefficient is defined as the averaged value of that of sub-surfaces. From this figure, one can observe that all relative friction values are smaller than 2 except for a point whose value reaches 4, hence this high value point will be removed. Figure 4.17 shows the corrected relative friction coefficients. At this time, the mean relative friction coefficient is 1.09 and the corresponding standard deviation is 0.071.

Operating condition effects

The method to predict friction variations for measured rough surfaces has already described in the previous subsection. In this subsection, effects of operating conditions on friction will be investigated. Except for the rolling speed u r which is varying, the values of all the other parameters are fixed as shown in Table 4 4.24 show the deformed surface patches and their pressure variations as functions of the wave vector q for all the operating conditions in table 4.5. The values at q = 0 in these figures are removed, because there are no deformations of sub-surfaces when q = 0. According to observation, when the rolling speed u r increases, the value of the deformed subsurface increases while the value of the pressure increase decreases.

Friction curves for measured surface roughness

Friction variations for two measured surface roughness are plotted in this subsection. For the roughness dART (shown in Figure 4.13), the variations of the Moes load parameter M and material parameter L using geometry and lubricant parameter in Table 4 shown in these two figures, minimum friction points do not occur. Reference [START_REF] Schipper | Transitions in the lubrication of concentrated contacts[END_REF] defined the transition position as the intersection between two tangent lines to the friction curve (shown as blue crosses in Figure 4. 27 and 4.31). For conformal contact, the transition from the mixed lubrication to full-film lubrication regime occurs around Λ = 3.0. For non-conformal contact cases considered in this subsection, the onset of the "mixed regime" at hc/R q varying from 4 to

Conclusion

This chapter proposes a rapid analytical method using the power spectral density to predict friction increase which is due to non-linear viscosity variations from pressure variations caused by surface roughness deformation. This method is validated for an artificial surface roughness by comparing the results with those of a full numerical simulation. A good agreement is found between the full numerical simulation and the PSD prediction. Then the PSD prediction method is also employed to analyse friction variations for a measured surface roughness. Prediction results show that as the "Λ ratio" increases, the friction decreases. The minimum friction point describing the transition from the mixed lubrication to the full-film lubrication regime does not occur for high pressure cases. The same result has also been found in experimental work [START_REF] Bonaventure | Influence of random surface roughness on friction in elastohydrodynamic, mixed and boundary lubrication[END_REF][START_REF] Schipper | Transitions in the lubrication of concentrated contacts[END_REF].

Chapter 5

Conclusion and perspective

Conclusion

Proper lubrication and low friction play essential roles in energy conservation. The Stribeck curve is a good tool to describe the friction variation throughout the entire lubrication regime. However, most of the work to obtain this friction curve was related to low pressure condition (conformal contact), and the work was mainly experimental. The motivation for conducting this work is to determine the frictional behaviour under high pressure conditions (nonconformal contact) when the surface roughness is taken into consideration. The main conclusions are listed below:

(i) Current MultiGrid codes show good efficiency to solve lubricated contact problems, but are not sufficiently robust to treat the rough surface problem in a general way. This is because for very rough surfaces, large variations of the coefficient ρh 3 /η in Reynolds equation occur on a small scale. Alcouffe et al [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] proposed an efficient way to restore the performance by constructing the coarse grid operator and the intergrid transfers. An new transient EHL code is modified by implementing the coarse grid construction method into the existing EHL MultiGrid code. The new code shows better performance in solving the rough contact problem.

(ii) The Stribeck curve shows that the "Λ ratio" is a suitable parameter to plot friction variation. According to the Amplitude Reduction Theory [START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF][START_REF] Venner | An engineering tool for the quantitative prediction of general roughness deformation in ehl contacts based on harmonic waviness attenuation[END_REF], under very high pressure, the surface roughness will deform and this deformation depends on the operating conditions as well as on the surface roughness parameters. This means that the old "Λ ratio" is not a proper parameter. In this work, using harmonic surface waviness, an elaborate scaling parameter θ 2 ff(r ) including the old "Λ ratio", operating conditions and surface anisotropy parameters is found. Using this new parameters, all simulation results can be unified into a single curve and a curve-fitting equation is obtained.

(iii) A rapid analytical prediction method using the power spectral density is proposed in this work to predict friction increase for a complex surface roughness. This friction increase is due to non-linear viscosity variations from pressure variations caused by surface roughness deformation. This method is validated by the comparison between predictions and full numerical simulations when an artificial surface roughness is employed. Then this rapid prediction method is applied to analyse measured surface roughness.

CHAPTER 5. CONCLUSION AND PERSPECTIVE

Perspective

The work carried out in this thesis offers a useful tool to detect friction variations under high pressure conditions. However, what we presented in this work is related to Newtonian lubricant, sufficient lubrication and isothermal model. Future work can be suggested as follows:

• In current work, a small slide-to-roll ratio is used to prediction frictional force and a unique curve is obtained. Big slide-to-roll ratio value should be applied to detect friction variations in a further study.

• The EHL model used in this work is isothermal and Newtonian. Thermal effects and non-Newtonian lubricant can be considered in future study.

• The work studied in this thesis is numerical, more relevant experiments should be conducted to give comparable results.

Appendix A Construction of the coarse grid operator

The coarse grid operator is:

L H = J H h L h I h H (A.1)
In order to get the coarse grid operator, orthogonal basis vectors are needed. The first step is to compute I h H 〈e H iC , jC ,k 〉:

I h H 〈e H iC , jC ,k 〉 = C h i F, j F,k e h i F, j F,k + NW h i F, j F,k e h i F -1, j F +1,k + SW h i F, j F,k e h i F -1, j F -1,k + SE h i F, j F,k e h i F +1, j F -1,k + N E h i F, j F,k e h i F +1, j F +1,k + N h i F, j F,k e h i F, j F +1,k + S h i F, j F,k e h i F, j F -1,k + W h i F, j F,k e h i F -1, j F,k + E h i F, j F,k e h i F +1, j F,k (A.2)
Then, substituting Equation A.2 in L h 〈I h H e H iC , jC ,k 〉 reads:

L h 〈I h H e H iC , jC ,k 〉 = L h [C h i F, j F,k e h i F, j F,k + NW h i F, j F,k e h i F -1, j F +1,k + SW h i F, j F,k e h i F -1, j F -1,k + SE h i F, j F,k e h i F +1, j F -1,k + N E h i F, j F,k e h i F +1, j F +1,k + N h i F, j F,k e h i F, j F +1,k + S h i F, j F,k e h i F, j F -1,k + W h i F, j F,k e h i F -1, j F,k + E h i F, j F,k e h i F +1, j F,k ] (A.3)
in which the influence of L h need to act on every term, as illustrated in

L h 〈C h i F, j F,k e h i F, j F,k 〉 = C h i F, j F,k L h 〈e h i F, j F,k 〉 = C h i F, j F,k [L h c (i F, j F, k)e h i F, j F,k + L h nw (i F, j F, k)e h i F -1, j F +1,k + L h sw (i F, j F, k)e h i F -1, j F -1,k + L h se (i F, j F, k)e h i F +1, j F -1,k + L h ne (i F, j F, k)e h i F +1, j F +1,k + L h n (i F, j F, k)e h i F, j F +1,k + L h w (i F, j F, k)e h i F -1, j F,k + L h s (i F, j F, k)e h i F, j F -1,k + L h e (i F, j F, k)e h i F +1, j F,k ]
• The influence of L h on the east point (i + 1, j ) (shown as From the three equations above one observes that the operator L h influences fifteen points around point (i , j ) in total. By that analogy, with respect to nine terms in Equation A.3 twenty five points around point (i , j ) are affected by the operator L h . Coefficients of 25 points around (i , j ) is denoted as:

L h 〈E h i F, j F,k e h i F +1, j F,k 〉 = E h i F, j F,k L h 〈e h i F +1, j F,k 〉 = E h i F, j F,k [L h c (i F + 1, j F, k)e h i F +1, j F,k + L h nw (i F + 1, j F, k)e h i F, j F +1,k + L h sw (i F + 1, j F, k)e h i F, j F -1,k + L h se (i F + 1, j F, k)e h i F +2, j F -1,k + L h ne (i F + 1, j F, k)e h i F +2, j F +1,k + L h n (i F + 1, j F, k)e h i F +1, j F +1,k + L h w (i F + 1, j F, k)e h i F, j F,k + L h s (i F + 1, j F, k)e h i F +1, j F -1,k + L h e (i F + 1, j F, k)e h i F +2, j F,
L h 〈N E h i F, j F,k e h i F +1, j F +1,k 〉 = N E h i F, j F,k L h 〈e h i F +1, j F +1,k 〉 = N E h i F, j F,k [L h c (i F + 1, j F + 1, k)e h i F +1, j F +1,k + L h nw (i F + 1, j F + 1, k)e h i F, j F +2,k + L h sw (i F + 1, j F + 1, k)e h i F, j F,k + L h se (i F + 1, j F + 1, k)e h i F +2, j F,k + L h ne (i F + 1, j F + 1, k)e h i F +2, j F +2,k + L h n (i F + 1, j F + 1, k)e h i F +1, j F +2,k + L h w (i F + 1, j F + 1, k)e h i F, j F +1,k + L h s (i F + 1, j F + 1, k)e h i F +1, j F,k + L h e (i F + 1, j F + 1, k)e h i F +2, j F +1,k ]
               IC (i F -2, j F + 2, k) IC (i F -1, j F + 2, k) IC (i F, j F + 2, k) IC (i F + 1, j F + 2, k) IC (i F + 2, j F + 2, k) IC (i F -2, j F + 1, k) IC (i F -1, j F + 1, k) IC (i F, j F + 1, k) IC (i F + 1, j F + 1, k) IC (i F + 2, j F + 1, k) IC (i F -2, j F, k) IC (i F -1, j F, k) IC (i F, j F, k) IC (i F + 1, j F, k) IC (i F + 2, j F, k) IC (i F -2, j F -1, k) IC (i F -1, j F -1, k) IC (i F, j F -1, k) IC (i F + 1, j F -1, k) IC (i F + 2, j F -1, k) IC (i F -2, j F -2, k) IC (i F -1, j F -2, k) IC (i F, j F -2, k) IC (i F + 1, j F -2, k) IC (i F + 2, j F -2, k)                (A.4)
where: Hence, the coarse grid operator whose stencil is:

IC (i F -2, j F + 2, k) = NW h i F, j F,k L h nw (i F -1, j F + 1, k) APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR IC (i F -2, j F + 1, k) = NW h i F, j F,k L h w (i F -1, j F + 1, k) IC (i F -2, j F, k) = NW h i F, j F,k L h sw (i F -1, j F + 1, k) + SW h i F, j F,k L h nw (i F -1, j F -1, k) + W h i F, j F,k L h sw (i F -1, j F, k) IC (i F -2, j F -1, k) = SW h i F, j F,k L h w (i F -1, j F -1, k) + W h i F, j F,k L h sw (i F -1, j F, k) IC (i F -2, j F -2, k) = SW h i F, j F,k L h sw (i F -1, j F -1, k) IC (i F -1, j F + 2, k) = NW h i F, j F,k L h n (i F -1, j F + 1, k) + N h i F, j F,k L h nw (i F, j F + 1, k) IC (i F -1, j F + 1, k) = C h i F, j F,k L h nw (i F, j F, k) + NW h i F, j F,k L h c (i F -1, j F + 1, k) + N h i F, j F,k L h w (i F, j F + 1, k) + W h i F, j F,k L h n (i F -1, j F, k) IC (i F -1, j F, k) = C h i F, j F,k L h w (i F, j F, k) + NW h i F, j F,k L h s (i F -1, j F + 1, k) + SW h i F, j F,k L h n (i F -1, j F -1, k) + N h i F, j F,k L h sw (i F, j F + 1, k) + S h i F, j F,k L h nw (i F, j F -1, k) + W h i F, j F,k L h c (i F -1, j F, k) IC (i F -1, j F -1, k) = C h i F, j F,k L h sw (i F, j F, k) + SW h i F, j F,k L h c (i F -1, j F -1, k) + S h i F, j F,k L h w (i F, j F -1, k) + W h i F, j F,k L h s (i F -1, j F, k) IC (i F -1, j F -2, k) = SW h i F, j F,k L h s (i F -1, j F -1, k) + S h i F, j F,k L h sw (i F, j F -1, k) IC (i F, j F + 2, k) = NW h i F, j F,k L h ne (i F -1, j F + 1, k) + N E h i F, j F,k L h w (i F, j F -1, k) + N h i F, j F,k L h n (i F, j F + 1, k) IC (i F, j F + 1, k) = C h i F, j F,k L h n (i F, j F, k) + NW h i F, j F,k L h e (i F -1, j F + 1, k) + N E h i F, j F,k L h w (i F + 1, j F + 1, k) + N h i F, j F,k L h c (i F, j F + 1, k) + W h i F, j F,k L h ne (i F -1, j F, k) + E h i F, j F,k L h nw (i F + 1, j F, k) IC (i F, j F, k) = C h i F, j F,k L h c (i F, j F, k) + NW h i F, j F,k L h se (i F -1, j F + 1, k) + SW h i F, j F,k L h ne (i F -1, j F -1, k) + SE h i F, j F,k L h nw (i F + 1, j F -1, k) + N E h i F, j F,k L h sw (i F + 1, j F + 1, k) + N h i F, j F,k L h s (i F, j F + 1, k) + S h i F, j F,k L h n (i F, j F -1, k) + W h i F, j F,k L h e (i F -1, j F, k) + E h i F, j F,k L h w (i F + 1, j F, k) APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR The last step is to compute J H h 〈L h I h H e H iC ,
L H =        L H nw (iC , jC , k) L H n (iC , jC , k) L H ne (iC , jC , k) L H w (iC , jC , k) L H c (iC , jC , k) L H e (iC , jC , k) L H sw (iC , jC , k) L H s (iC , jC , k) L H se (i c, jC , k)        (A.5)
in which:

L H nw (iC , jC , k) = IC (i F -2, j F + 2, k)c h i F -2, j F +2,k + IC (i F -1, j F + 2, k)w h i F -1, j F +2,k + IC (i F -2, j F + 1, k)n h i F -2, j F +1,k + IC (i F -1, j F + 1, k)nw h i F -1, j F +1,k L H sw (iC , jC , k) = IC (i F -2, j F -2, k)c h i F -2, j F -2,k + IC (i F -2, j F -1, k)s h i F -2, j F -1,k + IC (i F -1, j F -2, k)w h i F -1, j F -2,k + IC (i F -1, j F -1, k)sw h i F -1, j F -1,k L H se (iC , jC , k) = IC (i F + 2, j F -2, k)c h i F +2, j F -2,k + IC (i F + 1, j F -2, k)e h i F +1, j F -2,k + IC (i F + 2, j F -1, k)s h i F +2, j F -1,k + IC (i F + 1, j F -1, k)se h i F +1, j F -1,k L H ne (iC , jC , k) = IC (i F + 2, j F + 2, k)c h i F +2, j F +2,k + IC (i F + 1, j F + 2, k)e h i F +1, j F +2,k + IC (i F + 2, j F + 1, k)n h i F +2, j F +1,k + IC (i F + 1, j F + 1, k)ne h i F +1, j F +1,k APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR L H n (iC , jC , k) = IC (i F, j F + 2, k)c h i F, j F +2,k + IC (i F -1, j F + 2, k)e h i F -1, j F +2,k + IC (i F + 1, j F + 2, k)w h i F +1, j F +2,k + IC (i F, j F + 1, k)n h i F, j F +1,k + IC (i F -1, j F + 1, k)ne h i F -1, j F +1,k + IC (i F + 1, j F + 1, k)nw h i F +1, j F +1,k L H s (iC , jC , k) = IC (i F, j F -2, k)c h i F, j F -2,k + IC (i F -1, j F -2, k)e h i F -1, j F -2,k + IC (i F + 1, j F -2, k)w h i F +1, j F -2,k + IC (i F, j F -1, k)s h i F, j F -1,k + IC (i F -1, j F -1, k)se h i F -1, j F -1,k + IC (i F + 1, j F -1, k)sw h i F +1, j F -1,k L H w (iC , jC , k) = IC (i F -2, j F, k)c h i F -2, j F,k + IC (i F -2, j F + 1, k)s h i F -2, j F +1,k + IC (i F -2, j F -1, k)n h i F -2, j F -1,k + IC (i F -1, j F, k)w h i F -1, j F,k + IC (i F -1, j F + 1, k)sw h i F -1, j F +1,k + IC (i F -1, j F -1, k)nw h i F -1, j F -1,k L H e (iC , jC , k) = IC (i F + 2, j F, k)c h i F +2, j F,k + IC (i F + 2, j F + 1, k)s h i F +2, j F +1,k + IC (i F + 2, j F -1, k)n h i F +2, j F -1,k + IC (i F + 1, j F, k)e h i F +1, j F,k + IC (i F + 1, j F + 1, k)se h i F +1, j F +1,k + IC (i F + 1, j F -1, k)ne h i F +1, j F -1,k
L H c (iC , jC , k) = IC (i F, j F, k)c h i F, j F,k + IC (i F, j F + 1, k)s h i F, j F +1,k + IC (i F, j F -1, k)n h i F, j F -1,k + IC (i F -1, j F, k)e h i F -1, j F,k + IC (i F + 1, j F, k)w h i F +1, j F,k

+ IC (i F -1, j F + 1, k)se h i F -1, j F +1,k + IC (i F -1, j F -1, k)ne h i F -1, j F -1,k + IC (i F + 1, j F -1, k)nw h i F +1, j F -1,k + IC (i F + 1, j F + 1, k)sw h i F +1, j F +1,k

Appendix B Derivation of matrix A j for line relaxation

According to the discrete Reynolds equation in Chapter 2, the dynamic residual r l i , j ,k is defined as: r l i , j ,k = p f l i , j ,k -(L l c P l i , j ,k + L l n P l i , j +1,k + L l w P l i -1, j ,k + L l s P l i , j -1,k + L l e P l i +1, j ,k + L l nw P l i -1, j +1,k + L l sw P l i -1, j -1,k + L l se P l i +1, j -1,k + L l ne P l i +1, j +1,k ) + h l (1.5 ρl i , j ,k H l i , j ,k -2.0 ρl i -1, j ,k H l i -1, j ,k + 0.5 ρl i -2, j ,k H l i -2, j ,k )

+ (h l ) 2 ht (1.5 ρi,j,k H i , j ,k ) (B.1)
and the left hand side equation is:

L l i , j ,k 〈P l 〉 =(L l c P l i , j ,k + L l n P l i , j +1,k + L l w P l i -1, j ,k + L l s P l i , j -1,k + L l e P l i +1, j ,k + L l nw P l i -1, j +1,k + L l sw P l i -1, j -1,k + L l se P l i +1, j -1,k + L l ne P l i +1, j +1,k ) h l (1.5 ρl i , j ,k H l i , j ,k -2.0 ρl i -1, j ,k H l i -1, j ,k + 0.5 ρl i -2, j ,k H l i -2, j ,k )

- (h l ) 2
ht (1.5 ρi,j,k H i , j ,k ).

(B.2)

Reference [START_REF] Venner | Multi-Level Methods in Lubrication[END_REF] recommend that the switch parameter ξ l limit = 0.3 is a good choice to have an efficient smooth performance.

B.0.1 Gauss-Seidel line relaxation

When the local coefficient ξ l satisfies the below condition, the Gauss-Seidel line relaxation is applied.

ξ l i ±1/2, j ,k (h l ) 2 > ξ l limit and A j i ,i -1 = L l wh l (1.5 ρl i , j ,k K l l 1,0 -2.0 ρl i -1, j ,k K l l 0,0 + 0.5 ρl i -2, j ,k K l l 1,0 ) -(h l ) 2 ht (1.5 ρi,j,k K l l 1,0 ) (B.7)

ξ l i , j ±1/2,k (
• if i < n x -1:

A j i ,i +1 = L l eh l (1.5 ρl i , j ,k K l l 1,0 -2.0 ρl i -1, j ,k K l l 2,0 + 0.5 ρl i -2, j ,k K l l 3,0 ) h l (1.5 ρl i , j ,k ∆K l l 0,0 -2.0 ρl i -1, j ,k ∆K l l 1,0 + 0.5 ρl i -2, j ,k ∆K l l 2,0 )

- (h l ) 2
- (h l ) 2
ht (1.5 ρi,j,k ∆K l l 0,0 ) (B.12)

• for i > 2:

A j i ,i -2 = - 1 4
L l wh l (1.5 ρl i , j ,k ∆K l l 2,0 -2.0 ρl i -1, j ,k ∆K l l 1,0 + 0.5 ρl i -2, j ,k ∆K l l 0,0 )

- (h l ) 2
ht (1.5 ρi,j,k ∆K l l 2,0 ) (B.13)

• for i > 1: where the pressure increase ∆P is the following function of the deformation [START_REF] Johnson | Contact Mechanics[END_REF], i.e.:

A j i ,i -1 = L l w - 1
∆P = π 2 Ai 2λ/a h (1 - Ad Ai ).
Using a first order approximation of the dimensionless pressure increase ∆P , ᾱ•∆P reduces to:

ᾱ • ∆P ≈ ᾱ π 2 Ai 2λ/a h ≈ [ L π ( 3M 2 ) 1/3 ] π 2 H c 2λ/a h ( H c Ai ) -1 (C.2)
where ᾱ is expressed as ᾱ = (L/π)(3M /2) 1/3 . Defining H D c the dimensionless film thickness film thickness value using the well-known Hamrock-Dowson Equation [START_REF] Hamrock | Isothermal elastohydrodynamic lubrication of point contacts: Part 1-theoretical formulation[END_REF], H D c = 1.69G 0.53 U 0.67 W -0.067

2

(1-0.61exp(-0.73k)) with k = 1.03 for circle contact (i.e. H D c = 1.2L 0.53 U 0.49 M -0.067 ). Now the dimensionless central film thickness H c can be rewritten as: where Ad /Ai ≈ 1 -0.15∇ 2 ≈ 1 -0.15(λ/a h )(M /L) 0.5 .

H c = R 2
Observing Equation C.4 and Equation C.5, the exponent for the parameter M , L and λ/a h are summarized in Table C.1. Hence the expression of the θ 2 parameter using M 0.33 , L -1.1 and (λ/a h ) 0.67 obtained by curve-fitting, has coefficients that fall in the range outlined above.

Appendix D

The relation between the elastic deformation and corresponding pressure for 2D wavy surfaces

Reference [START_REF] Johnson | Contact Mechanics[END_REF] gives the relation between the elastic deformation and pressure for 1D wavy surfaces is:

δp = πE ∆ λ (D.1)
and for 2D wavy surfaces, which consists of two separate components:

δp = πE ∆ x λ x + πE ∆ y λ y . (D.2)
For the surface waviness expressed as Eq.3.1, there is no explicit equation to describe this relation. Hence, it is necessary to give such relation. Assuming a pressure distribution: Ensuite ce solveur modifié est étendu aux cas transitoires pour résoudre le problème de contact avec rugosité. La courbe de frottement est généralement représentée en fonction du « λ ratio », le rapport entre l'épaisseur du film d'huile et la valeur moyenne quadratique de la rugosité de la surface. Cependant, ce paramètre est moins approprié pour tracer les variations de frottement dans des conditions de haute pression (régime élasto piézovisqueux). Dans la deuxième partie de ce travail, le coefficient de frottement est calculé à l'aide du code EHL modifié pour de nombreuses conditions de fonctionnement ainsi que pour les paramètres d'ondulation de surface. Les résultats de la simulation montrent qu'il n'y a pas de courbe de frottement unique lorsque l'ancien paramètre « λ ratio » est utilisé. En se basant sur la théorie de la réduction d'amplitude, un nouveau paramètre de dimensionnement qui dépend des conditions de fonctionnement et des paramètres d'ondulation est trouvé, ce qui peut donner une courbe de frottement unique pour les situations de haute pression.

1. 1 69 D. 1 77 viii

 169177 World primary energy consumption (red column: Non-OECD, blue column: OECD). (Source: IEEJ Outlook 2019 and Scenario) . . . . . . . . . . . . . . . . . . . . . . . . 1.2 The variation of mean global surface temperature relative to 1880-2017. (Source: NASA/GISS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Total consumption by End-Use Sector, from 2000 to 2017. (Source: Data from the U.S. Energy Information Administration) . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Pressure flow factors. (Source: Reference [21]) . . . . . . . . . . . . . . . . . . . . . 1.5 Relative amplitude as a function of ∇ 2 under pure rolling, where ∇ 2 is dimensionless wavelength parameter, A i and A d are amplitude of surface roughness and deformed surface roughness respectively. (Source: Reference [58] ) . . . . . . . . 1.6 Results obtained from measurements compared with theoretical attenuation curve defined by [56]. (Source: Reference [61] ) . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Friction coefficient versus speed for different loads. (Source: Ref. [83]) . . . . . . 2.1 Mesh point (x i , y j ) and it's related mesh region r i , j . . . . . . . . . . . . . . . . . . 2.2 Mesh point (i , j ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Interpolation process (green points: coarse grid points, black dots: fine grid points, blue dots: middle points on the fine grid, red point: central point on the fine grid). 2.4 Weighting factors for the interpolation (blue points: coarse grid points, black dots: fine grid points). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Flow chart of the hybrid relaxation process . . . . . . . . . . . . . . . . . . . . . . . 2.6 Implementation of the Multi-Grid method with a two level "V" cycle. . . . . . . . 2.7 The time-dependent "V" cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Shear stress-shear rate relationship for the EHL contact. . . . . . . . . . . . . . . . 3.2 Comparison of the relative deformed amplitude (Ad /Ai ) as a function of f (r )∇ 2 for the current model (blue squares) with those on Reference [53] (solid line) . . 3.3 Top view of the surface waviness with λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic surface waviness r = 1, (b) the longitudinal surface waviness r = 2, (c) the transverse surface waviness r = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0) (blue lines) of isotropic surface waviness (r = 1) for M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c during a time period. The central pressure line (red line) for the smooth case is plotted as a reference. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0) (blue lines) of longitudinal surface waviness (r = 2) for M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c during a time period. The central pressure line (red line) for the smooth case is plotted as a reference. . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0) (blue lines) of transverse surface waviness (r = 0.5) for M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c during a time period. The central pressure line (red line) for the smooth case is plotted as a reference. . . . . . . . . . . . . . . . . . . . . . . . . v LIST OF FIGURES 3.7 The dimensionless central film thickness H c rough as a function of the dimensionless time T for: M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic surface wavy case, (b) the longitudinal surface wavy case r = 2, (c) the transverse surface wavy case r = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 The relative friction coefficient µ r /µ s as a function of the dimensionless time T for: M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) isotropic surface wavy case, (b) longitudinal surface wavy case r = 2, (c) transverse surface wavy case r = 0.5. (Blue dotted line: average value of the relative friction coefficient.) . . . . . . . . . 3.9 Relative friction coefficient as a function of H c/Ai for a specific operating condition 3.10 Effect of the load parameter M on the relative friction coefficient for L = 10 and λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient as a function of M 0.33 • (H c/Ai ) . . . . . . . . . . . . . . . . . . . . 3.11 Effect of material parameter L on the relative friction coefficient for M = 2000 and λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient as a function of L -1.1 • (H c/Ai ) . . . . . . . . . . . . . . . . . . . 3.12 Effect of wavelength λ/a h on the relative friction coefficient for M = 1000 and L = 10: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient as a function of (λ/a h ) 0.67 • (H c/Ai ) . . . . . . . . . . . . . . . . . . . . . 3.13 Relative friction coefficient as a function of the classical parameter "lambda ratio" i.e. H c/Ai for a large range of operating conditions. . . . . . . . . . . . . . . . . . 3.14 Relative friction coefficient as a function of the new parameter θ 2 , simulation results: black circles; fitted curve: the black dashed line. . . . . . . . . . . . . . . . . 3.15 Relative friction coefficient (µ r /µ s ) as a function of H c/Ai for different r (1 ≤ r ≤ 32) values for: M = 1000, L = 10 and λ y /a h = 0.5 (left), zoom from 2.3 -2.7 (right). 3.16 Relative friction coefficient (µ r /µ s ) as a function of ff(r ) • (H c/Ai ) for different r (1 ≤ r ≤ 32) values for: M = 1000, L = 10 and λ y /a h = 0.5. . . . . . . . . . . . . . . . 3.17 Relative friction coefficient (µ r /µ s ) as a function of (H c/Ai ) for different r (0 ≤ r ≤ 1) values for: M = 1000, L = 10 and λ x /a h = 0.5(upper), zoom from 2.4 -2.6 (lower). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Relative friction coefficient (µ r /µ s ) as a function of ff(r ) × (H c/Ai ) for different r (0 ≤ r ≤ 1) values for: M = 1000, L = 10 and λ x /a h = 0.5. . . . . . . . . . . . . . . . 3.19 ff (r ) as a function of r . Numerical results: red squares. Fitted curve: solid lines. . 3.20 Comparison between the transient relative friction coefficient and that of the stationary case. Transient results: black line. Stationary results: magenta dashdotted line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Relative friction coefficients for the purely longitudinal wavy case: (a) relative friction coefficient as a function of original "lambda ratio" H c/Ai parameter, (b) relative friction coefficient as a function of the new parameter θ * 2 . . . . . . . . . . . 4.1 Flow chart for the relative friction coefficient prediction . . . . . . . . . . . . . . . 4.2 Surface roughness (a) and its 2D power spectral density (b) . . . . . . . . . . . . . 4.3 Power spectral density C i so of the self-affine surface(Figure 4.2(a)) with H = 0.8. . 4.4 The selected artificial surface roughness (a), amplitude distribution of this surface roughness (b) and its power spectral density (c). . . . . . . . . . . . . . . . . . 4.5 The ratio of the deformed amplitude and the initial amplitude fitted as Equation (3.17) (a) and the deformed surface roughness in frequency domain (b). . . . . . 4.6 Comparison between the initial surface roughness (a) and the deformed surface roughness (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Pressure increase distribution in frequency (a) and space (b) domains, respectively. vi LIST OF FIGURES 4.8 Shear stress distributions for the smooth case (a) and for the rough case (b). . . . 4.9 The generated surface roughness patch (a), the roughness patch in the high pressure zone (b) and the periodical roughness pattern for full the numerical simulation (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Top view of the deformed surface roughness for a specific time step for a full numerical simulation (a) and for a PSD prediction (b). Central line r r d (x, 0) of the deformed surface roughness for the full numerical simulation (c) and for the PSD prediction (d). Central line p(x, 0) of the pressure distribution for the full numerical simulation (e) and for the PSD prediction (f). . . . . . . . . . . . . . . . . . . . 4.11 The relative friction as a function of dimensionless time employing the surface roughness pattern in Figure 4.9 (c) for the full numerical simulation method. . . 4.12 Top view of the twenty generated artificial random rough surfaces from N • 1 to N • 20 with a same standard divation value σ = 0.05µm and a same set of operating conditions listed in Table 4.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Measured surface roughness dART: (a) corrected surface roughness and (b) raw surface roughness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Effective prediction areas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 An extracted surface patch of dART: (a) surface roughness height of this surface patch; (b) deformed surface patch; (c) pressure fluctuation of the surface patch. . 4.16 Extracted 529 surface patch (left) and their relative friction coefficient values (right). 4.17 Corrected relative friction coefficient values for 528 surface patches (left) and its histogram (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Relative friction coefficient as a function of total number of surface patches under the operating condition in Table 4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.19 A extracted surface patch (left) and its average initial amplitude as a function as q. 4.20 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Moes parameters M and L as a function of u r for the roughness dART. . . . . . . 4.26 The relative friction coefficient as a function of u r for the roughness dART. . . . . 4.27 The relative friction coefficient as a function of "Λ ratio" for the roughness dART. 4.28 Surface roughness dARTEb (left) and the top view of this roughness (right). . . . 4.29 Moes parameters M and L as a function of the rolling speed for the surface roughness dARTEb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.30 The relative friction coefficient as a function of u r for the surface roughness dAR-TEb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.31 The relative friction coefficient as a function of the "Λ ratio" for the surface roughness dARTEb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.32 The relative friction coefficient as a function of the rolling speed u r for the surface roughness dART and dARTEb, respectively. . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF FIGURES A.1 Influences of the coarse grid operator L h on central point, east point and northeast point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 A.2 Influences of the injection operator J H h on nine coincidental points (blue points). Pressure distribution and the corresponding elastic deformation. . . . . . . . . . 75 D.2 Amplitude of the elastic deformation A D d as a function of initial pressure amplitude A p i for the following cases: (a) isotropic, (b) purely transverse, (c) purely longitudinal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 D.3 Amplitude of the elastic deformation A D d as a function of wavelength λ (λ = min(λ x , λ y )) for the following cases: (a) isotropic, (b) purely transverse, (c) purely longitudinal. 76 D.4 Amplitude of the elastic deformation A D d as a function of anisotropy parameter r = λ x /λ y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .List of Tables

Figure 1 . 1 :Figure 1 . 2 :Figure 1 . 3 :

 111213 Figure 1.1: World primary energy consumption (red column: Non-OECD, blue column: OECD). (Source: IEEJ Outlook 2019 and Scenario)

Figure 1 . 4 :

 14 Figure 1.4: Pressure flow factors. (Source: Reference [21])

Figure 1 . 5 :Figure 1 . 6 :

 1516 Figure 1.5: Relative amplitude as a function of ∇ 2 under pure rolling, where ∇ 2 is dimensionless wavelength parameter, A i and A d are amplitude of surface roughness and deformed surface roughness respectively. (Source: Reference [58] )

Figure 1 . 7 :

 17 Figure 1.7: Friction coefficient versus speed for different loads. (Source: Ref. [83])

Chapter 5 : 10 2. 2 11 2. 3 13 2. 4 21 2. 7 22 2. 8

 5102113134217228 The main results of current work are summarized and recommendations of future work are also made. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Transient EHL model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2 Dimensionless equations and parameters . . . . . . . . . . . . . . . . . . The finite difference scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Transfer operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.2 Injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5 Coarse grid operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Implementation of the Multi-Grid method . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

21 )

 21 Compared to this equation, the Reynolds equation is of the same type with U = P , σ = 0, D = -( ρH 3 )/(η λ) and f = ∂( ρH )/∂X + ∂( ρH )/∂T . Rearranging the Reynolds equation (2.11) yields:

Figure 2 . 2 :

 22 Figure 2.2: Mesh point (i , j ).

Figure 2 . 3 :

 23 Figure 2.3: Interpolation process (green points: coarse grid points, black dots: fine grid points, blue dots: middle points on the fine grid, red point: central point on the fine grid).

Figure 2 . 4 :

 24 Figure 2.4: Weighting factors for the interpolation (blue points: coarse grid points, black dots: fine grid points).

Figure 2 . 5 :

 25 Figure 2.5: Flow chart of the hybrid relaxation process

Figure 2 . 7 :

 27 Figure 2.7: The time-dependent "V" cycles

Contents 3 . 1 25 3. 2 25 3. 3 28 3. 4 29 3. 5 34 3. 6

 31252253284295346 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Lubricant rheological models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.1 Relative friction coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.2 Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Time-dependent solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Effect of operating conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Effect of surface anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.6.1 Longitudinal and transverse wavy cases . . . . . . . . . . . . . . . . . . . . 38 3.6.2 Purely longitudinal wavy case . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

Figure 3 . 1 :

 31 Figure 3.1: Shear stress-shear rate relationship for the EHL contact.

Figure 3 . 2 :

 32 Figure 3.2: Comparison of the relative deformed amplitude (Ad /Ai ) as a function of f (r )∇ 2 for the current model (blue squares) with those on Reference [53] (solid line)

Figure 3 . 3 :

 33 Figure 3.3: Top view of the surface waviness with λ/a h = 0.5 and Ai = 0.5H c: (a) the isotropic surface waviness r = 1, (b) the longitudinal surface waviness r = 2, (c) the transverse surface waviness r = 0.5.

   show dimensionless pressure distribution P (X , 0) and dimensionless film thickness distribution H (X , 0) for three surface waviness cases, in which H c is the central film thickness for the smooth case. For the operating condition of M = 1000 and L = 10, the dimensionless central film thickness for the smooth case is: H c = 0.0243. In this work, we only focus on small surface waviness amplitude values i.e. Ai /H c ∈[START_REF] Holmberg | Global energy consumption due to friction in passenger cars[END_REF][START_REF] Venner | Multi-Level Methods in Lubrication[END_REF].

CHAPTER 3 .Figure 3 . 4 :Figure 3 . 5 :Figure 3 . 6 :Figure 3 . 7 :

 334353637 Figure 3.4: Central line pressure P (X , 0) (black lines) and central line film thickness H (X , 0) (blue lines) of isotropic surface waviness (r = 1) for M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c during a time period. The central pressure line (red line) for the smooth case is plotted as a reference.

Figure 3 . 8 :

 38 Figure 3.8: The relative friction coefficient µ r /µ s as a function of the dimensionless time T for: M = 1000, L = 10, λ/a h = 0.5 and Ai = 0.5H c: (a) isotropic surface wavy case, (b) longitudinal surface wavy case r = 2, (c) transverse surface wavy case r = 0.5. (Blue dotted line: average value of the relative friction coefficient.)

Figure 3 . 9 :

 39 Figure 3.9: Relative friction coefficient as a function of H c/Ai for a specific operating condition

Figure 3 . 10 :Figure 3 .

 3103 Figure 3.10: Effect of the load parameter M on the relative friction coefficient for L = 10 and λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient as a function of M 0.33 • (H c/Ai )

  Figure 3.10 displays the influence of the load parameter M on the relative friction coefficient for L = 10 and λ/a h = 0.5. From Figure 3.10 (a), one can observe that for each M value, there is a single curve. As M increases, the relative friction coefficient decreases, meanwhile each single curve has a small left-shift. It seems that there is a scaling factor can be used to combine all curves together.Figure 3.10 (b) shows that M 0.33 is an appropriate factor to scale all curves.

Figure 3 .

 3 Figure 3.11 shows the influence of the material parameter L on the relative friction coefficient for M = 1000 and λ/a h = 0.5. From Figure 3.11 (a), one can observe that the parameter L affects the relative friction coefficient. As L increases, the relative friction coefficient increases, each single curve has a small right-shift. Once again, there exists a scaling parameter L -1.1 which can be applied to shift all curves together (shown in Figure 3.11 (b)).

Figure 3 . 11 :Figure 3 . 12 :

 311312 Figure 3.11: Effect of material parameter L on the relative friction coefficient for M = 2000 and λ/a h = 0.5: (a) relative friction coefficient as a function of H c/Ai , (b) relative friction coefficient as a function of L -1.1 • (H c/Ai )

CHAPTER 3 .

 3 FRICTION INFLUENCE OF HARMONIC SURFACE WAVINESS 3.12 (a)). The parameter λ/a 0.67 h is a suitable factor to unify those curves (shown in Figure3.12 (b)).

Figure 3 . 13 :

 313 Figure 3.13: Relative friction coefficient as a function of the classical parameter "lambda ratio" i.e. H c/Ai for a large range of operating conditions.

  -3.12, it is possible to combine all results obtained for different values of λ/a h , M , L as well as H c/Ai into a single curve using a dimensionless parameter θ 2 . Figure(3.14) shows the relative friction coefficient as a function of the new parameter θ 2 for 500 ≤ M ≤ 2000, 5 ≤ L ≤ 15, 0.25 ≤ λ/a h ≤ 1.0 and 1.0 ≤ H c/Ai ≤ 10. After curve-fitting, the single curve can be described by the following equation:

CHAPTER 3 .Figure 3 . 14 :

 3314 Figure 3.14: Relative friction coefficient as a function of the new parameter θ 2 , simulation results: black circles; fitted curve: the black dashed line.

Figure 3 .

 3 Figure 3.15 shows the relative friction coefficient (µ r /µ s ) as a function of H c/Ai for different r (1 ≤ r ≤ 32) values. It can be observed that the results of the longitudinal case show the same trend as those of the isotropic wavy case (r = 1). That is, a decreasing trend of the relative friction coefficient when H c/Ai increases. For each r value a single curve exists but it shifts to the left compared to the curve with r = 1. For the case considered, the relative friction coefficient monotonically decreases as r increases (shown in the right small figure in Figure 3.15).

Figure 3 . 15 :

 315 Figure 3.15: Relative friction coefficient (µ r /µ s ) as a function of H c/Ai for different r (1 ≤ r ≤ 32) values for: M = 1000, L = 10 and λ y /a h = 0.5 (left), zoom from 2.3 -2.7 (right).

Figure 3 . 16 :

 316 Figure 3.16: Relative friction coefficient (µ r /µ s ) as a function of ff(r ) • (H c/Ai ) for different r (1 ≤ r ≤ 32) values for: M = 1000, L = 10 and λ y /a h = 0.5.

Figure 3 . 17 :

 317 Figure 3.17: Relative friction coefficient (µ r /µ s ) as a function of (H c/Ai ) for different r (0 ≤ r ≤ 1) values for: M = 1000, L = 10 and λ x /a h = 0.5(upper), zoom from 2.4 -2.6 (lower).

Figure 3 . 18 : 40 CHAPTER 3 .

 318403 Figure 3.18: Relative friction coefficient (µ r /µ s ) as a function of ff(r ) × (H c/Ai ) for different r (0 ≤ r ≤ 1) values for: M = 1000, L = 10 and λ x /a h = 0.5.

Figure 3 .

 3 [START_REF] Christensen | Some aspects of the functional influence of surface roughness in lubrication[END_REF] presents the comparison between the transient results and the stationary results for M = 1000, L = 10, H c/Ai = 2, λ y /a h = 1.0 and λ x = ∞. In order to compare these results, the stationary result is also plotted as a function of the dimensionless time T in Figure3.20. Figure3.20 shows that the value of the relative friction coefficient is the same as the value of the stationary case, i.e. µ r /µ s = 1.673. This is due to the waviness term (Equation 3.1) reduces to RR(X , Y ) = Ai × cos(2π Y λ y /a h ), which is a time-independent term.

Figure 3 . 20 :

 320 Figure 3.20: Comparison between the transient relative friction coefficient and that of the stationary case. Transient results: black line. Stationary results: magenta dash-dotted line.

Figure 3 . 21 :

 321 Figure 3.21: Relative friction coefficients for the purely longitudinal wavy case: (a) relative friction coefficient as a function of original "lambda ratio" H c/Ai parameter, (b) relative friction coefficient as a function of the new parameter θ *2

Figure 3 .

 3 Figure 3.21 gives the results of the relative friction coefficient for the purely longitudinal wavy case for a relative large range of operating conditions: M ∈ [500, 2000], L ∈ [5, 15] and λ y /a h ∈ [0.25, 0.8]. As usual, in Figure 3.21 (a), a very smooth curve as a function of H c/Ai for each operating condition. While, the old scaling parameter θ 2 can not scale all cases together. For the purely longitudinal cases considered, there is indeed a new parameter θ * 2 expressed as Equation (3.22), and shown in Figure 3.21 (b) to unify all cases onto a single curve. The operating conditions are: 500 ≤ M ≤ 2000, 5 ≤ L ≤ 15 and 0.25 ≤ λ/a h ≤ 0.8.

Contents 4 . 1 44 4. 2 46 4. 3 59 4. 5

 41442463595 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Power spectral density friction method . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.1 PSD friction model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.2 Model validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The artificial surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.3.1 Surface roughness power spectrum . . . . . . . . . . . . . . . . . . . . . . 47 4.3.2 Friction increase prediction of a rough surface . . . . . . . . . . . . . . . . 48 4.3.3 Comparison between the EHL simulation and the PSD prediction . . . . 50 4.4 Measured surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.4.1 Friction prediction under a specific operating condition . . . . . . . . . . 54 4.4.2 Operating condition effects . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.4.3 Friction curves for measured surface roughness . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figure 4 . 1 .Figure 4 . 1 :

 4141 Figure 4.1: Flow chart for the relative friction coefficient prediction

  Figure 4.2 (b)).

Figure 4 . 2 :

 42 Figure 4.2: Surface roughness (a) and its 2D power spectral density (b)

Figure 4 . 3 :

 43 Figure 4.3: Power spectral density C i so of the self-affine surface(Figure 4.2(a)) with H = 0.8.

  .4.

Figure 4 . 4 :

 44 Figure 4.4: The selected artificial surface roughness (a), amplitude distribution of this surface roughness (b) and its power spectral density (c).

CHAPTER 4 .

 4 FRICTION OF COMPLEX ROUGH SURFACESThe first step is to compute the deformed surface roughness based on the Amplitude Reduction Theory[START_REF] Venner | Amplitude reduction of non-isotropic harmonic patterns in circular ehl contacts, under pure rolling[END_REF] with the operating condition in Table4.1.

  Figure 4.5 (a) shows the ratio of the deformed amplitude and the initial amplitude described by Equation 3.17. Combining the initial amplitude distribution in Figure 4.4, the deformed surface roughness distribution is depicted in Figure 4.5 (b).

Figure 4 . 5 :

 45 Figure 4.5: The ratio of the deformed amplitude and the initial amplitude fitted as Equation (3.17) (a) and the deformed surface roughness in frequency domain (b).

Figure 4 .

 4 Figure 4.6 compares the initial surface roughness with the deformed surface roughness. It can be observed that the initial surface roughness is deformed.

Figure 4 . 6 :

 46 Figure 4.6: Comparison between the initial surface roughness (a) and the deformed surface roughness (b).

Figure 4 . 7 :

 47 Figure 4.7: Pressure increase distribution in frequency (a) and space (b) domains, respectively.

Figure 4 . 8 :

 48 Figure 4.8: Shear stress distributions for the smooth case (a) and for the rough case (b).

Figure 4 . 9 :

 49 Figure 4.9: The generated surface roughness patch (a), the roughness patch in the high pressure zone (b) and the periodical roughness pattern for full the numerical simulation (c).

CHAPTER 4 .Figure 4 . 10 :

 4410 Figure 4.10: Top view of the deformed surface roughness for a specific time step for a full numerical simulation (a) and for a PSD prediction (b). Central line r r d (x, 0) of the deformed surface roughness for the full numerical simulation (c) and for the PSD prediction (d). Central line p(x, 0) of the pressure distribution for the full numerical simulation (e) and for the PSD prediction (f).

Figure 4 . 11 :

 411 Figure 4.11: The relative friction as a function of dimensionless time employing the surface roughness pattern in Figure 4.9 (c) for the full numerical simulation method.

Figure 4 .

 4 Figure 4.11 shows the periodical relative friction coefficient variation when the surface roughness pattern in Figure 4.9 (c) is used. As mentioned before, the relative friction coefficient value should be the averaged value of the relative friction coefficient in one time period. The PSD method only gives the value of the relative friction coefficient for one time step, hence one should predict this value for each time step within one time period.

CHAPTER 4 .

 4 Figure 4.12: Top view of the twenty generated artificial random rough surfaces from N • 1 to N • 20 with a same standard divation value σ = 0.05µm and a same set of operating conditions listed in Table 4.1.

Figure 4 . 13 :

 413 Figure 4.13: Measured surface roughness dART: (a) corrected surface roughness and (b) raw surface roughness.

CHAPTER 4 .

 4 Figure 4.14: Effective prediction areas.

Figure 4 . 15 :

 415 Figure 4.15: An extracted surface patch of dART: (a) surface roughness height of this surface patch; (b) deformed surface patch; (c) pressure fluctuation of the surface patch.

. 4 .

 4 Under this specific operating condition, Moes parameters M = 9.09 and L = 16.23, Figure4.15 shows a deformed surface geometry and pressure increase distribution of a surface patch, meanwhile its relative friction coefficient is 1.03.

Figure 4 . 16 :

 416 Figure 4.16: Extracted 529 surface patch (left) and their relative friction coefficient values (right).

Figure 4 .

 4 Figure 4.[START_REF] Lubrecht | Multigrid, an alternative method of solution for two-dimensional elastohydrodynamically lubricated point contact calculations[END_REF] shows the extracted 529 surface patches (left) and their values of the relative friction coefficient (right). For this measured surface roughness dART under this specific operating condition, the relative friction coefficient is defined as the averaged value of that of sub-surfaces. From this figure, one can observe that all relative friction values are smaller than 2 except for a point whose value reaches 4, hence this high value point will be removed. Figure4.17 shows the corrected relative friction coefficients. At this time, the mean relative friction coefficient is 1.09 and the corresponding standard deviation is 0.071.

Figure 4 . 17 :

 417 Figure 4.17: Corrected relative friction coefficient values for 528 surface patches (left) and its histogram (right).

Figure 4 .

 4 Figure 4.18 shows the mean relative friction coefficient as a function of total number of extracted surface patches with operating condition in Table4.4. It can be observed that the variation of mean relative friction for different the number of surface patch are small. For the sake of saving computation time, the number of surface patch for the rest study is selected as

Figure 4 . 18 :

 418 Figure 4.18: Relative friction coefficient as a function of total number of surface patches under the operating condition in Table 4.4.

  .4. Equation C.1 shows that the pressure fluctuation δp makes a significant contribution to the relative friction coefficient. Equation D.1 shows that the amplitude deformation r r d of surface roughness is also an important variable. Therefore, it is necessary to investigate δp and r r d . A surface patch shown extracted from the surface roughness dART in figure 4.19 is extracted from the rough surface dART to study the variations of δp and r r d when different operating conditions are employed (shown in table 4.5).

Figure 4 . 57 CHAPTER 4 .Figure 4 . 20 :

 4574420 Figure 4.19: A extracted surface patch (left) and its average initial amplitude as a function as q.

Figure 4 . 21 :

 421 Figure 4.21: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 2.

Figure 4 . 22 :

 422 Figure 4.22: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 3.

Figure 4 . 23 :

 423 Figure 4.23: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 4.

CHAPTER 4 .Figure 4 . 24 :

 4424 Figure 4.24: The deformed surface patch shown in Figure 4.19 (left) and its pressure variations for case 5.

Figures 4. 20

 20  show the deformed surface patches and their pressure variations as functions of the wave vector q for all the operating conditions in table 4.5. The values at q = 0 in these figures are removed, because there are no deformations of sub-surfaces when q = 0. According to observation, when the rolling speed u r increases, the value of the deformed subsurface increases while the value of the pressure increase decreases.

   show the deformed surface patches and their pressure variations as functions of the wave vector q for all the operating conditions in table 4.5. The values at q = 0 in these figures are removed, because there are no deformations of sub-surfaces when q = 0. According to observation, when the rolling speed u r increases, the value of the deformed subsurface increases while the value of the pressure increase decreases.

  .4 are shown in Figure 4.25.

Figure 4 . 25 :

 425 Figure 4.25: Moes parameters M and L as a function of u r for the roughness dART.

Figure 4 .

 4 Figure 4.26 shows variations of the relative friction coefficient employing the operating conditions shown in Figure 4.25. It can be observed from this figure that as rolling speed u r increases, the relative friction coefficient decreases. This decrease trend can also be found in Reference [140] Figure 6.1 (b). According to Hamrock-Dowson central film thickness equation [120], h c is a function of rolling speed u r , hence Figure 4.26 can be re-plotted as a function of the "Λ ratio" (shown in Figure 4.27).

Figure 4 . 26 :

 426 Figure 4.26: The relative friction coefficient as a function of u r for the roughness dART.

Figure 4 . 27 :

 427 Figure 4.27: The relative friction coefficient as a function of "Λ ratio" for the roughness dART.

Figure 4 . 28 :

 428 Figure 4.28: Surface roughness dARTEb (left) and the top view of this roughness (right).

Figure 4 . 29 :

 429 Figure 4.29: Moes parameters M and L as a function of the rolling speed for the surface roughness dARTEb.

Figure 4 . 30 :

 430 Figure 4.30: The relative friction coefficient as a function of u r for the surface roughness dAR-TEb.

5 .

 5 At the same time, two different surface roughness show different friction behaviour(shown in Figure 4.32).

Figure 4 . 31 :

 431 Figure 4.31: The relative friction coefficient as a function of the "Λ ratio" for the surface roughness dARTEb.

Figure 4 . 32 :

 432 Figure 4.32: The relative friction coefficient as a function of the rolling speed u r for the surface roughness dART and dARTEb, respectively.

Contents 5 . 1 64 5. 2

 51642 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

  Figure A.1 for 3 terms: • The influence of L h on the central point (i , j ) (shown as Figure A.1 (a)):

  Figure A.1 (b)):

  k ] APPENDIX A. CONSTRUCTION OF THE COARSE GRID OPERATOR • The influence of L h on the north-east point (i + 1, j + 1) (shown as Figure A.1 (c)):

Figure A. 1 :

 1 Figure A.1: Influences of the coarse grid operator L h on central point, east point and north-east point.

  jC ,k 〉. The way to extend the restriction coefficients on 25 points in Equation A.4 is the same as that of computing L h 〈I h H e H iC , jC ,k 〉. There is no need to do a restriction for 25 points. As was shown in Figure A.2, besides coincidental points (blue points) themselves, 3 points contribute to each corner point (shown as Figure A.2 (a)), 5 points influence each middle point (shown as Figure A.2 (b)) and 8 points impact on central point (shown as Figure A.2 (c)).

Figure A. 2 :

 2 Figure A.2: Influences of the injection operator J H h on nine coincidental points (blue points).

- (h l ) 2 ht ( 1 .

 21 h l ) 2 > ξ l limit .(B.3) Then the matrix A j i ,m is given by:A j i ,m = ( ∂L l i , j ,k 〈P l 〉 ∂P l m, j ,k ) P h = P h (B.4)for 0 < m < n x and 0 < i < n x .The matrixA j i ,m have different expressions for different conditions: • if |i -m| > 1: A j i ,m =h l (1.5 ρl i , j ,k K l l |i -m|,0 -2.0 ρl i -1, j ,k K l l |i -m-1|,0 + 0.5 ρl i -2, j ,k K l l |i -m-2|,0 ) 5 ρi,j,k K l l |i -m|,0 ) (B.5) • if i = m: A j i ,i = L l ch l (1.5 ρl i , j ,k K l l 0,0 -2.0 ρl i -1, j ,k K l l 1,0 + 0.5 ρl i -2, j ,k K l l 2,0 ) -(h l ) 2 ht (1.5 ρi,j,k K l l 0,0 ) (B.6) APPENDIX B. DERIVATION OF MATRIX A J FOR LINE RELAXATION • if i > 1:

2 4 (

 24 Jacobi distributive line relaxationWhen the local coefficient ξ l can not satisfy the condition Equation B.3, the Jacobi distributive line relaxation is used. Subsequently, the matrix A j i ,m is:A j i ,m = [ ∂L l i , j ,k 〈P l 〉 ∂P l m, j ,k -1 ∂L l i , j ,k 〈P l 〉 ∂P l m+1, j ,k + ∂L l i , j ,k 〈P l 〉 ∂P l m-1, j ,k + ∂L l i , j ,k 〈P l 〉 ∂P l m, j +1,k + ∂L l i , j ,k 〈P l 〉 ∂P l m, j -1,k )] P h = P h (B.9)for 0 < m < n x and 0 < i < n x . For convenience, introducing a parameter ∆K m, j to simply the n + K l l m+1,n + K l l m,n+1 + K l l m,n-1 ). (B.10)For different i values, the matrix of A j i ,m are represented as: • for |i -m| > 2: A j i ,m =h l (1.5 ρl i , j ,k ∆K l l |i -m|,0 -2.0 ρl i -1, j ,k ∆K l l |i -m-1|,0 + 0.5 ρl i -2, j ,k ∆K l l |i -m-2|,0 ) -(h l ) 2 ht (1.5 ρi,j,k ∆K l l |i -m|,0 )

-Derivation of the scaling parameter θ 2

 2 h l (1.5 ρl i , j ,k ∆K l l 1,0 -2.0 ρl i -1, j ,k ∆K l l 0,0 + 0.5 ρl i -2, j ,k ∆K l l 1,0 ) -(h l ) 2 ht (1.5 ρi,j,k ∆K l l 1,0 ) (B.14)Appendix CAccording to the Barus[START_REF] Barus | Isothermals, isopiestics and isometrics relative to viscosity[END_REF] viscosity-pressure equation, the shear stress ratio can be approximated as:τ r τ s ≈ e ᾱ•∆P = 1 + ᾱ • ∆P + ( ᾱ • ∆P )

R 2 x /a 2 h 2

 222 is expressed as R 2 x /a 2 h = (3/2) -2/3 M -2/3 U -1/2 . Substituting Equation C.3 into Equation C.2 gives: ᾱ • ∆P ≈ 1.6467[L -1.03 M -0.4 (λ/a h ) 1 (H c/Ai )] -1 . (C.4)Applying a second order approximation of ∆P , ᾱ∆P yields:ᾱ • ∆P ≈ π Ai 2λ/a h (1 -Ad Ai ) ≈ 0.24[L -1.03 M -0.1 (λ/a h ) 0 (H c/Ai )] -1 (C.5)

  Figure D.1: Pressure distribution and the corresponding elastic deformation.

Figure D. 1

 1 Figure D.1 shows that the elastic deformation due to the pressure shown as Figure D.1 (a) has the same harmonic waviness. Thus the elastic deformation Equation D.4 can be re-written as: D(x, y) = A D d cos( 2πx λ x )cos( 2πy λ y ) (D.5) in which A D d is the amplitude of the elastic deformation.

Figure D. 2 :

 2 Figure D.2: Amplitude of the elastic deformation A D d as a function of initial pressure amplitude A p i for the following cases: (a) isotropic, (b) purely transverse, (c) purely longitudinal.

Figure D. 2 Figure D. 3 :

 23 Figure D.2 shows the relation between the amplitude of the elastic deformation A D d and the initial pressure amplitude A p i for three r values. From this figure, it is can be found that the deformation amplitude A D d is proportional to the initial pressure amplitude A p i :

Figure D. 3

 3 Figure D.3 shows the relation between the amplitude of the elastic deformation A D d and the wavelength λ for three r values. From this figure, it is can be found that the deformation amplitude A D d is proportional to the the wavelength λ: A D d πE ∝ λ. (D.7) Figure D.4 shows the relation between the amplitude of the elastic deformation A D d and the anisotropy parameter r . From this figure, it is can be found that the deformation amplitude A D d

  1.1 U.S. CO 2 emissions from end-use sectors, 2008-2017. (Source: U.S. Energy Information Administration, August 2018 Monthly Energy Review) . . . . . . . . . . . . 3.1 Relative friction coefficient versus the number of mesh points for: M =1000, L=10, λ x /a h = 0.5, r=0.4 and H c/Ai = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Relative friction coefficient versus different surface anisotropy parameters for: M = 1000, L = 10, λ y /a h = 0.5 and H c/Ai = 2. . . . . . . . . . . . . . . . . . . . . . 4.1 Operating condition parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Relative friction coefficients as a function of the mesh points for two prediction schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 The relative friction coefficient obtained by EHL simulation and PSD prediction for 20 artificial random rough isotropic surfaces. . . . . . . . . . . . . . . . . . . . 4.4 Measured operating condition and lubricant parameters. . . . . . . . . . . . . . . 4.5 Operating conditions of selected cases. . . . . . . . . . . . . . . . . . . . . . . . . . C.1 Range of the exponent for each parameter. . . . . . . . . . . . . . . . . . . . . . . .

			LIST OF TABLES
	Nomenclature p pressure	Pa
	p s	pressure for smooth cases	Pa
	p h	Hertzian pressure	Pa
	A d δp	dimensionless deformed amplitude in the center of the contact pressure fluctuations	Pa
	A i ∆P	dimensionless initial amplitude dimensionless pressure fluctuations ∆P = δp/p h
	a h q x , q y wavenumbers in x and y directions respectively the radius of the contact area a h = 3 3wR x /(2E )	m 1/m
	C 2D q r	2d power spectral density roll-off wavenumber	m 4 1/m
	C i so R q	2d power spectral density of an isotropic surface root mean square of the surface roughness	m 4 m
	E 1 , E 2 elastic moduli of the two contact bodies R x reduced radius of curvature in x : 1/R x = 1/R 1x + 1/R 2x	Pa m
	E R	reduced modulus of elasticity 2/E = (1 -v 2 1 )/E 1 + (1 -v 2 2 )/E 2	Pa
	f	the friction force induced by the shearing of the lubricant	N
	F	the dimensionless friction force
	G	dimensionless material parameter G = αE
	G e	elastic shear modulus	Pa
	G ∞	the limiting elastic shear modulus	Pa
	h	film thickness	m
	H	dimensionless film thickness H = hR x /a 2 h
	h 0	mutual approach	m
	H 0	dimensionless mutual approach
	h c	central film thickness	m
	I h H	interpolation operator
	I H h	restriction operator
	L x , L y lengths of final topography	m
	L l	coarse grid operator on the l th level
		ix x

H c dimensionless central film thickness for a smooth case H c = h c R x /a 2 h H c rough dimensionless central film thickness for a rough case hx, h y dimensionless mesh sizes on the fine grid H x, H y dimensionless mesh sizes on the coarse grid L dimensionless material parameter (Moes) L = G(2U ) 0.25 M 2d dimensionless load parameter (Moes) M = W 2 (2U ) -0.75 y reduced radius of curvature in y : R y = R x m r wavelength ratio used to describe the surface anisotropy r = λ x /λ y r r undeformed surface roughness m r r d deformed surface roughness m RR dimensionless surface roughness RR = r r • R x /a 2 h SRR slide to roll ratio t time s T dimensionless time T = t ū/a h ∆T dimensionless time step u 1 , u 2 velocities of lower surface and upper surface respectively m/s u r mean velocity of contact surfaces u r

  Methods to solve the rough contact problem . . . . . . . . . . . . . . . . . 3 1.2.2 Friction in rough EHL contact problem . . . . . . . . . . . . . . . . . . . . 6

				LIST OF TABLES LIST OF TABLES
	z i so Chapter 1 pressure viscosity index isotropic
	α l Introduction pressure viscosity index the l th level		1/Pa
	ᾱ Subscripts dimensionless viscosity index ᾱ = αp h
	τ a, b	shear stress induced by the shearing of the lubricant inlet, outlet	Pa
	τ 0 i , d	the Eyring stress initial, deformed		Pa
	τ L r, s	the limiting shear stress rough, smooth		Pa
	∇ 2 st	dimensionless wavelength parameter ∇ 2 = (λ/a h ) M /L start
	λ x, y	dimensionless speed parameter space domain
	λ t q x , q y frequency domain the time constant for the fluid
	λ x , λ y wavelength in x, y direction		m
	η	viscosity		Pa • s
	η 0	the atmospheric viscosity		Pa • s
	η	dimensionless viscosity η = η/η 0
	ρ	density of the lubricant		Kg • m -3
	ρ 0	atmospheric density of the lubricant	Kg • m -3
	ρ	dimensionless density ρ = ρ/ρ 0
	θ 2	dimensionless new lambda ratio parameter
	θ * 2	dimensionless new lambda ratio parameter for purely longitudinal rough surfaces
	γ	strain rate		1/s
	µ	friction coefficient
	σ	standard deviation of surface roughness	m
	Abbreviations
	ART	Amplitude Reduction Theory
	EHL Elastohydrodynamic Lubrication
	FFT	Fast Fourier Transform
	PSD Power Spectral Density
	Superscripts
	d	deformed	a h
	y h	coordinate perpendicular to x the fine grid		m
	Y H	dimensionless coordinate Y = y/a h the coarse grid
			xi xii xiii

Contents 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 1.3 Research aims and Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.1 Research aims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.2 Outlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Table 1 .

 1 1: U.S. CO 2 emissions from end-use sectors, 2008-2017. (Source: U.S.

	Energy Informa-

Table 3 . 1 :

 31 Relative friction coefficient versus the number of mesh points for: M =1000, L=10, λ x /a h = 0.5, r=0.4 and H c/Ai = 2.

	Mesh points Relative friction coefficient
	129 × 129	1.489
	257 × 257	1.457
	513 × 513	1.462
	1025 × 1025	1.465

Table 3 .

 3 2: Relative friction coefficient versus different surface anisotropy parameters for: M = 1000, L = 10, λ y /a h = 0.5 and H c/Ai = 2.

	r	1	2	4	8	16	32	∞
	µ r /µ s 1.789 1.566 1.508 1.491 1.489 1.487 2.041

Table 4 .

 4 1: Operating condition parameters.

	Parameter	Value	Units
	w	600	N
	u r	0.84	m/s
	R x	0.018	m
	E	2.26 × 10 11 Pa
	α	2.2 × 10 -8	Pa -1
	η 0	4 × 10 -2	Pa • s
	h c	0.233	µm
	σ	5 × 10 -8	m
	L x = L y	8.29 × 10 -4 m
	q r	0	m -1
	Hurst exponent 0.8

Table 4 .

 4 2: Relative friction coefficients as a function of the mesh points for two prediction schemes.Mesh points N x × N y hx = h y µ r /µ s (PSD) µ r /µ s (EHL)

	257 × 257	1/64	1.534	1.487
	513 × 513	1/128	1.534	1.518
	1025 × 1025	1/256	1.534	1.510

Table 4 . 3 :

 43 The relative friction coefficient obtained by EHL simulation and PSD prediction for 20 artificial random rough isotropic surfaces.

	Surface number µ r /µ s (EHL) µ r /µ s (PSD) Deviation(%)
	1	1.40	1.30	7
	2	1.45	1.36	6
	3	1.12	0.96	14
	4	1.15	1.03	11
	5	1.42	1.34	6
	6	1.23	1.11	9
	7	1.37	1.29	6
	8	1.14	1.00	11
	(to be continue)
		53

Table 4 . 4 :

 44 Measured operating condition and lubricant parameters.

	Parameter Value	Units Parameter Value	Units
	w	10	N	η 0	0.7590 Pa • s
	u r	0.29	m/s	T r oom	22	• C
	R x	9.525	mm	SRR	25	%
	E	210	GPa	a h	87.95	µm
	α	2.0 × 10 -8 Pa -1	p h	617.23 MPa

Table 4 .

 4 5: Operating conditions of selected cases.

	Cases u r (m/s) M	L	µ r /µ s
	1	0.010	114.79 6.97	5.21
	2	0.046	36.40	10.22 2.02
	3	0.100	20.35	12.41 1.54
	4	0.154	14.73	13.82 1.38
	5	0.195	12.38	16.65 1.31

Table C .

 C 1: Range of the exponent for each parameter.

	Parameter Range
	M	-0.10 ∼ 0.40
	L	-1.03 ∼ -1.53
	λ/a h	0 ∼ 1
	H c/Ai	1
		74

  Le nouveau paramètre d'échelle et l'augmentation du frottement prédite par la méthode PSD montrent une bonne précision technique pour une utilisation pratique. MOTS-CLÉS : Lubrification élastohydrodynamique, Simulation numérique, Régime élasto piézo-visqueux, Théorie de réduction d'amplitude, Variation de frottement Laboratoire (s) de recherche : Laboratoire de Mécanique des Contacts et des Structures UMR CNRS 5259 -INSA de Lyon

	Bâtiment Sophie Germain
	Avenue Jean Capelle
	69621 Villeurbanne Cedex FRANCE
	Directeur de thèse: LUBRECHT Antonius. A
	Président de jury : CAYER-BARRIOZ Juliette
	Composition du jury : EVANS Pwt	KŘUPKA Ivan
	VENNER Cornelis. H	CAYER-BARRIOZ Juliette
	BIBOULET Nans	LUBRECHT Antonius. A

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2019LYSEI063/these.pdf © [Y. Zhang], [2019], INSA Lyon, tous droits réservés

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2019LYSEI063/these.pdf © [Y. Zhang], [2019], INSA Lyon, tous droits réservés Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2019LYSEI063/these.pdf © [Y. Zhang], [2019], INSA Lyon, tous droits réservés

Acknowledgments

To all those people who have supported me in carrying out my PhD study.

Firstly, I would like to express my sincere gratitude to my supervisor Professor Antonius Lubrecht for his guidance, encouragement and support of this research work. I would also like to thank my supervisor for opportunities I have been given to participate in scientific conferences as well as two months academic exchanges. Meanwhile, I am very grateful to my supervisor for caring my life in Lyon. During these three years, my supervisor has taught me more than this specific subject, he has taught me how to think more logically and how to solve problems more efficiently. Without his persistent help, this work would not have materialized.

I would also like to extend my thanks to Professor Cornelis Venner and his group for their warmly welcome and help in the MultiGrid method studying and the first manuscript writing during the period I have visited the University of Twente. I also want to thank my friends: Mr. Haichao Liu, Mr. Balan Ramani and Mrs. Xiaozhen He for their help when I was in the University of Twente. Special thanks to Mr. Binbin Zhang for taking care of me in Enschede.

I would like to show my appreciations to my colleagues: Dr. Romain Bugnicourt, Mr. Martin Denni, Ms. Simona Dahdah and Ms. Marjolaine Gonon-Caux. With them, we built a pleasant and friendly working environment. I would like to offer my special thanks to Dr. Marie-Pierre Noutary for teaching me the coarse grid construction method. I also thankful to Dr. Nans Biboulet for his time and advice of my work. I would like to thank the members of the secretary group of LaMCoS for their help over the last three years. Thanks to my following dear friends with whom I spent a happy time: Ms. Jing Lu, Ms. Pei Niu, Mr. Minhu Lv, Ms. Jingya Yuan, Mr. Minyang Lou et al. I would especially like to thank Dr. Hanfeng Gu and Dr. Teng Zhang, for their selfless help when I first came to Lyon. This research work is financially supported by the China Scholarship Council (CSC).

Finally, my deepest appreciation goes to my parents: Mr. Ming Zhang and Mrs. Shulan Zeng. Meanwhile, I also thank my boyfriend for his support and waiting.

Thank you for all of you again.

Appendix B Derivation of matrix

• for i < n x -1:

• for i < n x -2: