in a Radial Lip Seal Using Different Viscosity Models

M A S T E R’S T H E S I S

2006:242 CIV

EVA WETTERHOLM

Simulation of Flow

in a Radial Lip Seal Using Different Viscosity Models

MASTER OF SCIENCE PROGRAMME Engineering Physics

Luleå University of Technology

Department of Applied Physics and Mechanical Engineering

Division of Machine Elements

Preface

This thesis is the final project for the Master of Science in Engineering Physics and a part of the Research Trainee program. This work has been carried out at Luleå Univer- sity of Technology at the Division of Machine Elements, in cooperation with SKF Engi- neering Research Centre, The Netherlands.

I would like to thank my supervisor and examiner Elisabet Kassfeldt for her help and support during the work. I would also like to thank COMSOL support in Stockholm for their help with the software problems and all in the Research trainee group 2005/2006.

Finally I appreciate the time my supervisors at SKF Netherlands, Piet Lugt and Joop Vree have spent in guiding me.

Luleå, May 2006

Eva Wetterholm

Abstract

Radial lip seals for rotating shaft have been used for a long time to keep lubricant within bearings and air and dust particles outside. How they work is not completely under- stood. This study examines operation of a plain radial lip seal lubricated with grease and aims at developing a model that explains the lubrication and sealing behaviour of this type of seal.

One existing hypothesis how the seals function is the Weissenberg effect hypoth- esis. The hypothesis state that the lubricant is exposed to extreme shear load, which will contribute to the pumping effect and the maintenance of the lubricated film. An investigation is therefore made to examine if the non Newtonian characteristics of the grease prevents leakage.

Comparison between Newtonian and non Newtonian fluid calculated by using the finite element software program COMSOL Multiphysics has been performed. The results show that there are differences between Newtonian and non Newtonian flow.

The non Newtonian characteristics contribute to a larger flow in the axial direction as compared to the corresponding Newtonian fluid. This may contribute to a pumping effect and transportation of heat near the minimum film thickness. At the edge of the contact, the non Newtonian fluid becomes more viscous than the Newtonian fluid.

The magnitude of the velocity in the pumping direction for the non Newtonian fluid decreases, which may contribute to the maintenance of the lubricated film.

Contents

Nomenclature 1

1 Radial lip seals 3

1.1 Seal design . . . 4

1.1.1 Elastomeric materials for radial lip seals . . . 5

1.1.2 Garter spring . . . 5

1.2 Lubricating grease . . . 6

2 Non Newtonian fluids 7 2.1 The Weissenberg effect . . . 8

3 Method 9 3.1 Assumptions . . . 9

3.2 Course of action . . . 9

3.3 Swirl flow . . . 10

3.4 The geometry . . . 10

3.5 The boundary conditions . . . 10

3.6 The mesh . . . 11

3.7 Solving the model . . . 11

3.8 The accuracy of the results . . . 12

3.8.1 Numerical estimation of the shear rate . . . 12

4 Results 13 4.1 Cursory investigation of the flow . . . 13

4.1.1 Non Newtonian fluid . . . 13

4.1.2 Newtonian fluids . . . 13

4.2 Thorough investigation of the flow . . . 15

4.2.1 Axial velocity . . . 15

4.2.2 Net flow . . . 15

4.2.3 Effects of the centrifugal force . . . 15

5 Error analysis 19 6 Discussion 21 6.1 Accuracy of the result . . . 21

7 Future work 23

References 25

A The script 27

Nomenclature

F Body force [N]

h₀ Minimum film thickness [m]

m Constant in the power law

n Constant in the power law

p Pressure [Pa]

u (ux, uy, uz) or ur, uϕ, uz

[m/s]

(ux, uy, uz) Velocity vector in Cartesian coordinates [m/s]

u_r, u_ϕ, uz

Velocity vector in Cylindrical coordinates [m/s]

α lip angel air side [^◦]

β lip angel seald side [^◦]

γ˙ Total shear rate

η Dynamic viscosity [Pas]

ρ Density [kg/m³]

τyx Shear stress [kg/m³]

1 RADIAL LIP SEALS

1 Radial lip seals

The rotary lip seal, as shown in Figure 1, is the most common type of rotary shaft seals.

They are used throughout the industry to withstand differences in pressure, contain lubricant and exclude contaminants such as air, water and dust particles.

Figure 1: Radial lip seals are widely used to seal lubricant in many kinds of industrial compo- nents, for example rolling element bearings. The figure shows a schematic cross section of a radial lip seal.

The rotary lip seals has bean used for many decades but how they work is still not completely understood. Between the lip and the running shaft a 1-2 µm thin lubricated film is established [1]. Despite the presence of the film, a successful lip seal exhibits no leakage and a meniscus separates the sealed liquid from the surrounding atmosphere [2], see Figure 2. The meniscus is located on the air side of the lip, at the air side edge of the sealing zone or at some location within the sealing zone. If the air side of a successful seal is flooded with liquid, then the liquid will be pumped to the liquid side.

This reverse pumping is essential for the operation of the seal, if the reverse pumping is very low or zero, the seal will not perform well under normal operating condition [3]. A balanced pumping without net flow in any direction is most favourable.

Figure 2: A schematic figure of the contact area between the shaft and the lip. If the lip is successful a lubricated film will separate the lip from the shaft and a meniscus will be established.

The meniscus separates the lubricant from the surrounding media and the seal will not leak.

Several studies have been performed on oil lubricated radial lip seals. The most common approximation in these studies is that the oil behave like a Newtonian fluid because of the narrow gap between the lip and the shaft, see for example [4]. Other

1.1 Seal design 1 RADIAL LIP SEALS

studies have indicated that it is the non Newtonian behaviour of the lubricant due to the extreme shear load that contributes to the maintenance of the film and the pumping effect [5]. This hypothesis is called the Weissenberg effect.

This study examines the grease lubricated radial lip seal where the non Newto- nian characteristics of the lubricant can not be neglected. The investigation is made by using computer fluid dynamics (CFD) and the lubricating grease is assumed to be- have according to the power law. A comparison between non Newtonian fluids with power law and Newtonian fluids will be performed, to investigate if the non Newtonian behaviour prevents leakage.

A part of this project was to carry out a literature study and to determine how the model should be simplified to make it solvable. The project considered different directions, but the deformation of the lip was the main concern. The investigation also considered as to how the anglesα andβ in Figure 2 influenced the deformation of the lip due to the flow of grease, a fluid structure interaction model. The grease was then assumed to be a thick Newtonian fluid. After a literature survey about radial lip seals and a course on introduction to COMSOL Multiphysics, the project went on with simulation of the flow.

1.1 Seal design

The lip is a circular elastic part made of an elastomer bonded to a rigid case, see Figure 3. The inner diameter of the lip is smaller than the shaft’s outer diameter so that when the seal is installed the lip is stretched outward and creates a force between the lip and shaft. As the shaft rotates, oscillates or reciprocates, the lip will flex and follow the shafts motion to prevent leakage. A garter spring is often added to compensate for changes in the rubber properties that occur when the material is subjected to heat and different kinds of chemicals.

Figure 3: A schematic figure of the sealed bearing unit where the lip seal is bounded to a rigid case. A garter spring is added to compensate for the changes in the rubber properties that will occur when the lip is subjected to chemicals.

1 RADIAL LIP SEALS 1.1 Seal design

An important factor of the lip is the cross section, for an oil seal the angle relative to the shaft axis on the air side must be smaller than the corresponding angle on the liquid side, see Figure 2. The angles also must lie within a certain limit if the seal is to be successful [6]. If a successful seal is mounted by inverting the air and the liquid side, then the seal will leak excessively [7].

Since the seal behaviour is very complex, most of the development in both seal design and materials have been due to trial and error. That is why there are many types of radial lip seal designs and material to choose from.

1.1.1 Elastomeric materials for radial lip seals

The elastomeric material in the lip consists mostly of a base polymer mixed with fillers, antioxidants and curing agents. The material must be compatible with the sealed fluid, be abrasion and tear resistant, and maintain physical properties for an extended period of time. The elastomer usually tends to swell in contact with oil and chemicals, which can influence the lip quality. Choice of the right material is essential for the seals life and there are many different kind of materials to choose from. Two of them are Flouroelastomer (FKM) and Nitrile (NBR).

Fluoroelastomer

Fluoroelastomers is an expensive material but it is often used in rotating shaft seal applications because it has a well oil and chemical resistance. It has good high temper- ature properties and will not fracture in low temperature. Low temperature can make the elastomer stiff, and it may not follow the shaft wobble and vibrations until fric- tion heat warms up the elastomer. This could cause sporadic leakage if the lubricant remains at low temperatures. Fluoroelastomer compounds are mostly used for high speed, high temperature operations with spars lubrications.

Nitrile

Nitrile is a low cost material relative to other polymers, with good oil and abrasion resistance. The biggest disadvantage of the compound is the lack of heat resistance at elevated temperatures. The material will harden, crack and lose interference with the shaft, which will result in failure. Nitrile compounds are mostly used when the operating condition is mild.

1.1.2 Garter spring

The elastomer of the lip will absorb the lubricant when exposed to elevated tempera- tures. The seal lip swell and softens, which results in a decrease of lip interference and load between the lip and the shaft. In the worst scenario the lip expands completely of the shaft and results in a leakage. A garter spring is used to compensate for material changes and provide a uniform load which will extend the life of the seal. The garter spring is a helical coiled spring whose ends are connected so that the spring forms a circle. The seal life and lip wear width is related to spring tension when the spring is used.

1.2 Lubricating grease 1 RADIAL LIP SEALS

1.2 Lubricating grease

Lubrication is used to minimise wear and help bearings and other mechanical applica- tions to run smoothly. The most common way is to lubricate with oil or grease.

The main difference between these two lubricants is their difference in consistency.

Grease consistency is more semi-solid than oil and that is the key to many of its ad- vantages. For example, grease can prevent external contaminants from entering the bearing and it is less easily displaced from bearings surfaces than oil. It will therefore not drain from the bearing under gravity. However, unlike oil, grease provides a very poor means of transferring friction heat away from the sliding surface, which usually leads to higher friction due to viscous forces.

Lubricating grease consists of a mineral based oil and thickening agents. The ad- ditives can be of metallic soap, urea compound, carbon black or other materials. The components of a grease make it very hard to describe its rheological properties. The polymeric substance has non Newtonian characteristics hence its resistance to shear and elastic properties. The apparent viscosity varies with temperature, pressure, time and shear. At very low shear rate the grease acts as a solid or semisolid. When the shear rate increases, the viscosity approaches that of the base fluid.

2 NON NEWTONIAN FLUIDS

2 Non Newtonian fluids

The mechanical behaviour of a fluid when it is subjected to forces depends on its rheo- logical properties. One way to define these properties is how the shear strain responds to shear stress. For Newtonian fluids the shear stress is proportional to the relative rate of movement

τyx= µ˙γyx (2.1)

where µ is the dynamic viscosity, which is by definition only dependent on temperature and pressure. ˙γxyis the rate of shear strain or shear rate

˙γyx=dγyx

dt =du_x

dy. (2.2)

When the stress is removed from the fluid, the shear rate goes to zero, i.e. the motion stops, but there is no tendency of the fluid to return to any past state. If the properties of the fluid are such that the shear stress and the shear rate are not directly proportional but are instead related by some more complex function, the fluid is said to be non Newtonian. One approximation of this complex function is the power law where the dynamic viscosity depends on the shear rate

η= m˙γⁿ⁻¹ (2.3)

where m and n are scalars that characterize the non Newtonian fluid. The total shear rate ˙γare in Cartesian respectively in cylindrical coordinate systems defined as

γ˙= s

2

∂ux

∂x

2

+ 2

∂u_y

∂y

2

+ 2

∂uz

∂z

2

+

+

∂ux

∂y +∂uy

∂x

2

+

∂uy

∂z +∂uz

∂y

2

+

∂ux

∂z +∂uz

∂x

2

, (2.4)

γ˙= s

2

∂ur

∂r

2

+ 2 1 r

∂uϕ

∂ϕ +u_r r

2

+ 2

∂uz

∂z

2

+

+

 r∂

∂r

u_ϕ r

+1 r

∂ur

∂ϕ

2

+ 1 r

∂uz

∂ϕ +∂uϕ

∂z

2

+

∂ur

∂z +∂uz

∂r

2

. (2.5) To solve the velocity profile needed to calculate 2.4 and 2.5 Navier-Stokes equation is used, with the assumption that the fluid is incompressible

ρ∂u

∂t −∇η(∇u+ (∇u)^T) +ρ(u ·∇)u +∇p= F (2.6)

∇· u = 0. (2.7)

In the equations above is u the velocity vector,(ux, uy, uz) in Cartesian coordinates and u_r, uϕ, uz
in cylindrical coordinates,ρis the density of the fluid, p is the pressure and F is the body force term.

2.1 The Weissenberg effect 2 NON NEWTONIAN FLUIDS

2.1 The Weissenberg effect

There are many ways the human eye can see differences between the Newtonian and non Newtonian fluid. One example is the Weissenberg rod climbing effect, which will occur when a vertical rod is rotating in a cup of non Newtonian fluid. If the fluid is of Newtonian character the centrifugal force causes the fluid to move radially towards the cups wall, see Figure 4(a). For a non Newtonian or polymeric liquid the fluid moves to the rotating rod, see Figure 4(b). This is caused by the influence of normal stresses on flow properties. These normal stresses create tension along the circular lines of flow and generate a pressure towards the centre, which drives the fluid up the rod. This phenomenon was first described by Garner and Nissan [8], the experiment was then analyzed by Weissenberg [9].

In the radial lip seal should it therefore be, caused by the normal stresses, a flow in the axial direction. That flow could be a reason why the lubricated film is established and why it is remained.

(a) In Newtonian fluids, centrifugal forces generated by the rotation push the fluid away from the rod.

(b) In non Newtonian fluids, normal forces are stronger than centrifugal force and force the fluid inward toward the rod.

Figure 4: Weissenberg effect, also called the rod climbing effect.

3 METHOD

3 Method

Radial lip seals are very complex to describe mathematically; to make it solvable sev- eral assumptions and approximations had to be made. The equations, boundary- and surface-settings that should be used were also needed to be determined. The geometry was also modified to minimise the singularities to make it solvable.

3.1 Assumptions

• The model is axial symmetric, the gradient in the rotational direction is zero.

• Run-in lip, the examined lip is worn and has a smooth curvature.

• Smooth surfaces, there is no roughness on either the lip or the shaft.

• There will be no deformation on the lip due to shaft motion or flow of the lubri- cant.

• Fully flooded contacts, there is no contact between the lip and the shaft.

• No external forces. The gravitational force, the radial force from the garter spring and forces from the lip are assumed to be zero.

• Isothermal conditions, there is no changes in temperature.

• The Newtonian fluid is incompressible, with constant density and viscosity.

• The non Newtonian fluid is incompressible with constant density, and the dy- namic viscosity is according to the power lawη= m˙γⁿ⁻¹.

3.2 Course of action

In this investigation, a 2D axial symmetry model with swirl flow was studied. To make the comparison between Newtonian and non Newtonian, a non Newtonian model was calculated first. It had the properties m= 1000 and n = 0.2 in the power law described in chapter 2 on page 7 and the density was set to 960 kg/m³. The properties were chosen in collaboration with SKF ERC, the Netherlands.

Two Newtonian fluids were then chosen, based on the results from the non New- tonian model. The first one was chosen to have the same viscosity as the non Newtonian at the area of highest shear i.e. the area of the thinnest film thickness. The second one was chosen to have the same viscosity as the non Newtonian had at the edge of the area of interest. The velocity profiles in the z direction, see Figure 2 on page 3, for the different models were then compared to investigate if the non Newtonian characteristics prevents leakage.

The models were calculated with a radius of 0.05 m. The radius was then increased for the non Newtonian model. It was interesting to see how the flow in the z direction would change when the centrifugal forces decreased. The velocity on the shaft periph- ery was kept at the same speed, i.e. the angular velocity of the shaft was decreased so that the velocity on the boundary was the same for every case.

3.3 Swirl flow 3 METHOD

The models were solved with aid of COMSOL Multiphysics and the post process- ing was modified in an m-file in MATLAB. Comparisons between the maximum and minimum velocity in the z direction, net flow in the z direction and the velocity profile at the area of interests were made. See appendix A for the m-file that were used for the 2D axial symmetric model.

3.3 Swirl flow

Swirl flow is a swirling flow in 2D axial symmetric coordinate systems. This means that the rotational speed is important, although constant, so that the gradient in this direction can be considered to be zero. This application take the volume force in the rotational direction, F_ϕinto consideration. The rotational velocity u_ϕwill still remain to be solved in the Navier Stokes equation, see equation 2.6. The equation of motion can therefore be simplified as

ρ u_r∂ur

∂r −u_ϕ² r + uz

∂ur

∂z

 +∂p

∂r = η 1 r

∂

∂r

 r∂ur

∂r



−u_r r²+∂²u_r

∂r²



ρ u_r∂u_ϕ

∂r −u_ru_ϕ r + uz

∂ur

∂z

 +∂p

∂r = η 1 r

∂

∂r

 r∂u_ϕ

∂r



−u_ϕ r² +∂²u_ϕ

∂r²

 (3.1)

ρ u_r∂uz

∂r + uz

∂uz

∂z

 +∂p

∂r = η 1 r

∂

∂r

 r∂uz

∂r

 +∂²u_r

∂z²

 ,

which will be used when solving the models.

3.4 The geometry

The lubricated geometry was modelled in drawing mode in COMSOL Multiphysics.

The geometry was chosen to have a smooth curvature with a minimum film thickness of 4 µm, because the geometry should have as few discontinuities as possible. The contact width was 1.1 mm and extended a bit longer on the lubricated side and the position of minimum film thickness was located at z= 0. The curvature near the air side had a smaller angle than the lubricated side, all this because the geometry should be as similar to the reality as possible. For the same reason, the outflow on the air side was chosen to have smaller width than the outflow on the lubricated side, see Figure 5.

The outflows on both sides were then extended, so the boundary condition should have as little affect on the flow near the area of interest as possible. The shaft radius was chosen to be 0.05 m. All the parameters above were chosen in collaboration with SKF ERC, the Netherlands.

3.5 The boundary conditions

At the inlet and outlet the boundary conditions were set to pressure zero so the flow in the z direction should be affected as little as possible. The extension on the inlet and outlet were chosen to be as long as it were needed so the pressure distribution were set to normal when entering the lip sealing zone.

3 METHOD 3.6 The mesh

The lip side of the geometry had a no slip condition, i.e. there was no movement on the boundary. On the shaft side, a velocity field with a constant velocity of 5.24 m/s in theϕdirection were implemented. This corresponded to an angular velocity of 1000 rpm when the shaft had a radius of 0.05 m.

3.6 The mesh

For this application, a structured mesh was the best choice and received the best con- vergence but because of some problem in the post processing, caused by the large difference between the film thickness and the shaft radius, an unstructured mesh had to be chosen. After convergence tests of the mesh were made, a mesh consisted of element with the size of 0.5 µm at the area of contact and 4 µm at the extended inlet and outlet were implemented.

3.7 Solving the model

COMSOL Multiphysics solves all physics with finite element method (FEM). Math- ematically this method is used for finding approximations of partial differential equa- tions (PDE), where the domain is broken into a set of discrete volumes or finite ele- ments. The equations are multiplied by a weight function before they are integrated over the domain. In the simplest finite element method, the solution is approximated by a linear shape function within each element in a way that guarantees continuity of the solution across element boundaries. In these calculations, a second order Lagrange elements are used for the velocity components and linear elements for the pressure, which are the default elements for the incompressible Navier-Stokes equation mode.

The solvers in COMSOL Multiphysics are either direct or iterative. Both type of solvers break down the problem into the solution of one or several linear systems. The direct solvers solve linear systems with Gaussian elimination, which is a stable process well suited for ill-conditioned systems according to COMSOL Multiphysics guide. The problem with direct solvers is that they are memory consuming. The iterative solvers are less memory consuming but they are not as fast as the direct solvers. Both types

Figure 5: The geometry used in the calculation, it was chosen to be asymmetric to resemble the real seal lip. The whole geometry is lubricated.

3.8 The accuracy of the results 3 METHOD

of solvers were used to solve this model. Because the model had to be solved several times for obtaining accurate result, the mesh was refined each time.

For this particular problem a stationary non linear solver was needed. One of the parameters in this mode is the relative tolerance, which give the criterion for conver- gence. The software stops calculating when the relative error, err is less than the relative tolerance. The relative error is given by

err= 1 N

∑

i

(|Ei|/Wi)²

!1/2

. (3.2)

where N is the number of degrees of freedom, Wiis a weight factor and E is an estimated error in to a current approximation to the true vector.

3.8 The accuracy of the results

The software checked for convergence of iterations when the relative tolerance was set to 10⁻⁹. A convergence test of the mesh was also made to investigate if the solution at a certain point will converge to a specific value, when the number of elements in the geometry goes to infinity. A Richardson interpolation test was to be constructed for an estimation of the error.

3.8.1 Numerical estimation of the shear rate

One more way to check the accuracy of the model is to compare the simulations with a numerical estimation of the shear rate. Equation 2.5 on page 7 which describes the total shear rate was used. After terms that were zero or close to zero were cancelled, the shear rate could be approximated as

γ˙≈ s

 r∂

∂r

u_ϕ r

2

+

∂u_ϕ

∂z

2

. (3.3)

Because the contact width is much longer than the film thickness, equation 3.3 can be shortened and simplified as

˙γ≈ r∂

∂r

u_ϕ r

=∂uϕ

∂r −u_ϕ r ≈u_ϕ

h₀−u_ϕ r ≈u_ϕ

h₀. (3.4)

The numerical values of u_ϕ= 5.24 m/s and h0= 4 µm gives a shear rate of 1.31 Ms⁻¹.

4 RESULTS

4 Results

The results from the simulations are presented below. To investigate the pumping abil- ity of a seal, the flow in the axial direction is in the main concern.

The non Newtonian fluid was first calculated and the two Newtonian fluids were chosen based on that result. The simulations were computed for the geometry in Figure 5 on page 11 but only the contact area was analysed. The contact area is located between the z coordinates -50 µm and 60 µm. The extended inlet and outlet are only used to minimize the effects of the boundary conditions and are therefore not examined.

4.1 Cursory investigation of the flow

A cursory investigation of the non Newtonian flow and the Newtonian flows were first made to see if there are any distinct differences.

4.1.1 Non Newtonian fluid

The pressure distribution, shear rate, dynamic viscosity and the velocity profile of the non Newtonian fluid were calculated and results are shown in Figure 6. The Reynolds number is smaller than 0.7 for the whole domain, the flow is therefore laminar with no turbulent tendency.

Figure 6(a) shows the pressure distribution, which has fluctuations in the extended inlet and outlet. As desired, the distribution has stabilized when entering the contact area. Figure 6(b) shows the shear rate, which has the maximum value of 1.41 · 10⁶ located at z= 0. The shear rate is decreasing towards the edges.

The dynamic viscosity in Figure 6(c), has the smallest value of 0.012 Pas at the point of highest shear, i.e. at the minimum film thickness. Near z= 30 µm and at the edge of the contact area at the air side are the dynamic viscosity 0.03 Pas. Those values will be used in the following calculations of the Newtonian fluids.

In Figure 6(d), a pumping effect is visible, from both sides near the lip is the fluid pumping towards the point of minimum film thickness. But it is not certain if it is due to the non Newtonian characteristics together with the asymmetric geometry or due to the centrifugal force.

4.1.2 Newtonian fluids

Calculations of the Newtonian fluids with viscosity 0.012 Pas and 0.03 Pas shows sim- ilar surface velocity profile as the non Newtonian fluid in Figure 6(d). The normalized velocity arrows for the three fluids can not be distinguished from each other. The thicker fluid has the same magnitude of the z velocity as the non Newtonian fluid and the thinner fluid has a slightly larger.

4.1 Cursory investigation of the flow 4 RESULTS

(a) Surface plot of the pressure distribution. The lowest pressure occurs near the shaft.

(b) Surface plot of the shear rate. Maximum shear is at the thinnest film thickness, which is indicated.

(c) Surface plot of the dynamic viscosity. The minimum viscosity is 0.012, at the edge of area of interest is the viscosity 0.03.

(d) Surface plot of the velocities profile in the z direction, i.e. the axial direction. The arrows indicate the normalized velocity in the r and z direction.

Figure 6: Results from the calculation of the non Newtonian fluid

4 RESULTS 4.2 Thorough investigation of the flow

4.2 Thorough investigation of the flow

Because the velocity profiles of the non Newtonian and Newtonian fluids were alike a more thorough investigation was needed. Comparisons between the cross section of the flow in the axial direction and net flow in the axial direction were performed. An investigation of the centrifugal effects for the non Newtonian fluid was also made.

4.2.1 Axial velocity

Cross section plots in Figure 7 show that the non Newtonian fluid is more alike the thicker Newtonian fluid. However the non Newtonian fluid has larger axial velocity near z= 0, i.e. at the minimum film thickness. The curves of the Newtonian fluids are very similar to each other; the only difference is the magnitude of the z velocity. The cross sections are taken at the same positions for all three models, where the z was kept constant and r value is taken from the lip to the shaft.

The cross section plots, Figure 7 and the surface velocity plot in Figure 6(d) indicate that the flow in z direction will go towards the point of minimum film thickness at the lip side and to the other direction at the shaft side for the same z value. The maximum z velocity will therefore lie near the shaft for positive z and the minimum near the lip.

The opposite occurs for negative z positions. A simple summary of Figure 6(d) and 7 is therefore to plot the maximum and minimum z velocities as a function of the axial position. The results are shown in Figure 8.

The two Newtonian fluids in Figure 8(b) and 8(c) are comparable, the only differ- ence is the power of the velocity. Figure 8 also shows that the non Newtonian fluid is more alike to the thicker Newtonian fluid, regarding the magnitude of the axial velocity.

In Figure 8 the positions of minimum film thickness are marked. For the three fluids, the axial velocity is slightly increased at the position of minimum film thickness as compared to that in the nearest neighbourhoods. A stirring effect occurs when z<

−30 µm and z > 20 µm. The effect is much larger for the Newtonian fluids than for the non Newtonian fluid.

4.2.2 Net flow

The integral of the z velocity profiles in Figure 7, i.e. the net flows, are constant through the contact area for all three models. The non Newtonian has a calculated net flow of 8.2 · 10⁻¹²s⁻¹ throughout the contact in the z direction. The Newtonian fluids has a slightly larger net flow of 1.2 · 10⁻¹¹s⁻¹for the thinner fluid and a slightly smaller net flow of 5.0 · 10⁻¹²s⁻¹for the thicker.

4.2.3 Effects of the centrifugal force

The shaft radius was then increased while the velocity at the periphery of the shaft was kept at the same speed. The test was only performed for the non Newtonian fluid.

The results show similar velocity profiles for the lubricant in the axial direction for all different shaft radiuses. The only difference is the magnitudes of the axial velocities and they are exponentially decreased when the radius is increased by a factor 10. Figure 9 shows the maximum velocities in axial direction at z= 0 for different shaft radius.

4.2 Thorough investigation of the flow 4 RESULTS

(a) Velocity in the z direction at different cross section for the non Newtonian fluid.

(b) Velocity in the z direction at different cross section for the thicker Newtonian fluid.

(c) Velocity in the z direction at different cross section for the thinner Newtonian fluid.

Figure 7: Cross section plots of the velocity in the z direction for the three different models.

4 RESULTS 4.2 Thorough investigation of the flow

−5 −4 −3 −2 −1 0 1 2 3 4 55 x 10⁻⁵

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1 1x 10⁻⁵

z [m]

velocity z−direction [m/s]

max min

(a) Maximum and minimum velocities for the non Newtonian fluid.

(b) Maximum and minimum velocities for the thicker Newtonian fluid.

−5 −4 −3 −2 −1 0 1 2 3 4 55 x 10⁻⁵

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1 1x 10⁻⁵

z [m]

velocity z−direction [m/s]

max min

(c) Maximum and minimum velocities for the thin- ner Newtonian fluid.

Figure 8: Maximum and minimum velocities trough the contact area in the z direction. The position of minimum film thickness is indicated with a circle.

10⁻² 10⁻¹ 10⁰ 10¹

10⁻⁸ 10⁻⁷ 10⁻⁶ 10⁻⁵

Shaft radius [m]

Velocity [m/s]

Figure 9: Maximum axial velocity at z= 0 with different shaft radius when the periphery speed of the shaft is constant.

4.2 Thorough investigation of the flow 4 RESULTS

5 ERROR ANALYSIS

5 Error analysis

Structured mesh gave the best convergence but because of problem in the post pro- cessing an unstructured mesh had to be used. The problem was caused by the large difference between the shaft radius and the film thickness.

Convergence tests show that the structured mesh is in the asymptotic area, see Fig- ure 10(a), the unstructured is not, see Figure 10(b). The tests were made for different numbers of element with respective values for the shear rate at the point(0.05, 0). It is clear that the unstructured mesh need more elements, but because of the computer available was not as powerful, a finer mesh could not be constructed. A Richarsons extrapolation test could not therefore be performed to get an estimation of the error of discretization. However the values were in a small range of variation and compar- isons between the surface plots of the different element sizes of the triangular mesh and comparisons between the structured mesh do not show any visible differences.

The error of iteration should also be insignificant because the software is pro- grammed to check for convergence.

The estimated value of the shear rate on page 12 was calculated to 1.31 · 10⁶. A comparison with the calculated value in Figure 10 of 1.25 · 10⁶shows that calculations seems to be in order.

0 0.4 0.8 1.2 1.61.6

x 10⁻³ 1.253

1.2532 1.2534 1.2536 1.2538

1.254x 10⁶

Elements⁻¹

Shear rate

(a) Test of convergence for the structured mesh

0 1 2 3 4

x 10⁻⁴ 1.253

1.2532 1.2534 1.2536 1.2538

1.254x 10⁶

Elements⁻¹

Shear rate

(b) Test of convergence for the unstructured mesh Figure 10: Test of convergence, the shear rate is calculated at(0.05, 0)

5 ERROR ANALYSIS

6 DISCUSSION

6 Discussion

The pumping ability is essential for a radial lip seal; how the flow appears beneath the lip affect that phenomenon. The calculations show that there are differences between the axial flow for the non Newtonian and Newtonian fluids.

The velocity profiles in Figures 7 and 8 show that the non Newtonian fluid is almost similar to the thicker Newtonian fluid than to the thinner. But the non Newtonian fluid had a larger net flow through the contact compared to the thicker Newtonian fluid.

Even the maximum velocity of the non Newtonian fluid was larger at the position of minimum film thickness. This may indicate that the non Newtonian fluid has a larger pumping ability.

Lubricating grease provides a poor mean of transporting friction heat away. But at the small scale, near the location of minimum film thickness, the larger pumping may contribute to transportation of heat away from that area. This will influence the viscosity of the lubricant.

Near the edge of the contact, at about z= 40 µm and z = −40 µm there is a smaller difference between the maximum and minimum z velocity for the non Newtonian than for the Newtonian fluids, see Figure 8. This may contribute to the maintained lubricated film because the magnitudes of the velocity in the axial direction are smaller for the non Newtonian fluid. The non Newtonian fluid may be compared to as a plug, the movement is getting smaller where the shear rate is decreasing and the fluid is getting more viscous.

A Weissenberg effect can not be separated from the results. The velocity in the axial direction seems to go towards zero when the centrifugal effect is minimized. However, the Weissenberg effect can not be discarded as a hypothesis. The power law is a rough approximation of the grease rheology and it is only a shear thinning formula. The formula does not take the established normal forces into consideration, as described in chapter 2.1. A more complex approximation of the grease rheology is needed if the normal forces are to be considered.

6.1 Accuracy of the result

One essential question in this model is if the axial velocities arise only because of nu- merical errors. Test has been made when the speed of shaft have been decreased and increased. The investigation showed that the velocity in the axial direction decreases when the speed of the shaft decreased and increases when the speed of the shaft in- creased, which is a normal behaviour. The convergence test of the structured mesh, which was in the asymptotic area, also indicates that the velocity in the axial direction was due to the asymmetric geometry and due to the curvature of the shaft.

6.1 Accuracy of the result 6 DISCUSSION

7 FUTURE WORK

7 Future work

This is a simplified model of the radial lip seal and the power law is a rough approxima- tion of grease rheology. To have a more accurate investigation, a better approximation is needed to examine the importance of the non Newtonian characteristics. It is also interesting to take the meniscus into consideration, how it behaves and influences the flow of the lubrication.

Because the viscosity of the grease is strongly temperature dependent, it would also be interesting to have heat transfer within the model.

A different direction of the project would also to have a fluid structure interaction.

The lip is made of rubber and should be deformed under operating condition, but that would require a lot of computer capacity.

7 FUTURE WORK

REFERENCES REFERENCES

References

[1] E.T. Jagger. Study of the lubrication of synthetic rubber rotary shaft seal. Proc.

Conf. Lubric. Wear, pages 409–415, 1957.

[2] M. J. L. Stakenborg. On the sealing mechanism of radial lip seal. Tribol. Int., 21(6):335–340, Dec 1988.

[3] L. Horve. The correlation of rotary shaft radial lip seal service reliability and pump- ing ability to wear track roughness and microasperities formation. SAE, 100(sect 6):620–627, 1991.

[4] R. Salant. Modelling rotary lip seals. Wear, 207(1-2):92–99, Jun 1997.

[5] F. Schulz, K. Wiehler, V. M. Wollesen, and M. Voetter. A molecular-scale view on rotary lip sealing phenomena. Leed- Lyon Symposium on Tribology, (25):457–466, 1998.

[6] F. Hirano and H. Ischiwata. The lubricating condition of a lip seal. Proc. Inst.

Mech. Engrs., 180(Pt. 3B):138–147, 1965-66.

[7] Y. Kawahara, M. Abe, and H. Hirabayashi. An analysis of sealing characteristics of oil seals. ASLE Trans., 23:93–102, 1980.

[8] F. H. Garner and A. H. Nissan. Rheological properties of high viscosity solution of long molecules. Nature, 158:634–635, 1946.

[9] K. Weissenberg. A continuum theory of rheological phenomena. Nature, 159:310–

311, 1947.

REFERENCES REFERENCES

A THE SCRIPT

A The script

This is a shorten version of the script that was used to solve the model. In reality the model had to be solved repeatedly to get a fine mesh. The solution was saved for every time it was solved and the solution was used as an initial value for the next time. In this extraction the subdomain settings are not included.

% COMSOL Multiphysics Model M-file

% Generated by COMSOL 3.2 (COMSOL 3.2.0.222, $Date: 2005/09/01 18:02:30 $) flclear fem

% COMSOL version clear vrsn

vrsn.name = ’COMSOL 3.2’;

vrsn.ext = ’’;

vrsn.major = 0;

vrsn.build = 222;

vrsn.rcs = ’$Name: $’;

vrsn.date = ’$Date: 2005/09/01 18:02:30 $’;

fem.version = vrsn;

% Constants

r_shaft=0.05; % Radius of the shaft

r_nere=0.00003; % Constant to construct the geometry z_nere=-5e-5; % Constant to construct the geometry z_uppe=5e-5+1e-5; % Constant to construct the geometry r_uppe=(z_uppe+5e-6)/5.5; % Constant to construct the geometry z_langd_in=2e-5; % Length of the extended outlet air side

z_langd_ut=3e-5; % Length of the extended outlet lubricated side

% Boundary constants

P1=0; % Pressure at air side

Ps=0; % Pressure at lubricated side

omega_shaft=1000*2*pi/60; % Speed of shaft

% Properties for grease

rho_grease=960; % Density of grease m_grease=1000;

n_grease=0.2;

m_Newmin=0.0121;

m_Newmax=0.03;

n_New=1;

m=m_grease;

%m=m_Newmin;

%m=m_Newmax;

A THE SCRIPT

n=n_grease;

%n=n_New;

fem.const = {’r_shaft’,r_shaft,...

’P1’,P1,...

’Ps’,Ps,...

’omega_shaft’,omega_shaft,...

’m_grease’,m_grease,...

’n_grease’,n_grease,...

’rho_grease’,rho_grease,...

’m’,m,...

’n’,n};

% Descriptions clear descr

descr.const= {’r_shaft’,’radius of shaft [m]’,...

’P1’,’pressure [Pa]’,...

’Ps’,’pressure [Pa]’,...

’omega_shaft’,’speed of shaft [rad/s]’,...

’m_grease’,’m’,...

’n_grease’,’n’,...

’rho_grease’,’density of grease [Kg/m^3]’};

fem.descr = descr;

% Drawing the geometry

carr={curve2([0,0],[5.0E-5,z_nere],[1,1]), ...

curve2([0,r_nere],[z_nere,z_nere],[1,1]), ...

curve2([r_nere,0,0,1.0E-5],[z_nere,-1.5E-5,5.0E-6,5.0E-5],[1,1,1,1]),...

curve2([1.0E-5,0],[5.0E-5,5.0E-5],[1,1])};

g5=geomcoerce(’solid’,carr);

g1=rect2(r_nere,z_langd_ut,’base’,’corner’,’pos’,[0,z_nere-z_langd_ut]);

carr={curve2([0.00001,r_uppe],[5.0E-5,z_uppe],[1,1]), ...

curve2([r_uppe,0],[z_uppe,z_uppe],[1,1]), ...

curve2([0,0],[z_uppe,5.0E-5],[1,1]), ...

curve2([0,0.00001],[5.0E-5,5.0E-5],[1,1])};

g3=geomcoerce(’solid’,carr);

g2=geomcomp({g5,g3},’ns’,{’CO1’,’CO3’},’sf’,’CO1+CO3’,’edge’,’all’);

g6=rect2(r_uppe,z_langd_in,’base’,’corner’,’pos’,[0,z_uppe]);

g3=mirror(g1,[0,0],[0,1]);

g4=mirror(g2,[0,0],[0,1]);

g5=mirror(g6,[0,0],[0,1]);

A THE SCRIPT

clear g1 clear g2 clear g6

g3=move(g3,[r_shaft,8.7e-6]);

g4=move(g4,[r_shaft,8.7e-6]);

g5=move(g5,[r_shaft,8.7e-6]);

clear s

s.objs={g3,g4,g5};

s.name={’CO1’,’CO3’,’CO5’};

s.tags={’g3’,’g4’,’g5’};

fem.draw=struct(’s’,s);

fem.geom=geomcsg(fem);

% Initialize mesh

fem.mesh=meshinit(fem, ...

’hmax’,[], ...

’hmaxfact’,1, ...

’hgrad’,1.3, ...

’hcurve’,0.3, ...

’hcutoff’,0.001, ...

’hnarrow’,1, ...

’hpnt’,10, ...

’xscale’,1.0, ...

’yscale’,1.0, ...

’mlevel’,’sub’, ...

’hmaxsub’,[1,4e-6,2,0.5e-6,3,4e-6]);

% Application mode 1 clear appl

appl.mode.class = ’NonNewtonian’;

appl.mode.type = ’axi’;

appl.dim = {’u’,’v’,’w’,’p’};

appl.sdim = {’r’,’phi’,’z’};

appl.name = ’chnn’;

appl.module = ’CHEM’;

appl.shape = {’shlag(2,’’u’’)’,’shlag(2,’’v’’)’,...

’shlag(2,’’w’’)’,’shlag(1,’’p’’)’};

appl.gporder = {4,2};

appl.cporder = {2,1};

appl.sshape = 2;

appl.border = ’off’;

appl.assignsuffix = ’_chnn’;

clear prop

A THE SCRIPT

prop.elemdefault=’Lagp2p1’;

prop.analysis=’static’;

prop.stensor=’full’;

prop.nisot=’Off’;

prop.swirl=’On’;

prop.frame=’rz’;

clear weakconstr

weakconstr.value = ’off’;

weakconstr.dim = {’lm1’,’lm2’,’lm3’};

prop.weakconstr = weakconstr;

appl.prop = prop;

clear pnt pnt.pnton = 0;

pnt.p0 = 0;

pnt.ind = [1,1,1,1,1,1,1,1,1,1];

appl.pnt = pnt;

clear bnd bnd.u0 = 0;

bnd.type = {’noslip’,’neutral’,’uv’,’out’};

bnd.p0 = 0;

bnd.w0 = {0,0,’omega_shaft*0.05’,0};

bnd.v0 = 0;

bnd.name = ’’;

bnd.ind = [3,4,3,2,3,3,2,4,1,1,1,1];

appl.bnd = bnd;

clear equ

equ.eta_inf = 0;

equ.delcd = 0.35;

equ.init = 0;

equ.F_phi = 0;

equ.shape = [1;2;3;4];

equ.F_z = 0;

equ.cporder = {{1;1;1;2}};

equ.delid = 0.5;

equ.sdtype = ’pgc’;

equ.cdon = 0;

equ.usage = 1;

equ.type_visc = ’power’;

equ.delsd = 0.25;

equ.kappadv = 0;

equ.delps = 1;

equ.rho = ’rho_grease’;

equ.cdtype = ’sc’;

equ.eta0 = 1;

equ.gporder = {{1;1;1;2}};

equ.m = ’m’;

A THE SCRIPT

equ.sdon = 0;

equ.lambda = 0;

equ.pson = 0;

equ.n = ’n’;

equ.F_r = 0;

equ.idon = 0;

equ.ind = [1,1,1];

appl.equ = equ;

fem.appl{1} = appl;

fem.sdim = {’r’,’z’};

fem.frame = {’rz’};

% Shape functions

fem.shape = {’shlag(2,’’u’’)’,’shlag(2,’’v’’)’,...

’shlag(2,’’w’’)’,’shlag(1,’’p’’)’};

% Integration order fem.gporder = {4,2};

% Constraint order fem.cporder = {2,1};

% Geometry shape order fem.sshape = 2;

% Simplify expressions fem.simplify = ’on’;

fem.border = 1;

% Equation form fem.form = ’general’;

fem.units = ’SI’;

% Global expressions fem.expr = {};

% Functions clear fcns

fem.functions = {};

% Descriptions clear descr

descr.const= {’r_shaft’,’radius of shaft [m]’,...

’P1’,’pressure [Pa]’,...

’Ps’,’pressure [Pa]’,...

’omega_shaft’,’speed of shaft [rad/s]’,...

A THE SCRIPT

’m_grease’,’m’,...

’n_grease’,’n’,...

’rho_grease’,’density of grease [Kg/m^3]’};

fem.descr = descr;

% Solution form fem.solform = ’weak’;

% Multiphysics

fem=multiphysics(fem, ...

’sdl’,[]);

% Extend mesh

fem.xmesh=meshextend(fem,’geoms’,[1],’eqvars’,’on’,...

’cplbndeq’,’on’,’cplbndsh’,’off’);

% Solve problem

fem.sol=femnlin(fem, ...

’method’,’eliminate’, ...

’nullfun’,’auto’, ...

’blocksize’,5000, ...

’complexfun’,’off’, ...

’solfile’,’off’, ...

’conjugate’,’off’, ...

’symmetric’,’off’, ...

’solcomp’,{’w’,’u’,’p’,’v’}, ...

’outcomp’,{’w’,’u’,’p’,’v’}, ...

’rowscale’,’on’, ...

’ntol’,1.0E-6, ...

’maxiter’,25, ...

’hnlin’,’off’, ...

’linsolver’,’umfpack’, ...

’thresh’,0.1, ...

’umfalloc’,0.7, ...

’uscale’,’auto’, ...

’mcase’,0);

% Save current fem structure for restart purposes fem0=fem;

% Plot solution, surface plot of the z velocity figure(1)

postplot(fem, ...

’tridata’,{’v’,’cont’,’internal’}, ...

’triedgestyle’,’none’, ...

’trifacestyle’,’interp’, ...

A THE SCRIPT

’tribar’,’on’, ...

’trimap’,’jet(1024)’, ...

’solnum’,’end’, ...

’phase’,(0)*pi/180, ...

’title’,’Surface: Velocity field [m/s]’, ...

’refine’,3, ...

’geom’,’on’, ...

’geomnum’,[1], ...

’sdl’,{[1,2,3]}, ...

’axisvisible’,’on’, ...

’axisequal’,’on’, ...

’grid’,’off’);

% Cross-section plot figure(2)

postcrossplot(fem,1,[r_shaft (r_shaft+0.00003);0 0], ...

’lindata’,’v’, ...

’linstyle’,’-’, ...

’lincolor’,’cycle’, ...

’linmarker’,’none’, ...

’npoints’,1000, ...

’spacing’,[3e-5 2e-5 0 -2e-5], ...

’solnum’,’all’, ...

’phase’,(0)*pi/180, ...

’axislabel’,{’Arc-length’,’z-velocity [m/s]’}, ...

’axistype’,{’lin’,’lin’}, ...

’geomnum’,[1], ...

’transparency’,1.0);

legend(’z=3e-5’,’z=2e-5 m’,’z=1e-5 m’,’z=0 m’,’z=-1e-5 m’,...

’z=-2e-5 m’,’z=-1e-5 m’,’z=-2e-5 m’,’Location’,’SouthEast’);

title([’z-velocity [m/s] m=’,num2str(m),’, n=’,...

num2str(n) ’, r_{shaft}=’,num2str(r_shaft)])

% Construction of the max and min velocity in the z direction

% and the plot of the netflow thought the whole domain i=1;

zi=[];

zmax=[];

zmin=[];

I=[];

z_mini=-6e-5-z_langd_in;

z_maxi=6e-5+z_langd_ut;

for zii=z_mini:0.5e-6:z_maxi;

zi(i)=zii;

pd = postcrossplot(fem,1,[r_shaft (r_shaft+0.00003);zii zii], ...

A THE SCRIPT

’lindata’,’v’,’npoints’,1000,’outtype’,’postdata’);

I(i) = meshintegrate(pd.p);

zmax(i)=max(pd.p(2,:));

zmin(i)=min(pd.p(2,:));

i=i+1;

end figure(3) plot(zi,I)

title([’The integral over z-velocities m=’,num2str(m),...

’, n=’,num2str(n) ’, r_{shaft}=’,num2str(r_shaft)]) xlabel(’z [m]’)

ylabel(’netflow z-direction’) figure(4)

plot(zi,zmax,’b’, zi,zmin,’r’)

title([’max min velocities m=’,num2str(m),’, n=’,...

num2str(n) ’, r_{shaft}=’,num2str(r_shaft)]) xlabel(’z [m]’)

ylabel(’velocity z-direction [m/s]’) legend(’max’,’min’);