Development of a Rolling Chatter Model Considering Unsteady Lubrication

Abstract

Chatter, limiting the rolling speed for thin steel strips, has been identified as a serious limitation in increasing the efficiency of the cold rolling process. Many experimental investigations have led to the point that lubrication is one of the effective factors causing chatter. In this article, a new chatter model of the cold strip rolling with consideration of unsteady lubrication is proposed. The limiting shear stress concept and the variation of the lubricant viscosity with pressure are considered in this model. The results of the simulation are verified by comparing to experimental data. Also, a parametric study on the effect of some of the major parameters of rolling lubricant on the chatter critical speed is conducted. The results show that under the operating conditions of this paper, the limiting shear stress is an effective parameter.

1. Introduction

Chatter is the most destructive vibration problem for thin steel strips in the cold rolling processes. Rolling chatter has three different types: third octave, fifth octave, and torsional chatter.¹⁾ The first is the most important type because it can lead to large gauge fluctuations, surface finish imperfections, and damage to mills.²⁾ Also the lubricant flow between the work roll and the strip surface becomes extremely unsteady as a result of chatter. Many researchers have mentioned that lubrication is one of the basic factors incurring chatter.^3,4,5,6) Due to the complexity of the problem, however, most of them have used simple friction models in their analysis,^7,8,9,10) and only very few researches have described how chatter occurs through unsteady lubrication phenomenon.

Yarita et al.³⁾ concluded from their practical investigations that chatter arises from impaired lubrication. They suggested that lubricant film strength, emulsion stability, and particle size distribution should be taken into consideration in selecting the rolling lubricants, and also the use of a lubricant with good lubricity successfully prevent chatter. Johnson and Qi^11,12) and Chen et al.^13,14) considered friction at the roll-strip interface as a function of relative velocity between the rolls and strip. Their study shows the importance of friction in chatter predictions.

Tung¹⁵⁾ developed an unsteady lubrication model for vibrations in cold rolling. Kimura et al.⁴⁾ discovered the existence of an optimum range for friction coefficient in which the vibrations are damped and the rolling mill is stable despite the disturbances. Also, they showed that with increasing the rolling speed, the rolling mill tends to chatter.

Yun et al.^16,17,18,19) developed a fundamental understanding of chatter in rolling. Based on a constant friction factor, they presented a suitable dynamic model for studying chatter. Hu et al.^{2,20,21,22,23)} expanded Yan et al.’s work and pointed out that there are four principle mechanisms for third octave chatter: negative damping, mode coupling, regenerative, and model matching. They expressed that an appropriate friction model is necessary to perform accurate chatter simulation.

Meehan⁷⁾ investigated mode shapes and natural frequencies of a 4-high stand. It was seen from his research that third octave chatter and fifth octave chatter were vertically symmetrical about the stand pass line. Yang et al.⁸⁾ built a chatter model for cold rolling including the hydraulic servo system model. Also a mathematical model containing the driving system was offered by Kim et al.^9,24)

Wilson²⁵⁾ has defined four main lubrication regimes: the thick film and thin film regime (usually referred to as the full film regime), the mixed regime and the boundary regime. Most of models which considered lubrication have been developed under the assumption of steady flow of lubricant.^26,27,28) This assumption, however, is not appropriate for studying chatter.

In this paper, a new chatter model for the cold rolling process under the unsteady lubrication is presented. In order to verify the proposed model, experimental results obtained from the cold strip tandem mill unit of Mobarakeh Steel Company (MSC) are used. Based on this model, effect of rolling lubricant parameters such as the rolling lubricant viscosity, the viscosity-pressure coefficient, and the limiting shear stress on chatter critical speed is studied.

2. Methodology

Rolling chatter is discussed as a consequence of interactions between the rolling process and the mill structure.

2.1. Rolling Process Model

In this research, the lubrication regime is assumed to be the full film. Figure 1 illustrates the coordinate system and geometry of the rolling process. The roll bite is divided into inlet, work and outlet zone. The first zone is where the lubricant pressure increases from the atmospheric pressure to the pressure required incurring strip plastic deformation. The second zone is where the rolling force and torque are generated. The variation of pressure in the last zone is reverse of that in the first zone. The lubricant pressure in the first and last zone only cause strip elastic deformation.

Fig. 1.

Rolling process geometry.

2.1.1. Inlet Zone Analysis

As the first step in the rolling process model, an inlet zone analysis should be developed. The key output from this analysis is the inlet film thickness (h₀(t)).

The hydrodynamic pressure causes the separation between the strip and roll. This pressure (p), built up by wedge and squeeze action, is obtained from the Reynolds equation:   

$$\frac{\partial }{\partial x^{\prime }}\left[\frac{h_1{ }^3}{12\eta }\frac{\partial p}{\partial x^{\prime }}\right]=–{\overline{u}}_1\frac{\partial h_1}{\partial x^{\prime }}+\frac{\partial h_1}{\partial t}$$(1)

Where x' = x – x₁, x₁ is position of the boundary between the inlet and work zone, h₁ is the lubricant film thickness in the inlet zone, η is the lubricant viscosity, t is the time and ${\overline{u}}_1$ is the mean value of strip speed at entry and roll peripheral velocity. Since the lubricant viscosity depends on the pressure, the Barus equation²⁹⁾ is used to calculate the effective viscosity:   

$$\eta ={\eta }_0\hspace{0.17em}{e}^{\alpha p}$$(2)

In which η₀ is viscosity at the atmospheric pressure and α is the viscosity-pressure coefficient.

The boundary conditions in the inlet zone are:   

$$\begin{array}{l}{x}^{\prime }=0\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{h}_1=\hspace{0.17em}{h}_0\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}p={\sigma }_k–{\sigma }_{x,1}\\ x^{\prime }=\infty \hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{h}_1=\hspace{0.17em}\infty \hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}p=0\end{array}$$(3)

In the above equation, σ_k is the yield stress in the plane strain condition and σ_x,1 is the strip tensile stress at entry. Due to high ratio of strip width-to-strip thickness, it can be assumed that the rolling process operates in the plane strain condition.

In each time step increment, Eq. (1) is solved by applying the boundary conditions in order to determine the inlet film thickness (h₀(t)). The inlet film thickness shows the amount of the oil input into the work zone.

2.1.2. Work Zone Analysis

Having determined the variation of inlet film thickness from the inlet zone analysis, the work zone can now be analyzed. The output of this analysis is the work zone film thickness distribution (h(x, t)) which can be obtained from the continuity equation in fluid mechanics:   

$$h=\frac{v_r+u_1}{v_r+u}\hspace{0.17em}{h}_0$$(4)

Where v_r is the roll peripheral velocity and u₁ is the strip speed at entry.

2.1.3. Friction Model

Figure 2 illustrates the stresses acting on a vertical slab element inside the work zone. Assuming that the tensile stress (σ_x) and the pressure (p) are principal stresses and with eliminating the work hardening effect and the comparatively small terms, the force balance equation in the x direction can be represented as:   

$${\sigma }_k\frac{dy}{dx}–y\frac{dp}{dx}+2{\tau }_f=0$$(5)

In this equation, y is the local strip thickness and τ_f is the friction stress.

Fig. 2.

A slab element.

The friction stress in the full film regime can be expressed as:   

$${\tau }_f=\eta \frac{u–v_r}{h}$$(6)

Because of the high pressure in the work zone, the Roelands’ equation³⁰⁾ instead of the Barus equation should be applied. This formula is given as:   

$$\eta (p,T)={\eta }_{0}exp\left[\left\{ln({\eta }_0)+9.67\right\}\left\{{\left(1+\frac{p}{p_r}\right)}^z{\left(\frac{T_0–138}{T–138}\right)}^{s_0}–1\right\}\right]$$(7)

Where p_r is a constant value of 1.962×10⁸Pa, z is the viscosity-pressure coefficient, T₀ is the ambient temperature, s₀ is the viscosity temperature coefficient. Defining $\alpha =\frac{1}{\eta }{\left(\frac{\partial \eta }{\partial p}\right)}_{p=0,T=T_0}$ , the relation between the viscosity pressure coefficient in the Barus and Roelands equation is given by:   

$$\alpha =z\left[\frac{1}{p_r}\left(Ln{\eta }_0+9.67\right)\right]$$(8)

Meanwhile, the limiting shear stress (τ_L) concept is used in this study, which means that a viscoplastic behavior is assumed for the lubricant.^31,32) Thus, when the friction stress exceeds the limiting shear stress, the friction stress is set equal to the limiting shear stress.

In the friction model, the pressure in the work zone is computed based on the Eqs. (5), (6), (7), (8) and the film thickness distribution obtained from the work zone analysis. Integrating the pressure distribution gives an estimation of the rolling force. Also, considering the mass conservation and given the fact that the strip speed at neutral point is equal to the roll peripheral velocity, the strip speed at entry and exit have been obtained.^33,34)

2.2. Mill Structural Model

The focus of this research is on the third octave chatter, which is symmetric with respect to the center-plane of the strip and is only a function of the vertical vibration.⁷⁾ By neglecting the horizontal and torsional vibration, a unidirectional structure model is used as shown in Fig. 3. This model was used by many other researchers.^7,12,34)

Fig. 3.

Mill structural model.

The relationship between the displacements of the center of the work roll (y_w) and the center of the backup roll (y_b) and the rolling force variation (df_y) can be expressed by two second order differential equations:   

$$\{\begin{array}{l}{M}_w{\ddot{y}}_w+C_w({\dot{y}}_w–{\dot{y}}_b)+K_w(y_w–y_b)=wd{f}_y\\ M_b\ddot{y}{}_b+C_b{\dot{y}}_b+K_b{y}_b=C_w({\dot{y}}_w–{\dot{y}}_b)+K_w(y_w–y_b)\end{array}$$(9)

The subscripts ‘w’ and ‘b’ signify the work and backup roll respectively. M is the roll mass, C is the damping coefficient, K is the spring constant and w is the strip width. It is worth mentioning that the displacement of work roll center and the strip thickness variation at exit (dy₂) can be related as y_w = dy₂/2.

2.3. Chatter Model

The rolling vibration model is achieved by combining the rolling process model with the mill structural model. Figure 4 shows the chatter model for a single stand.

Fig. 4.

Chatter Model.

The variation of strip tensile stress at entry and exit of stand is proportional to the integral of the difference between the variations of the exit speed of a stand and the entry speed of next stand.   

$${\left(d{\sigma }_{x,1}\right)}_i=\frac{E}{L_{i–1}}{\int }_t(d{u}_{1,i}–d{u}_{2,i–1})dt$$(10)

  

$${\left(d{\sigma }_{x,\hspace{0.17em}2}\right)}_i=\frac{E}{L_i}\underset{t}{\int }(d{u}_{1,\hspace{0.17em}i+1}–d{u}_{2,i})dt$$(11)

Where i is indicative of the (i)’ th stand, E is the Young’s modulus of the strip and L_i is the distance between Stands i and i+1.

The inputs of the rolling process model are the strip thickness at entry, the strip tensile stress at entry and exit. The outputs of that are the strip thickness at exit, the strip speed at entry and exit and the rolling force. The output speeds are used to calculate the strip tensile stress. The dynamic rolling force enters to the mill structural model. It leads to variation of the roll gap spacing and this variation enters to the process model. So the relationship between the rolling process and the mill structural model continues.

3. Experimental Equipment

The experimental data of this research are from a two-stand tandem mill unit of Mobarakeh Steel Company (MSC). Figure 5 illustrates a view of this unit.

Fig. 5.

A view of two-stand tandem mill unit of Mobarakeh Steel Company.

Incoming coils to this unit are often 2 mm in thickness. Coils open up and pass through two 4-high stands two or three times. Usually thicknesses of the strip are reduced by 60 to 90 percent depending on the desired thickness. Outgoing strips from the second stand in any pass are coiled by a coiler machine (Pickup reel). The pickup reel of the previous pass plays role of the payoff reel for the current pass which is in the reverse direction of the previous pass. In practice in this mill, chatter usually occurs in second stand during the third pass.

Rolling control instrumentation in this unit consists of several analog and digital signals, all connected with a data collection system (IBA). All rolling conditions are recorded in this system and can be used online or offline. To evaluate the status of the rolling system in chatter conditions, rolling speed is increased for the moment. This is done by the operator via manual part of the control system which its online monitoring system is shown in Fig. 6.

Fig. 6.

Online monitoring system of the control system of the rolling process.

3.1. Mil Specifications

Mill modulus of a rolling mill is experimentally obtained using calibration test. In the calibration test, the deformation of the mill under rolling force is measured while the rolls are rotating but no strip material is passing through the roll gap. This deformation and the rolling force are used to calculate the mill modulus. Mill modulus is the equivalent stiffness of the contacting rolls and the stand. The deformation of the rolls is nonlinear. This nonlinearity is associated with the change in the contact area but can be linearized around the steady state rolling conditions. A finite element model is utilized to calculate the stiffness based on the contact between work roll-work roll and work roll-backup roll (K_w). Knowing both the mill modulus and the stiffness of the contacting rolls, the stand stiffness (K_b) can be calculated.

While no damping is assumed for the work rolls (C_w), a slight damping effect due to the lubricant film is present in the journal bearing of the backup rolls.^12,22) Mill defects such as eccentricity and/or ovality in the rolls enforce the mill to vibrate while strip rolling. A simulation program of the forced vibration strip rolling is utilized to calibrate the backup roll damping coefficient (C_b). Actually in this program the backup roll damping coefficient is chosen so that the result of the oscillation domain of the rolling force matches with the experiment. This value is then verified by simulating mill response in various rolling conditions. This coefficient is used in the current project for the chatter simulation. This method has been applied in the cold and hot rolling process for the first time by the authors.^1,35)

Specifications of the rolling mill under consideration are presented in Table 1. It also presents characteristics of the strip in the third pass. Inter-stand parameters are presented in Table 2. Lubricant parameters are reported in Table 3.

Table 1. Specifications of two-stand tandem mill unit.

Parameter	Stand 1	Stand 2
Work roll mass, Mw, (Kg)	14000	14000
Backup roll mass, Mb, (Kg)	38000	38000
Work roll damping coefficient, Cw, (N.s/m)	0	0
Backup roll damping coefficient, Cb, (N.s/m)	1.254e6	1.254e6
Work roll spring constant, Kw, (N/m)	4.97e10	4.97e10
Backup roll spring constant, Kb, (N/m)	1.22e10	1.22e10
Work roll radius, aw, (m)	0.245	0.245
Backup roll radius, ab, (m)	0.675	0.675
Strip width, w, (m)	0.762	0.762
Yield stress in the plane stress condition, σk, (MPa)	896	953
Entry thickness, y1, (mm)	0.409	0.317
Exit thickness, y2, (mm)	0.317	0.210

Table 2. Inter-stand parameters of the two-stand tandem mill unit.

Parameter	Payoff real-Stand 1	Stand 1-Stand 2	Stand 2-Pickup real
Distance, L, (m)	5.675	4.725	6.6
Tensile Stress, σx, (MPa)	95	147	79

Table 3. Parameters of the used lubricant.

Parameter	Value
Dynamic viscosity at 50°C, η50, (Pa.s)	0.0178
Viscosity pressure coefficient in the Barus equation, α, (GPa–1)	1.5
Viscosity pressure coefficient in the Roelands equation, z	0.5216
Limiting shear stress, τL, (MPa)	24.3

4. Simulation Algorithm

The chatter simulation model of the second stand of two-stand tandem mill unit of Mobarakeh Steel Company during the third pass has been constructed as shown in Fig. 7. The “stand 2” module represents the single-stand chatter model including the rolling process model, the mill structural model and the inter-stand model. The last deals with the transportation time delay and the inter stand tension effects. The “stand 1” and “pickup reel” module are applied in order to consider the tension variations before the entry and after the exit from the “stand 2”. Also, the “stand 1” module provides an initial excitation to the “stand 2”. The “clock” is used to trace the simulation time.

Fig. 7.

Chatter simulation program.

5. Results and Discussions

A positive equivalent damping coefficient tends to dampen vibration and to stabilize the rolling mill, while a negative one increases vibration and makes the mill unstable. On the other hand the equivalent damping coefficient is proportional to the inverse of the rolling speed.^4,22) So the fact that increasing rolling speed causes a considerable decrease of the equivalent damping coefficient partially explains why the rolling speed is the most important parameter in controlling vibration of the rolling mill.³³⁾ Also the main characteristic in a chatter simulation program is its dependency to rolling speed. To investigate this issue, fluctuation in the work roll of the second stand has been examined. The results in the two different velocities (9 and 13 m/s) are presented in Figs. 8 and 9. Figure 8 shows that vibrations are damped at low speed (9 m/s) and the system is stable in this case. With increasing speed to 13 m/s, vibrations increase with time and the system becomes unstable (Fig. 9). It is clear that the undamped vibrations make the system unstable. Thus, these figures clearly show that the program depends on the rolling speed.

Fig. 8.

Work roll fluctuation in the stable case for rolling speed of 9 m/s.

Fig. 9.

Work roll fluctuation in the unstable case for rolling speed of 13 m/s.

With trial and error, one can get the case where vibrations are undamped and vibrations amplitude versus time is unchanged. This represents the case in which the system is in the critical status. Accordingly, critical speed of the considered strip has been found to be 10.23 m/s. Work roll fluctuation of the second stand in this case is presented in Fig. 10.

Fig. 10.

Work roll fluctuation in the critical case.

Figure 11 shows the FFT spectrum of the vibration signals in the unstable case for rolling speed of 13 m/s. Dominant frequency is calculated 169.5 Hz, which is within the frequency range of third octave chatter.

Fig. 11.

FFT spectrum of the vibration signals in the unstable case.

The above procedure has been repeated for four strips rolled in the two-stand tandem mill unit of MSC. Then, the critical speeds obtained from these simulations and experiments have been compared (Fig. 12). The results show that the critical speeds of simulations are on average about 42 percent larger than the critical speeds of experiments. Much of this difference can be related to occurrence mechanism of rolling chatter in the simulations with the single stand rolling and in the experiments with the two stands rolling. In the single stand rolling, the negative damping mechanism is activated while the regenerative chatter also occurs in the two stands rolling. In the other hand, the critical speed is much higher in the negative damping mechanism than in the regenerative chatter.^33,34) So as expected, the rolling chatter in the experiments occurs much faster than that in the simulations.

Fig. 12.

Comparison of critical speeds obtained from simulations and experiments.

Using the simulation program, the effect of some main parameters of the rolling lubricant on the chatter critical speed has been identified. Figures 13, 14, and 15 show the effect of the rolling lubricant viscosity (η₅₀), the viscosity pressure coefficient in the Roelands equation (z), and the limiting shear stress (τ_L) on the chatter critical speed, respectively. The results indicate that excessive increasing of the first and the second parameter do not have considerable effect on the chatter critical speed; but the limiting shear stress variation does. This conclusion is consistent with the experimental results of Yarita et al.³⁾ This can be expressed by the fact that in the investigated conditions, the friction stress exceeds the limiting shear stress, hence the friction stress is set equal to the limiting shear stress. As a result, the effect of the rolling lubricant viscosity and viscosity pressure coefficient variations in the studied conditions on the friction stress and chatter critical speed are negligible.

Fig. 13.

Effect of rolling lubricant viscosity on chatter critical speed.

Fig. 14.

Effect of viscosity pressure coefficient on chatter critical speed.

Fig. 15.

Effect of limiting shear stress on chatter critical speed. (Online version in color.)

6. Conclusion

In this research, a new rolling chatter model by considering unsteady lubrication is presented. This chatter simulation program has been run for four strips rolled in the two-stand tandem mill unit of Mobarakeh Steel Company. The applied model has been verified by comparing the critical speeds obtained from the simulations to the experimental results. In addition, the effect of the rolling lubricant viscosity, the viscosity pressure coefficient, and the limiting shear stress on the chatter critical speed have been investigated quantitatively. The results show that under the reported operating conditions, the first and the second parameters do not have important effect on the chatter critical speed; but the last parameter is directly proportional to the chatter critical speed. Future work could focus on improving the lubrication model by formulating mixed lubrication regime, consideration the effect of surface roughness, and employing a more accurate relation for the limiting shear stress.

Acknowledgment

The authors are grateful for the assistance of the Mobarakeh Steel Company.

References

• 1)   M. R.  Niroomand,  M. R.  Forouzan,  M.  Salimi and  H.  Shojaei: ISIJ Int., 52 (2012), 2245.
• 2)   P.-A.  Hu,  H.  Zhao and  K. F.  Ehmann: Proc. Inst. Mech. Eng., Part B: J. Eng. Manuf., 220 (2006), 1267.
• 3)   I.  Yarita,  K.  Furukawa,  Y.  Seino,  T.  Takimoto,  Y.  Nakazato and  K.  Nakagawa: Trans. Iron Steel Inst. Jpn., 18 (1978), 1.
• 4)   Y.  Kimura,  Y.  Sodani,  N.  Nishiura,  N.  Ikeuchi and  Y.  Mihara: ISIJ Int., 43 (2003), 77.
• 5)   J.  Zhong,  H.  Yan,  J.  Duan,  L.  Xu,  W.  Wang and  P.  Chen: J. Mater. Process. Technol., 120 (2002), 275.
• 6)   R.  Arnken,  C.  Bilgen and  W.  Hennig: Iron Steel Technol., 4 (2007), 81.
• 7)   P. A.  Meehan: J. Vib. Acoust., Trans. ASME, 124 (2002), 221.
• 8)   X.  Yang,  C.-n.  Tong,  G.-f.  Yue and  J.-j.  Meng: J. Iron Steel Res. Int., 17 (2010), 30.
• 9)   Y.  Kim,  H.  Park,  S. S.  Lee and  C. W.  Kim: ISIJ Int., 52 (2012), 1135.
• 10)   A.  Heidari and  M. R.  Forouzan: J. Adv. Res., 4 (2013), 27.
• 11)   R. E.  Johnson: Proc. Royal Soc. London, Series A: Math. Phys. Sci., 445 (1994), 479.
• 12)   R. E.  Johnson and  Q.  Qi: Int. J. Mech. Sci., 36 (1994), 617.
• 13)   Y.  Chen,  S.  Liu,  T.  Shi,  S.  Yang and  G.  Liao: Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. Sci., 216 (2002), 1224.
• 14)   Y.  Chen,  T.  Shi and  S.  Yang: Key Eng. Mater., 245 (2003), 123.
• 15)   A. C.-P.  Tung: Ph.D. 9832705, Northwestern University, Ill, USA, (1998).
• 16)   I. S.  Yun,  W. R. D.  Wilson and  K. F.  Ehmann: J. Manuf. Sci. Eng., 120 (1998), 330.
• 17)   I. S.  Yun,  K. F.  Ehmann, and  W. R. D.  Wilson: J. Manuf. Sci. Eng., Trans. ASME, 120 (1998), 337.
• 18)   I. S.  Yun,  K. F.  Ehmann and  W. R. D.  Wilson: J. Manuf. Sci. Eng., Trans. ASME, 120 (1998), 343.
• 19)   I. S.  Yun,  W. R. D.  Wilson and  K. F.  Ehmann: Int. J. Mach. Tool. Manuf., 38, (1998), 1499.
• 20)   P.-H.  Hu and  K. F.  Ehmann: Int. J. Mach. Tool. Manuf., 40 (2000), 1.
• 21)   P.-H.  Hu and  K. F.  Ehmann: Int. J. Mach. Tool. Manuf., 40 (2000), 21.
• 22)   P.-H.  Hu,  H.  Zhao and  K. F.  Ehmann: Proc. Inst. Mech. Eng., Part B: J. Eng. Manuf., 220 (2006), 1279.
• 23)   P.-H.  Hu,  H.  Zhao and  K. F.  Ehmann: Proc. Inst. Mech. Eng., Part B: J. Eng. Manuf., 220 (2006), 1293.
• 24)   Y.  Kim,  C. W.  Kim,  S. J.  Lee and  H.  Park: ISIJ Int., 52 (2012), 2042.
• 25)   W. R. D.  Wilson: J. Appl. Metalwork., 1 (1979), 7.
• 26)   M.  Saniei and  M.  Salimi: J. Mater. Process. Technol., vol. 177 (2006), 575.
• 27)   P. B.  Kosasih and  A. K.  Tieu: Tribol. Int., 40 (2007), 709.
• 28)   S.-W.  Lo,  T.-C.  Yang,  Y.-A.  Cian and  K.-C.  Huang: J. Tribol., 132 (2010), 011801.
• 29)   C.  Barus: Am. J. Sci., 45 (1893), 87.
• 30)   C. J. A.  Roelands: PhD Thesis, Delft University of Technology, The Netherlands, (1966).
• 31)   Y.  Zhang: Ind. Lubr. Ttribol., 58 (2006), 72.
• 32)   S.  Bair: High-pressure Rheology for Quantitative Elastohydrodynamics, Elsevier Science, Amsterdam, (2007).
• 33)   H.  Zhao: PhD. 3303624, Northwestern University, Ill, USA, (2008).
• 34)   H.  Zhao and  K. F.  Ehmann: ASME Int. Manufacturing Science and Engineering Conf., New York, USA, (2009), 627.
• 35)   Y.  Jahangardy,  M. R.  Forouzan and  M.  Rostamipour: Usak Univ. J. Mater. Sci., 1 (2012), 173.