*1.2. Previous Works*

Several methods have been proposed to solve isothermal, non-Newtonian EHL problems. Kauzlarich and Greenwood [2] obtained the exact solution of line contact EHL problems considering the Herschel–Bulkley model. Conry et al. [3] obtained an exact solution of line contact EHL problems by considering Eyring's model. Specifically, the velocity *u*(*z*) was represented by a function containing cosh, indicating that, in general, *u*(*z*) and *v*(*z*) cannot be precisely represented using polynomials of *z*. Dong and Qian [4] obtained an approximate solution of line contact EHL problems considering Bauer's model and using the weighted residual method. Peiran and Shizhu [5] proposed a method to obtain an exact solution of point contact EHL problems for general rheological models by using an equally divided Z-direction mesh. Subsequently, the authors applied the method to a line contact EHL problem. Kim et al. [6] attempted to obtain an exact solution of point contact EHL problems by considering Eyring's model. Ehret et al. [7] considered the X-direction to align with the sliding direction and obtained an approximate solution of *u*(*z*) and *v*(*z*) represented by the 2nd order polynomial of *z* by using the perturbation method. Thus, researchers have obtained two effective viscosities: one along the sliding direction and another along the perpendicular direction.

Greenwood [8] focused on the considerable amount of computation time required to obtain an exact solution and compared two approximation methods. Sharif et al. [9] considered the X-direction to align with the sliding direction and developed a method to obtain an exact solution of point contact EHL problems for an arbitrary rheological model. Using this approach, the authors obtained two effective viscosities along the X- and Y-directions. Liu et al. [10] formulated a method to obtain an exact solution of point contact thermal EHL problems considering Eyring's model. Yang et al. [11] formulated a general Reynolds equation for point contact EHL problems by dividing the flow into Couette and Poiseuille flows. Subsequently, the authors obtained an exact solution of line contact EHL problems considering the power law model and demonstrated the effectiveness of their proposed method. Furthermore, Bordenet et al. [12] obtained an exact solution of pure rolling point contact EHL problems considering Bauer's model for *n* = 1/2 and applied the approach to grease.

#### **2. Overview of the Proposed Method**

In this work, an isothermal, non-Newtonian EHL formulation considering Method 2-B was developed. Although this approach does not yield the exact solution of *u*(*z*) and *v*(*z*), the calculation is simple and fast. As shown in Figure 2, the local Xc-direction is considered as the direction of the pressure gradient *d*, and the Yc-direction is considered to be perpendicular to Xc. The flow velocities in the Xc- and Yc-directions are denoted as *u*c and *v*c, respectively. As the pressure gradient toward Yc is zero, the flow *v*c is assumed to be a Couette flow, which can be represented using a linear function of *z*. As the pressure gradient toward Xc is generally non-zero, the flow *u*c is assumed to be a Poiseuille flow, which can be represented using a 4th-order polynomial of *z*. As *u*c and *v*c cannot be precisely represented by polynomials of *z*, they are expanded using polynomials of *z*. In such cases, the 6th-order or even higher order polynomials can be considered; however, in this work, a lower 4th-order polynomial was employed. A viscosity corresponding to the Newtonian flow was obtained and termed as the equivalent viscosity. To replace the viscosity of the Newtonian flow with the equivalent viscosity, which is a typical process when evaluating EHL problems, the method proposed by Venner and Lubrecht [13] can be used without any change for the isothermal, non-Newtonian EHL calculation. Given a rheological model, any non-Newtonian isothermal point contact EHL problem can be solved using the proposed method.

**Figure 2.** Local coordinate system Oc, Xc, Yc.

The calculation procedure consists of two steps. First, the Newtonian EHL calculation is executed for the base oil that yields the Newtonian pressure distribution. Second, using the pressure distribution as the initial value, the non-Newtonian EHL calculation is executed for the non-Newtonian flow. The local coordinate system, Oc, Xc, and Yc, at a particular point, depends on the pressure gradient and it is determined simultaneously in the process to obtain the pressure distribution. The equivalent viscosity is calculated based on the local coordinate system in each iteration loop to obtain the pressure distribution. Here, the method was applied to an experimental grease characterized by Bauer's model.

#### **3. Calculation of Velocity Distribution as a Function of** *z*

Figure 1 presents the XYZ coordinate system. The force balance of a fluid can be expressed as follows [10]:

$$\frac{\partial \pi\_{\chi}}{\partial z} = \frac{\partial P}{\partial x} \tag{1}$$

$$\frac{\partial \pi\_y}{\partial z} = \frac{\partial P}{\partial y} \tag{2}$$

where *P* is the pressure of the fluid, and *τx* and *τy* are the shear stresses in the X- and Y-directions, respectively. The parameters *τ* and . *γ* are defined as follows:

$$
\tau = \sqrt{\tau\_\mathbf{x}^2 + \tau\_\mathbf{y}^2} \tag{3}
$$

$$
\dot{\gamma} = \sqrt{\frac{\partial \mu^2}{\partial z} + \frac{\partial v^2}{\partial z}}\tag{4}
$$

In Bauer's model, *τ* is assumed to be represented as a function of .*γ*, as follows [1,14,15]:

$$
\tau = \left(\tau\_0 + k\_1 \cdot \dot{\gamma} + k\_2 \cdot \dot{\gamma}^T\right) \cdot \frac{\eta\_n}{\eta\_0} \tag{5}
$$

Here, *τ*0, *k*1, *k*2, and *n* are Bauer's rheological parameters, and *η*n and *η*0 represent the *P* dependent viscosity and ambient viscosity of the base oil, respectively. In this work, according to Dong and Qian [4], the parameters *τ*0, *k*1, *k*2, and *n* were assumed to be *P*-independent known values determined from the *τ*- .*γ* curve measured at the ambient pressure. In Eyring's model, according to Conry et al. [3] and Johnson and Tevaarwerk [16], the relationship between *τ* and .*γ* can be expressed as follows:

$$\dot{\gamma} = \frac{\pi\_0'}{\eta\_{\rm n}} \sin \mathbf{h} \left(\frac{\tau}{\tau\_0 \prime}\right)$$

Here, *<sup>τ</sup>*0 is Eyring's rheological parameter. Therefore, the relationship between *τ* and . *γ* can be rewritten as follows:

$$\tau = \tau\_0' \cdot \sin \mathbf{h}^{-1} \left(\frac{\eta\_\mathbf{n} \dot{\gamma}}{\tau\_0'}\right)$$

The effective viscosity *η*∗ can be defined as follows:

$$
\eta^\* = \frac{\pi}{\dot{\gamma}} \tag{6}
$$

The effective viscosity *η*∗ is a function of .*γ*, which in turn is a function of *z*. The shear stresses *τx* and *τy* are assumed to be represented as

$$
\pi\_{\mathfrak{X}} = \eta^\* \, \frac{\partial \mathfrak{u}}{\partial z} \tag{7}
$$

$$
\pi\_{\mathcal{Y}} = \eta^\* \frac{\partial v}{\partial z} \tag{8}
$$

Substituting Equations (7) and (8) into Equations (1) and (2), respectively, yields

$$\frac{\partial}{\partial z} \left( \eta^\* \frac{\partial \mu}{\partial z} \right) = \frac{\partial P}{\partial x} \tag{9}$$

$$\frac{\partial}{\partial z} \left( \eta^\* \frac{\partial \upsilon}{\partial z} \right) = \frac{\partial P}{\partial y} \tag{10}$$

The pressure gradient vector *d* is defined as follows:

$$d = \left(\frac{\partial P}{\partial \mathbf{x}'}, \frac{\partial P}{\partial y}\right) \tag{11}$$

Similar to the method employed by Yang et al. [11], this method involves the flow being divided into Couette and Poiseuille flows. As shown in Figure 2, the Xc-direction is considered to be along *d*, and its direction vector is *n*1. The Yc-direction is perpendicular to *d*, and its direction vector is *n*2. The local coordinate system, Oc, Xc, and Yc, at a given point depends on the pressure gradient and is determined simultaneously in the process to obtain the pressure distribution. In the Xc and Yc coordinate system, Equations (9) and (10) can be rewritten as follows: 

$$\frac{\partial}{\partial z} \left( \eta^\* \, \frac{\partial \mu\_c}{\partial z} \right) = d \tag{12}$$

$$\frac{\partial}{\partial z} \left( \eta^\* \frac{\partial v\_c}{\partial z} \right) = 0 \tag{13}$$

$$d = \sqrt{d \cdot d} \tag{14}$$

Furthermore, the velocities *u*c(*z*) and *v*c(*z*) satisfy the following boundary conditions:

$$u\_{\mathfrak{c}}(0) = \mathcal{U}^-,\ u\_{\mathfrak{c}}(h) = \mathcal{U}^+,\ v\_{\mathfrak{c}}(0) = \mathcal{V}^-,\ v\_{\mathfrak{c}}(h) = \mathcal{V}^+ \tag{15}$$

Here, *U*<sup>+</sup> and *U*− denote the velocities of the upper and lower surfaces in the Xcdirection, respectively; *V*<sup>+</sup> and *V*− denote the velocities of the upper and lower surfaces in the Yc-direction, respectively. Although the velocities *u*c(*z*) and *v*c(*z*) cannot be represented by polynomials exactly [3], here they are approximated and expanded using polynomials of *z* so that the velocities satisfy Equation (15), as follows. In this case, the variables *a*1 and *a*2 are unknown.

$$u\_{\mathbf{c}}(z) = \frac{\Delta l I}{h} \cdot z + lI^- - a\_1 \cdot z(h - z) - a\_2 \cdot z^2(h - z)^2 \tag{16}$$

$$v\_{\mathbb{C}}(z) = \frac{\Delta V}{h} \cdot z + V^- \tag{17}$$

Here, *h* is the fluid film thickness, and Δ*U* and Δ*V* denote the velocity differences; specifically, Δ*U= U*<sup>+</sup> − *U*− and Δ*V= V*<sup>+</sup> − *V*<sup>−</sup>. A higher order term of *z*, for example, *z*<sup>3</sup>(*<sup>h</sup>* − *z*)3, can also be considered; however, in this work, the lower-order approximation was chosen. As the pressure gradient toward the Yc direction is zero, *v*c was assumed to be a Couette flow and approximated considering a linear equation of *z*. Furthermore, as the pressure gradient toward the Xc-direction is generally non-zero, *u*c was assumed to be a Poiseuille flow and approximated using a 4th-order polynomial of *z*. If *d* = 0, then *u*c is also a Couette flow, and Equation (16) can be replaced with the following equation:

$$
\mu\_{\mathbf{c}}(z) = \frac{\Delta \mathcal{U}}{h} \cdot z + \mathcal{U}^- \tag{18}
$$

The following equations are derived from Equations (16) and (17):

$$\frac{\partial u\_{\mathbb{C}}}{\partial z} = u'(z) = \frac{\Delta U}{h} - a\_1 \cdot (h - 2z) - a\_2 \cdot 2z(h - z)(h - 2z) \tag{19}$$

$$\frac{\partial v\_{\text{c}}}{\partial z} = v'(z) = \frac{\Delta V}{h} \tag{20}$$

The following equations are derived from Equation (19):

$$u'(0) = \frac{\Delta U}{h} - a\_1 h \tag{21}$$

$$
\mu'(h) = \frac{\Delta U}{h} + a\_1 h \tag{22}
$$

$$
\mu' \left( \frac{h}{4} \right) = \frac{\Delta l I}{h} - a\_1 \frac{h}{2} - a\_2 \frac{3h^3}{16} \tag{23}
$$

$$
\mu' \left(\frac{3h}{4}\right) = \frac{\Delta l I}{h} + a\_1 \frac{h}{2} + a\_2 \frac{3h^3}{16} \tag{24}
$$

If rheological constitutive equations are given, the effective viscosity *η*<sup>∗</sup>(*z*) can be calculated using Equations (4)–(6), (16) and (17). The integration of Equation (12) from *z* = 0 to *z* = *h* yields Equation (25), and the integration of Equation (12) from *z* = *h*/4 to *z* = 3*h*/4 yields Equation (26). These equations are used to determine the values of *a*1 and *a*2.

$$\eta^\*(h) \cdot \left(\frac{\Delta lI}{h} + a\_1 h\right) - \eta^\*(0) \cdot \left(\frac{\Delta lI}{h} - a\_1 h\right) = dh \tag{25}$$

$$
\eta^\* \left( \frac{3h}{4} \right) \cdot \left( \frac{\Delta l I}{h} + a\_1 \frac{h}{2} + a\_2 \frac{3h^3}{16} \right) - \eta^\* \left( \frac{h}{4} \right) \cdot \left( \frac{\Delta l I}{h} - a\_1 \frac{h}{2} - a\_2 \frac{3h^3}{16} \right) = \frac{dh}{2} \tag{26}
$$

Equations (21) and (22) show that both *η*<sup>∗</sup>(*h*) and *η*<sup>∗</sup>(0) do not contain *a*2. Therefore, Equation (25) does not contain *a*2 and contains only the unknown variable *a*1. The equation can be solved using the one-variable Newton–Raphson method. Although Equation (26) contains both *a*1 and *a*2, *a*1 has been determined using Equation (25). Consequently, Equation (26) can be considered as an equation involving only the unknown variable *a*2. Thus, it can also be solved using the one-variable Newton–Raphson method. To determine *a*1, a non-dimensional variable *b*1 defined using Equation (27) and a function *f*1(*b*1) defined using Equation (28) are introduced. The value of *b*1 can be calculated considering *f*1(*b*1) = 0.

$$b\_1 = \log\left(\frac{2a\_1\eta\_n}{d}\right) \tag{27}$$

$$f\_1(b\_1) = \eta^\*(h) \cdot \left(\frac{\Delta l I}{h} + \mathbf{e}^{b\_1} \frac{dh}{2\eta\_{\rm in}}\right) - \eta^\*(0) \cdot \left(\frac{\Delta l I}{h} - \mathbf{e}^{b\_1} \frac{dh}{2\eta\_{\rm in}}\right) - d \cdot h \tag{28}$$

When *b*1 is near the solution, Δ*b*1 can be calculated using the following equation:

$$0 = f\_1(b\_1 + \Delta b\_1) \doteq f\_1(b\_1) + \frac{\text{d}f\_1}{\text{d}b\_1} \cdot \Delta b\_1 \doteq f\_1(b\_1) + \text{e}^{b\_1} \frac{\text{d}h}{2\eta\_{\text{n}}} \left[\eta^\*(h) + \eta^\*(0)\right] \cdot \Delta b\_1 \tag{29}$$

In other words, the new candidate *b*1new of *b*1 is calculated using the iterative process of Newton–Raphson's method, as follows:

$$b\_{1\text{new}} = b\_1 + \Delta b\_1 = b\_1 - \frac{f\_1(b\_1)}{a\_1 h \left[\eta^\*(h) + \eta^\*(0)\right]} \tag{30}$$

As *η*<sup>∗</sup>(*h*) and *η*<sup>∗</sup>(0) are originally functions of *b*1, d*f*1/d*b*1 includes d*η*<sup>∗</sup>/d*b*1; however, in this work, the dependency was ignored, and Δ*b*1 was approximated as in Equation (30). To determine *a*2, a non-dimensional variable *b*2 defined using Equation (31) and a function *f*2(*b*2) defined using Equation (32) are introduced. The value of *b*2 can be calculated considering *f*2(*b*2) = 0.

$$b\_2 = \frac{2a\_2h^2\eta\_{\rm in}}{5d} \tag{31}$$

$$f\_2(b\_2) = \eta^\* \left(\frac{3h}{4}\right) \cdot \left(\frac{\Delta lI}{h} + a\_1 \frac{h}{2} + a\_2 \frac{3h^3}{16}\right) - \eta^\* \left(\frac{h}{4}\right) \cdot \left(\frac{\Delta lI}{h} - a\_1 \frac{h}{2} - a\_2 \frac{3h^3}{16}\right) - \frac{dh}{2} \tag{32}$$

When *b*2 is near the solution, Δ*b*2 can be calculated using the following equation:

$$0 = f\_2(b\_2 + \Delta b\_2) \doteq f\_2(b\_2) + \frac{\text{d}f\_2}{\text{d}b\_2} \cdot \Delta b\_2 \doteq f\_2(b\_2) + \frac{15 \text{d}h}{32 \eta\_n} \left[ \eta^\* \left(\frac{3h}{4}\right) + \eta^\* \left(\frac{h}{4}\right) \right] \cdot \Delta b\_2 \tag{33}$$

In other words, the new candidate *b*2new of *b*2 is calculated using the iterative process of Newton–Raphson's method, as follows:

$$b\_{2\text{new}} = b\_2 + \Delta b\_2 = b\_2 - \frac{f\_2(b\_2)}{\left[\eta^\*(3h/4) + \eta^\*(h/4)\right]} \cdot \frac{32\eta\_{\text{n}}}{15dh} \tag{34}$$

As *η*<sup>∗</sup>(3*h*/4) and *η*<sup>∗</sup>(*h*/4) are originally functions of *b*2, d*f*2/d*b*2 includes d*η*<sup>∗</sup>/d*b*2; however, in this work, the dependency was ignored, and Δ*b*2 was approximated as in Equation (34). Subsequently, in the iteration process of Newton–Raphson's method, only *η*∗ depends on the rheological characteristics. Therefore, if the rheological equation corresponding to Equation (5) is incorporated, any isothermal, non-Newtonian EHL calculation can be performed. As per the Newton–Raphson method, the initial value for both *b*1 and *b*2 can be zero. As both variables *a*1 and *a*2 are solved using the one-variable Newton–Raphson method, the calculation can be performed within a reasonable time.

#### **4. Calculation of Equivalent Viscosity, Flow, and Surface Force**

The flow *q*1 along the *n*1 direction can be defined using Equation (16), as follows. The density *ρ* is assumed to be independent of *z*.

$$\begin{split} q\_1 = \int\_0^h \rho u\_\varsigma \, dz &= \int\_0^h \rho \left[ \frac{\Delta l I}{h} \cdot z + l I^- - a\_1 \cdot z (h - z) - a\_2 \cdot z^2 (h - z)^2 \right] \, dz \\ &= \rho l I h - \frac{\rho a\_1 h^3}{6} - \frac{\rho a\_2 h^5}{30} \end{split} \tag{35}$$

Here, *U* is the average velocity in the Xc-direction and can be expressed as follows:

$$
\mathcal{U} = \frac{\mathcal{U}^+ + \mathcal{U}^-}{2} \tag{36}
$$

The equivalent viscosity *ηeq* is defined as follows:

$$\eta\_{eq} = \frac{5d}{2(5a\_1 + a\_2h^2)}\tag{37}$$

Consequently, *q*1 can be represented as

$$q\_1 = \rho Ulh - \frac{\rho h^3}{12\eta\_{c\eta}} \cdot d \tag{38}$$

The flow *q*2 along the *n*2 direction can be derived from Equation (17), as follows:

$$q\_2 = \int\_0^h \rho v\_\mathbb{C} \, dz = \rho V h \tag{39}$$

Here, *V* is the average velocity in the Yc direction and can be expressed as follows:

$$V = \frac{V^{+} + V^{-}}{2} \tag{40}$$

Hence, in the XYZ coordinate system, the flow vector *q* can be expressed as

$$\begin{split} q &= \left(\rho \text{III} \mathbf{\hat{n}} - \frac{\rho \mathbf{h}^3}{12\eta\_{\text{eq}}} \cdot d \right) \cdot \mathbf{n}\_1 + \rho V \mathbf{h} \cdot \mathbf{n}\_2 \\ &= \rho \text{III} \mathbf{\hat{n}} \cdot \mathbf{n}\_1 + \rho V \mathbf{h} \cdot \mathbf{n}\_2 - \frac{\rho \mathbf{h}^3}{12\eta\_{\text{eq}}} \cdot d \cdot \mathbf{n}\_1 \\ &= \rho \text{III} \mathbf{\hat{n}} - \frac{\rho \mathbf{h}^3}{12\eta\_{\text{eq}}} \cdot \left(\frac{\partial P}{\partial \mathbf{x}}, \frac{\partial P}{\partial \mathbf{y}}\right) \end{split} \tag{41}$$

Here, *U* is the average velocity vector of the upper and lower surfaces, defined as follows:

$$\mathcal{U} = \mathcal{U} \cdot \mathfrak{n}\_1 + \mathcal{V} \cdot \mathfrak{n}\_2 = \begin{pmatrix} \mathcal{U}\_{\mathcal{X}} \ \mathcal{U}\_{\mathcal{Y}} \end{pmatrix} \tag{42}$$

*Ux* and *Uy* denote the average velocities in the upper and lower surfaces in the XYdirection, respectively. When the mass conservation law is applied to Equation (41), the following modified Reynolds equation is obtained.

$$\frac{\partial \rho L\_x h}{\partial x} + \frac{\partial \rho L\_y h}{\partial y} - \frac{\partial}{\partial x} \left( \frac{\rho h^3}{12 \eta\_{cq}} \cdot \frac{\partial P}{\partial x} \right) - \frac{\partial}{\partial y} \left( \frac{\rho h^3}{12 \eta\_{cq}} \cdot \frac{\partial P}{\partial y} \right) = 0 \tag{43}$$

The difference in the representation of Equations (41) and (43) and that of Newtonian flow only pertains to the viscosities *ηeq* and *η*n, respectively. In fact, the equivalent viscosity *ηeq* defined by Equation (37) was determined so that Equations (41) and (43) maintain the same form as that of Newtonian flow. Therefore, the EHL calculation procedure for Newtonian flows, such as the method proposed by Venner and Lubrecht [13], can be

applied to the current calculation by simply replacing *η*n with *ηeq*. The shear stress *τ*1 along the Xc-direction can be expressed as

$$\tau\_1(z) = \eta^\* \frac{\partial u\_\mathbb{C}}{\partial z} = \eta^\* \left[ \frac{\Delta \mathcal{U}}{h} - a\_1 \cdot (h - 2z) - a\_2 \cdot 2z(h - z)(h - 2z) \right] \tag{44}$$

Therefore, the surface forces *P*10 and *P*1h acting on the lower and upper surfaces along the Xc-direction, respectively, can be defined as

$$P\_{10} = \tau\_1(0) = \eta^\*(0) \cdot \left(\frac{\Delta U}{h} - a\_1 h\right) \tag{45}$$

$$P\_{\rm Ih} = -\tau\_{\rm I}(h) = -\eta^\*(h) \cdot \left(\frac{\Delta U}{h} + a\_1 h\right) \tag{46}$$

The shear stress *τ*2 along the Yc-direction is expressed as

$$
\pi\_2 = \eta^\* \frac{\partial \upsilon\_\mathbb{C}}{\partial z} = \eta^\* \frac{\Delta V}{h} \tag{47}
$$

Therefore, the surface forces *P*20 and *P*2h acting on the lower and upper surfaces along the Yc-direction, respectively, can be defined as

$$P\mathbb{2}0 = \tau\mathbb{2}(0) = \eta^\*(0)\cdot\frac{\Delta V}{h} \tag{48}$$

$$P\_{2\text{h}} = -\tau\_2(h) = -\eta^\*(h) \cdot \frac{\Delta V}{h} \tag{49}$$

In the XYZ coordinate system, the surface force vectors *P*0 and *P*h that act on the lower and upper surfaces, respectively, can be expressed as

$$P\_0 = P\_{10} \cdot \mathfrak{n}\_1 + P\_{20} \cdot \mathfrak{n}\_2 \tag{50}$$

$$P\_{\rm h} = P\_{\rm 1h} \cdot \mathfrak{m}\_1 + P\_{\rm 2h} \cdot \mathfrak{m}\_2 \tag{51}$$

#### **5. Application to a Grease**

As mentioned previously, Kochi et al. [1] conducted experiments on grease under soft EHL conditions and measured the film thickness and traction forces. In the present study, the proposed method was applied to one of the greases considered in the study by Kochi et al. [1] so as to validate the theoretical approach. Grease A in the literature [1] was chosen to test the proposed method. The modified Reynolds equation, as expressed in Equation (43), was solved using a multi-level method, as reported by Venner and Lubrecht [13]. The commercial program Tribology Engineering Dynamics Contact Problem Analyzer (TED/CPA) V852 was employed. Figure 3 illustrates the calculation conditions. The upper body was a steel ball, and the lower body was a disk composed of glass or polycarbonate (PC). The rheological properties of the grease were assumed to be represented by Bauer's model, according to existing literature [1]. Detailed properties of the steel, glass, PC, and grease are described in the previous study [1]. The pressure dependency of the density was defined using Dowson–Higginson's formula, as follows:

$$
\rho(P) = \rho\_0 \cdot \frac{p\_0 + \beta \cdot P}{p\_0 + P}, \; \rho\_0 = 1, \; p\_0 = 590 \text{ MPa}, \; \beta = 1.34 \tag{52}
$$

**Figure 3.** Calculation condition.

The pressure dependency of the base oil viscosity was assumed to be defined using Barus' formula:

$$
\eta\_{\rm n}(P) = \eta\_0 \cdot \exp(aP), \text{ } \eta\_0: \text{ 49.5 mPa} \cdot \text{s}, \text{ } \text{ } \text{ 14.} \,\text{GPa}^{-1} \tag{53}
$$

In particular, Equation (6) diverges when *γ* approaches zero. It was assumed that if . *γ* is lower than a certain value *cmin*=100.0 s<sup>−</sup>1, then *η*∗ varies linearly with the gradient d*η*<sup>∗</sup>/d .*γ* at *cmin*. Bauer's parameter *k*1 was assumed to be the base oil ambient viscosity *η*0. The values of Bauer's parameters *τ*0, *k*2, and *n* were determined from the apparent viscosity curve of grease A, as shown in Figure 9 of Kochi et al. [1]. When *P* = 0, the curve was assumed to pass through the following three points: 100 mm/s, 6.46475 Pa·s; 10,000 mm/s, 0.16712 Pa·s; and 1,000,000 mm/s, 0.06573 Pa·s.

.

The values of *τ*0, *k*2, and *n* were determined to ensure that Bauer's curve passes through the abovementioned three points, as follows:

$$
\pi\_0 = 0.000621839 \text{ MPa}, \ k\_2 = 6.99118 \cdot 10^{-\gamma}, n = 0.7248 \tag{54}
$$

Figure 4 presents the central film thickness of grease A and base oil as a function of the rolling velocity in the case of pure rolling and *Fz* = 10 N. The solid lines show the calculation results and the dotted lines show the experimental results. The experimental data were read using the caliper from Figure 4 of Kochi et al. [1]. Figure 4a,b shows the cases of a PC disk and glass disk, respectively. The calculation range was set as follows: −1.0 ≤ X ≤ 0.4 and −0.6 ≤ Y ≤ 0.6.

**Figure 4.** Comparison of film thickness between the calculation and the experiment.

However, in the case of grease A, a glass disk, and a velocity of less than or equal to 300 mm/s, the range was set as follows: −0.35 ≤ X ≤ 0.14 and −0.2 ≤ Y ≤ 0.2.

In these cases, the oil film was thin and the wide range calculation became hard to perform. In the case of the PC disk, the calculation results exhibited good agreemen<sup>t</sup> with the experimental results. In the case of the glass disk, the results of the base oil showed some difference but the other data showed reasonable agreement. Figure 5, which illustrates a sample calculation, shows the distribution of *P*, *h*, and *ηeq* for the case involving a PC disk, a pure rolling velocity of 1200 mm/s, and grease A. Typically, in the case of pure rolling velocity, .*γ* is small and *ηeq* is large. In addition, only the appearance of the distribution of *ηeq* differs from that of the base oil. Figure 5d shows the distribution of *ηeq* at section Y = 0. It can be seen that *ηeq* becomes extremely large at the center, where the pressure gradient is small and the flow volume is low.

**Figure 5.** Distribution of *P*, *h*, and *ηeq* of grease A at a rolling velocity of 1200 mm/s.

Figure 6 presents the traction coefficient as a function of the slide roll ratio when *Fz* = 20 N. The solid lines show the calculation results and the dotted lines show the experimental results. The experimental data were read using the caliper from Figure 8 of Kochi et al. [1]. The slide roll ratio is the difference between the upper and lower velocities divided by their average value. The traction coefficient was calculated according to the approaches proposed in the existing literature [1,17]:

$$TRC = \frac{TX0 - TXh}{2F\_z} \tag{55}$$

Here, *TX*0 and *TXh* denote the X-direction traction forces acting on the lower and upper surfaces, respectively; *Fz* is the load. The calculation range was −0.8 ≤ X ≤ 0.5, −0.5 ≤ Y ≤ 0.5. The calculation results were in fairly good agreemen<sup>t</sup> with the experimental results.

**Figure 6.** Comparison of the traction coefficient between the calculation and the experiment.

Figure 7, which illustrates a sample calculation, shows the distribution of *P*, *h*, and *ηeq* for the case involving grease A and a slide roll ratio of 10%. Only the appearance of the distribution of *ηeq* differs from that of the base oil, exhibiting a figure eight in the vicinity of the contact point. Figure 7c,d shows the same distribution of *ηeq* in different display ranges. It can be observed that *ηeq* reduces in the rapid flow region. The figure eight phenomenon is characteristic of non-Newtonian flow, in which the apparent viscosity becomes large at a low velocity gradient. This phenomenon can be explained as follows. Let the XY coordinates of points A, B, C, and D be (−0.1, <sup>0</sup>), (+0.1, <sup>0</sup>), (0, <sup>−</sup>0.1), and (0, <sup>+</sup>0.1), respectively. The velocity gradient vector *n*1 at these points are approximately (+1, <sup>0</sup>), (−1, <sup>0</sup>), (0, <sup>+</sup><sup>1</sup>), and (0, <sup>−</sup><sup>1</sup>), respectively, as shown in Figure 8a. When *a*2 in Equation (37) is neglected, *ηeq* can be expressed as follows:

$$
\eta\_{c\eta} = \frac{d}{2a\_1} \tag{56}
$$

Calculating *a*1 from Equation (25) and substituting it into Equation (56) yields:

$$\eta\_{eq} = \frac{\left[\eta^\*(h) + \eta^\*(0)\right]}{2} \cdot \frac{1}{1 + \frac{\Delta H}{dh^2} \left[\eta^\*(0) - \eta^\*(h)\right]} \tag{57}$$

At points C and D, the direction of flow and that of *n*1 is orthogonal, so Δ*U* = 0. Consequently, the equivalent viscosity parameters *ηeq*,*<sup>C</sup>* and *ηeq*,*<sup>D</sup>* at points C and D are given as follows:

$$
\eta\_{\text{eq},\text{C}} = \eta\_{\text{eq},\text{D}} = \frac{[\eta^\*(h) + \eta^\*(0)]}{2} \tag{58}
$$

At point A, where *n*1 directs towards +X and the velocity gradient at Z = *h* is greater than that at Z = 0 (as shown in Figure 8b), the following equation is satisfied:

$$
\eta^\*(0) > \eta^\*(h), \ \Delta \!\!/\!I = \!\!\!I\_{\rm x}^+ - \!\!\!I\_{\rm x}^- > 0\tag{59}
$$

At point B, where *n*1 directs towards −X, and the velocity gradient at Z = *h* is smaller than that at Z = 0 (as shown in Figure 8c), the following equation is satisfied:

$$
\eta^\*(0) < \eta^\*(h), \,\Delta l I = -\left(\mathcal{U}\_x^+ - \mathcal{U}\_x^-\right) < 0 \tag{60}
$$

In any case, the equivalent viscosity parameters *ηeq*,*<sup>A</sup>* and *ηeq*,*<sup>B</sup>* at points A and B become smaller than *ηeq*,*<sup>C</sup>* and *ηeq*,*<sup>D</sup>* by the effect of the second term of Equation (57). In the

pure rolling case, where Δ*U* is 0, Equation (57) results in Equation (58). It can be understood that in such a case, no figure-eight-shaped distribution appears as is shown in Figure 5c.

**Figure 7.** Calculation results of grease A at a 10% slide ratio.

**Figure 8.** Explanation of the figure eight phenomenon.
