Evaluation of Various Shear-Thinning Models for Squalane Using Traction Measurements, TEHD and NEMD Simulations

Abstract

The accurate prediction of friction in highly loaded concentrated contacts is one of the most challenging aspects of thermal elastohydrodynamic (TEHD) simulation. The correct modelling of fluid behaviour on the macroscale, in particular non-Newtonian flow behaviour, is an essential prerequisite. For many years, shear-thinning models have been developed and validated with different approaches and controversially discussed. In basic research, model fluids are often used in this context, which have a similar behaviour to practical lubricants. Accompanied by earlier research results, this paper carries out comprehensive investigations on the rheometric behaviour of the model fluid squalane. Based on traction measurements at four different tribometers, an overall parameter optimisation and performance evaluation of three different shear-thinning models is performed using numerical TEHD simulations. In order to additionally validate the theoretical viscosity behaviour, the optimised shear-thinning curves are then compared with comprehensive non-equilibrium molecular dynamics (NEMD) simulations. The key aspect of this paper is the simultaneous consideration of the shear-thinning models in terms of rheometric, experimental, and simulative investigations without changing the parameters. All investigations show that the Eyring model, despite its simplicity, provides the best agreement in both the numerical contact simulation and the NEMD simulations.

1. Introduction

When optimising machine elements with concentrated contacts, as found in gears, rolling bearings, or valve trains, different objectives are pursued depending on the application. On the one hand, the intention is to create a separating lubricating film to prevent wear or increase the lifetime; on the other hand, the focus is on minimising friction to achieve the highest possible efficiency of the tribological system. The objective is an optimum of high efficiency and long lifetime.

Often, lubricated systems are designed according to the trial-and-error method based on many years of experience and a large number of tests. This approach reaches practical limits because of the high costs involved and the increasingly shorter product development cycles. Further optimisations require a precise understanding of the contact processes in the lubrication gap. This knowledge can be developed using powerful sensor technology or numerical thermal elastohydrodynamic (TEHD) simulations. Numerical simulation tools enable the analysis of all variables in the contact area, such as hydrodynamic flow, lubrication gap height, elastic deformation of solids, pressure distribution, or temperature development in fluid and solids. These simulations are carried out either using CFD calculation tools [1,2] or software based on the generalized Reynolds differential equation, which additionally allows the consideration of mixed friction or the influence of rough surfaces [3,4]. The decisive advantage of simulations over measurements is that all variables and processes can be analysed locally and time-resolved, thus enabling targeted optimisation of systems.

For a valid friction calculation, the modelling of lubricant behaviour is an essential aspect. Lubricants usually have a complicated chemical composition, which leads to a complex behaviour. Especially in concentrated contacts, this must be considered precisely due to the occurring high pressures (up to 3 GPa) and shear rates (up to 1 × 10⁸ s⁻¹). In classical TEHD simulation tools, the lubricant behaviour in the gap must be described as a continuum for reasons of computing time. The objective and challenge of fluid modelling is, therefore, to develop a model that is generally applicable.

The viscosity η(ϑ, p, $\dot{\gamma }$) has the main influence on lubricant film thickness and friction caused by shearing. Furthermore, the density ρ(ϑ, p), as well as the specific heat capacity c_p(ϑ, p) and the thermal conductivity λ(ϑ, p), which influence the temperature distribution in the contact, must be taken into account.

1.1. Viscosity

The viscosity behaviour of a lubricant is highly non-linear. Various models have already been developed for the continuum-mechanical description of this behaviour, which are more or less suitable depending on the lubricant used. At this point, only a brief overview of such modelling will be given here. A detailed presentation can be found in [5]. Usually, the viscosity η of a lubricant is modelled by a combination of the so-called low-shear viscosity η₀, depending on temperature and pressure, and a function for shear rate dependency (1):

$$\eta =f\left({\eta }_0,\dot{\gamma }\right) with {\eta }_0=f\left(\vartheta ,p\right)$$

(1)

The low-shear viscosity η₀ of lubricants decreases with rising temperature. The temperature dependence at ambient pressure η₀(ϑ, p = 0) (Figure 1, left) can be easily measured on rheometers. Usually, rotational viscometers [6], oscillating viscometers [7] or falling-body viscometers [8] are used for this purpose. Starting from this value, the low-shear viscosity η₀ increases approximately exponentially with pressure (Figure 1, centre), which is shown in a single logarithmic plot as a straight line with the slope α_p. Actually, the isothermal curves in this representation often have a slightly degressive or degressive-progressive progression, which indicates a transition to the glassy state [9]. In addition, the curves usually spread, i.e., α_p = f(ϑ, p). This low-shear viscosity is often determined in rheology with oscillating viscometers or falling-body viscometers [8] and can be interpolated with various equations [5]. The results are well reproducible, but the maximum attainable pressures and measurable viscosities in a rheometer are limited, so that the relationships found must be extrapolated into the range of high pressures for application in TEHD contacts.

1.2. Shear-Thinning

The viscosity η₀ determined for a certain temperature and pressure can be considered constant at low shear rates (Newtonian behaviour). However, for most lubricants, the viscosity η decreases again at high shear rates (Figure 1, right). At ambient pressure, this shear-thinning only starts at very high shear rates, however, the transition shear rate decreases with rising pressure. It is, therefore, assumed that the shear stress τ occurring in the fluid is the decisive parameter for the transition to shear-thinning.

Since this effect is decisive for valid numerical friction calculations, it is the subject of current research. The cause and the best modelling have always been the subject of controversial discussions, as described in detail in the review paper [10] and its comments [11,12]. While Bair et al. [11] assume that the orientation of molecules under shear is responsible for the decrease in viscosity and suggests the application of the Carreau model from the findings of their rheometric experiments [13], Spikes et al. [12] argue that spherical molecules, such as cyclohexane, also show shear-thinning. According to Spikes et al., experiments on model test rigs can best be calculated using a hyperbolic sine function, which is given by the Eyring model [12,14]. Bair rejects this model for the description of shear-thinning [13].

Overall, the shear behaviour of a lubricant is very difficult to determine. There are three different approaches, as described below.

1.2.1. Option I—Rheometry

First of all, rheometric measurements can also be carried out for the analysis of shear-thinning. The most commonly used viscometer for measuring the so-called high-temperature high-shear viscosity at ambient pressure is the tapered bearing simulator according to ASTM D4683-20 [15]. The central problem of these investigations is that a high shear is always accompanied by a temperature increase in the liquid in the measuring gap, which means that the temperature distribution over the gap height is unknown [10]. Due to the shear heating, the viscosity in the measurement appears to decrease in addition to the shear-thinning. The influence of shear heating and shear-thinning cannot be clearly separated in the measurement data afterwards. Newer concepts (ultra-shear viscometer) try to reduce this problem by conditioning the solids, with extremely short measuring times and online temperature measurement [16].

In addition to the shear behaviour at ambient pressure, the behaviour at very high pressures is particularly relevant for concentrated contacts. For this purpose, Bair [17] developed a high-pressure high-shear stress viscometer, allowing the investigation of pressures up to approximately 1 GPa. The behaviour of fluids at very high shear rates and pressures was also investigated by Jacobsen [18] using high-pressure chamber measurements. However, in high-pressure investigations, the deformation of the measuring gap, the mechanical strength and sealing of the measuring chamber, and the related secondary friction represent an additional challenge. The maximum achievable shear rates are also limited by the rapid shear heating [10].

1.2.2. Option II—Traction Measurements

Another possibility to determine the flow behaviour of lubricants under real conditions is to carry out traction measurements on model tribometers with very precise friction measurement (Figure 2) [10,14,19]. The simple test specimens (ball, plane, or disc) are equipped with extremely smooth surfaces so that a very homogeneous lubrication gap is created and an influence of the flow by roughness or mixed friction is avoided. It is assumed that the measured friction is generated entirely by hydrodynamic shear in the lubrication gap. By varying the oil supply temperature, load, speed, and slide-to-roll ratio (SRR), the traction behaviour of the lubricant can be determined for different nominal operating conditions. The problem is that the measured frictional force only represents an integral measured variable, which results from the sum of the varying conditions at each point of the gap (pressure, temperature, and lubrication gap height).

The flow behaviour can, therefore, only be determined indirectly by simulating the contact conditions and comparing the simulated integral friction force resulting from the shear stress on the surfaces with the measurements. By varying the assumed flow behaviour in the simulation (model equation and parameters), the actual fluid behaviour can be estimated iteratively [20,21,22].

1.2.3. Option III—Molecular Dynamics Simulations

Both rheometric lubricant investigations and traction measurements deal with the behaviour of a lubricant on macroscopic scales and attempt to describe these with model equations. In order to better understand the causes for the behaviour of lubricants, molecular dynamics simulations of fluids under EHD-like conditions have increasingly been carried out in recent years. Synthetic fluids with a simple structure, mostly pure substances [23,24,25] or synthetic lubricants that are composed of a small number of molecules, e.g., PAO [26,27], are considered.

This method cannot yet be used for fully formulated oils due to their complexity (additives). The aim is, therefore, to build up an understanding of the basic shear behaviour and to transfer this knowledge to the description of more complex lubricants. A comprehensive overview of the current state of work is given in [28].

Among the various methods for molecular dynamic viscosity determination, non-equilibrium molecular dynamics (NEMD) simulations are often used for lubricants [23,25,26,27], in which a planar Couette shear is applied to a certain number of molecules (N) in a cell with a fixed volume (V). The walls of the cell are provided with periodic boundary conditions, so that a small section of a larger fluid volume is considered. The energy input by the shear is eliminated by thermostatting conditions (usually the Nosé-Hoover thermostatting [29] is used) so that constant temperature conditions (T) can be considered at extremely high shear rates. To apply the planar shear, the NVT-SLLOD equations of motion [30] are solved. The viscosity is eventually derived from the resulting stress tensor. Furthermore, molecular structural changes and orientation can be studied. To achieve the desired pressure (P) and density for a given system, NPT pre-simulations with varying cell size are carried out.

To investigate the flow profiles in small lubricating gaps, NEMD simulations can be extended by modelling solids. Here, the planar shear is applied by the moving wall, which enables the detailed investigation of molecular arrangement, molecular orientation and solid–fluid interaction [31].

Due to statistical fluctuations and the very small time span that can be simulated NEMD simulations can only be used reliably for very high shear rates of about $\dot{\gamma }$ > 1 × 10⁶ s⁻¹ [28]. Low shear rates, on the other hand, require increasingly longer calculation times or are subject to increasingly larger statistical uncertainties (signal-to-noise ratio increases). However, especially at high pressures, shear-thinning begins at much lower shear rates; thus, the transition from Newtonian to non-Newtonian flow behaviour and its functional characteristics cannot be evaluated sufficiently well with NEMD simulations. Coarse grain models would be a promising alternative to atomistic models in order to overcome the short time and length scale limitations. However, the parameterisation of these models to accurately predict lubricant rheology is very difficult [28].

Another possibility is to perform equilibrium simulations according to the Green–Kubo approach [27,32,33,34]. Here, the calculation cell is not sheared, but the viscosity at a certain pressure and temperature is related to the time integral of the autocorrelation function of the stress. This value, thus, represents the pressure viscosity at zero shear (Figure 3, intersection with the y-axis). However, this procedure is problematic for high pressures, as the time span to be simulated increases until a reliable statistical evaluation can be obtained [27]. This uncertainty can be reduced by a time decomposition approach [35].

Due to the limitations of all aforementioned methods, it can be stated that only the combination of rheometric measurements, molecular dynamics simulations, and the application of the models in TEHD simulations in comparison with traction measurements can improve the understanding of shear-thinning behaviour.

1.3. Limiting Shear Stress

The analysis of flow behaviour is directly related to the issue of the limiting shear stress τ_lim. The existence of such a limiting shear stress is controversially discussed in the literature [10]. The concept is that the maximum shear stress that can be transmitted by a fluid is limited, even if the shear rate is constantly increased. This hypothesis can be justified by observations of highly loaded traction curves, which always show a maximum value in the friction. If these maximum values are correlated with the Hertzian pressure, a linear relationship between pressure and limiting shear stress can be established. This increase is only noted above a certain pressure. A slight decrease is recorded with increasing temperature [20,21,22]. At this point, however, it should be pointed out once again that the derivation of a limiting shear stress from traction tests is fraught with uncertainties, since only an integral consideration of the fluid behaviour is possible and the boundary conditions differ locally (p = f(x,y); (ϑ, $\dot{\gamma }$) = f(x,y,z)). Traction tests are always accompanied by shear heating, which also influences the maximum value in the traction curve. Especially for highly loaded contacts, simulations can show considerable temperature developments; thus, curves with a significant maximum can also be calculated without applying a limiting shear stress. In recent investigations, Zhang and Spikes [14] could not find any evidence for limiting shear stress from their traction measurements. Additionally, molecular dynamics simulations, such as those found in [23] or [25], do not show any limitation of the shear stress for the cell, even for extreme shear rates.

In contrast, rheometric experiments with the Couette viscometer presented by Bair et al. [36,37] indicate a limiting shear stress. Furthermore, shear bands and intermittent shear localisation could be observed in a high-pressure flow visualisation cell for the fluids 5P4E and MCS1218, which are associated with a limiting shear stress [38]. Jacobsen [18] confirms a pressure-dependent shear stress increase in his investigations on various practical fluids.

When applied in simulations, the limiting shear stress is usually modelled as a simple linear [19,39] or bilinear [20,22,40] dependence on pressure, with the temperature dependence being neglected in most models.

1.4. Model Fluid Squalane as Object of Investigation

To get closer to the description of the flow behaviour of lubricants in the TEHD regime, squalane C₃₀H₆₂ has often been used in recent years as a synthetic high-purity substance with an extremely simple, defined structure. It is a short hydrocarbon with side branches and has a similar viscosity behaviour as fully formulated lubricants. The advantage of this substance is that all three described methods (options I to III) can be used to determine the lubricant behaviour. Rheometric measurements provide the temperature and pressure dependence of the viscosity, and with traction measurements, molecular dynamics simulations, and supplementary rheometric measurements, the flow behaviour can be analysed. Hence, there is a broad data base that can be relied on.

Besides many older works on the determination of viscosity at moderate pressures, the first high-pressure viscosity measurements with a falling-body viscometer up to approx. p = 1.2 GPa at 40 °C and 100 °C and the measurement results on shear rate dependence with a Couette viscometer were presented by Bair [17] in 2002. At the same time, Bair et al. [23] carried out extensive molecular dynamic calculations with a united atom model for many different conditions, despite the limited computing technology available in those days. The considered pressure range of p = 316 MPa was still very limited with regard to TEHD contacts. However, in combination with experimental data, a complete Shifted-Carreau flow model could be derived by reducing the curves found. In 2006 [8], Bair published further extensive low-shear viscosity measurements of squalane with multiple falling-body viscometers for pressures up to 1.2 GPa and improved the model. In addition, comparative data for temperature dependence at ambient pressure are given in [8], which were determined in [41]. In 2014, Bair et al. [42] investigated the relationship of viscosity of squalane measured stationary under pressure and oscillating under ambient pressure. In 2018, Bair et al. [43] investigated the specific heat capacity and thermal conductivity of squalane.

In 2013, Björling et al. [19] carried out traction measurements with squalane on a WAM11 tribometer (ball on disc). These measurements were then compared with numerical traction calculations using the Shifted-Carreau model derived by Bair et al. In this model, the shear stresses were limited to a pressure-dependent limiting shear stress. Due to an error in the calculation, which is corrected in [19], the results in [19] can only be evaluated qualitatively. It can be stated that the calculated friction is dominated in a wide range by the limiting shear stress.

Liu et al. [24] confirmed with NEMD simulations the shear-thinning curves calculated by Bair et al. [23], using a novel method to derive the viscosity from the simulations. Jadhao et al. [25] performed NEMD simulations up to p = 1 GPa in 2017 and were able to fit the reduced results very well with the Eyring model. In 2017, Ewen et al. [31] considered solid surfaces in their NEMD simulations and investigated in detail the flow profiles and mass distribution in the lubrication gap.

1.5. Conclusion and Aim of the Study

When analysing the individual publications, it is noticeable that the fitting of the respective data (rheometric measurements, MD simulations, and traction measurements) is mostly carried out under the assumption of a certain functional relationship (e.g., Eyring and Carreau). An evaluation of the best fit of different functions is hardly ever carried out. If additional TEHD simulations are performed to validate the derived fluid model, these are often based on a single tribometer or a small number of experiments. However, the investigations presented here show that it is relatively easy to find parameter sets for different shear models (Eyring, Carreau, etc.) that reproduce the behaviour of a single tribological system well in the simulation. However, if this fluid model is then applied to another system, the results are often worse.

In order to be able to assess the performance of individual shear-thinning models, it is essential not only to consider a single aspect but also to simultaneously take into account rheometry, MD simulations, and the application to traction calculation without changing the parameter sets. From our point of view, a central validation aspect is the representation of the fluid traction behaviour in arbitrary tribometers.

This holistic analysis is carried out in this paper. The data from different publications of previous years are combined and supplemented by our own measurements and calculations. The aim is a systematic comparison of different models for shear-thinning using different tribological systems. The procedure is shown in Figure 4.

First of all, equations for the description of the pressure and temperature dependence of the viscosity are derived on the basis of rheometric measurements from the literature. Using three selected shear-thinning models (Eyring, Eyring-Hyperbolic, and Carreau), numerical TEHD simulations are carried out for four different tribological systems. By comparing the simulation results with corresponding test rig measurements, a parameter optimisation of the models is carried out in such a way that all measured data are matched as well as possible. Subsequently, the applicability of the optimised shear-thinning models is evaluated. For theoretical evaluation, the viscosity behaviour is finally compared with comprehensive NEMD simulation results.

Furthermore, the rheometric measurements, the traction calculations, and the MD simulations are each compared with the Shifted-Carreau model for squalane developed by Bair [8].

2. Materials and Methods

2.1. Fluid Modelling

To model the low-shear viscosity as a function of temperature and pressure, the data from several sources were combined, which were determined using different measurement methods. For the temperature dependence of the viscosity at ambient pressure (Figure 5a), measurements were carried out using a rotational viscometer for ϑ = 20 °C … 120 °C [22]. The pressure dependence (Figure 5b) was determined using a high-pressure quartz viscometer for 20 °C, 40 °C, and 90 °C [22]. These data were compared with measurements from the literature [8,17,41].

Figure 5 shows that the measurement series of the different rheometers are very consistent with each other. However, it is also apparent that the available values are limited to pressures p < 1500 MPa. In order to represent the relevant pressure range up to 2.7 GPa in the TEHD contacts investigated, the curve must, therefore, be extrapolated into the range of high pressures using the selected equation.

On the basis of the available measurement values, an interpolation fitting was carried out so that all data are matched as well as possible. To describe the slightly degressive pressure–viscosity curve (in the semi-logarithmic representation), the less common equation according to Rodermund [44] (Equation (2)) was used. With the parameters A, B, and C, this contains the well-known Vogel equation. The parameters D, E, and p₀ describe the pressure dependence, whereby a curvature of the curves can be modelled.

$${\eta }_0\left(\vartheta ,p\right)=A·e^{\left[\frac{B}{C+\vartheta }{·\left(\frac{p}{p_0}+1\right)}^{\left(D+E·\frac{B}{C+\vartheta }\right)}\right]}$$

(2)

where $A=0.1014 \mathrm{m}\mathrm{P}\mathrm{a}\mathrm{s}$; $B=674.9375 °\mathrm{C}$; $C=95.825 °\mathrm{C}$; $D=0.6377$; $E=0.0004$; and $p_0=182.56 \mathrm{M}\mathrm{P}\mathrm{a}$.

For both the temperature and pressure dependency, good agreement with the measured data can be achieved with Equation (2). In the range of high pressure, a further degressive progression is assumed.

In comparison, Figure 5 also shows the curves of the interpolation fitting of the Tait-Doolittle equation as it can be found in [8] or has been applied in a similar form in [19]. At ambient pressure (Figure 5a), slightly higher values at low temperatures and slightly lower values at high temperatures are calculated than those measured with the given parameters. For the pressure dependence (Figure 5b), the curves are very similar to Equation (2) in the range of the measurement data, but at very high pressures p > 1500 MPa, they show an incipient glass transition and, thus, increasingly larger viscosity values.

Altogether, the question arises which pressure–viscosity relation is more appropriate for the application in TEHD simulations. While the measurement in rheometers indicates solidification for very high pressures, our investigations show that the assumption of a further degressive curve, in connection with the available shear-thinning models, leads to better agreements with the traction measurements than a progressive approach.

Based on the low-shear viscosity η₀ at a given temperature and pressure (Equation (2)), the shear behaviour is described functionally. Three different models (Eyring, Eyring-Hyperbolic, and Carreau) are applied, which equations are listed in

Table 1. Shear-thinning models.

Model	Equation		Ref.	Used Parameters for Squalane	Limiting Shear Stress
Eyring	$\dot{\gamma }=\frac{{\tau }_c}{{\eta }_0\left(\vartheta ,p\right)}sinh\left(\frac{\tau }{{\tau }_c}\right)$	(3)	[45]	$a=0.08\frac{\mathrm{M}\mathrm{P}\mathrm{a}}{\mathrm{K}}$ $b=0.003$ ${\vartheta }_0=0 °\mathrm{C}$ ${\tau }_{c0}=5 \mathrm{M}\mathrm{P}\mathrm{a}$ ${\eta }_0\left(\vartheta ,p\right) acc.Rodermund$	No
	${\tau }_c={\tau }_{c0}+a\left(\vartheta -{\vartheta }_0\right)+b·p$ ${\tau }_c>0$	(4)
Eyring- Hyperbolic	$\dot{\gamma }=\frac{{\tau }_c}{{\eta }_0\left(\vartheta ,p\right)}sinh\left(\frac{\tau }{{\tau }_c}\right)·\frac{1}{1-{\left(\frac{\tau }{{\tau }_{lim}}\right)}^t}$	(5)	[22,46]	${\tau }_c=20 \mathrm{M}\mathrm{P}\mathrm{a}$ $t=2$ ${\tau }_{lim,0}=40 \mathrm{M}\mathrm{P}\mathrm{a}$ $\beta =0.09$ $p_{lim,0}=400 \mathrm{M}\mathrm{P}\mathrm{a}$ ${\eta }_0\left(\vartheta ,p\right) acc.Rodermund$	Yes
	${\tau }_{lim}=\left\{\begin{array}{ll}{\tau }_{lim,0}·e^{{\beta }_{exp}·p}& for p\le p_{lim}\\ \beta ·\left(p-p_{lim,0}\right)& for p\ge p_{lim}\end{array}\right.$	(6)
Carreau	$\dot{\gamma }=\frac{\tau }{{\eta }_0\left(\vartheta ,p\right)·{\left(1+{\left(\frac{\tau }{{\tau }_c}\right)}^2\right)}^{\frac{1-\frac{1}{n}}{2}}}$	(7)	[47]	${\tau }_c=80 \mathrm{M}\mathrm{P}\mathrm{a}$ $n=0.085$ ${\eta }_0\left(\vartheta ,p\right) acc.Rodermund$	No
Shifted- Carreau	$\eta \left(\mu ,\dot{\gamma }\right)=\mu \left(T,p\right)·{\left[1+{\left(\frac{\mu \left(T,p\right)·\dot{\gamma }}{G}\right)}^2\right]}^{\frac{n-1}{2}}$ $G=G_R\frac{T}{T_R}\frac{V_R}{V}$	(8)	[8]	$n=0.463$ $G_R=6.94 \mathrm{M}\mathrm{P}\mathrm{a}$ $T_R=40 °\mathrm{C}$ $\mu \left(T,p\right) acc.Tait Doolittle$ $\Lambda =0.075$	Yes
	${\tau }_{lim}=\Lambda ·p$	(9)	[19,39]

, especially models with and without limiting shear stress are compared so that the effects on the traction curves can be verified. For the low-shear viscosity η₀, the curves modelled with the Rodermund equation (Equation (2)) are considered as defined and no longer changed.

Furthermore, the Shifted-Carreau model [8], which can be found in the literature, is considered. This remains completely unchanged in the following investigation, meaning that the low-shear viscosity μ corresponds to the Tait–Doolittle equation derived in [8] and shown in Figure 5.

Figure 6 shows the principal progressions of the individual models for the shear viscosity η(η₀, $\dot{\gamma }$) (top) and the shear stress (left, bottom) versus the shear rate and their degrees of freedom. Due to its simplicity, the Eyring model is probably the most widely used shear-thinning model. Its only parameter is the critical shear stress τ_c. This describes the shear stress at which Newtonian behaviour is left. The viscosity slope in the double logarithmic plot is predetermined by the hyperbolic sine function in the Eyring model and cannot be changed. Consequently, the shear stress is strongly degressive, but it is not mathematically limited for large shear rates. Investigations using the Eyring model by Johnson et al. [48] suggested that τ_c should be described linearly with increasing temperature and pressure. In order to have more degrees of freedom in the modelling, this approach is applied in this paper (Equation (4)).

In order to include a limiting shear stress in the Eyring model, the Eyring-Hyperbolic model, as shown in Equation (5), was developed [22] and has already been applied to various systems [40]. In this equation, a transition to the limiting shear stress is described mathematically by a scalable hyperbolic function following [46]. The result is a shear stress curve that corresponds to the Eyring curve for τ << τ_lim (described by τ_c) and is gradually restricted by the limiting shear stress τ_lim at higher shear rates. This moderate approximation to the limiting shear stress proves to be advantageous compared to other shear-thinning models with regard to the results and the numerical calculation stability. The pressure dependence of the limiting shear stress τ_lim is described by a continuous exponential-linear approach, Equation (6) (Figure 6, bottom right). The slope β of the linear region and the transition pressure p_lim are derived from the maxima of the traction measurements, and the limiting shear stress at ambient pressure τ_lim,0 is determined iteratively. Since the progression of the shear stress can be modelled very well by the limiting shear stress τ_lim, the critical shear stress τ_c is assumed to be constant. A total of 5 parameters must be found for the Eyring-Hyperbolic model.

Additionally, widely used are various forms of the Carreau model, which has a power-law behaviour. The version used here (Equation (7)) is presented in Bair et al. [47]. As shown in Figure 6, the transition to shear-thinning (τ_c) and the slope of the viscosity or the increase in shear stress at high shear rates (n) can be varied. Thus, only two parameters need to be determined. The shear stress in Equation (7) is not limited.

The Shifted-Carreau model basically has a similar behaviour. At this point, a discussion of the degrees of freedom will be omitted as the model is only used for comparison purposes. More details can be found in [8]. When the model was applied in TEHD simulations in [19,39], a limiting shear stress τ_lim linearly dependent on the pressure was added to the model. The coefficient Λ was derived from the traction measurements. This approach is adopted and evaluated here.

2.2. Traction Measurements

The tests on four different tribometers serve as the basis for deriving the coefficients and evaluating the shear-thinning models mentioned. Here, measurements from the literature are used and our own tests are supplemented. The aim is to use very different systems with large and small test specimens and with large and small Hertzian contact areas. In this way, a wide range of parameters is taken into account, and the most reliable statements possible can be made.

Figure 7 and

Table 2. Traction test rigs. Operating parameters in bold correspond to tests shown in the results Section 3.

	WAM11 (Data from [19])	Large Two-Disc Machine (Data from [22])	Small Two-Disc Machine (Own Measurements)	EHL2 Tribometer (Own Measurements)
Body 1	ball	disc	disc	ball
Dimensions of body 1	Ø 20.64 mm	Ø 120 mm crowning: R50 mm	Ø 20.4 mm crowning: R100 mm	Ø 19.05 mm
Body 2	disc	disc	disc	disc
Dimensions of body 2	raceway Ø 40 mm	Ø 120 mm cylindrical	Ø 20.4 mm crowning: R100 mm	raceway Ø 80 mm
Material	AISI52100 (100Cr6)	AISI52100 (100Cr6)	AISI52100 (100Cr6)	AISI52100 (100Cr6)
Roughness	ball: Rq = 25 nm disc: Rq = 35 nm comb.: Rq = 43 nm	body 1: Ra = 63 nm body 2: Ra = 26 nm comb.: Ra = 68 nm	Ra1,2 < 20 nm	Ra1,2 < 20 nm
Applied force	50 N, 300 N	1180 N, 1920 N 3750 N, 6500 N	97 N, 230 N, 396 N, 629 N, 938 N, 1333 N	8 N, 27 N, 50 N
Max. Hertzian Pressure	1070 MPa, 1940 MPa	1275 MPa, 1500 MPa 1875 MPa, 2250 MPa	1200 MPa, 1500 MPa 1800 MPa, 2100 MPa 2400 MPa, 2700 MPa	600 MPa, 900 MPa 1100 MPa
Oil feed temperature	40 °C	0 °C, 20 °C, 40 °C	40 °C, 60 °C	40 °C, 60 °C, 80 °C
Mean velocity	0.34 m/s … 9.6 m/s 2.51 m/s, 9.6 m/s	5 m/s, 10 m/s, 15 m/s	5 m/s	2 m/s, 3 m/s
SRR	0.2% … 49%	−15 … 15%	−10 … 10%	−50 … 50%

show an overview of the test rigs, the used test specimen geometry, and the operating parameters investigated. The operating parameters marked in bold are representative tests that are shown in the following figures. However, all the following statements also apply to the other operating points. All specimens are made of steel and have an extremely smooth surface.

First of all, the friction measurements of Björling et al. [19] on a Wedeven Associates Machine (WAM) 11 are used. The WAM11 is a ball-on-disc tribometer that allows very high pressures up to 2.91 GPa. The axis of the ball is positioned at an angle to avoid spin slip. The lubricant is fed at a controlled temperature in the centre of the disc and, thus, supplies the contact.

Furthermore, the results from [22] are used. A large two-disc machine with a disc diameter of 120 mm is used at the IMKT of the Leibniz University Hannover. This test rig allows very large contact areas and high velocities. Furthermore, the lubricant can be cooled so that oil temperatures of ϑ_oil = 0 °C are possible. The lubricant is supplied from above between the discs. The frictional torque is recorded via a rotatably mounted drive machine. The mass temperature of the discs is recorded and can be used as a boundary condition for the TEHD simulation.

Our own tests were carried out on a small two-disc machine. Here, both discs are provided with a crowning so that a very small contact area with extremely high pressure can be examined. The lubricant is also fed between the discs from above at a controlled temperature, and the disc temperature is measured. The friction torque is measured by a calibrated tilting balance.

The fourth test rig used is the EHL2 tribometer from PCS instruments. In this ball-on-disc tribometer, the ball is located underneath the disc and feeds the lubricant from a temperature-controlled oil bath into the contact. The Hertzian contact is similar to the WAM11, but the rotational axes of the specimens are at 90° to each other. The friction torque is measured via a torque sensor at the ball drive. Due to the very smooth surfaces, higher lubricant temperatures with small lubrication gap heights can be investigated. The temperature rise in the components due to friction is very low and can be neglected.

2.3. Iterative Optimisation Using TEHD Simulations

To optimise the model parameters, iterative TEHD simulations were carried out for the four tribometers. The accuracy of the simulated traction curves was evaluated simultaneously for all test rigs and operating points by using the method of least squares between the measurements and calculations. As the parameters of Equations (3)–(7) have certain effects on the curve progression, they can be specifically adjusted to reach a minimum error.

Only a brief overview is given of the numerical simulation model used. The calculations were carried out with the software Tribo-X from Tribo Technologies. Detailed theoretical information on the program can be found in [3,4,40,49,50]. The basis for the elastohydrodynamic calculation is the generalised Reynolds differential equation for compressible and laminar flow. The density of the fluid ρ = f(ϑ, p) and the viscosity η = f(ϑ, p, $\dot{\gamma }$) described by the shear-thinning models are taken into account variably in all spatial directions (x, y, and z) so that an equilibrium of acting shear stress and shear rate is always established in the gap height direction z, leading to a non-linear velocity profile. The deformation of the solid surfaces is represented by the elastic half-space. The resulting deformation slip is taken into account. To represent the shear heating and heat dissipation from the contact, the transient Fourier heat conduction equation is solved for the solids and the energy equation is solved for the fluid, taking into account convection, heat conduction, friction, compression, and expansion.

The thermal conductivity and specific heat capacity of the fluid are described as pressure and temperature dependent, analogous to the models of Björling et al. [19]. For the solids, a thermal conductivity of λ = 33 W/m∙K is specified. In the case of the two-disc machines, the mass temperature measured in the test specimens is specified as the boundary condition. Since all the systems considered have extremely smooth surfaces and the operating points lead to sufficiently large lubricating film, roughness is not considered, and the surfaces are regarded as ideally smooth.

3. Results and Discussion

The results of the parameter optimisation are shown in the following plots, and the values can be found in Table 1. It should be noted that it would be possible to find separate parameters for each test rig so that there is even better agreement between the measurement and the simulation. However, the aim of this paper is to find a common set of parameters for all tribological systems.

3.1. WAM11

Figure 8 shows the results for WAM11 [19]. As expected, the measurement shows that the friction coefficient increases with increasing pressure and decreases with increasing mean velocity. The greater the pressure and velocity (Figure 8, bottom right), the greater the drop in friction at high SRR. First of all, it can be stated that all shear-thinning models are able to represent the basic behaviour of the system. It is remarkable that even the models without limiting shear stress reproduce the pronounced maximum in the traction curve. This suggests that the friction drop results primarily from strong shear heating in the gap and does not necessarily indicate a limiting shear stress.

The Eyring model shows good agreement for all operating points. The shape of the curves corresponds with the measurement. With the Eyring-Hyperbolic model, friction is overestimated at low pressures and low SRR. At high pressures, the restriction by the limiting shear stress is visible (horizontal curve) so that the maximum measured friction is not reached. The Carreau model also shows good agreement with the measured values.

The results of the Shifted-Carreau model correspond qualitatively to those presented by Björling et al. in [19] and quantitatively to the corrected results in [19]. While a very good representation of the traction behaviour is given at low load, the model is visibly dominated by the limiting shear stress at the highly loaded contact. This hard limit of the shear-thinning curve leads to horizontal traction curves and kinks. This discontinuous behaviour is consistent with the results presented in [19] but is not realistic.

To illustrate the calculation, Figure 9 shows the results of the TEHD simulations of an operating point of the WAM11 marked in Figure 8, which are exemplary for the Eyring model. In addition to the known distributions for pressure and film thickness, it can be clearly seen that the viscosity in the contact area increases strongly due to the hydrodynamic pressure development. Since the selected point has a slip of SRR = 50%, there is considerable temperature development in the lubricant due to shearing of the fluid, whereby heat is dissipated by the fluid in the direction of the flow as well as on both surfaces of the solids.

The non-linear flow profile in the direction of the gap height results from the velocity boundary conditions of the surfaces (SRR) and at each point of the gap from the equilibrium of shear stress τ and shear rate $\dot{\gamma }$ according to Equation (3). This equilibrium (which corresponds to viscosity η) is strongly influenced by the temperature distribution in the gap; thus, the shear rate is greater in the centre of the gap than at the solid surfaces. In this case, shear rates of about $\dot{\gamma }$ = 1 × 10⁷ s⁻¹ are reached. In the area at the end of the gap, these rates are even higher due to strong pressure flows. The resulting shear stress at the surface represents the calculated frictional force.

3.2. Large Two-Disc Machine

Figure 10 shows the results for the large two-disc machine, which has a very large Hertzian contact and a high mean velocity. The strong decrease in friction with increasing SRR due to the strong shear heating is apparent. The decrease in friction with increasing temperature (from top to bottom) and the increase in friction with increasing pressure (from left to right) can also be seen clearly.

Additionally, in these tests, the Eyring model shows overall very good agreements and a well-formed, mild transition into the maximum friction at low SRR. In contrast, all other models show a steep increase and a harder transition into the maximum range. The Eyring-Hyperbolic model provides good results, but the maximum friction is again obviously limited by the limiting shear stress at high pressures and low temperatures (Figure 10, top right), which is not visible in the measurements.

The Carreau model, on the other hand, gives smooth results, but at low temperatures and low SRR, the friction is too high.

In these experiments with their high pressures, it can be noticed very well that the progressive progression of the Tait–Doolittle equation (representing the glass transition) leads to a very steep friction increase in the Shifted-Carreau model at low SRR. Again, the shear-thinning is not effective enough, resulting in a traction curve that runs strongly into the limiting shear stress and represents altogether too high friction.

3.3. Small Two-Disc Machine

The very small contact of the small two-disc machine produces very high friction values (Figure 11). The generated power is low and, thus, only a slight decrease in friction due to shear heating is visible.

The Eyring model is the only one that reacts correctly to changes in load. At a low load, the measured low friction is reproduced, and at a high load, the very high maximum friction is obtained. However, it seems to provide slightly too low friction in this system at very low SRR.

The Eyring-Hyperbolic model reacts too little to pressure variation. At low loads, the friction is too high, and at high loads, the traction curve is too limited. Increasing the slope of the limiting shear stress β could increase the spread, but this would lead to higher friction in the traction calculations of the other test rigs. While the Carreau model calculates too much friction for the large two-disc machine for higher pressures (Figure 10), it is clearly too low for the small two-disc machine (Figure 11). These results are, therefore, contradictory and a better approximation of both systems is not possible. The Shifted-Carreau model provides overall too high friction, but at the highest loaded contact, the curve is again limited and does not reach the measured values. The slope at low SRR is again very steep, and the transition to the maximum friction is too hard.

3.4. EHL2 Tribometer

Finally, the relatively lightly loaded contact of the EHL2 tribometer is considered. Figure 12 shows the results for three temperatures and two pressures. Overall, all shear-thinning models reproduce the contact behaviour similarly well. However, at high pressure and low temperatures (Figure 12, top right), the traction curve is again best reproduced by the Eyring model.

3.5. Comparison with Molecular Dynamics Simulations

In the previous section, the models are iteratively optimised in such a way that the traction behaviour on the tribometers is reproduced as best as possible. However, the modelling quality of the flow behaviour is only evaluated on the basis of the integral response. In order to evaluate how well the individual models represent the theoretical shear behaviour of squalane for very high shear rates, non-equilibrium molecular dynamics (NEMD) simulations were carried out at three temperatures and three pressures. The calculations were carried out with the programme LAMMPS [51].

For the structure of the squalane molecule C₃₀H₆₂, an automated force field topology builder (ATB) was used [52]. In order to reduce the computational effort, the squalane molecule is considered as a united-atom model, in which the structure is built from CH, CH₂, and CH₃ groups and the interactions are reduced to carbon atoms. Although systematic investigations have shown that united-atom models can lead to an underestimation of viscosity [34], valid viscosity calculations can be carried out for molecules containing multiple short methyl branches as it is the case for squalane [28]. According to former studies using NEMD simulations of lubricants [23,25,27], harmonic bonds, angles, dihedrals, and improper potentials are used, with index zero representing the equilibrium position (Equations (10)–(13)). The usual 1/2 term is included in the respective parameter K. The molecular structure used is shown in Figure 13 on the left.

$$E_{bond}=K_{bond}{\left(r-r_0\right)}^2$$

(10)

$$E_{angle}=K_{angle}{\left(\theta -{\theta }_0\right)}^2$$

(11)

$$E_{dihedral}=K_{dihedral}{\left(1+d·\cos \left(n\varphi \right)\right)}^2$$

(12)

$$E_{improper}=K_{improper}{\left(\chi -{\chi }_0\right)}^2$$

(13)

For a good balance between computational time and result quality, 216 squalane molecules are arranged in a cubic simulation cell with periodic boundary conditions. The Lennard-Jones potential Equation (14) is used for the interactions of the chemically unbound groups. The cutoff length is r_c = 14 Å. The neighbour list is updated every 10 iterations.

$$E=4·\epsilon ·\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right] for r<r_c$$

(14)

The GROMOS 54A7 force field [53] is used for the coefficients in Equations (10)–(14) given in Appendix A. The time-step size is defined as 1 fs for all calculations.

The simulation is carried out as follows. First of all, the molecules, which are arranged regularly at the beginning, are mixed statistically. This is performed as an NVE ensemble at a large volume (V), where a random initial velocity is assigned to the molecules for the respective temperature. This is followed by a compression step in which the volume of the simulation cell is varied using a thermo (T)- and barostatting (P) (NPT) algorithm until the desired pressure is reached. This is followed by a homogenisation step to achieve a random arrangement of the molecules. In this step, the simulation cell is already sheared planar with the target shear rate $\dot{\gamma }$. Since the pressure in the cell decreases again due to the redistribution and structural change of the molecules, this step is also carried out as NPT with iteratively varying volume.

For the actual viscosity simulation, the NVT-SLLOD algorithm is used, keeping the volume of the simulation cell constant. Nosé–Hover thermostatting [29] is used here, with the velocity component of the deforming box being subtracted and integrated into the SLLOD equations of motion. The viscosity is derived from the shear stress tensor according to Equation (15):

$$\eta =-\frac{〈{\sigma }_{xy}〉}{\dot{\gamma }}$$

(15)

The simulation is stopped when the viscosity reaches a statistically significant value (Figure 13, right). Due to the constant volume during the calculation (NVT), the pressure may continue to change slightly. If the pressure deviates too much from the target pressure after the end of the simulation, the calculation is repeated.

Figure 14 and Figure 15 show the MD simulation results in comparison with the curves of the four shear-thinning models. Three temperatures (ϑ = 0 °C, 40 °C, and 100 °C) are examined here, which are shown from top to bottom. For each temperature, the curves for pressures p = 0.1 MPa, 500 MPa, and 1000 MPa are determined at four points with the shear rates $\dot{\gamma }$ = 1 × 10⁸ … 1 × 10¹¹ s⁻¹. In order to validate the molecular dynamics model, the data from Bair et al. [23] and Liu et al. [24] for (40 °C; 0.1 MPa) and (100 °C; 316 MPa) are also presented. The comparison shows that the data from the literature can be reproduced well with the present model. Slightly lower values are obtained for very high shear rates. Furthermore, additional test simulations show exact agreements with the results of Jadhao et al. [25].

When considering the results, it becomes clear that the restriction of the simulations to high shear rates $\dot{\gamma }$ ≥ 1 × 10⁸ s⁻¹ prohibits statements about the functional transition from the Newtonian region to the shear-thinning region, as this already occurs at lower shear rates. Only at high temperatures and low pressures the Newtonian behaviour is reproducible. However, the marking of the maximum shear rates occurring in the TEHD simulations shows that the NEMD simulations are directly adjacent to the relevant range and that an evaluation of the models is, therefore, still meaningful.

It should be pointed out once again that the shear-thinning models (Eyring, Eyring-Hyperbolic, etc.) have not been adjusted for the MD simulations, but they continue to correspond to the curves that have been optimised with the help of the traction measurements (parameters from Table 1). Despite its simplicity, the Eyring model with the pressure- and temperature-dependent critical shear stress τ_c is able to meet all considered simulation results very well. Although the shear stress is not limited, the slope of the viscosity is well reproduced with the sinh function. Furthermore, the stronger viscosity decrease at a higher pressure (converging curves, clearly visible at ϑ = 100 °C) is implicitly included in the model. The transition curve’s progression from Newtonian behaviour to shear-thinning at low pressure corresponds to our own calculations and the data from the literature. The critical shear rate is correctly reproduced here.

The Eyring-Hyperbolic model (Figure 14, right) describes a slightly harder transition into the shear-thinning region compared to the Eyring model. Furthermore, it is evident that the limiting shear stress causes a steeper slope of the viscosity decrease than indicated by the MD simulations. This finding is consistent with the observation of an inadequate limitation of the traction curves at a high pressure (see Figure 11, p = 2400 MPa).

While the Carreau model (Figure 15, left) also delivers very good results for high temperatures, for low temperatures, larger values are reproduced than those determined with the simulations. This behaviour is also consistent with the findings from the traction calculations, as the friction at high pressures is overestimated, especially in the low-temperature tests (Figure 10). It could, therefore, be useful to describe the parameter n or τ_c as a function of temperature.

In addition, Figure 15 shows the curves of the Shifted-Carreau model. The results confirm the findings from [23] that this model is able to reproduce the molecular dynamics simulations well for comparatively low pressures (40 °C, 0.1 MPa) and (100 °C, 316 MPa). At higher pressures, however, the constant high value of n = 0.463 causes a strong parallel shift of the curves, and, thus, the MD simulation results carried out here are far overestimated. This is consistent with the findings in [25].

Figure 15 also shows the curves of the limiting shear stress according to Equation (9), which was used in the calculation according to the references [19,39]. It is obvious that the MD simulation results can only be reproduced to some extent with this measure for higher pressures. Since this limiting already takes effect at relatively small shear rates $\dot{\gamma }$ << 1 × 10⁸, it is absolutely necessary for the TEHD calculations, thus leading to the described strongly limited traction curves (see Figure 10). Furthermore, it can be seen that a simple, linear limiting shear stress modelling according to Equation (9) is not appropriate, since the shear stress for a low pressure p = 0.1 MPa is limited to a very low level (blue dashed line). According to the findings in the literature, a bilinear [20,22,40] or exponential-linear approach (Equation (6)) is more useful.

Overall, it can be concluded that the models without limiting shear stress represent the simulated flow behaviour at high shear better. This becomes particularly visible when the shear stress τ is plotted against the shear rate. Figure 16 shows this for an arbitrary calculation point. A limiting shear stress is not evident in the MD simulations even at extreme shear rates (far above the relevant shear rates in the TEHD contact simulations). While the Eyring-Hyperbolic model is already strongly constrained by the limiting shear stress, the shear stress of the Eyring and Carreau models continues to increase mathematically. However, even for these models, the shear stress in the practically relevant shear rate range will nevertheless not increase beyond all limits and will remain quasi-limited by the strongly degressive progression.

Although the traction curves and the NEMD simulations carried out can be reproduced very well even without a limiting shear stress, it must be pointed out that the statements only apply to the SLLOD simulations. Limiting shear stress phenomena can also be understood as wall slip on solids, which are not represented in the SLLOD simulations. Similarly, shear bands, as observed in the EHD contact, cannot be captured as a linear flow profile is applied. For both phenomena, complex moving-wall simulations with correct representation of the solid–fluid interactions are required [31]. Another limitation of the NEMD simulations carried out concerns the applicability of the thermostat algorithm. Although it has been shown that the established methods have the capability to maintain isothermal conditions [28], Bair [54] questions this for extreme high shear rates.

4. Conclusions

For the prediction of friction in concentrated contacts by TEHD simulations, the correct modelling of the lubricant behaviour is essential. Both the choice of the functional relationship and the implementation of a limiting shear stress are the subject of long-standing controversies.

In this paper, a comprehensive investigation of the flow behaviour of the model fluid squalane was carried out. For this purpose, the shear-thinning models Eyring, Eyring-Hyperbolic, and Carreau were used. For comparison, the Shifted–Carreau model for squalane presented in the literature was applied, which is based on the Tait–Doolittle equation for the pressure viscosity. Traction tests on four different tribometers and extensive NEMD simulations were used as the validation methods. The following statements can be made:

• For the modelling of low-shear viscosity, rheometric measurements from the literature were first collected, compared, and fitted with a degressive Rodermund model.
• In the measurement range, this modelling provides very similar values to the Tait–Doolittle equation used in the literature; however, in the extrapolated high-pressure range, the values are significantly lower due to the degressive curve.
• This modelling always leads to better agreements of the TEHD simulations with the traction curves than the degressive–progressive models.

Numerical TEHD contact calculations were carried out iteratively for each tribometer, varying the shear-thinning model parameters, so that the optimal parameters for each model were found as a result. The following statements can be made:

• The general behaviour of the traction measurements is reproduced by all models.
• The models with a limiting shear stress seem to limit the traction curve too much. Especially at high pressure, the traction behaviour is dominated by the limiting shear stress.
• The Carreau model only provides good agreement for certain cases. It can be assumed that this could be improved by variable parameters n and τ_c. Note: In subsequent calculations, no functional correlations for both parameters in terms of pressure and temperature could be found to improve the agreements.
• Overall, the best agreements are achieved with the Eyring model for the investigations carried out. This corresponds to the statements of Spikes and Zhang [10,12,14]. For good agreements, a pressure- and temperature-dependent critical shear stress τ_c = f(ϑ, p) is important.

Since only the integral system response could be evaluated using the traction curves, NEMD simulations were carried out to further validate the parameters found. In this case, the model parameters were not adjusted again, but the iteratively found curves were compared. The following statements can be made:

• It is possible to reproduce the NEMD simulation results from various sources with the described MD-model.
• The Eyring model again provides the best agreement when compared to the viscosity simulations for all temperatures and pressures.
• The Carreau model shows good agreement for high temperatures, but the modelled viscosity is too high for low temperatures. Variable parameters could help here.
• Analogous to the traction tests, a limiting shear stress acts too strongly. There is no evidence for a limiting shear stress from the SLLOD simulations carried out.

Despite the concerns, which are based on rheometric investigations, about the validity of the Eyring model, it can be stated that its application to the TEHD simulations leads to very good results for traction calculation. We cannot confirm the assessment that this is only the case with missing representation of the glass transition of the viscosity–pressure dependence. The use of progressive low-shear viscosity models rather leads to an unrealistic very steep increase in friction at low SRR, which cannot be compensated using other established shear-thinning models.

Of course, these statements only apply to the squalane examined here. However, in subsequent investigations with various practical oils, the use of the Eyring model in combination with degressive low-shear viscosity has also led to promising agreements between TEHD simulations and traction experiments. Due to the good agreement of the findings from the traction tests and MD simulations, it can also be assumed that it is possible to derive valid shear-thinning models exclusively through a combination of traction tests. This is particularly important with regard to TEHD simulations of practical lubricants, as it is not yet possible to carry out MD simulations here.