that thermal effects due to friction in the central zone of the contact play a role in ... there is no slip between the fluid and the solid boundaries ... Ω , it is assumed that considering the half-cylinder is sufficient to ensure accurate deformation ... Concerning kinematics conditions, the cylinder rolls with an angular velocity 2 Ï.
Computational Fluid Dynamics and Full Elasticity Model for Sliding Line Thermal Elastohydrodynamic Contacts DOI: 10.1016/j.triboint.2011.04.013
V. Bruyere a, N. Fillot a , G.E. Morales-Espejel b,a , P. Vergne a a
Université de Lyon, CNRS INSA-Lyon, LaMCoS UMR5259, F-69621, Villeurbanne, France b
SKF Engineering and Research Centre 3430 DT, Nieuwegein, The Netherlands
Abstract: Classically, the EHD problem is solved using the Reynolds assumptions to model the fluid behaviour, and the Boussinesq elastic deformation equation to model the solid response, both being coupled with the load balance equation. The development of an alternative approach is presented here in order to solve at once the Navier-Stokes equations (mass conservation and momentum equilibrium), the full elasticity and energy equations for the line EHD problem in a fluid-structure interaction approach. The Finite Element Method is used to solve the mathematical formulation in a fully coupled way, inspired from [1]. After linearization with a Newton procedure, all the physical quantities (pressure, velocity field, deformations, temperature), are solved together in a unique system. An important benefit of this approach is the possibility to implement in a simple manner the non-Newtonian and thermal effects; in fact all the quantities can vary through the film thickness. The extension to non-Newtonian rheology, pressure and temperature dependence for the viscosity and density are taken into account in a direct way to allow an acceptable prediction of the friction coefficient. Gradients across the film thickness and temperature fields in both the fluid and the two solids are naturally computed and analysed. As a case study, we focus firstly on the pure sliding cylinder-on-plane contact. It is shown that thermal effects due to friction in the central zone of the contact play a role in heating the lubricant at the inlet zone, via the heat conduction in the solids. By increasing the Slide-to-Roll Ratio (SRR), the occurrence of dimples and the subsequent effects in the different parts of the contact under Zero Entrainment Velocity conditions are then studied.
Keywords: Elastohydrodynamic Lubrication, Fluid mechanics, Thermal effects, Dimple 1. Introduction The solution of the EHD problem requires the coupled resolution of the equations that describe the lubricant flow and the elastic deformation of the contacting surfaces. Classically, the Reynolds equation is used to model the fluid behaviour. This equation gives proper results for the film thickness prediction, but requires many assumptions: 1
− the film thickness is very thin compared to other contact dimensions − the lubricant is Newtonian − the fluid flow is laminar − the lubricant boundary surfaces are parallel or at a small angle with respect to each other − there is no slip between the fluid and the solid boundaries − inertial and body forces are neglected − viscosity and density are constant across the film thickness In the 1960's, the first numerical line contact solution of Dowson and Higginson [2] was published and followed by many others. During the past decades, the EHD problem was solved using powerful methods based on finite differences, with multigrid techniques [3], or finite elements methods [4], [5]. The film thickness was accurately predicted, whereas under the assumption of Newtonian fluid, the friction coefficient was not. In fact, the shear rate viscosity dependence (non-Newtonian behaviour) needed to be considered, as suggested by Dowson [2]. Moreover under high speed or high sliding conditions, thermal effects become significant in EHD contacts [6], [7]. These effects can dramatically affect the film thickness and are essential for predicting the friction coefficient. Several "generalized" Reynolds equations have been developed [8], [9] to account for both thermal and non-Newtonian effects. This was achieved by carrying out different integrations at any point of the domain via the modified Reynolds equation. A different approach is to consider the original fluid flow formulation, expressed through the NavierStokes equations which are basically the fundamental fluid dynamics equations. This approach can be pursued now due to the progress made in Computational Fluid Dynamics (CFD) and the gradual increase in computational speed. Since the 1990's, CFD has been used to solve hydrodynamic problems [10], [11] in different types of bearings. In the last ten years, the EHD problem has been successfully solved by a CFD approach coupled together with the classical half-space elasticity formulation for the solid deformations, e.g. Schäfer et al. [12] and Yiping et al. [13] found agreement with classical Reynolds solutions. Almqvist et al. [14] and more recently Hartinger et al. [15] built thermal EHD line contact models. The latter used the finite volume method to model the fluid behaviour and studied the combined effects of temperature, pressure and shear thinning. They confirmed that these effects are essential to predict film thickness and friction coefficient and this can be related to the occurrence of significant gradients that appear across the contact film thickness. However in these previous works, classical Boussinesq elastic deformation and semi-infinite body temperature equations were still used. The purpose of the present paper is to describe a recently developed full fluid-structure model using Navier-Stokes, full solid mechanics, and the energy equations for the fluid and solids, solved together in a fully coupled way. Real lubricated contacts are defined by two solids of finite and given dimensions, separated by a lubricant fluid film. The three bodies are discretized with appropriate mesh size and 2
density according to the possible occurrence of gradients within the domain. Results are compared with previous literature results based on Reynolds and Navier-Stokes equations. Then, this new full fluid-structure formulation is used to study different cases where thermal effects are no longer negligible, i.e. under very high sliding conditions (pure sliding, zero entrainment velocity, etc.). 2. Physical Model and Fluid-Structure Equations In this section, the different physical phenomena and their implementation in the fully coupled model are described. Only the case of the line contact problem is considered in the present work. Both contacting solids are defined by the domains Ω1 and Ω 2 (Figure 1). The fluid is represented by the domain Ω 3 . It is assumed that the contact is fully flooded and that the fluid simulation domain is large enough to specify proper boundary conditions. The inlet and outlet abscissas are located at a distance from the contact centre of l = 6b with b is the Hertzian contact half-width along the rolling direction. This symmetry is required to compute cases where the fluid can come from the left or from the right of the contact (see section 4.2). Kaneta et al. [7] and Wang et al. [16] suggested that in most cases, a depth D of 3.15b for the contacting solids is enough to ensure zero temperature gradient in regions far away from the contact line. For the domain Ω 2 , it is assumed that considering the half-cylinder is sufficient to ensure accurate deformation and thermal results. Concerning kinematics conditions, the cylinder rolls with an angular velocity ω2 and the plane moves with a linear velocity u1 in the x direction. The velocity on the boundary ∂Γ3−2 at x = 0 is defined as u2 = ω2 R , the mean velocity is um = SRR =
u1 + u2 and the Slide-to-Roll Ratio is noted 2
u1 − u2 . The classical film thickness h is defined as the gap between the two solids, and ∆ is um
the displacement imposed to the solid Ω1 to obtain the normal load. The typical operating parameters used in this article are all given in Appendix A.
2.1. Fluid 2.1.1. Governing fluid dynamics equations Navier-Stokes equations, obtained from the mass conservation and momentum equilibrium equations apply in their stationary form on Ω 3 :
∇ ⋅ (ρ3u 3 ) = 0
(1)
ρ 3 (u 3 ⋅ ∇u 3 ) = ∇ ⋅ σ f (2) 3
where ρ3 is the fluid density, u 3 is the velocity vector and σ f is the stress tensor. The gravity term is neglected using Froude number analysis ( Fr =
are taken into consideration, the Reynolds number Re =
um >> 1 ). However, inertia effects gh
ρ3um h is at maximum equal to 100 , so the flow η
is considered as laminar. The constitutive law of the fluid behaviour writes:
(
σ f = − p + η∇ ⋅ u 3 I + η ∇u 3 + (∇u 3 ) 2 3
T
)
(3)
where p is the fluid pressure and η its viscosity. The viscosity can thus vary across the film thickness without additional requirements. It is also assumed here that the fluid sticks to the solid surfaces. Furthermore, the compressible form of the equations (eq.1) and (eq.2) is used. Consequently density may vary with pressure. Knowing that ∂Γ3−in and ∂Γ3−out are far enough from the contact, the pressure is set at the reference pressure and a zero normal velocity gradient is imposed on the inlet and outlet boundaries. For interface boundaries ∂Γ3−1 and ∂Γ3− 2 , velocities are imposed according to the SRR values while entrainment velocity remains constant (see Appendix A). 2.1.2. Rheology An essential issue concerning the EHD problem is the viscosity variation. A classical dependency of viscosity with pressure and temperature is the model inspired from Roelands [17] and further developed by Houpert [18]: z0 p T − 138 * (4) η N ( p, T ) = η 0 exp (ln(η 0 ) + 9.67 ) − 1 + 1 + δ ( T T ) − − 0 8 1 . 98 ⋅ 10 T − 138 0
S0 β (T0 − 138) α ⋅1.98 ⋅ 108 p with δ * = (ln(η 0 ) + 9.67 )1 + , z = , S0 = , with T the 0 8 ln (η 0 ) + 9.67 ln (η 0 ) + 9.67 1.98 ⋅10 T0 − 138 z0
temperature, α the pressure viscosity index, β the thermoviscosity constant and η 0 and T0 the viscosity and temperature references. Under severe conditions like those that can occur in EHL, the fluid behaviour can no longer be considered as Newtonian [19]. In fact, the apparent viscosity decreases as the shear stress increases, thus the corresponding film thickness decreases. In the literature, different constitutive laws were used, such as Carreau [20] or the non-linear Maxwell model [21], [22]. Bair and Winer [19] introduced a different type of non-Newtonian behaviour based on the limiting-shear-stress concept. The aim of the present work is to provide a general frame that allows the use of a wide variety of constitutive laws and not to compare quantitatively our results with experimental ones. For these reasons, the Ree-Eyring 4
sinh-law [23] is used here as a first example to model non-Newtonian effects: γ&eq =
τ eq τ0 sinh with ηN τ0
τ 0 the Ree-Eyring stress. To also take into account viscosity variations with pressure and temperature, the complete viscosity law is then defined as η ( p, T , γ&eq ) =
η N ⋅ γ&eq τ0 . sinh −1 γ&eq τ 0
In the Navier-Stokes approach, velocity components and obviously the shear rate tensor are directly
(
)
T solved. This tensor is then defined by γ& = ∇u 3 + (∇u 3 ) . According to [24], the viscous strain rate is
governed by the second tensor invariant. An equivalent shear rate can then be defined as γ&eq = as well as an equivalent viscosity ηeq =
1 γ& ⊗ γ& 2
f (γ&eq ,η N ) . The function f is suitable for any generalized γ&eq
rheological model. In this study, all the cases were computed with Ree-Eyring and modified Roelands equations because they were used in [15], i.e. in the reference chosen to validate our approach. 2.1.3. Density and cavitation model The model used for density variations with pressure and temperature is the Dowson-Higginson law [2]:
5.9 ⋅108 + 1.34 p − β DH (T − T0 ) 8 5.9 ⋅10 + p
ρ3 ( p, T ) = ρ0
(5)
with β DH the density-temperature coefficient and ρ0 the reference density. This law is obviously valid only for positive pressures. But it is well known that in the contact outlet, the pressure drops dramatically and thus cavitation phenomena may take place. In the Reynolds approach, this issue is handled by cancelling the negative pressures and the pressure gradient in the exit region. In the Navier-Stokes approach, different methods have been employed to eliminate these negative pressures. The most rigorous way would be to develop a multiphase solution, obtaining in this way information regarding what happens in the lubricant and the position of the free surface. This method is still in development and offers possibilities for future works. An alternative "hybrid" method was used by [15] to solve the cavitation area, applying specific operations. However if implemented in our approach, it would considerably delay the present model. Therefore the basic formulation used in this approach is finally inspired from [14]. A simple equation of state is defined and constrained by assuming the continuity of the density pressure-temperature law (eq.5) over the exit area, and neglecting temperature variation. If the pressure drops under the cavitation pressure pcav , the liquid is converted to vapour and the density of the resulting mixture falls significantly. In agreement with [14], a second order polynomial (eq.6) is then used to extrapolate the 5
density law to negative relative pressures. In order to ensure the convergence of this scheme, a further iterative computation is required. The value of p for which ρ3 = ρair , peps , starts from negative pressure and it gradually tends to the cavitation pressure.
(eq.5) if p ≥ 0 ρ − ρ 1 + 5.76 ⋅10 −10 p (6) ρ 3 ( p, T ) = air eps 0 2 −10 p + 5.76 ⋅10 p + 1 if p < 0 ρ 2 0 ρ 0 peps
(
)
2.2. Solid mechanics equations In order to simplify the problem, only the cylinder Ω 2 is considered as elastic, the solid Ω1 being considered as rigid. The linear elasticity problem consists in finding the displacement vector U on the domain Ω 2 such that:
∇ ⋅σ s = 0
(7)
(
T with σ s = Cε s (U ) , ε s = ∇U + (∇ U )
) the strain tensor,
C the stiffness tensor containing material
properties (Young's modulus and Poisson's ratio), U = {U , V }, U and V the X − and Y − displacements respectively. Since only the line contact case is studied here, the plane strain assumption applies and this simplifies the linear elastic problem. Concerning the boundary conditions, the upper surface ∂Γ2 , is considered as stationary ( U = V = 0 ) and the stress continuity is assumed at the fluid-solid interface ∂Γ3− 2 . The remaining boundaries are assumed to be free ( σ s ⋅ n = 0 ). Note that the linear elasticity assumption is used because the material and the applied forces studied in the present work lead to moderate strains and stresses. However, more complex solid behaviour laws can easily be implemented in this general method.
2.3. Energy conservation The need to incorporate thermal effects was recognized a long time ago, for instance by Cheng [6], and Dowson [2]. A first full numerical solution for the point contact problem was proposed by Zhu and Wen [25] followed by many authors, like Kim et al. [26] who included non-Newtonian effects. When entrainment velocity um and/or SRR increase(s), thermal effects play a major role on the viscosity evolution in the high pressure zone, hence on the friction coefficient. Under severe conditions, the film thickness can also dramatically be affected. The energy conservation principle is applied in the domains [Ωi ] ; it is expressed as a function of the temperature by considering that the heat capacities at constant pressure C p i and the thermal conductivities ki of domain i , are constant:
6
ρ i C p i (u i ⋅ ∇T ) = ∇ ⋅ (k i ∇T ) + Qi
(8)
with i = 1, 2 for the solids and 3 for the fluid, u i the velocity vectors of domain i and Qi the source terms. Concerning the source terms, there is no heat source in the solids, Q1 = Q2 = 0 . But in the fluid domain, the lubricant is subjected to high shear which causes viscous heating and to high compression which heats up the lubricant by the pressure work, Q3 = Qshear + Qcompr , expressed respectively by: 2 T 2 Qshear = η ∇u 3 ⊗ ∇u 3 + ∇u 3 ⊗ (∇u 3 ) − (∇ ⋅ u 3 ) 3 Qcompr = −
T ∂ρ3 ρ3 ∂T
(u3 ⋅ ∇p )
(9)
(10)
p
According to the Peclet number, Pei =
ρiui LiC p i ki
, ( Li being the characteristic length of domain i ) and
with the solids dimensions considered here, one has Pe1 ≈ 50u1 and Pe2 ≈ 750u2 . For classical velocity values, the convection mode clearly dominates the conduction one. But, for the pure sliding case, where one of the solids is stationary, the conduction effects are obviously important and must be taken into account, as it will be shown later. Temperature is imposed to a reference temperature T0 on the inlet and outlet fluid boundaries ∂Γ3−in and ∂Γ3− out (only for incoming fluid) and on the solid boundaries ∂Γ1 and ∂Γ2 . At the solid-liquid interfaces, ∂Γ3−1 and ∂Γ3− 2 , classical continuity conditions are imposed: the heat flux in the normal direction is assumed to be continuous across the boundary. At the remaining solid boundaries, the conductive heat fluxes are set to zero. Mathematically, the thermal boundary conditions are therefore defined as: T ∂Γ
= T0 if u > 0
T ∂Γ
= T0 if u < 0
3− in
3− out
T ∂Γ = T ∂Γ = T0 1
(11)
2
n ⋅ (− ki ∇T ) ∂Γ = n ⋅ (− k3∇T ) ∂Γ 3−i
3−i
for i = 1, 2
and for the remaining solid boundaries, n ⋅ (− ki ∇T ) = 0 for i = 1, 2
2.4. Load Balance To finalize the EHD problem formulation, the classical load balance equation must be included. With the present complete fluid-structure interaction approach, the mesh displacement ∆ enforces the load by pressing the solids against each other. To solve the imposed load problem and to ensure the real coupling, an extra equation is added to the problem: 7
∫σ
f
⋅ ndΓ − W = 0
(12)
∂Γ3−2
with W being the imposed normal load.
2.5. Moving Mesh After the description of the different physical mechanisms implemented in the model, the links are established by the use of a moving mesh to ensure the full fluid-structure coupling. The EHD problem is characterized by the fact that deformations of the solids are of the same order of magnitude (or larger) than the film thickness. Numerically, the solid deformations specify the fluid domain and the computational size of the problem. The classical method to express the fluid equations is by using the Eulerian description, which sets a fluid domain observing the moving fluid particles. Concerning the solid domain equations, the Lagrangian representation is classically used in which each node of the computational mesh follows the associated material particle during motion. In order to keep the classical form of each equation and to combine the best features of both approaches, the Arbitrary-Lagrangian-Eulerian (ALE) method is used to ensure the coupling [27]. Two coordinate systems are then defined: •
a fixed coordinate system, with the spatial coordinates ( X , Y )
•
a moving coordinate system, with the spatial coordinates ( x, y ) .
The fluid behaviour and energy conservation equations are written in the moving coordinate system and the solid equations in the fixed one. The gap is imposed by a vertical mesh displacement ∆ for all the nodes of the solid Ω1 . As for solid
Ω 2 , mesh displacements are set to the displacement values calculated with the solid deformation equations (eq.7). In the fluid domain, each node moves with a Laplace smoothing [27] driven by the imposed boundary conditions. For the lateral boundaries ∂Γ3−in and ∂Γ3− out , nodes move with a linear interpolation between the imposed displacement ∆ of solid Ω1 and the calculated ones for Ω 2 . To summarize the method: the cylinder is deformed by the high pressure in the fluid and its mesh follows these deformations. The fluid domain boundaries follow the solid deformation, implying therefore changes of the pressure and velocity fields. In order to properly take into account these fully coupled interactions, a specific numerical method is then required.
3. Numerical Procedure and Model Validation 3.1. Numerical Procedure The Finite Element Method is used to solve all the equations in their respective domains. Each variable is approximated by quadratic Lagrangian polynomials except for pressure for which linear ones are used 8
(one order lower than for the velocity field). Concerning the fluid domain, each node has then six degrees of freedom ( u, v, p, x, y and T ) and concerning the solids each node has three for the cylinder ( U , V , T ) and one for the plane ( T ) which is considered as rigid. Triangular elements were chosen to match the geometry and calculations (not detailed here) showed that at least six elements across the film thickness have to be considered. An advantage of this method is that all the equations are combined into a single non-linear system ensuring the full coupling (see Figure 2). The system is then linearized using a Newton-Raphson procedure. The variables are initialized from analytical formulae: Hertzian pressure distribution ( p ), analytical Reynolds formula for the horizontal velocity ( u ) and the Ertel-Grubin central film thickness [28] for the initial gap between the two solids ( ∆ ). The first resolution step (see Figure 3) is linear and consists to pre-load the upper solid (eq.7) and to initialize the fluid simulation domain with the moving mesh. All the equations are then solved in the same global loop, isothermally first and then with thermal effects. The load balance equation (eq.12), is then added to the global matrix and corresponds to the supplementary degree of freedom ∆ (Figure 2). From this solution, the density law is modified (eq.6), and iterations are accomplished to obtain the value of peps . Each iteration of the main loop takes approximately 10 minutes on a bi-processor Intel Xeon 2.27Ghz and 32Go RAM computer with typically 61000 elements for the fluid Ω3 , 33000 for the cylinder Ω 2 and 30000 for the plane Ω1 .
3.2. Validation with previous literature results In this section, a comparison is performed with some results obtained previously in the literature. Pressure and film thickness results of the present model are compared firstly with classical results obtained with a Reynolds model inspired from Habchi [29] (see Figure 4-a). A classical Newtonian isothermal case is computed with L = 5.6 and M = 21.9 (dimensionless Moes parameters), and a fairly good agreement is found. Secondly, in order to validate our approach with more representative cases, results are compared with those obtained by Hartinger et al. [15] using CFD to model the fluid behaviour (see Figure 4-b). Results are compared on a non-Newtonian thermal case computed with SRR = 1 , and again a satisfying agreement is found. The maximum film thickness deviation between the two series of results is lower than 2% (considering that the data were directly taken from figure 6-B in [15]). The model can now be used to study different sets of operating conditions.
9
4. Case Studies For classical EHD cases, the developed method can be used to check the validity of the Reynolds assumptions in the contact area. Moreover, there are some cases where this approach can provide interesting information for further understanding of highly-coupled EHD mechanisms. In the literature, little information was found concerning the thermal friction coefficient prediction at high SRR. Indeed, it is difficult to choose an appropriate fluid behaviour law to predict both friction and film thickness in the same model, especially under highly sliding operating conditions [30]. With this purpose, the evolution of the friction coefficient as a function of SRR (Figure 5) is firstly studied (input data are given in Appendix A). The friction coefficient is defined by C f (eq.13) calculated on the interface boundaries ∂Γ3−1 and ∂Γ3− 2 : l
1 C f = ∫ηγ& dx (13) W −l In Figure 5, the results are compared to those obtained by the "Reynolds approach" inspired from Habchi [29]. Considering isothermal behaviour, a good agreement is found and the friction coefficient increase is limited at high SRR by the fluid non-Newtonian behaviour. When thermal effects are taken into account in both models, the friction globally decreases with SRR, as normally expected [28]. It should be noticed that the two models predict a minimum friction for SRR = 2 (this particular behaviour is investigated in 4.1). Similar tendencies can be found in [31]. Three cases will be especially discussed: − the pure sliding case ( SRR = 2 ), which corresponds to the minimum of the friction curve (Figure 5) − an extreme case ( SRR = ∞ , with no entrainment velocity), where the surface velocities are equal but in opposite direction. − intermediate cases where SRR is high ( SRR > 2 ), to study the main variables evolution.
4.1. Pure sliding case ( SRR = 2 ) The pure sliding case is considered and solid Ω 2 is kept stationary. As shown in Figure 5, when thermal effects are not considered, the friction coefficient continuously increases for an increase of SRR. A different situation occurs when thermal effects are taken into account: a minimum in friction is found. For this case, the resulting temperature field in the solids and in the fluid is plotted in Figure 6. Temperature continuity between the fluid and the solids can be observed. The results also show that the stationary cylinder heats the fluid by conduction at the inlet of the contact. Heat is then transferred from the centre of the contact (where it is generated) towards the inlet via conduction in the solid. The hot stationary solid heats up the lubricant just before it enters the contact, modifying the lubricant viscosity in this area. Therefore, the quantity of lubricant dragged into the contact is reduced. 10
Different cases with SRR close to 2 are now analyzes to explain the particular friction behaviour around this value. For SRR = 1.9 , the surface velocities are oriented in the same direction and the contact "inlet" temperature is classically close to the reference one, set on the left boundary of the contact (Figure 7); advection drags the temperature out of the contact for both solids. For the case with SRR = 2 , the lubricant is heated all along the contact (Figure 7) by conduction provided by the heating of the stationary solid. The central film thickness and the friction coefficient decrease. Finally, for the case SRR = 2.1 , the surfaces move with opposite directions ( u1 = 5.125m ⋅ s −1 and
u2 = −0.125m ⋅ s −1 ). There are actually two inlets, one for each moving surface. On the left hand side of Figure 7 (main inlet) the temperature on the surface of Ω 2 is still higher due to the same conduction effects as those encountered in the case SRR = 2 . Even with more severe sliding conditions for
SRR = 2.1 than for SRR = 2 , the higher advection effect on the cylinder (since u2 increases) lowers the fluid temperature at the right hand side which increases the lubricant viscosity and produces higher friction than for the case SRR = 2 . This explains the local minimum of friction for SRR = 2 . After having studied the importance of thermal coupling between solids and fluid, an extreme case with velocities in opposite direction is then presented.
4.2. Zero Entrainment Velocity case ( SRR = ∞ , u1 = −u2 ) When the SRR increases, the typical film thickness shape with a flat plateau in the contact centre is not observed any more: Chiu [32] and Kaneta et al. [33, 34] found a local film thickness increase so-called "dimple" in this central area. An extreme case of nominal zero entrainment speed that can be encountered for example in cage-free ball bearings, is investigated here to better understand the dimple formation. With the same material for both solids, SRR is increased until u1 = −u2 (see Figure 8). The Reynolds equation predicts that no film thickness is generated since the entrainment speed is zero. But as observed by Dyson et al. [35], a film thickness is nevertheless well established for these conditions. Many explanations have been proposed and the most commonly used is “the viscosity wedge theory” originally postulated by Cameron [36]. In agreement with [37], the approach developed here gives the evidence of the existence of a temperature gradient through the film thickness (see Figure 9-a) and therefore a viscosity gradient. Indeed, due to advection in the opposite moving solids, the temperature at each inlet is different on both solid surfaces. In the case of equal but opposite solid velocities, the inlet and outlet definitions as we normally understand them in a Reynolds model in EHL are totally modified. But rather two inlets and two outlets
11
are created, one for each surface. This would certainly modify the classical film thickness and pressure generation. Using complete Navier-Stokes equations to obtain more precise information on the velocity field through the film thickness, the horizontal pressure gradient is obtained from:
∂u ∂u ∂p ∂ 2 ∂u ∂v ∂ ∂u ∂ ∂u ∂ ∂v + v = − − η + + 2η + η + η ∂y ∂x ∂x 3 ∂x ∂y ∂x ∂x ∂y ∂y ∂y ∂x ∂x
ρ 3 u
(14)
By analyzing the order of magnitude of each term, the one including viscosity variations across film thickness appears to be the main responsible of this pressure gradient: ∂ ∂u ∂ 2u ∂η ∂u η = η 2 + ⋅ ∂y ∂y ∂y ∂y ∂y The term
(15)
∂η ∂u , caused by the temperature difference at the moving walls, is represented in Figure ⋅ ∂y ∂y
8-a with the resulting pressure gradient in Figure 8-b. Superposed to the classical term η
∂ 2u , it creates ∂y 2
the pressure gradient and consequently, this particular film thickness profile. Using the present full fluid-structure approach, the flow streamlines are then plotted together with the temperature distribution in the fluid and in the solids in Figure 9-a and velocity field vectors in Figure 9-b. At each inlet or outlet ( x ≈ ±1⋅10−4 m ), the cold fluid is dragged into the contact by one solid and is also ejected out of the contact after heating by the other one. Between these two fluid layers, a stationary recirculation area is observed in the central part of the contact where the maximum temperature is reached (see Figure 9-a). Zero Entrainment Velocity conditions have been experimentally studied by Yagi et al. [38] on a ballon-disk contact with different material and lubricant properties that those used in this work. Nevertheless, a qualitative comparison is proposed showing the influence of the normal load and the sliding velocity ( u * = u1 ) on film thickness. When the load is increased (Figure 10-a), the dimple depth increases, the minimum film thickness remains constant and the contact width increases. When sliding velocity becomes greater (see Figure 10-b), the minimum film thickness grows and the dimple depth increases up to a critical value ( u* = 3m ⋅ s −1 ) and it then decreases. Indeed, as shown previously with the present operating conditions, the film thickness shape is dominated by thermal effects. By increasing the opposite solid velocities, the total input energy in the model is increased. Temperature difference on both fluid-solid interfaces at each inlet is then augmented due to differences in the convection rate between the surfaces produced by the higher difference of speeds, leading to a minimum film thickness increase. Furthermore, up to a critical velocity value, there is a competition between non-Newtonian and thermal effects. The shear rate shape is modified and a 12
high shear rate band is observed in the middle of the film modifying the viscosity wedge. The maximum pressure gradient is then reduced and headed outside of the contact (see Figure 10-b). The dimple depth consequently decreases. Note that similar tendencies have been experimentally found by Yagi et al. [38] concerning load and sliding velocity influence on film thickness profiles. Once the importance of the viscosity wedge on this extreme case is shown, highly SRR cases are computed to understand the dimple behaviour in less extreme cases.
4.3. Intermediate cases ( SRR > 2 ) As described previously, it is believed that the main explanation for the dimple creation is the "viscosity wedge" mechanism (e.g. Cameron [37]). With the present approach, it is possible to investigate thermal effects in the fluid and in the solids. Starting from pure sliding case, SRR is increased and dimples are observed (Figure 11). Indeed, as shown for the case SRR = 2.1 , on the secondary right hand side inlet a temperature gradient is observed and it increases with SRR. The quantity of lubricant entering the contact by the right-hand side is then increased, and the temperature gradient too. The viscosity wedge effect is amplified and the film thickness at x = 1 ⋅ 10−4 m is increased (Figure 11). The dimple moves towards the contact centre and it enlarges. Classically with non-Newtonian and thermal effects, the pressure peak is reduced when the SRR is increased. As shown in the previous extreme case, due to the viscosity wedge effect a pressure gradient is generated, modifying the pressure distribution and increasing the maximum pressure value that moves towards the contact centre. When SRR still increases, pressure and film thickness profiles tend to those found in the case of Zero Entrainment Velocity (see 4.2). Viscosity and temperature variations through film thickness are now investigated. Viscosity variations for SRR = 4 are plotted in Figure 12-a with flow streamlines. The viscosity maximum is reached where the dimple depth is maximum, near the cylinder boundary ∂Γ3−2 . The temperature (Figure 12-c) is low in this region due to higher advection of the cylinder; the pressure is maximum (Figure 11) and a resulting low shear rate band is observed in the fluid in Figure 12-b. Maximum temperature is reached in the middle of the contact, leading to lower viscosity values and to the occurrence of a high shear rate band. As in the previous extreme case, a recirculation area, where the fluid is trapped inside the contact, is generated by the opposite velocities (Figure 12-a). The recirculation centre is located at the maximum dimple depth abscissa. Horizontal velocity variations at different sections are then investigated Figure 13. Where viscosity variations are strong, velocity profiles are obviously affected. Viscosity is high near the lower boundary
∂Γ3−1 , and thus the lubricant is clearly more easily driven by this moving boundary. The particular 13
velocity profile in
x = 0.5 is thus explained by the high viscosity values near each one of the fluid-solid b
interfaces.
5. Conclusion An original approach has been developed in order to give physical insight for different cases in TEHL of line sliding contacts. Navier-Stokes equations, full elasticity formulation and full energy conservation have been used in order to reduce the number of assumptions in the formulation. The results have been successfully compared with classical Reynolds results and with CFD methods from literature applied to EHL. From the above, the following conclusions can be outlined: − By adding the fluid thickness dimension to the EHL problem and by modelling solids with finite dimension, our classical vision of the contact has been enlarged. Flow streamlines and heat transfer analysis provide more information on the EHL couplings understanding. − The thermal effects on friction have been investigated and the effect of conduction in the stationary solid explained the local minimum friction found in a pure sliding case. − High SRR cases have been computed, showing thermal effects preponderance. The "dimple" formation and shape were explained with the viscosity wedge theory by analyzing the gradients of viscosity, pressure and temperature in both the fluid thickness and in the solids. − Results obtained under Zero Entrainment Velocity conditions have been qualitatively compared with experimental ones published in the literature and similar trends were found. The present fluid-structure interaction model appears to be an appropriate tool to explain intriguing experimental results. Without too much extra work it could be possible to implement slippage effects, as studied by Ehret et al. [39]. With the accurate velocity field description offered by the model described in this work, it should be interesting to study the interaction between slippage and thermal effects and their consequences on dimple generation.
Acknowledgements The authors want to thank Mr. Thomas Doki-Thonon, for his help in the implementation of the Reynolds model. The third author wants to thank Mr. Alexander de Vries, Director SKF Group Product Development, for his kind permission to publish this article.
14
Appendix A Symbol [unit]
Value
R [m]
0.01
η0 [Pa ⋅ s ]
0.01
ρ 0 [kg ⋅ m −3 ]
850
p0 [Pa]
0
pcav [Pa]
− 1 ⋅ 105
um =
[
u1 + u 2 m ⋅ s −1 2
]
E [Pa]
2 ⋅ 1011
υ
0.3
α [Pa −1 ]
1.76 ⋅ 10−8
[
W N ⋅ m−1
]
1 ⋅ 105
τ 0 [Pa]
5 ⋅ 106
T0 [K ]
353
β [K −1 ]
0.0476
β DH [K −1 ]
[
k 3 W ⋅ (m ⋅ K )
−1
[
]
C p3 J ⋅ (kg ⋅ K )
[
−1
k i W ⋅ (m ⋅ K )
−1
[
C pi J ⋅ (kg ⋅ K )
]
7 ⋅ 10−4 0.15
]
−1
]
Nomenclature Roman Symbols
b [m]
Hertzian contact half-width
C [Pa]
Stiffness tensor −1
]
450 7850
Table 1 - Parameters properties used in the model
[
2300 47
ρ i [kg ⋅ m −3 ]
C pi J ⋅ (kg ⋅ K )
2.5
Heat Capacity of domain Ω i 15
D [m]
Depth of solid Ω1
Ei [Pa]
Young's modulus of solid Ω i
um gh
Fr =
Froude number
h [m]
Film thickness
hEG [m]
Ertel and Grubin analytical film thickness value
hc [m]
Central film thickness
[
k i W ⋅ (m ⋅ K )
−1
]
Thermal Conductivity of domain Ω i
L
Dimensionless Moes parameter
M
Dimensionless Moes parameter
n
Normal vector
p [Pa]
Pressure
pcav [Pa]
Cavitation pressure
peps [Pa ]
Pressure parameter tending to pcav
pRoelands [Pa]
Initial Roelands Pressure
p0 [Pa]
Reference pressure
[ ] Q [W ⋅ m ] qi W ⋅ m −2
−3
i
R [m] Re =
Heat source term in domain Ω i Cylinder Radius
ρ3um h η
SRR =
Heat flux of solid Ω i
u1 − u2 um
Reynolds number Slide-to-roll ratio
T [K ]
Temperature
T0 [K ]
Reference Temperature
[
u i m ⋅ s −1
[
u m ⋅ s −1
[
]
u i m ⋅ s −1
]
]
Velocity vector of domain Ω i Fluid Velocity along the x-axis Linear velocity of solid Ω i
16
um =
[
[
u1 + u 2 m ⋅ s −1 2
u * m ⋅ s −1
U [m]
[
v m ⋅ s −1
]
]
Entrainment velocity Sliding velocity Horizontal solid displacement
]
V [m]
[
W N ⋅ m−1
Fluid Velocity along the y-axis Vertical solid displacement
]
Load
x [m]
Spatial coordinate along the direction of flow, in moving frame
X [m]
Spatial coordinate along the direction of flow, in fixed frame
y [m]
Spatial coordinate across the film thickness, in moving frame
Y [m]
Spatial coordinate across the film thickness, in fixed frame
Greek Symbols
α [Pa −1 ] β [K −1 ]
β DH [K −1 ]
Pressure viscosity index Thermoviscosity constant Dowson & Higginson density-temperature coefficient
γ& [s −1 ]
Shear rate tensor
δ Hertz [m]
Analytical Hertzian displacement
∆ [m]
Imposed gap to ensure the load
εs
Solid strain tensor
η [Pa ⋅ s ]
Fluid Viscosity
η0 [Pa ⋅ s ]
Reference Viscosity
η N [Pa ⋅ s ]
Newtonian Viscosity
υi
Poisson ratio of solid Ω i
ρ air [kg ⋅ m −3 ]
Air Density
ρ 0 [kg ⋅ m −3 ]
Reference Density
ρ i [kg ⋅ m −3 ]
Density of domain Ω i
σ f [Pa ]
Fluid stress tensor
σ s [Pa]
Solid stress tensor
τ 0 [Pa]
Ree-Eyring stress 17
ω2 [rad ⋅ s −1 ]
Angular velocity of solid Ω 2
References [1]
W. Habchi, D. Eyheramendy, P. Vergne, and G. Morales-Espejel. A Full-System Approach of
the Elastohydrodynamic Line / Point Contact Problem. ASME Journal of Tribology, 130 (2): 1–10, 2008. [2]
D. Dowson and G.R. Higginson. Elastohydrodynamic lubrication, the fundamentals of roller and
gear lubrication. Pergamon Press, 1966. [3]
C.H. Venner and A.A. Lubrecht. Multilevel methods in lubrication. Tribology Series, 37, 2000.
[4]
H.-S. S. Hsiao, B. J. Hamrock, and J. H. Tripp. Finite element system approach to EHL of
elliptical contacts: Part i—isothermal circular non-Newtonian formulation. ASME Journal of Tribology., 120 (4): 695–704, 1998. [5]
H.-S. S. Hsiao, B. J. Hamrock, and J. H. Tripp. Finite element system approach to EHL of
elliptical contacts: Part ii—isothermal results and performance formulas. ASME Journal of Tribology, 121 (4): 711–720, 1999. [6]
H. S. Cheng. A refined solution to the thermal-elastohydrodynamic lubrication of rolling and
sliding cylinders. ASLE Trans., 8: 397–410, 1965. [7]
M. Kaneta, T. Shigeta, and P. Yang. Film pressure distributions in point contacts predicted by
thermal EHL analyses. Tribology International, 39 (8): 812 – 819, 2006. 4th AIMETA International Tribology Conference. [8]
B. Najji, B. Bou-Said, and D. Berthe. New formulation for lubrication with non-newtonian
fluids. ASME Journal of Tribology, 111 (1): 29–34, 1989. [9]
P. Yang and S. Wen. A generalized reynolds equation for non-Newtonian thermal
elastohydrodynamic lubrication. ASME Journal of Tribology, 112 (4): 631–636, 1990. [10]
P. G. Tucker and P. S. Keogh. On the dynamic thermal state in a hydrodynamic bearing with a
whirling journal using CFD techniques. ASME Journal of Tribology, 118 (2): 356–363, 1996. [11]
P. S. Keogh, R. Gomiciaga, and M. M. Khonsari. Cfd based design techniques for thermal
prediction in a generic two-axial groove hydrodynamic journal bearing. ASME Journal of Tribology, 119 (3): 428–435, 1997. [12]
C.T. Schäfer, P. Giese, W.B. Rowe, and N.H. Woolley. Elastohydrodynamically lubricated line
contact based on the Navier-Stokes equations. Tribology and Interface Engineering Series, 38: 57 – 69, 2000. [13]
H. Yiping, C. Darong, K. Xianmei, and W. Jiadao. Model of fluid-structure interaction and its
application to elastohydrodynamic lubrication. Computer Methods in Applied Mechanics and Engineering, 191 (37-38): 4231 – 4240, 2002. 18
[14]
T. Almqvist and R. Larsson. The Navier-Stokes approach for thermal EHL line contact
solutions. Tribology International, 35 (3): 163 – 170, 2002. [15]
M. Hartinger, M.-L. Dumont, S. Ioannides, D. Gosman, and H. Spikes. CFD modeling of a
thermal and shear-thinning elastohydrodynamic line contact. ASME Journal of Tribology, 130 (4): 041503, 2008. [16]
Y. Wang, H. Li, J. Tong, and P. Yang. Transient thermoelastohydrodynamic lubrication
analysis of an involute spur gear. Tribology International, 37 (10): 773 – 782, 2004. [17]
C. Roelands, J. Vluger, and H. Waterman. Correlational aspects of the viscosity-temperature-
pressure relationship of lubricating oils and in correlation with chemical constitution. Trans. ASME., J. Basic Eng., 11: 601–619, 1963. [18]
L. Houpert. New results of traction force calculations in EHD contacts. ASME Journal of
Tribology, 107: 241–248, 1985. [19]
S. Bair and W. O. Winer. A rheological model for elastohydrodynamic contacts based on
primary laboratory data. ASME J. Lubr. Techn, 101: 258–265, 1979. [20]
P. J. Carreau. Rheological equations from molecular network theories. Trans. Soc. Rheol., 16
(1): 99–127, 1972. [21]
J. L. Tevaarrwerk and K. L. Johnson. A simple non-linear constitutive equation for
elastohydrodynamic oil films. Wear, 35: 345–356, 1975. [22]
J. L. Tevaarrwerk and K. L. Johnson. The influence of fluid rheology on the performance of
traction drives. ASME J. of Lubr. Techn, 101: 266–274, 1979. [23]
H. Eyring. Viscosity, plasticity and diffusion as exalokes of absolute reaction rates. J. Chem.
Phys., 4: 283–291, 1936. [24]
J. L. Tevaarwerk. Traction drive performance prediction for the Johnson and Tevaarwerk
traction model. Technical report, NASA, October 1979. [25]
D. Zhu and S Wen. A full numerical solution for the thermo-elastohydrodynamic problem in
elliptical contacts. ASME Journal of Tribology 106: 246–254, 1984. [26]
H. J. Kim, P. Ehret, D. Dowson, and C. M. Taylor. Thermal elastohydrodynamic analysis of
circular contacts part 2: non-Newtonian model. Proceedings of the IMechE, Part J: Journal of Engineering Tribology, 215-4: 353–362, 2001. [27]
J. Donea, A. Huerta, J.-Ph. Ponthot, and A. Rodriguez-Ferran. Arbitrary lagrangian-eulerian
methods. Encyclopedia of Computational Mechanics., 2004. [28]
R. Gohar. Elastohydrodynamics. Imperial College Press, 2001.
[29]
W. Habchi. A full system finite element approach to elastohydrodynamic lubrication problems:
application to ultra-low-viscosity fluids. PhD thesis, INSA de Lyon, 2008.
19
[30]
M. Carli, K.J.Sharif, E.Ciulli, H.P.Evans, and R.W.Snidle. Thermal point contact EHL analysis
of rolling/sliding contacts with experimental comparison showing anomalous film shapes. Tribology International, 42: 517–525, 2009. [31]
B. J. Hamrock. Fundamentals of Fluid Film Lubrication. McGraw-Hill Companies, 1994.
[32]
Y.P. Chiu and L.B. Sibley. Contact shape and non-Newtonian effects in elastohydrodynamic
point contacts. Lub Eng, 28: 48–60, 1972. [33]
M. Kaneta, H. Nishikawa, K. Kameishi, T. Sakai, and N. Ohno. Effects of elastic moduli of
contact surfaces in elastohydrodynamic lubrication. ASME Journal of Tribology, 114 (1): 75–80, 1992. [34]
M. Kaneta, H. Nishikawa, T. Kanada, and K. Matsuda. Abnormal phenomena appearing in EHL
contacts. ASME Journal of Tribology, 118 (4): 886–892, 1996. [35]
A. Dyson and A. R. Wilson. Film thickness in elastohydrodynamic lubrication at high slide/roll
ratios. Proc. Inst. Mech. Eng., 183 (3P): 81–97, 1968-1969. [36]
A. Cameron. The viscosity wedge. ASLE Trans., 1: 248–253, 1958.
[37]
A. Cameron. Hydrodynamic lubrication of rotating disks in pure sliding. a new type of oil film
formation. J. of Institution of Petroleum, 37: 471–486, 1951. [38]
Kazuyuki Yagi, Keiji Kyogoku, and Tsunamitsu Nakahara. Relationship between temperature
distribution in EHL film and dimple formation. ASME Journal of Tribology, 127 (3): 658–665, 2005. [39]
P. Ehret,
D. Dowson,
and
C. M.
Taylor.
On
lubricant
transport
elastohydrodynamic conjunctions. Proc. R. Soc. A, 1998, vol. 454, pp. 763-787
20
conditions
in
Figures / Legends
Figure 1 - Model Geometry
Figure 2 - Fully coupled resolution principle
21
Figure 3 - Numerical algorithm
22
Figure 4 - Film thickness and pressure comparisons between present model and (a) Newtonian isothermal results with
L = 5.6 and M = 21.9 obtained from a Reynolds model [29] and (b) non-Newtonian thermal results at
SRR = 1 , obtained by [15] with a CFD approach
23
Figure 5 - Effect of SRR on friction coefficient: comparison with Reynolds model with isothermal and thermal conditions
24
Figure 6 - Temperature distribution (a) in the two solid domains in the cylinder
Ω1 and Ω 2 , and (b) in the contact inlet Ω 3 and
Ω 2 in the pure sliding case (white lines represent streamlines)
25
Figure 7 - Temperature comparison on the line
∂Γ3−2 for cases near pure sliding
26
Figure 8 - (a) Viscosity wedge term
∂η ∂u ⋅ N ⋅ m −3 responsible of the non-classical film thickness and (b) ∂y ∂y
resulting horizontal pressure gradient
∂p N ⋅ m −3 (streamlines in white) ∂x
[
[
]
]
Ω3
27
Figure 9 - (a) Temperature distribution in the solids and in the fluid, for with red arrows) through the film thickness
28
SRR = ∞ , (b) velocity field vectors (plotted
Figure 10 - (a) Film thickness versus normal load, (b) sliding velocity variations: on the left, horizontal pressure gradient generated by the viscosity wedge, on the right film thickness, at results)
29
um = 0 and SRR = ∞ (symmetrical
Figure 11 - Film thickness and pressure profiles on
∂Γ3−2 for SRR = [2 4 8]
30
Figure 12 - (a) Viscosity variations with flow streamlines (in white), (b) shear rate across film thickness and (c) temperature distribution in each domain for
SRR = 4
31
u h across film thickness h * = at different sections, for SRR = 4 u h local 1
Figure 13 - Horizontal velocity profiles
32
