DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

Numerical Approach of Elastohydrodynamically Lubricated Line Contact with Seed Oils Dedi Rosa Putra Cupu

Jamaluddin Md. Sheriff, Kahar Osman

Department of Mechanical Engineering University of Riau Pekanbaru, Indonesia e-mail: [email protected]

Faculty of Mechanical Engineering Universiti Teknologi Malaysia Skudai, Malaysia e-mail: [email protected]

Abstract — This paper proposes a numerical analysis for elastohydrodynamic lubrication (EHL) in line contact with seed oils as lubricants in order to estimate the pressure profiles and film thickness of EHL. The analysis is applied the contact between outer ring and cylindrical roller of bearing. Numerical solution was improved from Houpert and Hamrock’s fast approach calculation. The Newton-Raphson method was employed to solve two-dimensional Reynolds equation and Barus’ model was used to calculate the viscosity involving the pressure. Seed oils were used as lubricants. Results show that sunflower and soybean oil have relatively similar EHL pressure profile and film shape. All tested seed oils were able to replace mineral oil as lubricant due to their good pressure profiles and film thickness distributions. Increasing load, velocity and lubricant’s temperature affected the maximum spike pressure and minimum film thickness, where applied load and velocity of element bearing were increased as increasing load and velocity while increasing the lubricant’s temperature caused the pressure spike and film thickness decreased. Keywords – component, line contact of elastohydrodynamic lubrication, seed oils, Newton-Raphson method

I.

deformations and pressure gradient dP/dX, especially in the inlet region and near the pressure spike, where dP/dX is large. This method will be considered in this study because it offers more reliable results due to large pressure fluctuation. This fluctuation then will cause drastic changes in the film thickness.

INTRODUCTION

Elastohydrodynamic Lubrication is one of the hydrodynamic lubrication which involving physical interaction between the contacting bodies and the liquid of lubricant will cause these contacting surfaces will be deformed elastically and the changes of viscosity with pressure play fundamental roles. The contacting surfaces in many engineering applications, for example, rolling contact bearings, gears, cams, seals, etc., are non-conformal, therefore the resulting contacts areas are very small and the resulting pressures are very high [1].

For modeling of this simulation, the line contact of elastohydrodynamic lubrication is represented using the contact between outer ring and cylindrical roller of roller element bearing. Because the radius of outer ring is larger than the radius of cylindrical roller, it can be assumed that this contact is modeled as the contact between roller and flat plane, as shown in Fig. 1. A roller bearing is normally used to support a rotating shaft to the bearing housing and at the same time used to reduce the friction between those contacting surfaces.

Based on their solid contacted bodies, EHL is generally consists of two types problems, line contact and point contact. Contact between two spherical balls and contact between spherical ball and flat surface are represented as point contact problems. Cylindrical roller bearing is represented as the line contact problem. In this type, the rolling and load zones are angularly centered and rolling zone is smaller than the load zone [2]. In this paper, the simulation is proposed to solve the problem occurring on the real roller element bearing to calculate the generated pressure between roller element and outer raceway. Using calculation by numerically is very important before conducting experiments to measure the pressure and film thickness between the two contact bodies is aimed to reduce operational time and material destruction. This paper is extended from the paper of [3] and the numerical calculation is based on fast approach calculation of film thickness and pressure by Houpert and Hamrock [4]. Houpert and Hamrock proposed on how to solve EHL problem at high loads by calculating the elastic

DOI 10.5013/IJSSST.a.13.3A.01

Analysis will be conducted for two dimensional Reynolds only; an infinite length of steel roller will be simulated. Material roller and disc of bearing are steel AISI 52100 with elastic modulus is 210 GPa and poisons’ ratio is 0.3. Detail physical parameters of roller bearing that will be involved in the calculation are given in Table 1. A number of researchers have been attracted to investigate the influence of thermal effects on EHL experimentally [5] and also numerically [6]. Dow et al. [5] used a twin-disk machine which consists of a pair of 35 mm diameter disks. Based on their experimental results, it is clearly seen that the slip percentage of rolling/sliding had the major effects on the peak surface temperature which the higher the slip percentage, the higher the temperature generated during experiments. The increase of load caused small effects on the increase of peak surface temperature.

1

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

Seed oils will be used as lubricants in this analysis due to its high viscosity index and excellent coefficient of friction. L. A. Quinchia et al. [7] have tested the physical properties of seed oils, such as sunflower oil, soybean oil and castor oil. Densities of these oils were measured in a capillary densimeter in temperature range between 25 and 120oC. Dynamic viscosities were measured by using a rotational controlled-strain rheometer in the same temperature range. Kinematic viscosity was calculated by dividing the dynamic viscosity to the density. ASTM D 2270 was used as a standard to get the value of viscosity index (VI). The temperature for this simulation is assumed to be constant; Table 2 shows only dynamic viscosity (η), kinematic viscosity (υ) in temperature 40oC and viscosity index for seed oils studied which were based on works by L. A. Quinchia et al. [7]. II.

Figure 1. The contact between cylindrical roller and outer ring of roller element bearing.

Assuming the fluid flow is laminar and compressibility of the fluid is negligible, the Reynolds equation for an EHL line contact as shown in Fig. 1 can be written as:

Property

Dynamic viscosity (Pa.s) Kinematic viscosity (cSt) Viscosity Index

  3 p  h h   12uη  x  x  x where the entrainment speed, u 

u1  u 2 2

(2)

These dimensionless parameters will be used for performing magnitude analysis on Reynolds equation:

η ρ ;ρ  η0 ρ0



ρi  1 

  3 dP  3 π 2 U η H  H  X  dX  4 W 2 X

Elastic modulus (roller and disc)

Er = Ed = 210 GPa

Poison’s ratio

υr = υd = 0.3

Radius of roller

R = 0.0114 m

3.1 x 10-2

2.3 x 10-1

32.9 ± 2.3

33.6 ± 0.9

242.5 ± 21.7

248

250

116

0.6  10 9 p H p i 1  1.7  10 9 p H p i

 

The EHL film thickness at any point x on the contacted surface is:

PHYSICAL PARAMETERS OF ROLLER BEARING

Steel AISI 52100

2.9 x 10-2

Dimensionless density distribution is given by Dowson and Higginson [8] as follow:



Material (roller and disc)

Seed Oils Soybean Castor oil oil

η  η 0 e αp 

The Reynolds equation in dimensionless quantities can be written as:

TABLE I.

Sunflow er oil

Viscosity and density of lubricants are functions of pressure. Barus equation will be used for calculation of viscosity.

η0 u n w hR x W ;H  2 ;U  ; G  αE' ; X  n E' R b E' R b η

PHYSICAL PROPERTIES OF TESTED OIL AT 40OC

TABLE II.

BASIC EQUATION

Hi  H0 

X i2   δi 2 

Where the elastic deformation of solid is:

i   

1 2

X end

X

in





dP  X  X ' ln X  X '2  2 dX ' dX '

1  2 8W   ln R  4    

The equation for applied load in dimensionless form can be expressed as:

DOI 10.5013/IJSSST.a.13.3A.01

2

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

X end

X

PdX 

in

π  2 

 III.

NUMERICAL SOLUTION

Newton-Raphson iterative method as explained by [4], will be used in order to solve the non-linear of Reynold equation where the unknowns in this problem are: Xout - outlet boundary N - number of nodes - central film thickness H0 ρ e H e - the value of density times film thickness, ρH at P X  0 Pj - pressure at node j (j = 2, 3, …, N)

Figure 2. The geometry of EHL line contact.

Where subscript n is meaning new and o is old. The Jacobian factors are defined by [4] as:

Dimensionless Reynolds equation (4) at each node can be written as function of fi,

f i  (ρ e H e )

 ρ H   P  f i  H 3i     K η i H i  e e   0 ρ i   X  i   where K 



at X = Xmin, P = 0 and P at X = Xout, P  0 X



e

e

e









n

 

From the Newton-Raphson algorithm definition, the function of f i can be written as: o



 f   i H   0



n

o

 f     i  ΔPj   j 2  P j  N

 

n





  Hi  ρeHe  ρi 

   

  





f1 P2 f2 P2   f N P2 C2



    

f1 PN f2 PN   f N PN CN

  H 0  f2  H 0      f N   H 0  0  f1

o

Δ ρ e H e   Δ P2          Δ P    N   Δ H 0  

n

  f1    f2          f   ΔWN   

o



where ΔPj and ΔW are additional constant loads condition as required to complete equation (17) to be above matrix (18). This load condition can be written:

o

  ΔH 0 n   f i o   

DOI 10.5013/IJSSST.a.13.3A.01

H in  2

 f1   ρ H   ef e  2   ρ e H e       fN   ρ H   e e   0

H 0 n  H 0 o  ΔH 0 n 

 f i     Δ ρeHe   (ρ e H e ) 



A linear system of N + 1 equation is solved as below matrix:

 Δ ρ e H e    

(P j ) n  (P j ) o  Δ(P j ) n

(γγ n

 dP        dX  i   K  η i P j P j



successive iterations represent the new and old value after iteration. e

1 n

  1 ρi  K η i  D ij  D oj  ρ e H e P j 

The unknowns  e H e , Pj , H 0 are achieved by two

ρ H   ρ H 



f i n  2 1n n 1  dP   γ  H   D ij  D oj  P j n  dX  i

The boundary conditions for equation (10) are:

o

ρi

f i n  2 1 n n 1  dP  γ  H    K ηi   P j n  dX  i

3 2 U π and Hi as introduced in equation (7). 4 W2

n

Kηi



N  j2 C j ΔPj  n  ΔWn 

3

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

The convergence criteria used to stop the equation iteration for this study is given below:

Pj n  Pj o Pj n IV.

 ε  0.00001

pressure was used as pressure initial condition and Barus equation was used for relation of the viscosity and pressure. The results of these simulations for seed lubricants are shown in Figs. 3 to 6. Convergence of each simulation was taking about 4 until 7 seconds by using personal computer with 6 until 14 iterations.

 

Fig. 3 shows the EHL pressure profile and film thickness for soybean oil. This plot using Nmax = 321 nodes and Barus pressure-viscosity expression. In this simulation, starting point for dimensionless central film thickness was H0 = 0.6. Convergence iteration was got at the sixth iteration with absolute relative error was ε = 0.00001 and the outlet boundary at Xend = 1.16238. EHL pressure increased from the inlet flow and reached the maximum pressure about Pmax = 0.95908 the center roller, as shown in Fig. 2(a). Then, this pressure decreased from the center of roller onwards and rose again until reaching to the second maximum pressure, namely the pressure spike with Ps = 0.97086. The shape of elastohydrodynamic film thickness at contacted surface area was almost flat with H0 = 0.36479 but there was a dimple shape at the end of outlet flow. The minimum EHL film thickness was Hmin = 0.30910 (Fig. 2(b)).

RESULTS AND DISCUSSIONS

Pressure-viscosity coefficient was obtained from the following equation involving absolute viscosity and Roelands’ pressure-viscosity index as mentioned in [9].





α  Z 5.1  10 9 lnη 0  9.67 

η0 is the absolute viscosity and Z is a constant characteristic of the liquid (pressure-viscosity index) as below formulation, Z  7.81H 40  H100 1.5 F40   where

F40  0.885  0.864H 40  

Fig. 4 shows pressure profile for all seed oils. In this plot, parameters were identical to Fig. 2. Pressure profile and film thickness of sunflower oil almost coincided with soybean oil due to their low-temperature behavior, as reported by [5] that sunflower and soybean oil have relatively similar fatty acids compositions.

H100  loglogη 40   1.200  H100  loglogη100   1.200  Material parameter, G was calculated by involving pressure-viscosity index as mentioned in equation (21) and equivalent Young’s modulus ( E ' ), as following formula:

Fig. 5 shows the effect of loading variation on EHL pressure profile and film thickness with soybean oil as the lubricant. In this simulation, dimensionless speed is set to be U = 1.0 x 10-11.

G  αE' 

where, 1  1  ν 2r 1  ν d2 E'    2  E r Ed

At lower load (W = 1.0x10-5), as shown in Fig. 4(a), it needs 9 iterations to converge the problem. The value of spike pressure is higher than Hertzian pressure and other load variations where the pressure maximum is Pmax = 1.06166. As the load is increased up to W=10.0x10-5, the value of spike pressure is decreased and almost similar with Hertzian pressure, and Pmax is obtained about 0.99261.

     

subscript ‘r’ is roller and subscript ‘d’ is disc of bearing. In order to investigate the thermal effects on the characteristics of elastohydrodynamic lubrication, such as pressure and minimum film thickness, the program code is running at different temperatures between 0oC and 100oC with all vegetable oils. It should be noted that the dimensionless load and dimensionless speed are exactly the same as the parameters running at the first stage. However, in this study, the viscosity-temperature relationship proposed by Khonsari and Booser [9] is used to get the kinematic viscosity of oil (ν) at such temperature.

Referring to the Hertzian film distribution (Fig. 5 (d)), when the simulation runs at the dimensionless load of W = 1.5x10-05, the shape of the film thickness distribution is almost coincident to the Hertzian film. It means that this condition is ideal for the application of roller element bearing lubrication. TABLE III.

Complete calculated physical properties for this simulation are given in Table 3. The dimensionless speed and dimensionless load used in this paper were fixed as U = 1.0x10-11 and W = 2.0452x10-5. The inlet boundary was set at X = -4 and the outlet boundary was at X = 1.5. Hertzian

DOI 10.5013/IJSSST.a.13.3A.01

COMPLETE PARAMETERS AND PROPERTIES OF LUBRICANTS

Oil Name Sunflower oil Soybean oil Castor oil

4

µ0 at 40oC (Pa.s) 0.029 0.031 0.230

α (m2/N) 1.3008E-08 1.3229E-08 2.3283E-08

G 3002 3053 5373

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

(a) EHL Pressure and Hertzian pressure profile

(b) EHL film thickness

Figure 3. EHL pressure profile and film thickness shape for soybean oil. U = 1.0x10-11, and W = 2.0452x10-5.

(a) Pressure profiles for all tested oils

(c) Film thickness profiles for all tested oils

(b) Blowup pressure profiles for all tested oils at contact area

(d) Blowup film thickness profiles for all tested oils at contact area

Figure 4. EHL pressure profile and film thickness shape for all tested oils. U = 1.0x10-11, and W = 2.0452x10-5.

DOI 10.5013/IJSSST.a.13.3A.01

5

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

(b) Blowup pressure profile at contact area

(a) EHL pressure profile

(c) EHL film thickness distribution

(d) Blowup film thickness distribution at contact area

Figure 5. Effect of loading variation for soybean oil. U = 1.0x10-11.

elements while the higher the dimensionless speed the higher the pressure and the thicker the fim thickness shaped.

Fig. 6 shows the effect of velocity on pressure profile and film thickness. Material is still using soybean oil with dimensionless load is held fixed at 2.0452x10-5. At the lowest speed (U = 0.5x10-11), as shown in Fig. 6(a), converged simulation is finished at seventh iteration, and obtain the pressure spike value is 0.96899 (smaller than Hertzian pressure and other speed variation). Increasing speed cause raising the value of pressure spike, and reach the maximum pressure spike at 1.20721 (higher than Hertzian pressure).

Fig. 7 shows the effect of temperature on the EHL film thickness with soybean oil is used as the lubricant. At the lower temperature, the value of the film thickness formed is much higher (hmin = 55.3546 nm). As increasing temperature of lubricant, the value of minimum film thickness and central film thickness decreased, and the location of minimum film thickness shifted to the centre of roller. This can be explained because the value of the pressure-viscosity coefficient of α is reduced at higher temperature (see Table 5), and hence the film thickness is decreased. This result of the temperature effect on the film thickness is in accordance with the EHL theory [8] and the EHL experimental result [10].

Table 4 summarizes the variation of dimensionless pressure spike and minimum film thickness for the effects of various dimensionless load and dimensionless speed. In this simulation, soybean oil with material parameters of G = 3053 is used as the lubricant. It can be seen that the higher the dimensionless load the thinner the film thickness shaped and the lower the pressure working on the two contacted

DOI 10.5013/IJSSST.a.13.3A.01

6

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

(c) EHL film thickness distribution

(a) EHL pressure profile

(b) Blowup pressure profile at contact area

(d) Blowup film thickness distribution at contact area

Figure 6. Effect of velocity variation for soybean oil. W = 2.0452x10-5.

Figure 7. Thermal effect on EHL line contact problem using soybean oil.

DOI 10.5013/IJSSST.a.13.3A.01

7

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . . TABLE IV.

EFFECT OF DIMENSIONLESS LOAD AND SPEED ON DIMENSIONLESS PRESSURE AND FILM THICKNESS

Dimensionless speed parameter

Dimensionless load parameter

w W E' R

U

1.0 x 10-5 1.5 x 10-5 2.0 x 10-5 4.0 x 10-5 1.0 x 10-4 2.0452 x 10-5 2.0452 x 10-5 2.0452 x 10-5 2.0452 x 10-5

u

Pressure spike (Ps)

Minimum film thickness (Hmin)

1.06166 1.03413 0.97428 0.78821 0.49336 0.85334 0.97086 1.09415 1.20721

0.67714 0.43496 0.31680 0.14662 0.05309 0.19272 0.30910 0.49282 0.77793

Result

E' R

1.0 x 10-11 1.0 x 10-11 1.0 x 10-11 1.0 x 10-11 1.0 x 10-11 0.5 x 10-11 1.0 x 10-11 2.0 x 10-11 4.0 x 10-11

Load

Speed

 TABLE V.



PARAMETERS OF EHL FOR VARYING TEMPERATURE USING SOYBEAN OIL

T (oC)

Ps

ps (GPa)

H0

h0 (nanometer)

Hmin

hmin (nanometer)

0 20 40 60 80 100

1.2714 1.0087 0.9686 0.9442 0.9302 0.9223

0.5293 0.4200 0.4033 0.3931 0.3873 0.3840

0.4067 0.3828 0.3637 0.3479 0.3348 0.3236

64.5962 60.8093 57.7690 55.2672 53.1732 51.3961

0.3485 0.3252 0.3081 0.2942 0.2827 0.2729

55.3546 51.6519 48.9452 46.7388 44.9009 43.3490

V.

Proceedings of 5th Asia Modelling Symposium, Kuala Lumpur, Malaysia, 2011. [4] L. G. Houpert and B. J. Hamrock, "Fast approach for calculating film thicknesses and pressures in elastohydrodynamically lubricated contacts at high loads," ASME Journal of Tribology, vol. 108, 1986, pp. 411-420. [5] Dow, T. A., Stockwell, R. D. and Kannel, J. W. "Thermal Effects in Rolling/Sliding EHD Contacts: Part I – Experimental Measurements of Surface Pressure and Temperature," Journal of Tribology, Transaction of ASME, vol. 109(3), 1987, pp.503-510. [6] Sadeghi, F. and Dow T. A. "Thermal Effects in Rolling/Sliding Contacts: Part II – Analysis of Thermal Effects in Fluid Film," Journal of Tribology, Transaction of ASME, vol. 109(3), 1987, pp.512-517. [7] L. A. Quinchia, M. A. Delgado, C. Valencia, J. M. Franco, and C. Gallegos, “Viscosity modification of different vegetable oils with EVA copolymer for lubricant applications,” Industrial Crops and Products, vol. 32, 2010, pp. 607-612. [8] D. Dowson and G. R. Higginson, Elastohydrodynamic Lubrication. Pergamon Press: Oxford, 1977. [9] M. M. Khonsari and E. R. Booser, Applied Tribology, Second Edition. Chichester, England: John Wiley and Sons Ltd, ch. 2, pp. 2362. [10] Biresaw, G. "Elastohydrodynamic Properties of Seed Oils," Journal of the American Oil Chemists’ Society (JAOCS), vol. 83(6), 2006, pp. 559-566.

CONCLUSSION

The EHL pressures and film thicknesses in line contact have been simulated by using fast approach method, where Newton-Raphson method has been used to solve the Reynolds equation. The results indicate that at constant temperature 40oC, sunflower oil, soybean oil and castor oil were able to use to replace the current mineral oil and synthetic oil because these seed oils have good pressure profile and film shape during simulation. The spike pressure increases as increasing load and velocity, and film thickness is also increased. Dimple shape is occurred from center roller to the outlet region as increasing pressure spike. Increasing lubricant’s temperature decreased pressure and film thickness. REFERENCES [1]

G. W. Stachowiak and A. W. Batchelor, Engineering Tribology, Third Edition. Burlington: Elsevier ButterworthHeinemann, 2008, ch. 7, pp. 287-362.

[2]

E. Laniado-Jacome, J. Meneses-Alonso, V. Diaz-Lopez, "A study of sliding between rollers and races in a roller bearing with a numerical model for mechanical event simulations," Tribology International, vol. 43, 2010, pp. 2175-2182. J. MdSheriff, K. Osman and D., R., P. Cupu, "Numerical analysis of line contact elastohydrodynamic lubrication with biobased oils,"

[3]

DOI 10.5013/IJSSST.a.13.3A.01

8

ISSN: 1473-804x online, 1473-8031 print

DEDI ROSA PUTRA CUPU et al: NUMERICAL APPROACH OF ELASTO-HYDRO-DYNAMICALLY . .

NOMENCLATURE b Cj Dij E E G H He H0 h i,j N Nmax p P pH R u U W w x X   0

  0  

Hertzian contact radius (m) Weighting factor used in linear equation Influence coefficient Elastic Modulus (Pa) Effective Elastic Modulus Material Parameter Dimensionless film thickness Dimensionless film thickness where dP/dX=0 Dimensionless Central Film Thickness Film thickness (m) Nodes Number of nodes used in linear equation Maximum Number of nodes Pressure (Pa) Dimensionless Pressure Hertzian Pressure (Pa) Radius of roller (m) Average entrainment rolling speed (m/s) Dimensionless Speed Dimensionless Load Load per unit length (N/m) Distance along rolling direction (m) Dimensionless Distance Piezoviscous coefficient (m2/N) Viscosity of lubricant (Ns/m2) Viscosity at ambient pressure (Ns/m2) Dimensionless Viscosity Density of lubricant (kg/m3) Density at ambient pressure (kg/m3) Dimensionless Density Poisson Ratio

DOI 10.5013/IJSSST.a.13.3A.01

9

ISSN: 1473-804x online, 1473-8031 print