On Waviness and Two-Sided Surface Features in Thermal Elastohydrodynamically Lubricated Line Contacts

Machine components are designed to endure increasingly severe operating conditions due to the strive for improved energy efficiency of mechanical systems. Consequently, lubricated non-conformal contacts must rely on thin lubricant films where the influence of surface topography on the lubricating conditions becomes significant. Due to the complexity of the multiphysical problem, approximate assumptions are often employed to facilitate numerical studies of elastohydrodynamically lubricated (EHL) contacts. In this work, the rough, time dependent, thermal EHL problem is solved with focus on two main analyses. The first analysis focuses on the influence of sinusoidal roughness and the difference between a thermal non-Newtonian approach and an isothermal Newtonian approach. The second analysis is focused on the lubricating mechanisms taking place when two-sided surface features overtake within the thermal EHL contact. The results indicate that the film thickness in the outlet of the contact may be significantly overestimated by an isothermal Newtonian approach and that differences in the high-pressure region may also occur due to viscosity variations in the inlet of the contact. Moreover, for the studied two-sided surface features, it became evident that not only the surface feature combination but also the overtaking position influence the film thickness and pressure variations significantly.


Introduction
The current strive to improve energy efficiency of mechanical systems has led to machine components being designed to endure increasingly severe operating conditions. This generally implies an increased power density and reduced lubricant viscosity, which further results in reduced lubricant film thicknesses separating the contacting surfaces in highly loaded non-conformal contacts, typically found in for example, bearings, gears and cam-roller follower mechanisms. Thereby, a deeper understanding of how surface topography influences the lubricating performance becomes essential for improving the design of machine components and to achieve desired durability. Consequently, the general understanding of how surface roughness influences the lubricating performance of elastohydrodynamically lubricated (EHL) contacts has been significantly improved over the past years, see for example, Morales-Espejel [1].
Many important contributions to the understanding of non-smooth EHL contacts include those of a numerical nature, which were enabled with reasonable convergence rates and sufficient accuracy after the multi-grid technique was applied to the EHL problem by Lubrecht et al. [2,3] and Lubrecht [4]. This technique was further extended to the transient EHL problem and was used to study the influence ratio may again reduce due to a decay of the complementary function. This further implies that for sufficiently small wavelengths, the full solution may be described by the particular integral alone.
Apart from the aforementioned approach used to solve the EHL problem [6], a full-system finite element (FE) approach was more recently presented by Habchi et al. [43,44] and Habchi [45] that allowed the EHL problem to be solved in a fully coupled manner. This approach was extended to solve the thermal EHL (TEHL) problem by for example, Habchi et al. [46], Habchi [47,48], Lohner et al. [49] and Ziegltrum et al. [50]. While later also utilised and further extended to solve the transient, non-smooth, TEHL problem by Hultqvist et al. [51].
It is further known that under rolling-sliding conditions, significant shear heating takes place within the EHL contact due to shear of the relatively solid lubricant film, see for example, the works by Dowson and Whitaker [52], Kim and Sadeghi [53] and Hultqvist et al. [51]. Furthermore, it has been seen that due to the pressure sensitivity of viscosity, the entrainment flow still dominates within the high pressure region in EHL contacts. This has been used as the main justification for the often employed isothermal approaches to non-smooth rolling-sliding EHL contacts while studying film thickness variations. However, if the viscosity reduces sufficiently much for example, due to a sufficiently reduced viscosity as a result of shear thinning and temperature, especially in the inlet of the contact, the amplitude reduction theory may not prove accurate due to local viscosity variations.
In this work, thermal and non-Newtonian effects for both sinusoidal roughness and two-sided surface features are analysed, and the influence investigated. This is a problem that due to its complexity only recently has been studied in detail and still only to a small extent, see for example, Kaneta et al. [54]. Here, a full system FE approach [45] for the transient non-smooth problem [51] is further extended to enable a comparison between the thermal non-Newtonian solution and the isothermal Newtonian solution for sinusoidal roughness. Thereafter, the model is used to study the four possible scenarios that may occur in a rolling-sliding infinite line contact related to two-sided surface features as highlighted by Almqvist [55], that is, the cases where an asperity overtakes an asperity, an asperity overtakes a dent, a dent overtakes an asperity and a dent overtakes a dent.
The novelty of this work relates to the study of surface roughness and two-sided surface feature behaviour in transient EHL contacts with consideration to thermal and non-Newtonian effects. The purpose of this study is two-fold, first to establish a better understanding of the qualitative behaviour of surface roughness wavelength on the film thickness, but also of the thermal non-Newtonian behaviour in comparison to an isothermal Newtonian approach on a phenomenological level. Second, to understand the behaviour of two-sided surface features in EHL contacts, particularly in the context of the particular integral and the complementary function, and their influence.

Methodology
The modelling approach used in this study is based on the work done by Habchi [45,56], Lohner et al. [49] and Hultqvist et al. [51]. The theory, numerical approach and solution procedure for the transient TEHL problem including roughness features utilised in this work has previously been presented in detail by the authors, see Hultqvist et al. [51]. However, for the sake of completeness, the mathematical model that has been used in this work to solve the thermal non-Newtonian EHL problem is summarised in this section.
The problem at hand is assumed to be described as a typical lubricated non-conformal contact in between two infinitely long cylinders with different surface velocities, u 1 and u 2 , that are separated by a lubricant film. The difference in surface velocity between the surfaces defines the slide-to-roll ratio according to: The problem may be reduced to one cylinder by assuming an equivalent curvature in contact with a flat plane, as illustrated in Figure 1, according to: For detailed explanations of the aforementioned and upcoming parameters, the reader is referred to the nomenclature. Below, the theories describing the elastohydrodynamic model, the thermal model as well as the lubricant characterisation models are summarised in dimensional form. Moreover, the numerical approach is briefly summarised and important changes made to the aforementioned model highlighted.

Elastohydrodynamic Model
A generalised Reynolds equation that enables the study of a varying lubricant behaviour across the film thickness, related to non-Newtonian behaviour and thermal effects, was presented by Dowson [57] and improved by Yang and Wen [58]. This equation governs the lubricant pressure distribution within the EHL contact and reads: with the parameters describing the effective density and viscosity, based on the variations throughout the lubricant film, being written as follows: Moreover, the lubricant velocity throughout the film in the entrainment direction is written as: To solve the generalised Reynolds equation, knowledge about the viscosity and density behaviour over the film thickness becomes essential, which are affected not only by the pressure and the shear stress in the lubricant, but also by the heat distribution. Thereby, the energy dissipation and heat distribution in the lubricant film have to be solved properly, which is detailed in Section 2.2. Additionally, the lubricant film thickness in-between the contacting surfaces, h(x, t), has to be calculated and can be described according to: where the deformation in z-direction, w, may be found by utilising linear elasticity in combination with a plain strain approach due to the assumed infinite line contact. Moreover, a quasi-static response of the solids may be employed due to the difference in characteristic times between the EHL problem and elastic problem [59] and, in order to reduce the problem to only one deforming body, the equivalent elasticity approach may be utilised [45]. This further results in the stress equilibrium within the solid bodies to be formulated according to: where the stress tensor based on the strains in the material reads: and the parameters c 1 , c 2 and c 3 defines the equivalent material properties. For more details, the reader is referred to [45,51]. The surface features were in this work described in a manner that results in smooth profiles with no sharp corners [6,51]. The surface feature being located on the surface, i, was described according to: with the position of the surface feature changing in time according to: In the case of studying sinusoidal roughness, the roughness was only applied to one surface and hence described according to: again with the roughness moving through space in time according to Equation (10). The load balance between the applied load and the lubricant pressure was assumed to behave quasi-statically due to the stiffness of the lubricant film. This implies relatively small variations in film thickness also for relatively large rigid body movements [60]. These variations were assumed to be small enough to not affect the inlet of the contact significantly and thereby the film thickness negligibly. The load balance was therefore written according to: Finally, the cavitation model proposed by Wu [61] and employed by Habchi [45] was utilised in this work. This approach was employed due to its relatively simple implementation and fulfilment of the Reynolds boundary conditions. It is also assumed sufficiently accurate since lubricant film reformation does not need to be considered in this work.

Thermal Model
The heat generation and dissipation within the lubricant film and the solids were computed using the the energy equation. Assumptions reducing the equations based on the thin lubricant film as well as the infinitely long contact were employed [45]. In the lubricating film, convection and conduction in combination with heat generation from compression and shear governs the heat distribution, which are written according to: with the heat source from compression and shear following: The heat dissipation within the solids can be described by the convection and conduction alone since no heat is generated within the solids themselves, resulting in the energy equation being written, for each solid, i, according to: For more details about boundary conditions the reader is referred to, for example, References [48,49,51].

Lubricant Properties
The lubricant pressure, temperature and shear response was described in the same manner as explained in Reference [51], where the Tait equation of state (EoS) was utilised to describe the compressibility-pressure-temperature behaviour. A Vogel-like approach was employed for the viscosity-pressure-temperature dependency and a shifted Carreau model for the shear thinning behaviour of the lubricant. These equations with the corresponding input data are presented below and has reportedly been used to describe the pressure-shear-temperature behaviour of Squalane relatively successfully [62,63].
The Tait EoS describes the compressibility-pressure behaviour of the fluid according to: The compressibility-temperature variations are described by the ratio between the volume at ambient pressure and the volume at reference temperature, which is assumed to follow a linear relationship according to: By using K 0 = 11.74, a v = 8.36 · 10 −4 K −1 , K 00 = 8.658 GPa and β K = 6.332 · 10 −3 K −1 , the density variations of Squalane [63] may be calculated according to: The viscosity-pressure-temperature relationship describing the Newtonian viscosity behaviour, assumed to follow a Vogel-like relationship [62], reads: where the scaling parameter, ϕ, is dependent on the temperature and compressibility variations according to: Using Equations (19) and (20), the Newtonian viscosity-pressure-temperature response of Squalane may be found by using g = 3.921, ϕ ∞ = 0.1743, B F = 24.50 and µ ∞ = 0.9506 · 10 −4 Pa s.
The non-Newtonian shear thinning behaviour of the lubricant that was assumed to follow a shifted Carreau model may be written according to [62]: where the shear rate is calculated according to: and the use of µ r = 15.6 mPas, λ r = 2.26 · 10 −9 s, and n = 0.463 describes the shear thinning response of Squalane [62,64]. Note, however, that the shear stress was not limited by the function Λp in this work. Lastly, the thermal conductivity and volumetric heat capacity variations due to pressure and temperature were described according to [62]: and respectively. By using q = 2, the variations experienced by Squalane may be found [62].

Numerical Approach
The finite element (FE) numerical approach, including boundary conditions, weak formulations and implementation employed in this work, is thoroughly described in Hultqvist et al. [51] and therefore not repeated here. However, the main features are briefly described and modifications made in this work are highlighted below.
The transient TEHL problem was solved in a weakly coupled manner between a fully coupled EHL FE model and a thermal FE model using commercial software [65,66]. The two models were set up based on the work by Habchi [45] and the transient approach having features based on the work by Lohner et al. [49]. However, a major difference in solving the transient problem is employed where the full time dependent solution was solved at every iteration in the solution procedure, allowing for utilisation of numerical solvers in the commercial FE software [65], see Hultqvist et al. [51].
Moreover, the thermal model was in this work stabilised by the streamline upwind Petrov-Galerkin (SUPG) approach, see for example, Habchi [56]. The mesh was allowed to vary in element size along the entrainment direction in both the EHL model and thermal model, allowing for a similar accuracy to be achieved for a reduced number of degrees of freedom, leading to a reduced solution time compared to Reference [51]. Expressed in the dimensionless parameter X = x/a, the element size outside the high pressure region was set to 0.01, that is, for −4.5 ≤ X ≤ −1.5 and 1.5 ≤ X ≤ 2.5, while being set to 0.005 for the area comprising the high pressure region, that is, −1.5 ≤ X ≤ 1.5. In the thermal model, the film thickness height was discretised using an element size of 0.05, expressed in the dimensionless parameter Z = z/h. The transient solver used a backward differential formulation with an order 2 for the EHL model and an order varied between 1-5, adjusted to achieve the specified tolerances, for the thermal model.

Results and Discussion
The input data related to material properties of the contacting solids, contact parameters and operating conditions were chosen to represent a typical non-conformal contact found in for example, a rolling element bearing or mating gear teeth. However, relatively high sliding conditions were assumed in order to enlarge the effects related to sliding in the contact. For rolling conditions, shear thinning and thermal effects may be neglected [51], meaning that due to the small effect of compression on heat generation in EHL contacts [67], an isothermal Newtonian approach is a reasonable assumption and therefore such conditions are not studied in this work. The data used for this study are presented in Table 1 and were used in all simulations for straight forward comparison between the results. Table 1. Input data related to material properties and operating conditions used in the simulations.

Parameter Value
Unit Note that the material data in Table 1 describes a typical steel and may therefore not accurately describe the thermal properties of a hardened steel, which is typically used in the aforementioned applications. However, the influence of hardening on for example, the thermal conductivity was out of scope for this study and assumed to not influence the qualitative behaviour of the studied phenomena significantly.
The results are presented in two parts, first studying how sinusoidal roughness behaves within a thermal EHL contact, where focus lies on the differences between a thermal non-Newtonian approach and an isothermal Newtonian approach for various roughness wavelengths. Thereafter, the thermal EHL model is used to study different surface features that overtakes within the contact as well as different overtaking positions for one specific surface feature combination. For comparison, the smooth contact solutions for both the thermal non-Newtonian approach and the isothermal Newtonian approach are shown in Figure S1. Moreover, the results and input parameters are in this section expressed in dimensionless form using normalisation based on the following reference parameters and Hertz theory: where a and p h are the Hertzian dry contact parameters [68].

Sinusoidal Roughness
A comparison between the developed thermal non-Newtonian approach and an isothermal Newtonian approach was conducted for three different types of roughness. The roughness was sinusoidal according to Equation (11) and located on the slower surface, that is moving with speed u 1 = 0.75 m s −1 . For all sinusoidal waves, a fixed dimensionless amplitude ofĀ R,1 = 0.015 was used, however, with the dimensionless wavelength being varied between Λ R,1 = 1/2, Λ R,1 = 1/4 and Λ R,1 = 1/6, shown in Figure 2. This further implies that approximately 4, 8 and 12 waves were located inside the contact at one moment, respectively. A real surface is, however, made up of a complete spectrum of wavelengths that could possibly include all the aforementioned wavelengths and others, at once. However, in order to further understand the influence of wavelength on the lubricating mechanisms and film formation, this study was limited to the behaviour of the aforementioned separate wavelengths at a fixed amplitude. Initially, the case with the longest wavelength was studied, that is, Λ R,1 = 1/2. It was seen that the long wavelength roughness was generating a relatively smooth pressure profile with pressure ripples being directly connected to the peaks and valleys of the sinusoidal roughness, see Figure 3. The amplitude of the roughness within the contact was seen to deform significantly, as shown by the steady state solution used as initial conditions to the time dependent problem, that is, at Θ = 0. The steady state solution shows the film thickness and pressure associated with the deformed roughness within the contact due to the exclusion of time dependency and thereby also any inlet generated effects. This is commonly referred to as the particular integral [21] and dominates the pressure in the contact, however, with a small influence on the film thickness due to the significant roughness deformation within the contact.
Considering the transient solution, the roughness can be seen to be deformed noticeably less in the inlet of the contact and does, due to the viscosity magnitude, generate inlet perturbations that are formed somewhere in between −1 ≤ X ≤ −0.5 in what Morales-Espejel [23] refers to as the Couhier region based on the work by Couhier [69]. The amplitude of these inlet perturbations is noticeably smaller than the initial roughness amplitude and travels through the contact at the mean entrainment speed, u e . Moreover, these perturbations are wider than the wavelength of the roughness profile, which is explained by the negative slide-to-roll ratio in this case. Further meaning that the rough surface moves slower than the mean speed of the lubricant and thereby generates inlet perturbations with a wavelength of Λ = Λ R u e /u 1 due to the flow into the contact being limited by the roughness [21,38]. These inlet perturbations are commonly referred to as the complementary function and are, as opposed to the particular integral, related to small pressure variations but dominates the film thickness fluctuations. The temperature is seen to vary in direct relation to the roughness profile and is explained by the viscosity variations of the sheared lubricant due to the pressure variations, which thereby influences the shear heating [51]. A very good agreement related to both pressure and film thickness between the isothermal Newtonian approach and the thermal non-Newtonian approach can, for the case Λ R,1 = 1/2, be observed. However, with a significantly reduced film thickness during exiting of deformed roughness peaks for the thermal non-Newtonian approach, seen at for example, Θ = 2.25 in Figure 3. This is related to the reduced viscosity of the lubricant due to inlet shear thinning and shear heating, which leads to a further reduced viscosity and hence an increased pressure driven flow in the contact outlet (see Figure S2).
A smaller wavelength of Λ R = 1/4 shows significantly increased pressure ripple magnitudes compared to the previously studied wavelength. Here a reduced, but still almost complete, roughness deformation occurs within the contact as again seen by the steady state solution, see Figure 4. The amplitude of the complementary function is, however, slightly larger due to a reduced deformation of the roughness in the inlet. Furthermore, an increased amount of heat is generated due to the increased number of pressure peaks in combination with their increased magnitude. As can be seen in the temperature plots, the thermal effects also affect the film formation within the Couhier region and thereby the formation of the complementary function, which is especially visible for Θ = 2.25 in Figure 4.
A discrepancy between the thermal non-Newtonian solution and the isothermal Newtonian solution is starting to appear, which is not only visible in the outlet of the contact but also slightly within the high pressure region. This is explained by a reduction of the complementary function due to a reduced viscosity in the inlet region. A decay of the complementary function throughout the contact has previously been highlighted by, for example, Hooke [30] and Hooke and Morales [32,33]. This decay does, however, only seem to take place in the Couhier region for the case of Λ R,1 = 1/4, where the complementary function reduces as it passes through low viscosity regions connected to valleys on the surface roughness in this case, seen by the amplitude difference between the waves in-between approximately −0.8 ≤ X ≤ −0.4 at Θ = 2.5 in Figure 4. This is related to a reduced viscosity from shear thinning and shear heating in the case of the thermal non-Newtonian approach, which further leads to a difference between the two solutions. However, when the complementary function has passed further into the high pressure region, it moves through the rest of the contact relatively unaffected.
Moreover, it is interesting to note that the complementary function slightly increases in amplitude when passing through localised high pressure regions related to the surface roughness. This is explained by the density increase that further leads to a locally reduced film thickness as discussed by Venner and Bos [70]. The aforementioned effects become increasingly visible with a reduced wavelength of the waviness.
With a further reduced wavelength of the sinusoidal roughness using Λ R = 1/6, the pressure ripple magnitude as well as the number of pressure spikes within the contact increases further. This implies that the contact inlet becomes increasingly affected by fluctuations and thereby the formation of the complementary function is influenced to an increased extent, as can be seen in Figure 5. Due to the increased number of waves, the complementary function has time to pass through more roughness valleys with low viscosity already in the Couhier region, see Figure 6, which is seen to reduce the amplitude of the complementary function. The initially formed complementary function hence reduces in amplitude each time it passes through a low viscosity region and does basically vanish in the first half of the contact.
This occurs even though the lubricant flow is dominated by Couette flow throughout the contact in each of the studied cases (see Figure S3) and becomes increasingly pronounced for the thermal non-Newtonian approach because of the reduced viscosity connected to the roughness peaks in the inlet region following shear thinning and shear heating (see Figure S2). Moreover, In the second half of the contact, the film thickness reductions are seen to be directly connected to the temperature and viscosity variations, and thereby also the pressure and roughness itself, as can be seen in-between approximately 0.5 ≤ X ≤ 0.8 in Figure 5 at Θ = 2.25 and Θ = 2.5. This indicates that the complementary function effect has been fully reduced. To further study the differences between the thermal non-Newtonian solution and the isothermal Newtonian solution, the central and minimum film thickness values over time were extracted for each case analysed, see Figure 7. The relative difference between the two approaches is also visualised and was calculated according to: It can be seen that in the case of the complementary function defining the film formation throughout the contact, for example, in the case of Λ R,1 = 1/2, the difference between the thermal non-Newtonian and isothermal Newtonian approaches is negligible in the central region of the contact. However, in the outlet the solutions differ due to a reduced viscosity from the inlet shear thinning and shear heating. This does not seem to affect the central region due to the solidified behaviour of the entrapped lubricant, but does lead to an increased pressure driven flow when the pressure again drops in the outlet. With a reduced wavelength, however, the reduced amplitude of the complementary function leads to a discrepancy between the two approaches also in the central region due to viscosity differences in the Couhier region that influence the formation of the complementary function, which becomes increasingly pronounced for shorter wavelengths.
For the shortest wavelength studied in this work, the relative difference between central film thickness fluctuates between approximately 20 % and −14 % as a consequence of both the decaying complementary function and the compressibility-pressure-temperature dependency of the lubricant. Moreover, the minimum film thickness is overestimated by approximately 20 % during exiting of the deformed roughness as a result of reduced lubricant viscosity due to an increased temperature and the non-Newtonian lubricant behaviour.
This highlights the importance of applying a thermal non-Newtonian approach to accurately predict the film thickness behaviour within the high pressure region for short wavelength roughness, whereas for longer wavelength roughness the difference in the high pressure region is not as significant. Moreover, it is interesting to note that the largest difference in minimum film thickness is approximately the same for all studied wavelengths, and only slightly larger compared to the smooth contact case (see Figure S1). This implies that the minimum film thickness is not significantly affected by the local variations within the contact, but instead mainly dependent on the variations within the Couhier region, which further leads to a reduced outlet film thickness due to an increased outlet flow for the thermal non-Newtonian solution.
For the sake of completeness, the maximum pressure within the contact was compared between the three cases, see Figure 8. Due to the relatively small pressure fluctuations related to the complementary function [12,21], a less significant difference is seen compared to the film thickness variations. Nevertheless, as an effect of the reduced complementary function in the Couhier region, the pressure ripples reduce in magnitude within the contact, especially while considering thermal non-Newtonian effects, which further implies that the combined pressure component from the particular integral and the complementary function reduces for cases with shorter wavelengths. The results from this section show that for long wavelength roughness, the temperature and non-Newtonian lubricant behaviour have an effect on minimum film thickness in particular. It is thus important to make use of a thermal model in the study of two-sided surface features, which is presented in the next section.

Two-Sided Surface Features
The study of two-sided surface features was conducted using surface features with twice the amplitude compared to the aforementioned sinusoidal roughness, that is,Ā R,i = ±0.03. Moreover, a wavelength of Λ R,1 = Λ R,2 = 0.5 was employed, resulting in the surface features shown in Figure 9. . Asperity geometry (left) and dent geometry (right) using Equation (9) with the values A R,i = ±0.03, Λ R,1 = Λ R,2 = 0.5 and X 0 = 0.
The surface features were at first chosen to meet and completely overlap in the middle of the contact, that is the point of interference occurred at X = 0. The point of interference was later changed to X = −1 and X = 1, which is covered later in this section. The upper surface was moving slower at a speed of u 1 = 0.75 m s −1 and the lower surface moving faster with u 2 = 1.25 m s −1 . An overlap at X = 0 was achieved by setting the initial positions of the surface features to X 0,1 = −2.5 and X 0,2 ≈ −4.17. Four different combinations of overtaking situations were studied based on the work by Almqvist [55], which were defined as follows: • an asperity on the slower surface overtaken by an asperity on the faster surface (A + A), • a dent on the slower surface overtaken by a dent on the faster surface (D + D), • an asperity on the slower surface overtaken by a dent on the faster surface (A + D), • a dent on the slower surface overtaken by an asperity on the faster surface (D + A).
The case of two asperities overtaking, that is, the case A + A, is shown in Figure 10. Note that the split film thickness values in the left column are defined as half the combined geometry and half the combined deformation with the corresponding surface feature added, as defined by Morales-Espejel [23]. In the case of the contacting bodies having the same material properties, the split film thickness is defined according to: This improves visualisation of the problem and implies that the distance between S upper and S lower represents the film thickness H, which is shown in the right column in the figure together with the temperature field.
It can be seen that a complementary function is generated in the inlet based on the deformation of the surface feature in the Couhier region, similar to that of long wavelengths in Section 3.1. The complementary function moves with the mean entrainment speed and hence start to de-attach from the corresponding surface feature, which itself is significantly deformed. Moments thereafter, the other, faster, surface feature enters and generates another complementary function, smaller in size due to the increased speed, which implies a shorter time of inlet disturbance. This complementary function again moves with the entrainment speed and hence moves parallel in time with the first complementary function, indicating that the film thickness disturbances do not experience any interference.
Moreover, with the pressure and temperature being directly connected to the position of the surface features, as opposed to the complementary function, these values vary noticeably during the overtaking event and reaches their maximum values at X R1 = X R2 = 0. A similar behaviour was also seen for the case of two dents overtaking at X = 0, that is, the case D + D. However, the dents generate a complementary function in the form of a film thickness increase, as shown in Figure S4. Additionally, a significant pressure drop rather than a pressure spike occurs when the two dents overlap at X R1 = X R2 = 0.
In the case of a dent overtaking an asperity within the contact, that is, the case A + D, a similar behaviour of the complementary functions compared to the previously studied cases is seen, shown in Figure 11. Here, again, the formation of the complementary functions in the inlet moves through the contact with the mean entrainment speed. However, the pressure is experiencing negative interference as opposed to the cases A + A and D + D, meaning that at X R1 = X R2 = 0 an almost completely smooth pressure profile is seen and only small pressure ripples related to the complementary functions remain.  This further highlights that the influence on pressure from the complementary function is small in comparison to the particular integral. Consequently, the temperature is not as affected in this case due to the reduced amount of shear heating that is generated. A similar behaviour is seen for the case of a dent being overtaken by an asperity, that is, the case D + A, shown in Figure S5.
The film thickness values for all four cases are shown in Figure 12. It becomes clear that the two complementary functions formed in the inlet travel with the same speed throughout the contact in each of the studied cases, meaning that the film thickness variations are basically independent of each other. The width difference related to the time spent by the surface feature travelling through the contact inlet also becomes visible, which due to the limited flow into the contact results in an increasingly wide complementary function for the slower moving asperity. Moreover, only relatively small film thickness fluctuations are seen in-between the complementary functions that are related to the roughness deformation and hence the particular integral. It may also be seen that the thinnest film occurs in the outlet of the contact, which is especially pronounced for the case when an asperity exits the contact. Further studying the minimum film thickness over time between the cases reveals that the case A + A results in the thinnest film, followed by the case A + D, see Figure 13. The differences in minimum film during exiting of the surface features for the individual cases studied are explained by the combined exit of the faster surface feature with the complementary function of the slower surface feature, which occurs when the slower surface feature is located at X R,1 ≈ 0.5. The second significant film reduction taking place at X R,1 ≈ 0.9 is related to the exiting of the sole, slower, asperity. A similar behaviour can be seen for all four cases. Furthermore, it is interesting to note that the exiting of a dent in combination with the complementary function of an asperity gives a thinner minimum film thickness than that of an asperity exiting together with the complementary function of a dent, as shown by the differences between the cases A + D at X R,1 ≈ 0.5 and D + A at X R,1 ≈ 0.6.
Moreover, as noted before, the pressure is to a great extent related to the particular integral solution and thereby the surface feature itself. Hence, the two particular integrals interfere within the contact as opposed to the complementary waves, see Figure 14. This means that in the cases A + A and D + D, constructive interference between the particular integrals occurs and leads to a maximum and minimum pressure, respectively, when the surface features are located at X R1 = X R2 = 0. The same also applies to the cases A + D and D + A, however, with negative interference leading to an almost smooth pressure profile at X R1 = X R2 = 0. Only small pressure perturbations occurs away from the main pressure fluctuations, which are related to the complementary functions. The maximum pressure values over time are shown in Figure S6.  After studying the four different surface feature combinations overtaking at X = 0, a change of position was considered and applied to the case A + A. Now, the overtaking between the two asperities was set to occur at the outlet of the contact, implying that the complete overlap of the asperities takes place at X = 1. This reveals similar mechanisms as discussed above for the case of overlap at X = 0, that is, the complementary functions are formed in the inlet and moves parallel to each other throughout the contact, as shown in Figure 15. Additionally, due to the significant influence of the surface features on the outlet film thickness, the minimum film thickness becomes noticeably reduced during exiting of the two asperities together compared to the previously studied case. This is explained by the increased pressure driven flow related to the combined asperity exiting event. An additional change of overtaking position was also studied, that is, considering a complete overlap of the asperities at the inlet at X = −1. This revealed that in the case of overtaking occurring in the inlet, a combined inlet restriction from the two overlapping asperities occurs. Further implying that the complementary functions overlap and results in a single, combined, complementary function that moves through the contact at the mean entrainment speed, shown in Figure S7. The pressure does in this case, however, never experience interference due to the constant separation of the surface features throughout the contact.

Conclusions
In this work, an investigation of the rough thermal EHL problem was conducted using a previously developed thermal EHL model. The investigation was made up by two main analyses, one with focus on sinusoidal roughness and the differences between a thermal non-Newtonian approach and an isothermal Newtonian approach in the context of film formation with respect to the deformed roughness component and the generated inlet perturbations, also referred to as the particular integral and the complementary function, respectively. The other was focused on the transient event of two-sided surface features overtaking within the thermal EHL contact.
The general idea of the study was to learn more about rough EHL and the mechanisms of film formation. It was expected that thermal effects would make the minimum film thickness smaller than in the corresponding isothermal analysis. The results show that thermal effects must be taken into account when accurate values of minimum film thickness are sought for. The study can also be seen as a way to further investigate the mechanisms of film breakdown. Local thermal effects may lead to a reduction of viscosity that allow for asperities to grow out of their deformed state and thus causing metal-metal contacts. It is, however, not possible to estimate breakdown from the results presented here. Future studies of shorter wavelength roughness and a more pressure sensitive lubricant may potentially provide more information on mechanisms affecting film breakdown.
The salient conclusions from this study can be summarised as follows: (i) In case of sinusoidal roughness, differences in outlet film thickness between the thermal non-Newtonian and isothermal Newtonian approaches were for all studied wavelengths observed. These differences were explained by a reduced viscosity in the contact inlet due to shear thinning and shear heating, and further resulted in the isothermal Newtonian approach to overestimate the minimum film thickness by up to 20 %. (ii) Differences in the high pressure region were only found to be significant in case of short wavelength roughness. These differences were explained by a reduced viscosity in the contact inlet due to shear thinning and shear heating, further disturbing the formation of the complementary function. This resulted in central film thickness estimations to differ between 20 % and −14 % by using the thermal non-Newtonian approach compared to the isothermal Newtonian approach. (iii) Following the relatively small influence on pressure by the complementary function, only small differences in maximum pressure between the thermal non-Newtonian and isothermal Newtonian approaches were noticed. (iv) In case of two-sided surface features overtaking within the contact, it was found that an interference between the film thickness perturbations did not occur due to the complementary function being generated in the inlet of the contact and travelling with the speed of the lubricant. On the other hand, due to the particular integral being directly connected to the roughness feature, an interference between the pressure and temperature variations were obvious.
In case of the overtaking event occurring in the inlet of the contact, one single complementary function was formed with larger amplitude than the other cases. Moreover, if overtaking instead took place in the outlet of the contact, a significantly reduced film thickness was noted. (vi) For the studied cases, it can be concluded that a thermal non-Newtonian approach is required for quantitatively accurate outlet film thickness predictions and may also be necessary for accurate predictions of film thickness within the contact if the complementary function is affected during formation. Otherwise, in the case of long wavelength roughness or low sliding conditions, an isothermal approach may be sufficient due to the inlet dominated film formation.

Funding:
This research is part of the FFI-project P41215-1 funded by the Swedish Energy Agency (Energimyndigheten).

Acknowledgments:
The authors wish to thank Tomas Johannesson at Volvo Car Corporation for all the fruitful discussions and valuable input to this work.

Conflicts of Interest:
The authors declare no conflict of interest. Volume of lubricant at P 0 and T r , respectively m 3 x, y, z Spatial coordinates (m) x R,1 , x 0 Current-and initial position of roughness, respectively (m) X, Y, Z Dimensionless spatial coordinates (−) X R , X 0 Current-and initial dimensionless position of asperity, respectively (m)