- Email: [email protected]
Numerical simulation of piston ring lubrication
ARTICLE IN PRESS
Tribology International 41 (2008) 914–919 www.elsevier.com/locate/triboint
Numerical simulation of piston ring lubrication Christian Lotz Felter MAN B & W Diesel A/S, Technical University of Denmark, Nils Koppels Alle´, Bygning 404, DK - 2800 Kgs. Lyngby, Denmark Received 4 December 2006; received in revised form 20 November 2007; accepted 23 November 2007 Available online 22 January 2008
Abstract This paper describes a numerical method that can be used to model the lubrication of piston rings. Classical lubrication theory is based on the Reynolds equation which is applicable to confined geometries and open geometries where the flooding conditions are known. Lubrication of piston rings, however, fall outside this category of problems since the piston rings might suffer from starved running conditions. This means that the computational domain where the Reynolds equation is applicable (including a cavitation criteria) is unknown. In order to overcome this problem the computational domain is extended to include also the oil film outside the piston rings. The numerical model consists of a 2D free surface code that solves the time dependent compressible Navier–Stokes equations. The equations are cast in Lagrangian form and discretized by a meshfree moving least squares method using the primitive variables u, v, r for the velocity components and density, respectively. Time integration is performed by a third order Runge–Kutta method. The set of equations is closed by the Dowson–Higginson equation for the relation between density and pressure. Boundary conditions are the nonslip condition on solids and the equilibrium of stresses on the free surface. It is assumed that the surrounding gas phase has zero viscosity. Surface tension can be included in the model if necessary. The contact point where the three phases solid, liquid, and gas intersect is updated based on the velocity of the solid and the angle between the normals of the solid and the free surface. The numerical model is compared with the results from an analytical solution of the Reynolds equation for a fixed incline slider bearing. Then results from a more complicated simulation of piston ring lubrication are given and discussed. r 2007 Elsevier Ltd. All rights reserved. Keywords: Piston ring; Reynolds equation; Navier–Stokes equations; Free surface; Moving least squares
1. Introduction The performance of piston rings in combustion engines has been a topic of research for many years. Piston rings act as sealing between the liner and the piston. The piston rings are lubricated with oil and can thus be considered as slider bearings. This paper is concerned with theoretical predictions of the performance based on numerical simulations. Classical theory of lubrication is based on the Reynolds equation, which can be derived from control volume analysis under certain simplifying assumptions [1]. This equation calculates the oil film pressure given the film thickness, squeeze velocity, and the pressure at the boundaries (typically ambient pressure). It is well known that piston rings can have greatly changing running conditions. At some points the piston ring might experience E-mail address: [email protected] 0301-679X/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.triboint.2007.11.018
fully flooded conditions, at others it might suffer from starvation. This means that knowledge of the amount of available oil becomes very important for successful simulations. However, the Reynolds equation does not apply outside the piston ring and therefore cannot model e.g. buildup of oil in front of the moving piston ring. Therefore it becomes difficult to relate the undisturbed oil film thickness on the liner and the immediate thickness in front of the piston ring. The literature on simulation of piston rings is vast. Since the problem is very complex each publication typically does not treat all effects at once. However, the problem of the unknown inlet film thickness has been addressed by other workers. In the paper by Esfahanian et al. [2] a single ring is treated assuming fully flooded conditions at all times. Dowson [3] analyzes a single ring and a complete ring pack. In that paper the undisturbed oil film thickness is assumed to be one-half of the film thickness at the location under the piston ring where qp=qx ¼ 0 (thus the effect of squeeze
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919
action on the oil flow is neglected). A more complicated model is given in [4], where oil build up in front of the piston ring is modeled by a trapezoidal geometry. The same idea is followed in [5] expect the build up is assumed to have a parabolic shape. In this paper a 2D free surface model based on the Navier–Stokes equations is developed. The following section describes the governing equations and boundary conditions. Then follows a section on the numerical method. Finally we compare analytical results from the fixed incline slider bearing with the numerical results, and a more general problem including a free surface is solved.
2. Method 2.1. Governing equations The fluid flow is governed by the 2D compressible Navier–Stokes equations 2 du 1 q u q2 u 1 q2 u q2 v qp ¼ m þ þ þ , dt r qx2 qy2 3 qx2 qxqy qx
915
2.2. Boundary conditions Three types of boundary conditions are considered in this work. On solid walls the non-slip condition is used, which means that ufluid ¼ uwall ;
vfluid ¼ vwall .
(4)
In order to find the pressure on solid walls we write the following equation: ! ! nx du=dt RHSx ¼ 0, (5) ny dv=dt RHSy where nx is the x-component of solid wall normal, ny the ycomponent of solid wall normal, denotes the scalar product of two vectors, RHSx the right-hand side of first equation in (1) and RHSy the right-hand side of second equation in (1) which after rearrangement gives a Neumann condition for pressure qp=qn ¼ . Note that the velocity at the wall is known so one can evaluate du=dt and dv=dt. Note also that p and r are not independent since they are coupled through the equation of state. Thus one can either eliminate r using r ¼ f ðpÞ
2 dv 1 q v q2 v 1 q2 u q2 v qp ¼ m þ þ þ , dt r qx2 qy2 3 qyqx qy2 qy
(1)
where d=dt is the total derivative with respect to time t, u the x-component of velocity, v the y-component of velocity, p the pressure, r the density and m the dynamic viscosity. As can be seen from the expressions above viscosity is assumed to be constant, with a value given by an estimated average. One possibility is to evaluate viscosity at the mean temperature along the liner (based on measurements) and neglect the dependency from pressure. The continuity equation is stated as dr qu qv ¼ r þ . (2) dt qx qy The system is closed using the Dowson–Higginson equation of state [1, p. 71] 0:6p r ¼ r0 1 þ , (3) 1 þ 1:7p where r0 is the density when p ¼ 0, kg=m3 and p the gauge pressure, GPa. Piston rings in general have a barrel shaped running surface, which gives a converging/diverging geometry in connection with the liner. When the piston ring is moving pressure will build up on the converging part while cavitation might appear on the diverging part, due to the associated pressure drop. Currently the effect of cavitation is not included in the model, meaning that simulations are only valid when cavitation is not present. A cavitation model will be added in a future work.
or eliminate p by writing qp qp qr ¼ , qx qr qx qp qp qr ¼ qy qr qy and then solve (5) for just one unknown quantity. We use the second possibility solving for r. Because the equation of state for this application is nonlinear the solution procedure for r uses iteration. The boundary condition for a free surface is the balance of normal and shear stresses on the interface separating the two fluids. This can be stated as ðTij ni nj Þliquid ¼ ðTij ni nj Þgas , ðTij ni sj Þliquid ¼ ðTij ni sj Þgas , where T is the stress tensor, n the normal vector and s the tangent vector. It is assumed that the gas phase has constant pressure and zero viscosity. The effect of surface tension is neglected. Substituting the expression for the stress tensor [6] and using the fact that the normal and tangent vectors are perpendicular we get 2 qu qv qu 2 þ n pþ m 2m 3 qx qy qx x qv qu qv þ þ ¼ ½pgas , nx ny þ n2y qx qy qy liquid qv qu qu qv 2 2 þ ¼ ½0gas . ð6Þ ðnx ny Þ 2 nx ny qx qy qx qy liquid
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919
916
The pressure on the left-hand side (pressure in the liquid phase at the interface) can be calculated using (5) and (3). This gives three nonlinear equations for the three unknowns u, v, and p. Finally a special boundary condition for the triple point where the free surface touches the solid walls is needed. Clearly this point must be able to move and therefore the non-slip condition does not apply. The physic of the triple point acts at the level of molecules and cannot be simulated directly by a continuum model. Instead we require no flow through the wall, and adjust the triple point velocity such that the contact angle tends to p=2 at all times. This can be achieved by
vtriplepoint ¼ vwall
(7)
freesurface freesurface þ nwall ÞC, q ¼ ½cosðfss Þ ðnwall x nx y ny
where fss is the prescribed steady state contact angle and C the tuning parameter. See [7] for details and more advanced models. Pressure is not evaluated at the triple points since the normal is not well defined. 3. Numerical scheme In this section we describe how the expressions from the previous section are discretized. We employ the semidiscretization technique. For the spacial derivatives we use the method of moving least squares. The method is outlined below [8]. Consider the field f ¼ f ðxÞ ¼ f ðx; yÞ. If f is smooth within some distance from the point x it can be approximated by a polynomial surface. Define the functional Z Jðx Þ ¼ wðx x Þ½aðx ÞT pðx x Þ f ðxÞ2 dx (8) O
(11)
k ¼ 1; 2; . . . ; N; I index of points within smoothing length.
(12)
After rearrangement the equations can be put in matrix form 9 P 8 > > I wðxI x Þðf I f Þ > > > P > > > > > wðx x Þp ðx x Þðf f Þ = < I I 1 I I , (13) Ma ¼ .. > > . > > > > > > > ; : P wðxI x Þp ðxI x Þðf f Þ > 5 I I where M is called the matrix of moments. Note that the equation corresponding to the constant term in the polynomial basis has been eliminated by moving f ¼ f ðx Þ to the right-hand side. This ensures that the approximating surface passes through the field value at the point of expansion (x must coincide with one of the data points xI ). This property is crucial for the method. After determining the coefficients the field f can be approximated in the vicinity of x by (with a1 ¼ f ðx Þ) f ðxÞ ¼
N X
pi ðx x Þai .
(14)
i¼1
It is seen that the units of the ai ’s depend on the units of the field function f. Similarly the partial derivatives at the point of expansion x are approximated by
where x point of expansion; a coefficients to be determined; p ¼ f1 x y x2 xy g polynomial basis:
i¼1
N number of basis functions;
where
w kernel function with compact support;
I
where
utriplepoint ¼ uwall þ nwall y q, nwall x q,
functional (8) local. Thus, in order to evaluate J we do not need to sum over all points, but only those within the smoothing length of x . Furthermore the kernel has a belllike shape, with qw=qqjq¼2 ¼ 0. This means that a point can enter or leave the support of J is a smooth manner. The coefficients a are found by minimizing J using the stationarity condition qJ=qai ¼ 0. In a discrete setting this becomes " # N X X wðxI x Þ pi ðxI x Þai ðx Þ f I pk ðxI x Þ ¼ 0,
(9)
(We distinguish between pressure p and basis functions p. One element of p is denoted by the usual notation pi .) We use the following kernel expression: 8 2 3 > < 1 3=2q þ 3=4q ; qo1; 3 1pqo2; wðh; xÞ ¼ 1=4ð2 qÞ ; (10) > : 0; 2pq; where h is the smoothing length and q ¼ 2kxk=h. It is seen that the kernel has compact support, which makes the
N X qðjÞ f qðjÞ pi ð0Þ ¼ ai , qxðkÞ qyðlÞ qxðkÞ qyðlÞ i¼1
(15)
where k þ l ¼ j. These expressions are substituted directly into the governing equations. The time integration is performed by an explicit third order Runge–Kutta method. As indicated in (1) the computation is performed in the moving (Lagrangian) frame, which means that the position of the computational nodes must be updated according to dx ¼ u; dt
dy ¼ v. dt
(16)
Due to the distortion of the flow field it is necessary to redistribute the computational nodes at some time interval.
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919
4. Results First we compare three solutions for pressure from the fixed incline slider bearing. Consider the following situation: hi ¼ 25 mm,
inlet outlet
ho ¼ 15 mm,
width
l ¼ 1 mm,
speed
u ¼ 1 m=s,
viscosity
m ¼ 0:05 Pa s,
density r0 ¼ 900 kg=m3 . In Fig. 1 the analytical solution of the Reynolds equation is plotted together with the pressure obtained from the simulation program. The effect of compressibility is examined by solving the Reynolds equation for an incompressible fluid and also using the Dowson–Higginson equation of state. It is seen that compressibility leads to a slightly higher maximum pressure and that the location is shifted downstream. The solution of the Navier–Stokes equations gives results similar to the Reynolds equation. It was, however, not possible to achieve a fully converged solution of the Navier–Stokes equations because of an x 104 Reynolds (incompr) Reynolds (compr) Navier Stokes
gauge pressure [Pa]
5
instability at the inlet boundary condition. It is expected that this explains the somewhat smaller maximum pressure. A better implementation of the inlet boundary condition is a subject for future work. Fig. 2 shows the density distribution in the fixed incline slider bearing at the final time of the simulation, that is, at steady state. It is seen that the density (from a graphical view) is constant across the film. This feature is one of the assumptions under which the Reynolds equation is derived. The graph on Fig. 2 indicates that this assumption is perfectly valid for thin films at steady state. Second we simulate a more complicated situation including a free surface. Initially the piston ring is at rest and the oil film has uniform thickness outside the piston ring, see Fig. 3. The horizontal and vertical position of the piston ring is prescribed by sine functions. In a real application the vertical position will be controlled by the forces acting from the surrounding gas and the oil film. Modeling of that effect is omitted in this presentation. Here are the simulation parameters: ring shape
parabolic,
ring height
sh ¼ 1 mm,
ring width
l ¼ 2 mm,
ring horizontal speed ring vertical speed
6
917
ambient pressure
u ¼ 0:2p cosð200pt p=2Þ m=s,
v ¼ 0:05p cosð400pt p=2Þ m=s, pamb ¼ 0 Pa,
initial oil film thickness hoil ¼ 850 mm.
4
3
2
1
0 0
0.2
0.4
0.6
0.8
x [m] Fig. 1. Pressure curves for the fixed incline slider bearing.
1 x 10−3
It should be noted that these settings are unrealistic for any internal combustion engine. However, we use these values to illustrate the performance of the free surface code. It should also be noted that cavitation is not included in the model, so negative pressure can appear during the simulation. Gravity forces are not included. The first plot shows the velocity field after the first 100 time steps. It is seen that the triple points are accelerated in order to produce a contact angle of p=2 (velocity vectors are scaled using the maximum speed). The following plots show the build up of a vortex that develops under the piston ring. Also the deformation of the oil film is clearly seen.
Fig. 2. Density plot for the fixed incline slider bearing. Red color means high density and blue color means low density.
ARTICLE IN PRESS 918
C. Lotz Felter / Tribology International 41 (2008) 914–919
Fig. 3. Snapshots from piston ring simulation. From top to bottom t ¼ 3 105 s, t ¼ 8:1 104 s, t ¼ 1:7 103 s, t ¼ 2:8 103 s, t ¼ 3:9 103 s.
5. Discussion
6. Conclusion
Two simulation problems have been considered. For the fixed incline slider bearing the method agrees well when compared with the results using the Reynolds equation. For the case of piston ring simulation the model is not fully developed yet. From the modeling point of view the effect of cavitation needs to be considered. Currently, the model allows the generation of negative pressure, which is not physically admissible. The computational work required by the method must be investigated further. The time required for the calculation shown in this paper is approximately 1 h on a standard 3 GHz desktop computer. However, when solving real life problems the spacial resolution (and consequently the time step) must be improved. Despite these issues the presented method appears to be a promising alternative, when the Reynolds equation cannot be applied.
A new approach for the numerical simulation of piston ring lubrication has been presented. The main idea is to simulate also the free surface of the oil film outside the piston rings. The method has been compared with results for a confined geometry using the Reynolds equation with well agreement. The current state of the model indicates that more work is needed in order to include the effect of cavitation and also examine the computational cost induced by problems from real applications. References [1] Hamrock BJ. Fundamentals of fluid film lubrication. New York: McGraw-Hill, Inc.; 1994. p. 71. [2] Esfahanian M, Hamrock BJ, Elsharkawy AA. On the hydrodynamic lubrication analysis of piston rings. Lubr Sci; 1998;10(4):265–86.
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919 [3] Dowson D, Economou PN, Ruddy BL, Strachan PJ, Baker AJS. Piston ring lubrication—part II: theoretical analysis of a single ring and a complete ring pack. Energy conservation through fluid film lubrication technology: frontiers in research and design. Proceedings of the ASME winter annual meeting, 1979. p. 23–52. [4] Han D-C, Lee J-S. Analysis of the piston ring lubrication with a new boundary condition. Tribol Int 1998;31(12):753–60. [5] Gamble RJ, Priest M, Taylor CM. Detailed analysis of oil transport in the piston assembly of a gasoline engine. Tribol Lett 2003;14(2):147–56.
919
[6] Chung TJ. Computational fluid dynamics. Cambridge: 2002. p. 33. [7] Baer TA, Cairncross RA, Schunka PR, Rao RR, Sackinger PA. A finite element method for free surface flows of incompressible fluids in three dimensions. Part II. Dynamic wetting lines. Int J Numer Methods Fluids 2000;33(3):405–27. [8] Liu GR. Mesh free methods moving beyond the finite element method. Boca Raton, FL: CRC Press; 2003. p. 70–87.
Tribology International 41 (2008) 914–919 www.elsevier.com/locate/triboint
Numerical simulation of piston ring lubrication Christian Lotz Felter MAN B & W Diesel A/S, Technical University of Denmark, Nils Koppels Alle´, Bygning 404, DK - 2800 Kgs. Lyngby, Denmark Received 4 December 2006; received in revised form 20 November 2007; accepted 23 November 2007 Available online 22 January 2008
Abstract This paper describes a numerical method that can be used to model the lubrication of piston rings. Classical lubrication theory is based on the Reynolds equation which is applicable to confined geometries and open geometries where the flooding conditions are known. Lubrication of piston rings, however, fall outside this category of problems since the piston rings might suffer from starved running conditions. This means that the computational domain where the Reynolds equation is applicable (including a cavitation criteria) is unknown. In order to overcome this problem the computational domain is extended to include also the oil film outside the piston rings. The numerical model consists of a 2D free surface code that solves the time dependent compressible Navier–Stokes equations. The equations are cast in Lagrangian form and discretized by a meshfree moving least squares method using the primitive variables u, v, r for the velocity components and density, respectively. Time integration is performed by a third order Runge–Kutta method. The set of equations is closed by the Dowson–Higginson equation for the relation between density and pressure. Boundary conditions are the nonslip condition on solids and the equilibrium of stresses on the free surface. It is assumed that the surrounding gas phase has zero viscosity. Surface tension can be included in the model if necessary. The contact point where the three phases solid, liquid, and gas intersect is updated based on the velocity of the solid and the angle between the normals of the solid and the free surface. The numerical model is compared with the results from an analytical solution of the Reynolds equation for a fixed incline slider bearing. Then results from a more complicated simulation of piston ring lubrication are given and discussed. r 2007 Elsevier Ltd. All rights reserved. Keywords: Piston ring; Reynolds equation; Navier–Stokes equations; Free surface; Moving least squares
1. Introduction The performance of piston rings in combustion engines has been a topic of research for many years. Piston rings act as sealing between the liner and the piston. The piston rings are lubricated with oil and can thus be considered as slider bearings. This paper is concerned with theoretical predictions of the performance based on numerical simulations. Classical theory of lubrication is based on the Reynolds equation, which can be derived from control volume analysis under certain simplifying assumptions [1]. This equation calculates the oil film pressure given the film thickness, squeeze velocity, and the pressure at the boundaries (typically ambient pressure). It is well known that piston rings can have greatly changing running conditions. At some points the piston ring might experience E-mail address: [email protected] 0301-679X/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.triboint.2007.11.018
fully flooded conditions, at others it might suffer from starvation. This means that knowledge of the amount of available oil becomes very important for successful simulations. However, the Reynolds equation does not apply outside the piston ring and therefore cannot model e.g. buildup of oil in front of the moving piston ring. Therefore it becomes difficult to relate the undisturbed oil film thickness on the liner and the immediate thickness in front of the piston ring. The literature on simulation of piston rings is vast. Since the problem is very complex each publication typically does not treat all effects at once. However, the problem of the unknown inlet film thickness has been addressed by other workers. In the paper by Esfahanian et al. [2] a single ring is treated assuming fully flooded conditions at all times. Dowson [3] analyzes a single ring and a complete ring pack. In that paper the undisturbed oil film thickness is assumed to be one-half of the film thickness at the location under the piston ring where qp=qx ¼ 0 (thus the effect of squeeze
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919
action on the oil flow is neglected). A more complicated model is given in [4], where oil build up in front of the piston ring is modeled by a trapezoidal geometry. The same idea is followed in [5] expect the build up is assumed to have a parabolic shape. In this paper a 2D free surface model based on the Navier–Stokes equations is developed. The following section describes the governing equations and boundary conditions. Then follows a section on the numerical method. Finally we compare analytical results from the fixed incline slider bearing with the numerical results, and a more general problem including a free surface is solved.
2. Method 2.1. Governing equations The fluid flow is governed by the 2D compressible Navier–Stokes equations 2 du 1 q u q2 u 1 q2 u q2 v qp ¼ m þ þ þ , dt r qx2 qy2 3 qx2 qxqy qx
915
2.2. Boundary conditions Three types of boundary conditions are considered in this work. On solid walls the non-slip condition is used, which means that ufluid ¼ uwall ;
vfluid ¼ vwall .
(4)
In order to find the pressure on solid walls we write the following equation: ! ! nx du=dt RHSx ¼ 0, (5) ny dv=dt RHSy where nx is the x-component of solid wall normal, ny the ycomponent of solid wall normal, denotes the scalar product of two vectors, RHSx the right-hand side of first equation in (1) and RHSy the right-hand side of second equation in (1) which after rearrangement gives a Neumann condition for pressure qp=qn ¼ . Note that the velocity at the wall is known so one can evaluate du=dt and dv=dt. Note also that p and r are not independent since they are coupled through the equation of state. Thus one can either eliminate r using r ¼ f ðpÞ
2 dv 1 q v q2 v 1 q2 u q2 v qp ¼ m þ þ þ , dt r qx2 qy2 3 qyqx qy2 qy
(1)
where d=dt is the total derivative with respect to time t, u the x-component of velocity, v the y-component of velocity, p the pressure, r the density and m the dynamic viscosity. As can be seen from the expressions above viscosity is assumed to be constant, with a value given by an estimated average. One possibility is to evaluate viscosity at the mean temperature along the liner (based on measurements) and neglect the dependency from pressure. The continuity equation is stated as dr qu qv ¼ r þ . (2) dt qx qy The system is closed using the Dowson–Higginson equation of state [1, p. 71] 0:6p r ¼ r0 1 þ , (3) 1 þ 1:7p where r0 is the density when p ¼ 0, kg=m3 and p the gauge pressure, GPa. Piston rings in general have a barrel shaped running surface, which gives a converging/diverging geometry in connection with the liner. When the piston ring is moving pressure will build up on the converging part while cavitation might appear on the diverging part, due to the associated pressure drop. Currently the effect of cavitation is not included in the model, meaning that simulations are only valid when cavitation is not present. A cavitation model will be added in a future work.
or eliminate p by writing qp qp qr ¼ , qx qr qx qp qp qr ¼ qy qr qy and then solve (5) for just one unknown quantity. We use the second possibility solving for r. Because the equation of state for this application is nonlinear the solution procedure for r uses iteration. The boundary condition for a free surface is the balance of normal and shear stresses on the interface separating the two fluids. This can be stated as ðTij ni nj Þliquid ¼ ðTij ni nj Þgas , ðTij ni sj Þliquid ¼ ðTij ni sj Þgas , where T is the stress tensor, n the normal vector and s the tangent vector. It is assumed that the gas phase has constant pressure and zero viscosity. The effect of surface tension is neglected. Substituting the expression for the stress tensor [6] and using the fact that the normal and tangent vectors are perpendicular we get 2 qu qv qu 2 þ n pþ m 2m 3 qx qy qx x qv qu qv þ þ ¼ ½pgas , nx ny þ n2y qx qy qy liquid qv qu qu qv 2 2 þ ¼ ½0gas . ð6Þ ðnx ny Þ 2 nx ny qx qy qx qy liquid
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919
916
The pressure on the left-hand side (pressure in the liquid phase at the interface) can be calculated using (5) and (3). This gives three nonlinear equations for the three unknowns u, v, and p. Finally a special boundary condition for the triple point where the free surface touches the solid walls is needed. Clearly this point must be able to move and therefore the non-slip condition does not apply. The physic of the triple point acts at the level of molecules and cannot be simulated directly by a continuum model. Instead we require no flow through the wall, and adjust the triple point velocity such that the contact angle tends to p=2 at all times. This can be achieved by
vtriplepoint ¼ vwall
(7)
freesurface freesurface þ nwall ÞC, q ¼ ½cosðfss Þ ðnwall x nx y ny
where fss is the prescribed steady state contact angle and C the tuning parameter. See [7] for details and more advanced models. Pressure is not evaluated at the triple points since the normal is not well defined. 3. Numerical scheme In this section we describe how the expressions from the previous section are discretized. We employ the semidiscretization technique. For the spacial derivatives we use the method of moving least squares. The method is outlined below [8]. Consider the field f ¼ f ðxÞ ¼ f ðx; yÞ. If f is smooth within some distance from the point x it can be approximated by a polynomial surface. Define the functional Z Jðx Þ ¼ wðx x Þ½aðx ÞT pðx x Þ f ðxÞ2 dx (8) O
(11)
k ¼ 1; 2; . . . ; N; I index of points within smoothing length.
(12)
After rearrangement the equations can be put in matrix form 9 P 8 > > I wðxI x Þðf I f Þ > > > P > > > > > wðx x Þp ðx x Þðf f Þ = < I I 1 I I , (13) Ma ¼ .. > > . > > > > > > > ; : P wðxI x Þp ðxI x Þðf f Þ > 5 I I where M is called the matrix of moments. Note that the equation corresponding to the constant term in the polynomial basis has been eliminated by moving f ¼ f ðx Þ to the right-hand side. This ensures that the approximating surface passes through the field value at the point of expansion (x must coincide with one of the data points xI ). This property is crucial for the method. After determining the coefficients the field f can be approximated in the vicinity of x by (with a1 ¼ f ðx Þ) f ðxÞ ¼
N X
pi ðx x Þai .
(14)
i¼1
It is seen that the units of the ai ’s depend on the units of the field function f. Similarly the partial derivatives at the point of expansion x are approximated by
where x point of expansion; a coefficients to be determined; p ¼ f1 x y x2 xy g polynomial basis:
i¼1
N number of basis functions;
where
w kernel function with compact support;
I
where
utriplepoint ¼ uwall þ nwall y q, nwall x q,
functional (8) local. Thus, in order to evaluate J we do not need to sum over all points, but only those within the smoothing length of x . Furthermore the kernel has a belllike shape, with qw=qqjq¼2 ¼ 0. This means that a point can enter or leave the support of J is a smooth manner. The coefficients a are found by minimizing J using the stationarity condition qJ=qai ¼ 0. In a discrete setting this becomes " # N X X wðxI x Þ pi ðxI x Þai ðx Þ f I pk ðxI x Þ ¼ 0,
(9)
(We distinguish between pressure p and basis functions p. One element of p is denoted by the usual notation pi .) We use the following kernel expression: 8 2 3 > < 1 3=2q þ 3=4q ; qo1; 3 1pqo2; wðh; xÞ ¼ 1=4ð2 qÞ ; (10) > : 0; 2pq; where h is the smoothing length and q ¼ 2kxk=h. It is seen that the kernel has compact support, which makes the
N X qðjÞ f qðjÞ pi ð0Þ ¼ ai , qxðkÞ qyðlÞ qxðkÞ qyðlÞ i¼1
(15)
where k þ l ¼ j. These expressions are substituted directly into the governing equations. The time integration is performed by an explicit third order Runge–Kutta method. As indicated in (1) the computation is performed in the moving (Lagrangian) frame, which means that the position of the computational nodes must be updated according to dx ¼ u; dt
dy ¼ v. dt
(16)
Due to the distortion of the flow field it is necessary to redistribute the computational nodes at some time interval.
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919
4. Results First we compare three solutions for pressure from the fixed incline slider bearing. Consider the following situation: hi ¼ 25 mm,
inlet outlet
ho ¼ 15 mm,
width
l ¼ 1 mm,
speed
u ¼ 1 m=s,
viscosity
m ¼ 0:05 Pa s,
density r0 ¼ 900 kg=m3 . In Fig. 1 the analytical solution of the Reynolds equation is plotted together with the pressure obtained from the simulation program. The effect of compressibility is examined by solving the Reynolds equation for an incompressible fluid and also using the Dowson–Higginson equation of state. It is seen that compressibility leads to a slightly higher maximum pressure and that the location is shifted downstream. The solution of the Navier–Stokes equations gives results similar to the Reynolds equation. It was, however, not possible to achieve a fully converged solution of the Navier–Stokes equations because of an x 104 Reynolds (incompr) Reynolds (compr) Navier Stokes
gauge pressure [Pa]
5
instability at the inlet boundary condition. It is expected that this explains the somewhat smaller maximum pressure. A better implementation of the inlet boundary condition is a subject for future work. Fig. 2 shows the density distribution in the fixed incline slider bearing at the final time of the simulation, that is, at steady state. It is seen that the density (from a graphical view) is constant across the film. This feature is one of the assumptions under which the Reynolds equation is derived. The graph on Fig. 2 indicates that this assumption is perfectly valid for thin films at steady state. Second we simulate a more complicated situation including a free surface. Initially the piston ring is at rest and the oil film has uniform thickness outside the piston ring, see Fig. 3. The horizontal and vertical position of the piston ring is prescribed by sine functions. In a real application the vertical position will be controlled by the forces acting from the surrounding gas and the oil film. Modeling of that effect is omitted in this presentation. Here are the simulation parameters: ring shape
parabolic,
ring height
sh ¼ 1 mm,
ring width
l ¼ 2 mm,
ring horizontal speed ring vertical speed
6
917
ambient pressure
u ¼ 0:2p cosð200pt p=2Þ m=s,
v ¼ 0:05p cosð400pt p=2Þ m=s, pamb ¼ 0 Pa,
initial oil film thickness hoil ¼ 850 mm.
4
3
2
1
0 0
0.2
0.4
0.6
0.8
x [m] Fig. 1. Pressure curves for the fixed incline slider bearing.
1 x 10−3
It should be noted that these settings are unrealistic for any internal combustion engine. However, we use these values to illustrate the performance of the free surface code. It should also be noted that cavitation is not included in the model, so negative pressure can appear during the simulation. Gravity forces are not included. The first plot shows the velocity field after the first 100 time steps. It is seen that the triple points are accelerated in order to produce a contact angle of p=2 (velocity vectors are scaled using the maximum speed). The following plots show the build up of a vortex that develops under the piston ring. Also the deformation of the oil film is clearly seen.
Fig. 2. Density plot for the fixed incline slider bearing. Red color means high density and blue color means low density.
ARTICLE IN PRESS 918
C. Lotz Felter / Tribology International 41 (2008) 914–919
Fig. 3. Snapshots from piston ring simulation. From top to bottom t ¼ 3 105 s, t ¼ 8:1 104 s, t ¼ 1:7 103 s, t ¼ 2:8 103 s, t ¼ 3:9 103 s.
5. Discussion
6. Conclusion
Two simulation problems have been considered. For the fixed incline slider bearing the method agrees well when compared with the results using the Reynolds equation. For the case of piston ring simulation the model is not fully developed yet. From the modeling point of view the effect of cavitation needs to be considered. Currently, the model allows the generation of negative pressure, which is not physically admissible. The computational work required by the method must be investigated further. The time required for the calculation shown in this paper is approximately 1 h on a standard 3 GHz desktop computer. However, when solving real life problems the spacial resolution (and consequently the time step) must be improved. Despite these issues the presented method appears to be a promising alternative, when the Reynolds equation cannot be applied.
A new approach for the numerical simulation of piston ring lubrication has been presented. The main idea is to simulate also the free surface of the oil film outside the piston rings. The method has been compared with results for a confined geometry using the Reynolds equation with well agreement. The current state of the model indicates that more work is needed in order to include the effect of cavitation and also examine the computational cost induced by problems from real applications. References [1] Hamrock BJ. Fundamentals of fluid film lubrication. New York: McGraw-Hill, Inc.; 1994. p. 71. [2] Esfahanian M, Hamrock BJ, Elsharkawy AA. On the hydrodynamic lubrication analysis of piston rings. Lubr Sci; 1998;10(4):265–86.
ARTICLE IN PRESS C. Lotz Felter / Tribology International 41 (2008) 914–919 [3] Dowson D, Economou PN, Ruddy BL, Strachan PJ, Baker AJS. Piston ring lubrication—part II: theoretical analysis of a single ring and a complete ring pack. Energy conservation through fluid film lubrication technology: frontiers in research and design. Proceedings of the ASME winter annual meeting, 1979. p. 23–52. [4] Han D-C, Lee J-S. Analysis of the piston ring lubrication with a new boundary condition. Tribol Int 1998;31(12):753–60. [5] Gamble RJ, Priest M, Taylor CM. Detailed analysis of oil transport in the piston assembly of a gasoline engine. Tribol Lett 2003;14(2):147–56.
919
[6] Chung TJ. Computational fluid dynamics. Cambridge: 2002. p. 33. [7] Baer TA, Cairncross RA, Schunka PR, Rao RR, Sackinger PA. A finite element method for free surface flows of incompressible fluids in three dimensions. Part II. Dynamic wetting lines. Int J Numer Methods Fluids 2000;33(3):405–27. [8] Liu GR. Mesh free methods moving beyond the finite element method. Boca Raton, FL: CRC Press; 2003. p. 70–87.