Multigrid method for the solution of EHL line contact with bio-based oils as lubricants

The Elastohydrodynamic lubrication (EHL) is one of the important topics in Tribology, EHL problems have been studied by many researchers over the last three decades both theoretically and experimentally. EHL is commonly found in many machines components such as gear teeth contacts, cams and rolling element bearings during variations in load etc. There are generally two types of problems in EHL, line contact and point contact. In line contact EHL problem, the contact region of the elements are assumed to be infinitely long in one direction, that takes place as an infinitely long strips. In point contact EHL problem, the contact takes place within a finite elliptical region [1]. When two metals are in contact, the amount of asperities and interaction within the contact area increases, causing frictions which insist wear of the metal surfaces and generation of excessive heat. These friction, wear and excessive heat caused by the interaction between the surfaces are controlled by the inclusion of lubricant whose function is to reduce friction and wear to prevent oxidation and corrosion, acting as a coolant, facilitating heat dissipation from the engine.

The Bio-based oils are usually the extracts of various plants, animals, seeds, fruits etc. Bio based product plays an important role in lubrication [2, 3, 4, 5]. Certain physical and chemical property of these oils enables it to be a suitable lubricant due to their liquidity at room temperature. In spite of these advantages it has several drawbacks for their use as lubricant. Chemical modification can overcome these drawbacks by adding additives in these oils. Among these, vegetable oils consists of long chain molecule structure which interacts strongly with metal surface reducing friction and wear, due to the availability of their small range of viscosities. These oils have been discussed by many researchers; Quinchia et al. [6] have successfully modified the viscosity of vegetable oils by using Ethylene-Vinyl Acetate copolymer (EVA) as additive material. There is a strong concern of awareness in biodegradable low-toxicity lubricants and biodegradability is obtained by using suitable bio-degradable base fluid, by adding additives that is environmentally friendly [7]. There is a need to source out lubricants other than the computational ones obtained from biobased oils [8, 9], as an increasing world energy crisis.

Multigrid, Multilevel Multi-Integration (MLMI) with Full Approximation Scheme (FAS) are been discussed in [10, 11]. Multigrid technique has received immense attention from many investigators in EHL and is proved to be very efficient for solving EHL problems. Venner [12] examined FAS for the Reynolds equation, film thickness equation and the force balance equation by employing multigrid technique for the solution of EHL problems. Lubrecht et al. [13] presented its application in EHL and employed the reduced pressure in deriving the Reynolds equation. Brandt and Lubrecht [14] developed a MLMI algorithm for the fast calculation of elastic deformation. In the present work only the Gauss-Seidel relaxation method is adopted. The advantage of multigrid methods incorporate naturally and overcome the drawbacks of Newton-Raphson scheme, it enables to analyze EHL problems for higher loads required in non-conformal conjunctions such as rolling element bearings, gears etc., employing the Gauss-Seidel relaxation with small under-relaxation factors which illustrates that the Gauss-Seidal relaxation is a stable relaxation scheme with good smoothing properties. EHL problems have also been solved by many other different methods as in [15, 16], etc. Although many mathematical models of ordinary differential equations and partial differential equations [17, 18], etc., have been solved for its exact solution.

In this paper, we employ FAS for the solution of EHL line contact problem. This paper is organized as follows. Section 1 describes the introduction. Section 2 describes governing equation of the physical problem. Discretization of the governing equation is given in Section 3. In Section 4, Multigrid (FAS) method is presented. In Section 5, solution procedure is given in detail. Section 6 contains the results and discussions. Finally, the conclusion of the proposed work is given in Section 7.

Let us consider the flow is laminar and the compressibility of the fluid is negligible. The isothermal EHL line contact problem models the lubricant flow contact between outer ring and cylindrical roller of roller element bearing under an applied load which consists of three equations: the Reynolds equation, the film thickness equation and the force balance equation. All these equations are presented in a dimensionless form as [19]:

∂<!-- ? -->∂<!-- ? -->X(ε<!-- ? -->∂<!-- ? -->P∂<!-- ? -->X)−<!-- - -->∂<!-- ? -->∂<!-- ? -->X(ρ<!-- ? -->¯H),$$\begin{array}{} \displaystyle \frac{\partial}{\partial X}\left(\varepsilon \frac{\partial P}{\partial X} \right)-\frac{\partial}{\partial X}\left(\overline{\rho}H \right), \end{array} $$(1)

where ε<!-- ? -->=ρ<!-- ? -->¯H3η<!-- ? -->¯KandK=34(π<!-- ? -->W)2U.$\begin{array}{} \varepsilon =\frac{\overline{\rho}H^3}{\overline{\eta}K}\,\,\text{and}\,\,K=\frac{3}{4}\left(\frac{\pi}{W}\right)^2 U. \end{array} $

The non-dimensional parameters for magnitude analysis of Reynolds equation are

W=wE′R,H=hRb2,U=η<!-- ? -->0uE′R,G=α<!-- ? -->E′,X=xb,η<!-- ? -->¯=η<!-- ? -->η<!-- ? -->0,ρ<!-- ? -->¯=ρ<!-- ? -->ρ<!-- ? -->0,P=ppH.$$\begin{array}{} \displaystyle W=\frac{\, w}{E'R},\quad H=\frac{\, hR}{b^2},\quad U=\frac{\eta_0u}{E'R},\quad G=\alpha E',\quad X=\frac{x}{b},\quad \overline{\eta}=\frac{\eta}{\eta_0},\quad \overline{\rho}=\frac{\rho}{\rho_0},\quad P=\frac{p}{p_H}. \end{array} $$

The film thickness equation at any point on the contacted surface is [20]:

H(X)=H0+X22−<!-- - -->D¯,$$\begin{array}{} \displaystyle H(X)=H_0+\frac{X^2}{2}-\overline{D}, \end{array} $$

where H₀ is the central offset film thickness, X22$\begin{array}{} \frac{X^2}{2} \end{array} $ defines the un-deformed contact shape, and

D¯=−<!-- - -->12π<!-- ? -->∫<!-- ? -->XinXoutdPdX′(X−<!-- - -->X′)[log⁡<!-- -->(X−<!-- - -->X′)2]dX′−<!-- - -->14(R28Wπ<!-- ? -->)$$\begin{array}{} \displaystyle \overline{D}=-\frac{1}{2\pi} \int_{X_{in}}^{X_{out}}\frac{dP}{dX'}(X-X')[\log(X-X')^2]dX'-\frac{1}{4}\left(R^2\frac{8W}{\pi}\right) \end{array} $$

is the elastic deformation of the contacting surfaces. The equation for applied load in dimensionless form can be expressed as

∫<!-- ? -->XinXoutPdx=π<!-- ? -->2.$$\begin{array}{} \displaystyle \int_{X_{in}}^{X_{out}}Pdx=\frac{\pi}{2}. \end{array} $$

The non-dimensional viscosity-pressure relation is given by [21] and pressure viscosity index as mentioned in [22]:

η<!-- ? -->¯=exp{[log⁡<!-- -->(η<!-- ? -->0)+9.67][−<!-- - -->1+(1+pHPp0)z]}α<!-- ? -->=z[5.1×10−<!-- - -->9(log⁡<!-- -->η<!-- ? -->0+9.67)],$$\begin{array}{} \displaystyle \overline{\eta}=exp\left\{[\log(\eta_0)+9.67]\left[-1+\left(1+\frac{p_HP}{p_0}\right)^z\right]\right\}\\\displaystyle \alpha=z[5.1\times10^-9\,(\log \eta_0+9.67)], \end{array} $$

where η₀ is the absolute viscosity, p₀ is the ambient pressure and z is a constant characteristic of the liquid (pressure-viscosity index) given by

z=[7.81(H40−<!-- - -->H100]1.5F40,$$\begin{array}{} \displaystyle z=[7.81(H_{40}-H_{100}]^{1.5}F_{40}, \end{array} $$

where

F40=(0.885−<!-- - -->0.864),H40=log⁡<!-- -->(log⁡<!-- -->(η<!-- ? -->40)+1.200),H100=log⁡<!-- -->(log⁡<!-- -->(η<!-- ? -->100)+1.200).$$\begin{array}{} \displaystyle F_{40}=(0.885-0.864),\quad H_{40}=\log(\log(\eta_{40})+1.200),\quad H_{100}=\log(\log(\eta_{100})+1.200). \end{array} $$

The non-dimensional density-pressure relation is given by [1]:

ρ<!-- ? -->¯=(1+0.6p1+1.7p).$$\begin{array}{} \displaystyle \overline{\rho}=\left(1+\frac{0.6p}{1+1.7p}\right). \end{array} $$

The corresponding boundary conditions are

Inlet boundary condition at: X = X_in, P = 0.

Outlet boundary condition at: X = X_out, P=∂<!-- ? -->P∂<!-- ? -->X=0.$\begin{array}{} P=\frac{\partial P}{\partial X}=0. \end{array} $

Eq. (1) is discretized using second order finite differences with the grid size of N = 65 and the domain of interest initially is [X_in,Xout] = [−4, 2.2]. The discretized form of Reynolds equation is

ε<!-- ? -->i−<!-- - -->12Pi−<!-- - -->1−<!-- - -->(ε<!-- ? -->i−<!-- - -->12+ε<!-- ? -->1+12)Pi+ε<!-- ? -->i+12Pi+1Δ<!-- ? -->X2=ρ<!-- ? -->¯iHi−<!-- - -->ρ<!-- ? -->¯iHi−<!-- - -->1Δ<!-- ? -->X,$$\begin{array}{} \displaystyle \frac{\varepsilon_{i-\frac{1}{2}}P_{i-1}-(\varepsilon_{i-\frac{1}{2}}+\varepsilon_{1+\frac{1}{2}})P_i+ \varepsilon_{i+\frac{1}{2}}P_{i+1}}{\Delta X^2}=\frac{\overline{\rho}_{i}H_i-\overline{\rho}_{i}H_{i-1}}{\Delta X}, \end{array} $$

where Δ X² = X_i − X_i−1. The film thickness equation in discretized form is [16]:

Hi=H0+xi22−<!-- - -->1π<!-- ? -->∑<!-- ? -->i=1NDi,jPj$$\begin{array}{} \displaystyle H_i=H_0+\frac{x_{i}^{2}}{2}-\frac{1}{\pi}\sum_{i=1}^{N}{D_{i,j}P_j} \end{array} $$

where

Di,j=−<!-- - -->(i−<!-- - -->j+12)Δ<!-- ? -->X[log⁡<!-- -->(|i−<!-- - -->j+12|Δ<!-- ? -->X)−<!-- - -->1]+(i−<!-- - -->j+12)Δ<!-- ? -->X[log⁡<!-- -->(|i−<!-- - -->j−<!-- - -->12|Δ<!-- ? -->X)−<!-- - -->1]$$\begin{array}{} \displaystyle D_{i,j}=-\left(i-j+\frac{1}{2}\right)\Delta X\left[\log\left(|i-j+\frac{1}{2}|\Delta X\right)-1\right]+\left(i-j+\frac{1}{2}\right)\Delta X\left[\log\left(|i-j-\frac{1}{2}|\Delta X\right)-1\right] \end{array} $$

for i = 0,1,2, …, n, j = 0,1,2, …, n and the force balance equation in discrete form is

Δ<!-- ? -->X∑<!-- ? -->i=1∞<!-- ? -->Pi+Pi+12−<!-- - -->π<!-- ? -->2=0.$$\begin{array}{} \displaystyle \Delta X\sum_{i=1}^{\infty}{\frac{P_i+P_{i+1}}{2}}-\frac{\pi}{2}=0. \end{array} $$

The boundary condition in the inlet P(X_in) = 0 and the boundary in the outlet can be obtained by setting the negative pressure equal to zero. The discrete form is

P(Xout)=0andP(Xout)−<!-- - -->P(Xout−<!-- - -->1)Δ<!-- ? -->X=0.$$\begin{array}{} \displaystyle P(X_{out})=0\quad \text{and} \quad \frac{P(X_{out})-P(X_{out}-1)}{\Delta X}=0. \end{array} $$

The convergence criteria for pressure relation is:

∑<!-- ? -->|Pik+1−<!-- - -->Pik|∑<!-- ? -->Pik+1≤<!-- ? -->ε<!-- ? -->1,$$\begin{array}{} \displaystyle \frac{\sum|P_{i}^{k+1}-P_{i}^{k}|}{\sum P_{i}^{k+1}}\leq \varepsilon_1, \end{array} $$

where ε₁ is the given error.

The problem of solving the system of linear equations

Au=f,$$\begin{array}{} \displaystyle Au=f, \end{array} $$

where A is a square matrix and f is a column matrix, arising from discretization of a differential equation on grid Ω^h, where h is the step size. Ih2h$\begin{array}{} I_{h}^{2h} \end{array} $ interpolation operator, the superscript denotes to the finer grid and the subscripts denotes to the coarser grid, and I2hh$\begin{array}{} I_{2h}^{h} \end{array} $ restriction operator. We define a multigrid method recursively, for the two-level V-cycle method, first relax a few steps on the finer grid Ω^h to get u^h, then compute the residual r^h = f^h−A^hu^h; restrict the residual to the coarser grid Ω^2h, r^2h = Ih2h$\begin{array}{} I_{h}^{2h} \end{array} $r^h; and solve the error equation A^2he^2h = r^2h on the coarser grid. Finally, set u^h = u^h+Ih2h$\begin{array}{} I_{h}^{2h} \end{array} $e^2h and again relax a few steps on the finer grid. Based on this two level method, the V-cycle multigrid scheme is defined recursively, which are given in [12, 23, 24].

For the calculation of deformation surfaces, we employ MLMI. MLMI is a multiscale technique designed to speed up the evaluation of integrals. Apart from its use in EHL problems, this method has many applications in integral differentiation equation, integral equation, elasticity and acoustic problems, etc. The general formulation to solve the following integral is given by

v(x)=∫<!-- ? -->Ω<!-- ? -->K(x,y)u(y)dy,$$\begin{array}{} \displaystyle v(x)=\int_{\Omega}K(x,y)u(y)dy, \end{array} $$(2)

where u is a given function and K is the kernel on a domain Ω = (a, b). The discrete form of the Eq. (2) is descretized on a regular mesh of Nxh$\begin{array}{} N_{x}^{h} \end{array} $ points, Δ<!-- ? -->xh=(b−<!-- - -->a)Nxh−<!-- - -->1,$\begin{array}{} \Delta x^h=\frac{(b-a)}{N_{x}^{h}-1}, \end{array} $ where h is the grid level. The single grid routine for calculating multi summation at each node xih$\begin{array}{} x_{i}^{h} \end{array} $ = a+(i−1)Δ x^h for i = 1, …, N_x is

vih=vh(xih)=∑<!-- ? -->j=1NxhKi,jhhujh,$$\begin{array}{} \displaystyle v_{i}^{h}=v^h(x_{i}^{h})=\sum_{j=1}^{N_{x}^{h}}K_{i,j}^{hh}u_{j}^{h}, \end{array} $$

where Ki,jhh$\begin{array}{} K_{i,j}^{hh} \end{array} $ is the discretized kernel and ujh$\begin{array}{} u_{j}^{h} \end{array} $ is the approximation to u^h(x_i). The discrete form of Eq. (2) represents a full matrix-vector multiplication having certain smoothness property, wherein the kernel is dense rather than a sparse matrix, the method is possible only when the kernel itself has sufficient smoothing properties. The discretization of the elastic deformation integrals for both EHL line and point contact is discussed by Venner and Luberchet [11].

The process of multigrid includes the pressure correction and load balance by adjusting the rigid film thickness, which are carried out on the same grid. At any instant, the iteration procedure consists of the pressure iteration and the adjustment for H₀. The pressure must be relaxed on each grid level. If we denote (LP)_i by L_i(P) and ignoring superscripts which represents the grid level, the algebraic equation for pressure on each level can be expressed as L_i(P) = F_i (i = 1,2, …, n−1). For an under-relaxation factor, the process for the pressure iteration can be described as:

P~<!-- ~ -->i=P¯i+ciδ<!-- ? -->i,$$\begin{array}{} \displaystyle \widetilde{P}_i=\overline{P}_i+c_i\delta_i, \end{array} $$

where

δ<!-- ? -->i=(∂<!-- ? -->Li/∂<!-- ? -->Pi)−<!-- - -->1ri,ri=Fi[ε<!-- ? -->i−<!-- - -->1/2P~<!-- ~ -->i−<!-- - -->1−<!-- - -->(ε<!-- ? -->i−<!-- - -->1/2+ε<!-- ? -->i+1/2)P¯i+ε<!-- ? -->i+1/2P¯i+1]/(Δ<!-- ? -->X)2+(ρ<!-- ? -->¯iH¯i−<!-- - -->ρ<!-- ? -->¯i−<!-- - -->1H¯i−<!-- - -->1)/(Δ<!-- ? -->X),∂<!-- ? -->Li/∂<!-- ? -->Pi=(ε<!-- ? -->i−<!-- - -->1/2+ε<!-- ? -->i+1/2)/(Δ<!-- ? -->X)2+1π<!-- ? -->(ρ<!-- ? -->¯iK¯i,j−<!-- - -->ρ<!-- ? -->¯i−<!-- - -->1K¯i−<!-- - -->1,j)/Δ<!-- ? -->X,$$\begin{array}{} \displaystyle \delta_i=(\partial L_i/\partial P_i)^{-1}r_i,\\\\\displaystyle r_i=F_i[\varepsilon_{i-1/2}\widetilde{P}_{i-1}-(\varepsilon_{i-1/2}+\varepsilon_{i+1/2})\overline{P}_i+\varepsilon_{i+1/2}\overline{P}_{i+1}]/(\Delta X)^2+(\overline{\rho}_i\overline{H}_i-\overline{\rho}_{i-1}\overline{H}_{i-1})/(\Delta X),\\\\\displaystyle \partial L_i/\partial P_i=(\varepsilon_{i-1/2}+\varepsilon_{i+1/2})/(\Delta X)^2+\frac{1}{\pi}(\overline{\rho}_i\overline{K}_{i,j}-\overline{\rho}_{i-1}\overline{K}_{i-1,j})/\Delta X, \end{array} $$

and P_i, P_i+1 are initial values of pressure in an iterative process and P͂_i−1 is the new value of pressure obtained in the iteration. The load balance condition can be achieved by modifying the rigid film thickness H₀ as follows.

H¯0=H~<!-- ~ -->0+c2[GΔ<!-- ? -->−<!-- - -->Δ<!-- ? -->π<!-- ? -->∑<!-- ? -->j=1N−<!-- - -->1(Pj+Pj+1)],$$\begin{array}{} \displaystyle \overline{H}_0=\widetilde{H}_0+c_2\left[G^\Delta-\frac{\Delta}{\pi}\sum_{j=1}^{N-1}(P_j+P_{j+1})\right], \end{array} $$

where c₂ is the relaxation factor and G^Δ is the nondimensional load on the coarsest grid.

Following procedure is followed for the determination of pressure cavitation point [25]. Let X_c be a point in (X_in,X_out). If pressure gradient dPdX|Xc=0,$\begin{array}{} \frac{dP}{dX}|_{X_c}=0, \end{array} $ then X_c is the required pressure cavitation point. Otherwise, shift the point X_c on next interval to the left or right depending on dPdX|Xc><!-- ? -->0ordPdX|Xc<<!-- ? -->0$\begin{array}{} \frac{dP}{dX}|_{X_c} \gt 0\,\,\text{or}\,\,\frac{dP}{dX}|{X_c} \lt 0 \end{array} $ respectively. Next approximate another cavitation point Xc′$\begin{array}{} X'_{c} \end{array} $ and solve EHL problem. If pressure gradient dPdX|Xc′$\begin{array}{} \frac{dP}{dX}|_{X'_c} \end{array} $ is zero, we get converged solution, otherwise verify the condition

dPdX|Xc×dPdX|Xc′<<!-- ? -->0.$$\begin{array}{} \displaystyle \frac{dP}{dX}|_{X_c}\times\frac{dP}{dX}|_{X'_c} \lt 0. \end{array} $$(3)

If the condition Eq. (3) is not satisfied, repeat above process and find two approximate values Xc′andXc″$\begin{array}{} X'_c\,\,\text{and}\,\,X''_c \end{array} $ for which Eq. (3) is satisfied. A new cavitation point X_new is obtained by interpolating between the two points Xc′andXc″$\begin{array}{} X'_c\,\,\text{and}\,\,X''_c \end{array} $, which ensures pressure gradient nearly equal to zero at X_new and results into the required smooth solution of EHL problem.

The EHL line contact problem is studied with bio-based oil as lubricant, for an isothermal case with the influence of different loads and speed. The problem has been calculated to investigate the substantial variations in various parameters including its calculated physical properties as shown in Tables 3 and 4. Multigrid, multilevel techniques with FAS scheme is employed. The cavitation condition in the present analysis ranges from X = −4 to 2.2 with grid size of N = 65. The present computations are confined to the numerical parameters c₁ = 0.037 and c₂ = 0.037 used for the under-relaxation of the iterations P and H₀. The result of the simulation for the lubrication is as shown in Figs. 1-3 and are carried out for a different value of dimensionless speed (U) and load parameter (W) as specified in both Tables 3 and 4. Fig. 1 shows the effects of varying load for Sunflower, Soyabean, and Castor oils with constant speed (U = 1.0×10⁻¹¹), where the load decreases as the value of pressure spike and film thickness increases and almost similar to hertzian pressure, film thickness. Fig. 2 shows the effects of varying speed at constant load (W = 2.045×10⁻⁵) as the speed increases the value of pressure spike and film thickness decreases. Fig. 3 shows the pressure and film thickness profiles of both Sunflower and Soyabean oils. They are almost the same due to their low temperature behavior as reported by Dow T. A. Dow et al. [26].

Fig. 1

Fig. 2

Fig. 3

Table 1

Physical parameters of roller bearing.

Table 2

Physical properties of tested oil at 40°C.

Table 3

Effect of dimensional load and speed on dimensional pressure and filmthickness.

Table 4

Effect of dimensional load and speed on dimensional pressure and filmthickness.

In this paper, the EHL line contact problem with seed oils as lubricant is analyzed by using multigrid, multilevel techniques in order to clarify the best performance of lubricants. The characteristics of EHL under operating condition shows that, at constant temperature 40°C, sunflower, soyabean and castor oil are able to replace other lubricants, as they have good pressure and film profile during simulation.