lubrication Piston-Ring

(1)

Master's Degree Thesis ISRN: BTH-AMT-EX--2012/D-12--SE

Supervisors: Professor Roland Larsson, LTU

Department of Mechanical Engineering Blekinge Institute of Technology

Karlskrona, Sweden 2012

Mohammad Shirzadegan

Simulation of rough contact

lubrication Piston-Ring

(2)

(3)

Simulation of rough contact lubrication Piston-Ring

Mohammad Shirzadegan

Department of Mechanical Engineering Blekinge Institute of Technology

Karlskrona, Sweden Spring 2012

(4)

1

Abstract

In everyday life machine have a strong role, therefore investigation in machine elements are desired for many researchers. Piston ring is one of important component of engine. In order to understand the lubricant oil consumption, friction, wear, the first step is to analyze a Piston-Ring behavior. For this reason a theoretical models are considered and behavior of pressure in several situations are studied. This study is aimed to survey steady state and transient analysis condition among the Piston Ring. Also behavior of surface feature (Texture) on the lower surface is be modeled in some time step.

In this thesis by using MATLAB, a suitable cavitation algorithm with an iterative method is established.

Keywords:

Piston-Ring, Cavitation, theoretical model, MATLAB, dent, steady state, Transient

(5)

2

Acknowledgement

I would like to thank Professor Roland Larsson who gave me this invaluable opportunity and supporting me all the time. It is really my pleasure to work under his supervision. Studying in machine element division of LTU is one of the moments that I will never forget in my life.

Special thanks to Andreas Almqvist, Assistance professor on machine element of Luleå Tekniska Universitet, who helps me in MATLAB programing and mathematical and physical concepts.

I also thank Dr.Ansel Berghuvud from Blekinge Tekniska Högskola (BTH) for his valuable support and guidance through my work. I have learned many things from him during two years studying at BTH.

Finally, I would like to express my thanks to my family (Amir, Jamileh, Shahin, Marjan) that supporting me without any hesitation.

Mohammad Shirzadegan Spring 2012, Luleå, Sweden

(6)

3

1 Notation

Real area of contact Acceleration

Integration constant Integration constant Elastic modulus

Elastic modulus Elastic modulus Axial friction force

Asperity contact function

Axial friction force based on asperity contact Axial Hydrodynamic friction force

Radial friction force at ring groove pivot Gravity constant, switch function

Height function

Minimum film thickness Integral kernel

Length

Piston-Ring mass Number of node

(9)

6 Hydrodynamic pressure

Cavitation pressure

Axial applied gas pressure force Radial gas pressure relief force

Gas pressure behind ring

Radial applied gas pressure relief force Radial curvature

Axial component of Hydrodynamic force Radial for upper surface

Radial for lower surface Crank radius

Crank radius Right hand side

Engine Speed Velocity Angular velocity

External load Pressure viscosity index Pressure viscosity coefficient Bulk modulus

Indexes Fractional contact

(10)

7 Deformation

Viscosity

Viscosity at ambient pressure Density

Interval

(11)

8

2 Introduction

The main function of lubricant is to decrease wear, friction, heat and to facilitate load support of moving surfaces. There are five types of lubrication namely, Hydrodynamic, Hydrostatic, Elastohydrodynamic, Boundary and Solid film.

When two surfaces separated by a thick film of lubricant , pressure and load capacity of the system can be calculated based on fluid dynamic’s law (No metal contact), the lubrication is Hydrodynamic. But if the system does not have any motion and lubricant is inserted into the system with sufficient pressure, the lubrication is Hydrostatic (air and water are the most common lubricants).

Lubrication is Elastohydrodynamic when lubricant between the surfaces reduces to the specific amount and contact deformation of surfaces is considerable. In this situation combination of contact mechanic and fluid dynamics law are considerable. More details will be discussed in next sections.

Whenever load increases or velocity decreases, lacking of surface area, increase in lubricant temperature that lead to reducing of viscosity, any of these may put full film thickness in critical situation. Then highest

asperities may be detached by few drops of lubricant. This is called Boundary lubrication.

The selves solid lubricant, such as graphite, is necessary for operation if the system supposed to work on extremely high temperature and the load carried by the asperities. Mineral oil couldn’t have effective efficiency in high temperature. The following figure shows the difference between Hydrodynamic area and boundary condition. Viscosity (µ) and velocity (N) have a direct relation with friction. This graph has been obtained by Mckee brothers during friction test for bearing.

(12)

9

2.1 Piston-Ring

The Piston is one of the main parts of the engine that transfer power from combustion chamber to crank shaft. When gases burnt at top of cylinder (combustion chamber) pressure pushes the piston to the downside. The Piston reciprocal action move crank shaft and crank shaft moves the other engine components. This mechanism is repeated until the engine turned off.

Piston is surrounded by cylinder. For some important reasons space between the cylinder and piston should be sealed. This couldn’t be

happened without ring. Main duties of rings are divided in three main parts:

a) Sealing combustion chamber to prevent the gas leak from piston circumstance

b) Heat transfer from high temperature (piston) to low temperature (cylinder wall)

c) Adjust the oil

Figure 2.1. “The variation of the coefficient of friction with 𝜇𝑁 𝑃” [46].

(13)

10

All rings generally can be classified into two types namely compression rings and oil control ring. Compression rings are used for sealing of

compressed gases and they sit at the top of the piston. Cross section of these types is rectangular, barrel or tapered shape. Gas pressure moves to the back side of the ring and forces it towards the cylinder liner.

Oil control ring is set at the bottom of piston. Their main duty is controlling the oil and distributes it on to the cylinder wall. “The scraped oil is

collected in the oil control ring groove and transported through the piston back to the crankcase.” [1]

Many factors can fail a system high pressure and high temperature system causes a critical environment for Piston-Ring. If the system operates under such condition it could lead to wear, friction, deformation and other bad effects on the system.

Figure 2.2. Schematic piston ring.

(14)

11

Material have important role in running condition. As was said one function of ring is heat transfer between piston and liner. Normally applicable

material for ring production is called Grey cast iron. It shows good

performance under starved and dried conditions. Also chromium coating is used to prevent corrosion and abrasion. Due to the engine performance different type of coating such as aluminum-titania, tungsten carbide is used for piston. For further information see Andersson et al [1].

Figure 2.3. Rectangular and tapered ring [43].

(15)

12

The main object of this work is to simulate a Piston-Ring and study the behavior of pressure distribution under different conditions. Also it is a desire to understand the effect of surface roughness by adding a small dent in a lower surface, and finally, to study if elastohydrodynamic effects play any significant role.

Elastohydrodynamic effects are not normally taken into account in studies of the piston ring problem. But that deformation can cause load and film thickness improvements.

This study is based on the theoretical model for this reason some numerical methods like Jacobi relaxation and Gauss-Seidel relaxation and SOR are introduced. Most of references reported that multigrid method is the ideal method for lubricant analysis but it has some sophistication in programing and needs more expertise.

An introduction to Elastohydrodynamic parameter such as viscosity-

pressure, density-pressure and deformation will be explained. Also step-by- step approaches for cavitation algorithm based on modify Elrod’s work is implemented. Pressure distribution of piston-Ring under two different boundary conditions (fully flooded and starved) is studied.

A piston ring is modeled with two different geometries. Model one does not have any groove and model two has small groove in lower surfaces. Steady state and transient condition will be surveyed for these two models.

Friction wear and oil consumption are always occurs in reality but these concepts are beyond the aim of current survey.

(16)

13

3 Literature Review

Back in 1880s an English railway engineer Beauchamp Tower was first to discover high pressure in full film regime [2]. He did some experimental test on railroad bearing and founded unexpectedly high pressure. After that, on 1886 professor Osborn Reynolds [3] published his famous theorem. This theory was extracted from the Navier-Stoks equation and determined the pressure distribution across the arbitrary geometry. Nowadays his work is well known as “Reynolds equation”.

In 1916 Martin [4] modeled the meshing of gear teeth. He considered the lubricant as isoviscos. He extracted the relation between operation

condition and film thickness which was far from the roughness of gear in reality. On that period it was difficult to calculate the elastic deformation and viscosity-pressure of lubricant at the same time [5].

The viscosity- pressure relation is obtained from the effort of Barus [6]. He deliberated the various aspect of viscosity of marine glue in various

temperature and pressure. He observed that in any temperature the rate of viscosity increases with pressure and he deduced the exponential relation between pressure and viscosity.

Petrusevich [7] in 1951 solved the Reynolds and deformation equations together. He considered highly loaded effects on EHL contact and established the pressure spike for the first time in his study.

EHL line contact was solved with different methods. Dowson and Higginson [8] solved the problem with inverse method and Newton- Raphson was used by Okamura [9]. Besides Hamrock and Jacobson [10]

using the Gauss-Seidel relaxation method in order to solve the low load EHL contact.

Lubrecht [11]solved the EHL line and point contact in 1984 with

combination of nonlinear Gauss-Seidel and multigrid technic. Lubrecht and

(17)

14

Brandt [12] solved the line contact EHL for the low loaded contact. He used the multi-integration algorithm technique and increased the solvation process for deformation [13].

Almqvist [14] applied the block Jacobi method as iterative solution for EHL line contact problem. He investigated multigrid technique and Jacobi method in his work.

One important issue of lubricant analysis is cavitation. Famous solutions for cavitation were started by Sommerfeld [15]. He presented a pressure distribution for a full film lubricant in journal bearing. As the film rupture was not taken into consideration the negative pressure distribution was shown in his work. Then in 1914 Gumble [16] by changing the boundary conditions (Half-Sommerfeld) of Reynolds equation gained better pressure distribution but the conservation of mass was not fulfilled in his work.

Then, Swift in 1932 [17] and Stieber [18]in 1933 defined other boundary conditions for journal bearing which were known as Reynolds boundary conditions. Based on that, the pressure is started to build up from the beginning of the dominant and will be disappeared where the pressure gradient is equal to zero. These conditions also have a mass conservation problem because the film reformation was not considered. This problem was discussed by Brewe et al [19] work.

First experiment that presented the vapor cavitation was occurred at Luleå University of technology in machine element laboratory. They did

experiment on bearing in the motion of PMMA tube and shaft. After that NASA Lewis search center did the same experiment and gained better result [19]. Floberg [20]also published his experimental work on cavitation that was occurred on bearing. Mass continuity was totally considered by work of Floberg, Jakobsson and Olsson for moving boundary conditions.

This worked is known as JFO method and the boundary condition are varied due to the time depended load. This method was difficult to implement.

In 1981 Elrod [21] introduced his famous algorithm that computes the cavitation in the bearing. He defined the switch function that made the

(18)

15

Reynolds equation valid among the bearing. There are two advantageous for this method. The first one is consideration of mass conservation among the cavitation region and the second is easiest programing due to cavitation complexity.

After that Vijarvaraghavan and Keith [22]developed the Elrod algorithm.

For the cavitation boundary the second part of Reynolds equation (shear induced flow) automatically changes from central to upwind difference.

This happened by applying an artificial viscosity function to shear induced term. They analyzed their algorithm on slider and journal bearing under a heavy loaded.

Ausas et al [23] prepared another numerical algorithm due to the mass conversation and Elrod’s work. They used the Newark’s scheme and relaxation process to update the fluid fraction and pressure.

Dowson et al [24] studied the tribological behavior of a piston ring in eight four- stroke and six two-stroke diesel engines. The effects of squeeze film and elastohydrodynamics at top dead center were established too.

Priest et al [25] presented a free body diagram for the compression ring. He computed hydrodynamic pressure for four different cavitation theories.

Also minimum film thickness behavior was discussed for these models.

The dynamic behavior of the piston ring in a diesel engine was investigated by Tian [26]. He studied the performance and effect of lubricant between top two rings for different crank angles. This experimental investigation shows the transportation of oil in three main regions of piston and cylinder.

Rahnejat et al [27] investigated the gas force action behind the piston ring.

Variation of film thickness and load at different crank angle was discussed.

Rahnejat believed that beside deformation, asperity adhesion should be added to model to understand the behavior of piston at TDC position.

Dellis and Arcoumanis [28] developed an experimental reciprocating rig test in order to inspect cavitation and behavior of film thickness. They

(19)

16

discussed about the shape of pressure based on film thickness.

Development of cavitation in different load is shown in details.

Spencer et al [29] studied the surface texture for combustion engine experimentally and numerically. He determined the pressure distribution under the ring with Homogenized technique. He modeled an artificial texture on a piston ring that could be used in for real model. To deal with cavitation, Vijarvaraghavan algorithm was implemented in his work.

Andersson [1] compiled useful information about design, material, product, wear and friction in the Piston Ring-Cylinder liner. In his literature survive nearly 150 references were used.

(20)

17

4 Lubrication-Governing Equations

When the pressure distribution in lubricant reaches high enough level like 0.1 to 3 GPa the shape of film thickness changes due to the surface deformation. Then we can speak about Elastohydrodynamic lubrication (EHL). If we want to mention two major characters of EHL, they can be elastic deformation and piezo viscous effects. The elastic deformation was studied by Hertz in 1881. He considers a contact between two spherical bodies.

Mostly, EHL consider in non-conformal contacts which are line contact (cylinder-cylinder), elliptical contact (ball-cylinder) and circular contact (ball-ball).

The EHL problem is governed by various equations. Viscosity and density are function of pressure. As was said deformation equation must be added to the height function in Reynolds equation. Force balance is a parameter that should not be forgotten that balance an applied load and pressure.

In order to simplify the contact for example between two disks, one of the disks is replaced with plane in the x coordinate. The figure 4.1 shows this procedure [30].

Figure 4.1. Schematic simulation of contact [30].

(21)

18

For this system a reduced radius R is defined as below

( 4.1)

In order to calculate the parabolic shape, if considered as minimum film thickness an approximation for film thickness can be estimated as below

( ( ))( )

( 4.2) Which can be simplified to:

( 4.3) If

4.1 Film thickness equation

By adding deformation to the film thickness equation the final relation can be determined. A parabola carve could be plotted by this equation.

( ) ( 4.4)

Where deformation can be found from:

( )

∫ ( )

Extraction of deflection will be described briefly at next section.

(22)

19 4.2 Pressure-viscosity

In EHL two famous pressure-viscosity relationships are applicable. The first one is given by Barus [6]. He introduced a simple equation that describes behavior of viscosity under high pressure.

( 4.5)

Where is atmospheric viscosity and is a coefficient for pressure viscosity which is varies between for oils.

The second equation that is more complicated than Barus was introduced by Roelands in 1966 [31]

( ) ( ( ) )( ( ) )) ( 4.6) Where are constant. Comparison between these two equations as a function of pressure is plotted on figure 4.2.

Figure 4.2. Comparison between Barus and Roelands method [13].

(23)

20

Neither Barus nor Reolands are very accurate above 500 MPa. Free volume models should be applied in such cases. The Barus equation has been used in this study since pressure level is low.

4.3 Density-Pressure

Dowson and Higginson [8] defined a relation for density that varies due to the pressure. It reads

( ) (

) ( 4.7)

Where is atmospheric density. In this equation pressure is given in Pascal unit.

Figure 4.3. Relation between density and pressure due to Dowson [13].

(24)

21 4.4 Force Balance

The fluid pressure must balance externally applied load:

∫ ( )

( 4.8)

To satisfy relation 4.9 some numerical methods are suggested. The famous one is Secant method that is easy to implement. An algorithm in Matlab can be found in appendix 13.2. A bisection method is another mathematical root finding that could be applicable.

(25)

22

5 Numerical Approaches

In this chapter finite difference method which is used to solve the Reynolds equation and three algebraic techniques will be explained. In a boundary value problem suppose we have a second order equation which is generally describe by following formula

( ) ( ) ( ) ( ) ( )

( 5.1) That

( ) ( ) ( ) Let divide the interval [a,b] to N equal space therefore

For higher order equation it is essential to use this approximation ( ) (

)( ( ) ( ) ( ) ( 5.2) And using central, backward or forward difference for first order equation ( ) (

)( ( ) ( )) ( 5.3)

( ) ( )( ( ) ( )) ( 5.4)

( ) ( )( ( ) ( )) ( 5.5)

Now we can discretize Reynolds equation due to above procedure. It reads

(

) ( )

( 5.6)

(26)

23

Define

( )

and

( )

Substitute (e) and (RHS) into Reynolds equation

(

) ( )

( 5.7)

And extracting in a finite difference form

( )

( ) (

) ( ) (

)

( 5.8)

Central difference for right hand side of Reynolds equation

( 5.9)

This equation should only use for interior point.

( 5.10) For i=1

( ) ( ) ( )

( 5.11) For i=N

( ) ( ) ( )

( 5.12)

(27)

24

Rewriting the equation in a matrix form leads to following diagonal system:

Originally, Reynolds equation has a squeeze term that could estimate the film thickness. Since the film history adds to the simulation value of film thickness, location of cavitation and pressure distribution are getting changed.

(

) ( )

( )

To fulfill initial amount of film thickness and pressure are guessed and also initial elastic deformation is estimated. After that, the pressure of Piston- Ring is calculated by implementation of cavitation algorithm. Convergence criteria should have satisfied. Next second, bisection and inverse quadratic interpolation methods are applied to balance the load behind the ring [29].

The processes are repeated for every time increment. Results for each time step kept as history for next one. The summery of implementation of the model into software for transient condition is shown figure (5.2)

Figure 5.1. System of equation.

(28)

25 Start

Initial guess of height function

Apply pressure to calculate the

Estimate the Elastic deformation

Implement the cavitation algorithm and determine pressure distribution

Use h function of every cycle for Does P converge?

Force Balance

Adjust

𝑜

First cycle Use h for

previous time step & go to next time step

Yes

No

Print pressure distribution for every time step

End Yes

Figure 5.2. Algorithm for Transient analysis.

(29)

26

To solve this algebraic equation an iterative solution is necessary. Jacobi and Gauss-Seidel relaxation and successive over relaxation method are discussed briefly.

5.1 Jacobi method

This method needs an initial guess value for start the iteration loop. Then it produces several sequences that converge to initial guess. If we have equation system [A]{x}={b}, this technique changes the system to the

.

An algorithm for solving algebraic system with Jacobi method implement as follow.

Suppose A is formed by three diagonal U (upper triangular), L (lower triangular), D (diagonal A).

( 5.13) Each loop is identical to solve for every variable one. Now equation Ax=b can be changed to

( )

( 5.14)

( )

( 5.15) That ( )

This method can program to Matlab by following orders. This method is converged vary slowly.

(30)

27

Due to Venner and Lubrect [32] Jacobi relaxation in matrix form is:

( 5.16) ( )

5.2 Gauss-Seidel Method

This method is popular technique for solving an algebraic system of

equation. In contrast to Jacobi method that updates values at end of iteration this method uses the new approximation for solving the later one. This iteration method is diverged faster than Jacobi.

An algorithm for matrix form is:

( ) ( )

( 5.17) Figure 5.3. Matlab algorithm for Jacobi method.

(31)

28

It can be programed to Matlab by following orders.

5.3 Successive over relaxation

SOR help us to gain convergence faster than Gauss-Seidel and extracted by extrapolating Gauss-Seidel method.

( ) ( 5.18)

In a matrix form it is rewritten as:

( )( ( ) ( )

( 5.20)

If w is set to one the equation changes to the Gauss-Seidel method and If w<1 the system is under relaxed.

Figure 5.4. Matlab algorithm for G-S method

(32)

29 5.4 Deformation scheme

There are different approaches to calculate the deflection and film

thickness. The elastic deformation for line contact problem is obtained from the following integral.

( )

∫ ( )

(5.21)

Okamara [33] discretized the deformation equation by the following formula in dimensionless form. He used this equation throw to numerical analysis of isothermal elastohydrodynamic lubrication.

̌

∑ (|

| ) (|

|)

(5.22)

Houpert and Hamrock [34] proposed a faster approach to deformation.

They improved Okumara’s method and determined the integral inside the equation analytically. The baseline in this assumption were pressure distinguished by a polynomial of second degree in the interval [ ]

̌

∑ (

) ( 5.19)

Almqvist [14] used alternative approach which differentiate

̌ ∑ with respect to . Due to his work deformation could be calculated as

(33)

30 ̌( ) ∑

(5.24)

Where

-

Evans and Hughes [35] introduced differential deflection method that calculated the deformation of the system for any arbitrary pressure distribution among the lubricant. They assumed that for calculation of deflection numerically, it could be written in a quadrature form.

( ) ∑

( 5.25)

Where is a weighting function based on pressure under the area of integration.

By two times differentiation respect to x, the equation (4.14) can be written in the following form. (For the mathematical procedure see appendix 13.1)

∑( )

( 5.20) Which rewritten to the form

∑

( 5.21) Where

If we want to rearrange the equation (4.16) for any mesh point

(34)

31 ( )

∑

( 5.22) Equation 4.17 could calculate the deformation. This equation needs two boundary conditions that are important to implement correctly.

5.4.1 Direct integral

Due to singular kernel of evaluation of this integral seems difficult. Integral term of deflection equation can be determined by

∫

(5.23)

In the interval

Analytical evaluation of kernel (k) - mostly used in direct integral method - is extracted in the following equation.

( ) ( | | ) (

)( | | ) ( 5.30)

Therefore deformation can be approximated as

( ) ∑ ( )

( )

( 5.31) Where P is approximation for the pressure and K is called kernel

coefficient.

(35)

32

Figure 5.5. Effect of deformation on Height function.

(36)

33

6 Cavitation

If the pressure inside a lubricant becomes lower than the gage/ambient value, the liquid film couldn’t withstand the high tensile stress. Then the fluid film will break-up and cavitation occurs.

Four major category of cavitation can be summarized as hydrodynamic, Acoustic, optic and particle cavitation. Hydrodynamic cavitation takes place due to the variation of velocity. It can be subdivided into four groups namely, travelling, fix, vortex and vibratory cavitation.

- Travelling cavitation happens when the bubble starts to grow from inside the liquid, become larger and afterward collapse.

- Fix cavitation occurs when the liquid flow separates from the rigid boundary of dunk body and the cavity remain fix to the boundary.

- Vortex cavitation take place in cores of vortices which form in areas of high shear.(i.e. blade of ship)

- At vibratory cavitation due to the low velocity of lubricant a given element of liquid is exposed to many cycles of cavitation instead of one. [36]

6.1 Cavitation algorithm

In 1981 Elrod [21] introduced a switch function that computes the

cavitation automatically among the lubricant regime. The regime is divided into two parts, full film and cavity. Reynolds equation is valid only in the complete film and in the cavity regime the equation needs certain changes.

These two regimes are named as coquette flow and Shear flow. Due to Elrod, “universal” differential equation is contained fractional content ( ) which help to calculate the pressure distribution among the lubricant.

(37)

34

Fractional content or nondimensional density (in the full film regime) is defined as and the pressure can be obtained from ( ) where is cavitation pressure and is a bulk modulus.

Based on the variation of through the regime cavitation index g is defined. If cavitation index is zero which means the flow are in cavity zone. So if cavitation index will set as zero.

On X direction the mass conservation for the upwind stream and downwind stream when the velocity set positive is written as below

( ̇ ) ( ( )

( ) ( )) ( 6.1)

( ̇ ) (

) ( ( )

( ( )) ( 6.2)

Due to above equations (6.1-2), in the full film zone g =1 for parabolic system second order –central is applied. In the cavitation zone the pressure equation is omitted because g is equal to zero and upwind difference is used for connective term.

In order to discretize the Elrod algorithm Finite Difference is adopted for convective and pressure term.

Switch function

(38)

35 { [ ( ) ]

[

]}

{

[

]}

{ [( ) ]

[ (

) ]}

[ ( ) ( ) ( ) ]

[

(

)

]

( 6.3)

Vijarvaraghavan and Keith [22] modified the Elrod algorithm. This made it less time-consuming and the accuracy of the result is improved with a coarser grid. In the numerical formulation an artificial viscosity term is add to the shear flow term where is equal to zero in the full film zone and in the cavity zone it is equal to one. Artificial viscosity has the following relation with switch function g:

The one dimensional steady state Reynolds equation should be changed as follow

(

)

( 6.4) To discretize the equation in X direction, we can rewrite it in a following order. This helps us to fill in matrices easily. For two dimensional transient Reynolds equations see for example Spencer et al [29]

( 6.5) Where

[ ]

(39)

36 [ ]

[ ] [ ]

[( ) ]

[ ( ) ]

And coefficients are introduced as

(

)

To study cavitation algorithm a parabolic geometry is modeled. A slider parabolic bearing which has been modeled by Vijarvaraghavan [22]and Sahlin et al [37] is applied. The geometry of the slider bearing is shown below. The minimum film thickness is constant and boundary conditions are considered as fully flooded film.

The height function is defined with following equation Figure 6.1. Parabolic shape [37].

(40)

37 ( ) ( )

( ) ( 6.6)

Input parameters that used for the model is shown on table 2.

Input parameters

( m) 2.54e-5

( m) 7.62e-2

( m/s) 4.57

( pa .s) 0.039

(N/ ) 6.9e7

( ) 0

Figure 6.2. Pressure distribution for fully flooded boundary condition.

Table 6.1. Input parameter for slider bearing.

Figure 6.3. Pressure distribution for starved boundary condition.

(41)

38

7 Theoretical solutions

Here just the theoretical equation for the piston ring relation is extracted.

But for analyzing the system numerical method that explained in chapter 6 is used. Simplified mathematical model for piston ring is calculated.

Several assumptions have been made. In comparison with gas pressure, ring weight and friction are small and they set to zero. Also at z direction

friction force which made by groove is neglected. System assumed to be fully flooded. is Axial component of Hydrodynamic force and is Axial applied gas pressure force [25]

( 7.1)

( 7.2)

( 7.3)

( 7.4)

Where is considered as radial gas pressure relief force and is radial applied gas pressure relief force. is radial force due to asperity contact

(42)

39

and is radial component of hydrodynamic force. All the parameter and their direction are shown in following figure.

A model that presented by Greenwood and Tripp [38]was applied to determined . Approach proposed by Ruddy et al should be implemented.

Then [25]

√

( ) √ ∫ ( )

( 7.5)

Figure 7.1. Piston Ring [25].

(43)

40

Which is composite root mean square surface roughness of the piston ring and cylinder bore and is asperity radius of curvature

When the hydrodynamic pressure integrated among the inlet and outlet, radial hydrodynamic load will be obtained.

∫

( 7.6)

And hydrodynamic pressure can be obtained from Reynolds equation (by two times integration)

∫

∫ ∫ _{( 7.7)}

Which c1 and c2 will be calculated based on boundary conations. For more details about calculating cavitation and friction see Priest [25]

(44)

41

8 The model problem

To implement all information into our differential equation, two models are introduced. The first one is piston-ring in a steady state condition. Pressure distributions for three different crank angles -20, -10,-2 are extracted. It is assumed that the surface is totally smooth and it doesn’t have any

roughness. Boundary conditions are variable during the analysis due to the gas pressure behind the piston-ring and they are always set as fully flooded.

Effects of temperature on viscosity and other parameters are neglected.

All the geometry are defined in table 8.1

The force that applied behind the ring is obtained from ring tension (T) and combustion gas pressure. It has different values on different crank angle.

The maximum value of gas pressure is near TDC. According to the figure (8.2) this amount is close to 1.53 Mpa. Relation between ring tension and gas pressure with force can be determined by equation (8.1) [39].

( )

( 8.1) Figure 8.1. Schematic Piston-Ring, first model.

(45)

42

Parameter value

X-start (m) -7.375×1

X-end (m) 7.375×1

Ring tension T (MPa) 0.341

Bulk modulus (Pa) 1.72×1

Radius R (m) 0.0183

Cavitation pressure (MPa) 0.02

Bore diameter D (mm) 0.889

Viscosity coefficient (Pa-1) 1.8

Due to the cavitation algorithm two boundary conditions are demanded as function of fraction content( ). Both of them can be estimated by substituting pressure from top ring ( ) and pressure from the blow-by( ).

According to Yang et al [40]

( 8.2)

( 8.3)

and are constant and equal to . Figure 8.2. Combustion gas pressure [39].

Table 8.1. Input parameter for first model.

(46)

43

Whenever the piston move in every crank angle velocity of piston is changed therefore a relation needed to define this movement. It can be estimated by following equation [40].

( √ ( )

)

( 8.4) Where and n is an engine speed, R is a crank radius and L is connecting rod length.

The second model has the same geometry as the first one but lower surface has a small dent on it. Contrary to reality it is assumed that only one surface of the Piston-Ring considered rough. In order to study the behavior of surface feature, a mathematical statement should add to height function which is [41]

( ( ( )

) ) ( ( ( )

)) (‎8.5)

Where is the amplitude, is the center of dent which is varies in each time step and is the dent wavelength.

(47)

44

In the second model input data for gas pressure in every crank angle is needed here the load and velocity are assumed to be constant.

Characteristics value

Load (Mpa) 5.5

Amplitude (µm) 0.8

Wavelength (mm) 0.1

Time (s) 0.000926

Velocity(m/s) 2

Figure 8.3. Schematic view of model 2.

Table 8.2. Input parameter for second model.

(48)

45

9 Results

9.1 Pressure distribution among Piston-Ring model 1

Pressure distribution for different crank angle is calculated. At the velocity is set to and maximum pressure reaches to near . When engine goes to gas pressure behind the ring make the contact thinner also velocity reduces to and maximum pressure

reaches . Before TDC ( ) maximum pressure increases rapidly.

The amount reaches to

Whenever crank angle get close to TDC an applied load behind the ring increased, velocity decreased to and deformation influenced more on pressure build up. Maximum pressure value is close to

Figure 9.1. Pressure distribution- 20 and -10 crank angle.

(49)

46

It can be seen that position of the cavitation is shifted to the center. Also reformation in a divergent part is started to grow. See figure (9.1 and 9.2) Film thickness for three different crank angles is plotted. See figure (9.3).

When the Piston moves to the TDC position the value of film thickness reaches its minimum value.

Figure 9.2. Pressure distribution -2 and TDC.

(50)

47

Figure 9.3. Film thickness -20,-10 and TDC.

(51)

48

9.2 Pressure distribution among piston-ring model 2

In the analysis of small dent the load and velocity are kept constant.

Behavior of dent during in four different steps is shown in figure 9.4.

There is no effect of dent on figure 9.4.a. Then it moves through to system.

Because of the parabolic shape of height function in the beginning the pressure decreases but after a while it starts to build up and after that when the dent reaches to the middle of height function once again the pressure drops and gains its minimum values.

Finally, when the dent passes the middle the second peak pressure reduces until it disappears from sight. See figure 9.3.c and 9.3.d. behavior of

Figure 9.4. Pressure distribution for model 2 in different time step.

(a) (b)

(c) (d)

(52)

49

pressure during this action doesn’t change even when the amount of load varies.

9.3 Effect of surface features

In this section maximum pressure and minimum height function in every time step are plotted. Effects of surface feature that introduced by equation (8.5) in details has been plotted. In order to get converged results the load behind of Piston-Ring is set to 6.5×1 and the velocity 2 m/s.

Figure 9.5. Maximum pressure in each step.

Figure 9.6. Minimum height in each step.

Figure 9.7. Parabola shape when the dent is in the middle.

Figure 9.8. Pressure distribution when the dent is in the middle.

(53)

50

When the dent reaches to the middle of parabola shape as was expected the pressure started to decrease. While center of dent reaches to diverge part of parabola once again the pressure start to build up. This is shown in figure (9.7).

(54)

51 9.4 Effect of deformation

In order to see what will happen if the deformation is omitted from the analysis of the system a high constant load (5 MPa) is applied to the system. After one complete engine cycle the results show how the amount of pressure varies. These changes are become more important whenever the system have roughness (dent) see figure (9.9, 9.10).

Figure 9.10. Effect of deformation on dent in high pressure.

Figure 9.9. Pressure distribution with and without deformation.

(a) (b)

(55)

52

10 Discussions and Conclusions

- The results shows EHL have effects on film thickness and perhaps pressure when the load is high. These functions such as viscosity, deformation, density improve the value of load capacity and therefore it could be considerable for further calculation like

friction, wear or power loss. Elastohydrodynamic effect changes the peak pressure. Comparison between pressure for system with and without deformation shows that the value of pressure peak is decreases due to implementation of deformation. By increasing amount of applied pressure (load) to the system deformation have active role in results. As can be seen in figures (9.9) for constant load 5.5e6 MPa and velocity 2 m/s the peak of pressure is changed.

Also for the same position of dent in the second model the shape of pressure peak is different.

- Pressure reformation is started to grow up as long as the piston ring get close to the TDC. As was seen on fig (8.1) in it is small but near the TDC as the velocity decreases and gas pressure increases, the reformation started from . Another reason for this phoneme is elastohydrodynamic effects.

- It is necessary to satisfy the force balance. It has direct influence on squeeze film-term (). On the steady state condition, due to

elimination of squeeze term pressure peak is goes high near the top dead center.

- Transient analysis of the system may differ from steady state.

Position of cavitation, peak pressure and minimum film thickness are varying due to the squeeze-term.

- Investigation of the surface features shows that the effects of texture may not be neglected and even a small dent into system cause changes in output. This texture may increase the friction and wear and the life time of the machine will be reduced.

(56)

53

- Due to the iterative procedure for Piston-Ring analysis the initial value of some parameter like height function should be guessed. For gaining suitable results it is necessary to run the system at least for one engine cycle.

- Behavior of iteration method like Jacobi and G-S shows that for the high loads system these methods are totally unstable and they never converged.

- It can be mention that in higher load, the pressure spike will be added to the pressure distribution among the system.

- In order to analyze the extracted results from mathematic view coefficient matrix A is considered. A Is a square matrix which is sparse and band matrix with a constant bandwidth of three (see fig (10.1)).

At first the eigenvalue of the matrix plotted. In the spectrum of the

eigenvalues it is observed that most of the values are clustered around zero.

Furthermore, the large distance between largest and smallest eigenvalues concludes that the condition number of the analyzed matrix should to be large. This reveals that the coefficient matrix is an ill-conditioned matrix which means that it is close to be singular. Therefore, conventional direct strategies will not guarantee to maintain the desired accuracy of the results due to the sensitivity of the system because of unavoidable rounding errors.

Alternatively, an iterative solver has been employed in order to solve the equation of this ill-condition system, i.e. condition number=1.547e9.

Figure 10.1. Diagonal of matrix A. this figure is zoomed

(57)

54

In the further analysis the singular values of the coefficient matrix is studied. Fig. 9.3 demonstrates the singular values of the coefficient matrix.

Obviously the values decade slowly as expected, and no large and detectable gaps are observed. Hence, the evaluated system is not rank deficient and uniqueness of the answer of the system is ensured.

Figure 10.2. Spectrum of the eigenvalues of the system.

Figure 10.3. Spectrum of the singular values of the system.

(58)

55

11 Future works

Effect of wear and friction will be discussed. Investigation in this area seems challenging. The iteration solution shows that they are too slow and the rate of convergence for high pressure is too low. Multigrid techniques will be studied in order to spread solution to the high pressure area and could be able to analyze the pressure spike and surface texture at same time. This method also has more advantageous simulation speed and accuracy.

It is interesting to adopt another cavitation algorithm and compare with modified Elrod’s algorithm. FFT method for determining the deformation will be investigated. This method is commonly used in LTU.

And finally, In order to solve Reynolds equation and evaluate film thickness simultaneity coupled method will be investigated.

(59)

56

12 Bibliography

[1] Andersson,P. ,Tamminen, J., Sandström, C.E, Piston ring tribology A literature survey, Otamedia Oy, Espoo: JULKAISIJA – UTGIVARE – PUBLISHER, VTT 2002.

[2] Tower,B, "First report on friction experiment (friction of lubricated bearing)," proceeding of Inistitute of Mechanical Engineering (London), pp. 632-659, 1883.

[3] Reynolds,O., "On the theory of lubrication and its application to Mr Beauchamp Tower’s experiments, including an experimental determination of the viscosity of olive oil.," Philosophical

Transactions of the Royal Society of London, vol. 177, pp. 157-234, 1886.

[4] Martin,H.M., "Lubrication of gear teeth," Engineering (London), vol.

102, p. 199, 1916.

[5] Goodyer,C.E., "Adaptive Numerical Methods for Elastohydrodynamic Lubrication," The University of Leeds School of Computing, Leeds, 2001.

[6] Barus,C., "Isothermals, isopiestics and isometrics relative to viscosity,"

American Journal of Science, vol. 45, p. 87–96, 1893.

[7] Petrusevich, A.I., "fundamental conclusion from contact hydrodynamic theory of lubrication," Izu. Akad. Nauk SSR (OTN), vol. 2, pp. 209-223, 1951.

[8] Dowson, D. and Higginson, G.R., "A numerical solution of the elastohydrodynamic problem," Journal of Mechanical Engineering Science, vol. 1, no. 1, pp. 6–15,, 1959.

(60)

57

[9] Okamura,H, "A contribution to the numerical analysis of isothermal elastohydrodynamic lubrication," Proceedings of the 9th Leeds-Lyon Symposium on tribology, vol. UK, p. leeds, 1982.

[10] Hamrock, B.J and Jacobson, B.O., "Elastohydrodynamic lubrication of line contacts," ASLE Transactions,, vol. 4, no. 27, p. 275–287, 1984..

[11] Lubrecht,A.A., Napel,W., Bosma,R., "An alternative method for calculating film thickness and pressure profiles in

elastohydrodynamically lubricated line contacts," Journal of Tribology.

ASME, vol. 108, pp. 551-556, 1986.

[12] Brandt,A., Lubrecht,A.A.,, "Multilevel matrix multiplication and fast solution of integral equations," Journal of Computational Physics, vol.

90(2), pp. 348-370, 1990.

[13] Zargari, E. A., "Computational Analysis of Integral and Differential Formulations of the Elastohydrodynamic Lubrication Film Thickness Equation," University of Leeds school of computing, Leeds, 2007.

[14] Almqvist, A., A coupled solution method to the EHL line contact problem, luleå: universitetryckeriet,luleå,ISSN:1402-1617, 2001.

[15] Sommerfeld, A., "Zur hydrodynamische theorie der

schmiermittelreibung," Zeit Math. Phys, vol. 50, p. 97–155, 1904.

[16] Gumbel, L., "Das Problem der Lagerreibung," Mon.

Berl.Bezirksverein.,V.D.I., p. 109–120, 1914.

[17] Swift, H. W., "The stability of lubricating films in journal bearings,"

Proc. Inst. Civil Engnrs (Lond.), vol. 233, p. 267–288, 1932.

[18] Stieber, W., "Das schwimmlager:Hydrodynamische Theorie," des Gleitlagers, no. (V.D.I. Verlag GMBH, Berlin), p. 106, 1933.

[19] Brewe, D. E., Ball, J. H. and Khonsari, M. M., "introduction Part 2:

(61)

58

current research in cavitating fluid films.Proceedings of the Cavitation Symposium," STLE Annual Meeting, NASA, Vols. TM-103184, p. 25–

26, 1988.

[20] Floberg, L., "Sub-cavity pressures and number of oil streamers in cavitation regions with special reference to the infinite journal bearing," Acta Polytech. Scand. Mech.Engng Ser., vol. 37, p. 1–31., 1968.

[21] Elrod, H. G., "A Cawitation Algorithm," Journal of Lubrication Technology, vol. 103, pp. 350-354, 1981.

[22] Vijayaraghavan, D., and Keith, T. G., "Development and Evaluation of Cavitation Algorithm," Tribol. Trans, vol. 32, p. 225–233., 1989.

[23] Ausas, R.F., Jai, M. and Buscaglia, G.C., "A Mass-Conserving Algorithm for Dynamical Lubrication Problems With Cavitation,"

Journal of Tribology ASME, vol. 131, 2009.

[24] Dowson,D. Ruddy,B.L., and Econimou,P.L., "The

ElastoHydrodynamic lubrication of Piston Ring," Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 386, pp. 409-430, 1983.

[25] Priest,M.,Dowson,D. and Taylor,C.M., "Theoritical modeling of cavitation in piston ring lubrication," Proceeding of the institution of mechanical engineers, vol. 214, no. part c, pp. 435-447, 2000.

[26] Tian,T, "Dynamic behaviours of piston rings and their practical impact.

Part 2: oil transport, friction and wear of ring/ liner interface and the effects of piston and ring dynamics," Proc Instn Mech Engrs

Engineering Tribology, vol. 216, no. part J, pp. 229-246, 2002.

[27] Rahnejat,H.,Mishra,P.C, and Balakrishnan,S, "Tribology of compression ring-to-cylinder contact at reversal," Proc. IMechE,

(62)

59

Engineering Tribology, vol. 222, no. part J, pp. 815-826, 2008.

[28] Dellis,P. and Arcoumanis, C., "Cavitation development in the lubricant film of a reciprocating piston–ring assembly," Proc. Instn Mech, Journal of Engineering Tribology, vol. Vol. 218, no. Part J, pp. 157- 171, 2004.

[29] Spencer,A., Almqvist,A. and Larsson,R., "Effect of honing angle on hydrodynamic lubrication," Proceedings of the Institution of

Mechanical Engineers, Part J: Journal of Engineering Tribology, vol.

225, pp. 683-689, 2011.

[30] Gohar,R, Elastohydrodynamics, Chichester: Ellis Horwood Limited, 1988.

[31] Roelands,C.J.A., corelation aspectof viscosity-temperature-pressure relationship of lubrication oils, Groning,Netherlands: Ph.D Thesis , Technical university Delft,Delft,Netherlands, 1966.

[32] Venner, C.H., Lubrecht, A.A., Multi level method in lubrication, Amsterdam: Elsevier Science B.V, 2000.

[33] Okamura,H., "A contribution to the numerical analysis of isothermal elastohydrodynamic lubrication," In Tribology of Reciprocating Engines, Proceedings of the 9th Leeds-Lyon Symposium on Tribology, Leeds, Butterworth, London, p. 313–320, 1982.

[34] Houpert L.G., Hamrock B.J., "Fast Approach for Calculating Film Thicknesses and Pressures in Elastohydrodynamically Lubricated Contacts at High Loads," Journal of Tribology ASME, vol. 108/411, pp. 411-419, 1986.

[35] Evans,H.P. and Hughes,T.G., "Evaluation of deflection in semi-infinite bodies by a differential method," Proceeding of the institution of mechanical engineers, vol. 214, pp. 563-584, 2000.

(63)

60

[36] Shah, Y.T., Pandit ,A.B., Moholkar, V.S., Cavitation Reaction Engineering (The Plenum Chemical Engineering Series), New york:

Kluwer Academic/Plenum publisher, 1999.

[37] Sahlin, F., Almqvist, A.,Larsson, R. and Glavatskih R, "A cavitation algorithm for arbitrary lubricant compressibility," Tribology

International 40, pp. 1294-1300, 2007.

[38] Greenwood,J.A. and Tripp,J.H., "The contact of two nominally flat rough surface.," in proceeding of 9th Leeds-Lyon symposium on tibology: Tribology of reciprocating engines, vol. V(i), pp. 109-121, 1982.

[39] Chong,W.W.F, Teodorescu,M., Vaughan,N.D, "Cavitation induced starvation for piston-ring/liner trobological conjunction," Tribology International 44, pp. 483-497, 2011.

[40] Yang,Q.,Theo,G. and Keith,J.R, "Elastohydrodynamic cavitation algorithm for piston ring lubrication," Tribology transcations, vol. 38, no. 1, pp. 97-107, 1995.

[41] Almqvist,T.,Almqvist,A. and Larsson,R, "A comparision between computational fluid dynamic and Reynolds approch dor simulating transient EHL line contact," Tribology international, vol. 37, pp. 61-69, 2004.

[42] Maass, H. and Klier, H. , Kräfte,Momente und deren Ausgleich in der Verbrennungskraftmaschine (Forces, moments and their equilisation in internalcombustion engines / in German)., Vienna, Austria, Springer- Verlag, 422 p. ISBN 3-211-81677-1.: Die Verbrennungskraftmaschine, Neue Folge, Band, 1981.

[43] "Indo German industries (diesel engine)," Indo German industry, 2012.

[Online]. Available: http://indogermanind.com/compressionrings.html.

[Accessed 6 6 2012].

(64)

61

[44] Dowson, D. and Hamrock,B. J., "Numerical evaluation of the surface deformation of elastic solids subjected to a hertzian contact stress,"

International Lubrication Conference cosponsored by the Japanese Society of Lubrication Engineers and the American Society of Lubrication Engineers, 1975.

[45] Huttner, T., "Patentstorm, US Patent Application 20110132079 - FIXTURE FOR MEASURING THE DRAG OF A PISTON RING IN AN ENGINE CYLINDER AND METHODS OF USE THEREOF,"

PatentStorm LLC, 9 june 2011. [Online]. Available:

http://www.patentstorm.us/applications/20110132079/description.html.

[Accessed 6 june 2012].

[46] Shigley, J.E., Shigley's Mechanical Engineering Design, Newyork:

McGraw Hill Series in Mechanical Engineering, 2008.

[47] Woods, C.M., Brewe, E.M., "The Solution of the Elrod Algorithm for a Dynamically Loaded Journal Bearing Using Multigrid Techniques," in Tribology Conference cosponsored by the ASME and STLE,

Baltimore,Maryland, 1988.

(65)

62

13 Appendix

13.1 Evans and Hughes coefficient

Quadrature coefficient f split to three separate equations which is

[

( ) (

)]

[

( ) (

)]

[

( ) (

)]

To solve this equation a finite difference method is applied and extracted the left hand side.

( )

Two boundary conditions are needed to evaluate the above equation. There are two methods available for this reason. It can be substituted by direct integration method or by taking an arbitrary constant and use the following equation. [35]

( ) ( ) ∑( )

(66)

63 13.2 Secant Method

One of the numerical ways to estimate the force balance condition is secant method. In this method the root of function is estimated by line tangent of carve with two points. It is not necessary for the start and end points have different signs.

The secant method is define by the following equation For i=1,2,…max iteration

( )[

( ) ( )]

In Matlab the following file can calculate the root of equation with arbitrary tolerance and iteration. It just need too input function and interval.

function secant(f,x1,x2,tol,j) itr=0;

g1=feval(f,x1);

g2=feval(f,x2);

err=abs(x2-x1);

disp('______________________________________________________________

')

disp('itr xj f(xj) f(xj+1)-f(xj) |xj+1-xj|') disp('______________________________________________________________

')

fprintf('%2.0f %12.6f %12.6f\n',iter,x0,u)

fprintf('%2.0f %12.6f %12.6f %12.6f %12.6f\n',itr,x1,g2,g2- g1,err)

while (err>tol)&(itr<=j)&((g2-g1)~=0) x=x2-g2*(x2-x1)/(g2-g1);

x1=x2;

g1=g2;

x2=x;

g2=feval(f,x2);

err=abs(x2-x1);

itr=itr+1;

fprintf('%2.0f %12.6f %12.6f %12.6f

%12.6f\n',itr,x2,g2,g2-g1,err) end

if ((g2-g1)==0) disp(' NAN') end

if (itr>j)

disp(' Convergency problem ') end

(67)

lubrication Piston-Ring

Mohammad Shirzadegan

Simulation of rough contact

lubrication Piston-Ring

Simulation of rough contact lubrication Piston-Ring

Mohammad Shirzadegan

Abstract

Keywords:

Acknowledgement

Contents

1 Notation

2 Introduction

2.1 Piston-Ring

3 Literature Review

4 Lubrication-Governing Equations

5 Numerical Approaches

( )

( )

-

6 Cavitation

7 Theoretical solutions

8 The model problem

9 Results

10 Discussions and Conclusions

11 Future works

12 Bibliography

13 Appendix