Formulation to Calculate Isothermal, Non-Newtonian Elastohydrodynamic Lubrication Problems Using a Pressure Gradient Coordinate System and Its Verification by an Experimental Grease

: This paper presents a formulation of point contact elastohydrodynamic lubrication analysis for an isothermal, non-Newtonian flow. A coordinate system of the pressure gradient was employed herein. A Couette flow and a Poiseuille flow were considered along the directions of the zero and non-zero pressure gradients, respectively. The Poiseuille flow velocity was assumed to be represented by a 4th-order polynomial of z along the film thickness direction. The Couette flow velocity was assumed to be represented by a linear function of z . Subsequently, the modified Reynolds equation, which contains an equivalent viscosity, was obtained. Using Bauer’s rheological model, the formulation presented in this study was applied to a grease that has been previously experimented upon. The results of previous studies were compared with those of the present study and a reasonable agreement was noted. The distribution of the equivalent viscosity showed a no-table difference from that of Newtonian flow. The formulation can be incorporated easily to the usual elastohydrodynamic lubrication calculation procedure for Newtonian flow. The method can be easily applied to other non-Newtonian rheological models. The equivalent viscosity can be calculated using the one-parameter Newton-Raphson’s method; as a result, the calculation can be performed within a reasonable time.


Introduction
Performing experiments on non-Newtonian flows, including grease flows, is considerably time-consuming and costly. Therefore, it is important to numerically analyze the phenomena corresponding to non-Newtonian flows. Numerical approaches can help obtain a variety of data that cannot be determined experimentally. Grease flows can be well defined using Bauer's model; however, owing to the extreme complexity of this model, it is difficult to determine the exact solution as well as approximate solutions for point contact, isothermal, non-Newtonian elastohydrodynamic lubrication (EHL) analyses. Therefore, it is convenient that the non-Newtonian EHL calculation can be executed within a reasonable calculation time and without large modification to the usual Newtonian EHL calculation procedure. Kochi et al. [1] performed experiments on grease under soft EHL conditions and measured the film thickness and traction forces. The method proposed in the present study can be applied to the grease considered in Kochi et al. [1] to validate this theoretical approach.

Classification of Calculation Methods
As shown in Figure 1, the Z-direction is considered to be the film thickness direction. The flow velocities along the X-and Y-directions are denoted by and , respectively, which are functions of x, y, and z; however, when considering only z dependency, the velocities can be expressed as (z) and (z), respectively. The numerical methods for isothermal, non-Newtonian EHL analyses can be classified in terms of the accuracy of (z) and (z), as follows: Method 1: Exact solution of (z) and (z) is obtained. Method 2: Approximate solution of (z) and (z) is obtained.
Method 2 can be further classified in terms of the employed coordinate system.

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 (z) was represented by a function containing cosh, indicating that, in general, (z) and (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 (z) and (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 Xand 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 = 1/2 and applied the approach to grease.

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 (z) and (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 , and the Yc-direction is considered to be perpendicular to Xc. The flow velocities in the Xc-and Yc-directions are denoted as c and c , respectively. As the pressure gradient toward Yc is zero, the flow 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 c is assumed to be a Poiseuille flow, which can be represented using a 4th-order polynomial of z. As c and 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.
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.  Figure 1 presents the XYZ coordinate system. The force balance of a fluid can be expressed as follows [10]:

Calculation of Velocity Distribution as a Function of z
where is the pressure of the fluid, and and are the shear stresses in the X-and Y-directions, respectively. The parameters and γ̇ are defined as follows: In Bauer's model, is assumed to be represented as a function of γ̇, as follows [1,14,15]: Here, 0 , 1 , 2 , and are Bauer's rheological parameters, and n and 0 represent the dependent viscosity and ambient viscosity of the base oil, respectively. In this work, according to Dong and Qian [4], the parameters 0 , 1 , 2 , and were assumed to be -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: Here, 0 ′ is Eyring's rheological parameter. Therefore, the relationship between and γ̇ can be rewritten as follows: The effective viscosity * can be defined as follows: * = γ̇ The effective viscosity * is a function of γ̇, which in turn is a function of z. The shear stresses and are assumed to be represented as Substituting Equations (7) and (8) into Equations (1) and (2), respectively, yields The pressure gradient vector is defined as follows: 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 , and its direction vector is 1 . The Yc-direction is perpendicular to , and its direction vector is 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: Furthermore, the velocities c (z) and c (z) satisfy the following boundary conditions: Here, + and − denote the velocities of the upper and lower surfaces in the Xcdirection, respectively; + and − denote the velocities of the upper and lower surfaces in the Yc-direction, respectively. Although the velocities c (z) and 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 1 and 2 are unknown.
Here, ℎ is the fluid film thickness, and Δ and Δ denote the velocity differences; specifically, Δ = + − − and Δ = + − − . A higher order term of , for example, 3 (ℎ − ) 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, 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, c was assumed to be a Poiseuille flow and approximated using a 4th-order polynomial of z. If = 0, then c is also a Couette flow, and Equation (16) can be replaced with the following equation: The following equations are derived from Equations (16) and (17): The following equations are derived from Equation (19): If rheological constitutive equations are given, the effective viscosity * (z) can be calculated using Equations (4)-(6), (16), and (17). The integration of Equation (12) Equations (21) and (22) show that both * (ℎ) and * (0) do not contain 2 . Therefore, Equation (25) does not contain 2 and contains only the unknown variable 1 . The equation can be solved using the one-variable Newton-Raphson method. Although Equation (26) contains both 1 and 2 , 1 has been determined using Equation (25). Consequently, Equation (26) can be considered as an equation involving only the unknown variable 2 . Thus, it can also be solved using the one-variable Newton-Raphson method. To determine 1 , a non-dimensional variable 1 defined using Equation (27) and a function 1 ( 1 ) defined using Equation (28) are introduced. The value of 1 can be calculated considering 1 ( 1 ) = 0.
When 1 is near the solution, Δ 1 can be calculated using the following equation: In other words, the new candidate 1new of 1 is calculated using the iterative process of Newton-Raphson's method, as follows: As * (ℎ) and * (0) are originally functions of 1 , d 1 /d 1 includes d * /d 1 ; however, in this work, the dependency was ignored, and Δ 1 was approximated as in Equation (30). To determine 2 , a non-dimensional variable 2 defined using Equation (31) and a function 2 ( 2 ) defined using Equation (32) are introduced. The value of 2 can be calculated considering 2 ( 2 ) = 0. (31) When 2 is near the solution, Δ 2 can be calculated using the following equation: In other words, the new candidate 2new of 2 is calculated using the iterative process of Newton-Raphson's method, as follows: As * (3ℎ/4) and * (ℎ/4) are originally functions of 2 , d 2 /d 2 includes d * / d 2 ; however, in this work, the dependency was ignored, and Δ 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 1 and 2 can be zero. As both variables 1 and 2 are solved using the one-variable Newton-Raphson method, the calculation can be performed within a reasonable time.

Calculation of Equivalent Viscosity, Flow, and Surface Force
The flow 1 along the 1 direction can be defined using Equation (16), as follows. The density is assumed to be independent of z.
Here, is the average velocity in the Xc-direction and can be expressed as follows: The equivalent viscosity is defined as follows: Consequently, 1 can be represented as The flow 2 along the 2 direction can be derived from Equation (17), as follows: Here, is the average velocity in the Yc direction and can be expressed as follows: Hence, in the XYZ coordinate system, the flow vector can be expressed as Here, is the average velocity vector of the upper and lower surfaces, defined as follows: and denote the average velocities in the upper and lower surfaces in the XY-direction, respectively. When the mass conservation law is applied to Equation (41), the following modified Reynolds equation is obtained.
The difference in the representation of Equations (41) and (43) and that of Newtonian flow only pertains to the viscosities and n , respectively. In fact, the equivalent viscosity 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 . The shear stress 1 along the Xc-direction can be expressed as Therefore, the surface forces 10 and 1h acting on the lower and upper surfaces along the Xc-direction, respectively, can be defined as The shear stress 2 along the Yc-direction is expressed as Therefore, the surface forces 20 and 2h acting on the lower and upper surfaces along the Yc-direction, respectively, can be defined as In the XYZ coordinate system, the surface force vectors 0 and h that act on the lower and upper surfaces, respectively, can be expressed as

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: The pressure dependency of the base oil viscosity was assumed to be defined using Barus' formula: n ( ) = 0 · exp( )， 0 : 49.5 mPa · s， : 14 GPa −1 (53) In particular, Equation (6) diverges when γ̇ approaches zero. It was assumed that if γ̇ is lower than a certain value cmin=100.0 s −1 , then * varies linearly with the gradient d * dγ⁄ at cmin. Bauer's parameter 1 was assumed to be the base oil ambient viscosity 0 . The values of Bauer's parameters 0 , 2 , and were determined from the apparent viscosity curve of grease A, as shown in Figure 9 of Kochi et al. [1]. When = 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 , 2 , and were determined to ensure that Bauer's curve passes through the abovementioned three points, as follows: 0 = 0.000621839 MPa, 2 = 6.99118 · 10 −7 , = 0.7248 (54)   [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. 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 agreement 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 , ℎ, and 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 is large. In addition, only the appearance of the distribution of differs from that of the base oil. Figure 5d shows the distribution of at section Y = 0. It can be seen that becomes extremely large at the center, where the pressure gradient is small and the flow volume is low.  [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]: Here, 0 and ℎ denote the X-direction traction forces acting on the lower and upper surfaces, respectively; is the load. The calculation range was -0.8 ≤ X ≤ 0.5, -0.5 ≤ Y ≤ 0.5. The calculation results were in fairly good agreement with the experimental results.
At points C and D, the direction of flow and that of 1 is orthogonal, so Δ = 0. Consequently, the equivalent viscosity parameters , and , at points C and D are given as follows: At point A, where 1 directs towards +X and the velocity gradient at Z = ℎ is greater than that at Z = 0 (as shown in Figure 8b), the following equation is satisfied: At point B, where 1 directs towards −X, and the velocity gradient at Z = ℎ is smaller than that at Z = 0 (as shown in Figure 8c), the following equation is satisfied: * (0) < * (ℎ)，Δ = −( + − − ) < 0 In any case, the equivalent viscosity parameters , and , at points A and B become smaller than , and , by the effect of the second term of Equation (57). In the pure rolling case, where Δ 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.

Conclusions
In this study, an isothermal, non-Newtonian EHL formulation of Bauer's model was performed using the local coordinate system of the pressure gradient. The flow toward the pressure gradient was assumed to be a Poiseuille flow and was approximated using a 4 th -order polynomial of z. The flow along the direction of the zero pressure gradient was assumed to be a Couette flow and was approximated using a linear function of z. The following results were obtained.
(1) A modified Reynolds equation, which contains an equivalent viscosity, was obtained. (2) The EHL calculation procedure for Newtonian flows can be applied to non-Newtonian flows by simply replacing the Newtonian viscosity with the equivalent viscosity. (3) If rheological equations are incorporated, any isothermal, non-Newtonian EHL calculation can be performed easily. (4) As the equivalent viscosity is calculated using the one-variable Newton-Raphson method, the EHL calculation can be performed within a reasonable calculation time. (5) Using Bauer's model, the formulation was applied to a grease that was evaluated experimentally by Kochi et al. [1]. The results obtained using the proposed method and the experimental results were compared, and reasonable agreement was noted. (6) In the case of sliding velocity, the equivalent viscosity shows a figure-eight-shaped distribution in the vicinity of the contact point.
However, the proposed method yields an approximate solution. If the Poiseuille flow and Couette flow cannot be approximated using a 4th-order polynomial of z and a linear function of z, respectively, the obtained results may be inaccurate. The application limits of the current formulation are not clear. Therefore, future work should be focused on determining these limitations.
Funding: This research was funded by the TriboLogics Corporation.

Acknowledgments:
The author sincerely thanks TriboLogics Corporation (www.tribology.co.jp/indexEng.htm) for permitting the use of TED/CPA. The author is also very grateful to Hiroyuki OHTA of the Nagaoka University of Technology, Japan, for his constant guidance, encouragement, and valuable suggestions.

Conflicts of Interest:
The author declares no conflict of interest.