 Open access peer-reviewed chapter

# Modelling Cavitation in (Elasto)Hydrodynamic Lubrication

Written By

Andreas Almqvist and Peter Wall

Submitted: November 10th, 2015 Reviewed: April 8th, 2016 Published: October 26th, 2016

DOI: 10.5772/63533

From the Edited Volume

Edited by Pranav H. Darji

Chapter metrics overview

View Full Metrics

## Abstract

In this chapter we will present a derivation of a mathematical model describing how cavitation influences the pressure distribution in a thin lubricant film between two moving surfaces. The main idea in the derivation is to first describe the influence of cavitation on the mass flow and thereafter using a conservation law for the mass. This leads to a nonlinear system with two complementary variables: one is the pressure distribution and the other is related to the density, i.e. a nonlinear complementarity problem (NLCP). The proposed approach is used to derive a mass conserving cavitation model considering that density, viscosity and film thickness of the lubricant depend on the pressure. To demonstrate the applicability and evaluate the proposed model and the suggested numerical implementation, a few model problems are analysed and presented.

### Keywords

• Cavitation
• Reynolds equation
• Thin film flow
• Complementarity problem
• Hydrodynamic lubrication

## 1. Introduction

A central problem in hydrodynamic lubrication is to model the pressure in a lubricant (fluid) between two surfaces which are in relative motion. In this chapter we consider the full film regime, i.e. when the surfaces are fully separated by the lubricant. The pressure and the velocity field can be modelled using the Navier-Stokes equations and the continuity equation. However, in real applications the distance between the surfaces, h, is extremely small in relation to the size of the surfaces. This fact can be used in order to derive a simplified two-dimensional equation for the hydrodynamic pressure, p. Indeed, if the velocity of the upper surface is vuand the velocity of the lower surface is vl, then under the assumptions outlined in e.g. , the following Reynolds equation is obtained

t(ρh)=(ρh312μp12ρh(vu+vl)),E1

where μand ρare the viscosity and density of the lubricant. A generalization of this thin film approximation which considers the effects of non-Newtonian fluids can e.g. be found in the Encyclopaedia of Tribology, see the entry  by Larsson.

A fluid cannot sustain large tensile stress and it is known that when the pressure becomes too low the continuous film will rupture and air bubbles will be formed. This phenomenon is known as cavitation and has a huge impact on the hydrodynamic performance. In areas where cavitation takes place the pressure is usually treated as constant and is assumed to have the same value as the saturated vapour pressure, i.e., the pressure at which cavitation starts. This is commonly referred to as the cavitation pressure, here denoted by pc.

The common practice in the field is to build mathematical models that consider hydrodynamic cavitation based on the Reynolds equation (1). These models rests on the assumption that the pressure can be regarded as constant in areas where cavitation takes place and this constant level of pressure is typically also assumed to be the saturated vapour pressure, i.e., the pressure at which cavitation starts. The most important pioneering works are the papers presenting the Jakobsson-Floberg-Olsson (JFO) cavitation boundary conditions , Elrod’s work  comprising these boundary conditions in one universal equation for the unknown pressure (or saturation) and its corresponding finite difference scheme and then Vijayaraghavan’s generalization  of Elrod’s results. These results have been frequently used to study the effects of cavitation in real applications and they have also been subject for many further generalizations.

In Elrod’s approach the lubricant is treated as incompressible while Vijayaraghavan and Keith assumed that the relation between the density and the pressure is of the form

ρ=ρcexp((ppc)/β),E2

where ρcis the density at cavitation pressure pc, and βis the bulk modulus of the lubricant. One example of a generalization of the results in  and  is the paper  by Sahlin et al., who developed a cavitation algorithm for arbitrary pressure-density relationships. The case considering incompressible lubricants was further developed by Bayada et al. [9, 10]. Their results put the finite difference scheme in  into a framework which makes further mathematical analysis possible. Moreover, these works consider both cavitation and the effects of surface roughness, where the surface roughness is analyzed by homogenization. In [10, 11] these ideas were extended to also include elastohydrodynamic effects, i.e. elastic deformation of the surfaces caused by the hydrodynamic pressure. In particular this means that the film thickness, h, depends on the pressure.

A significant progress in mathematical modeling of cavitation in hydrodynamic lubrication was recently presented by Giacopini et al. . They reformulated the model presented in [9, 10] for incompressible fluids. The model in [9, 10] includes two unknowns, namely the pressure and the saturation (or density) while the unknowns in the reformulation are the pressure and a new unknown variable related to the saturation. The major advantage with the reformulation is that the two unknowns are complementary, i.e. their product is zero, in the whole domain. This implies that the discretized system of equations becomes a linear complementarity problem (LCP) which can be readily solved by standard numerical methods as e.g. the Lemke’s pivoting algorithm [13, 14]. The idea of using complementarity was further developed by Bertocci et al. . They present a very comprehensive cavitation model, which discretized assumes the form of a nonlinear complementarity problem (NLCP).

In the work by Almqvist et al. , a new approach for modeling cavitation was presented. The model was derived by first considering how the mass flow is influenced by cavitation and thereafter using the law of conservation of mass. This immediately leads to a linear complementarity problem formulation. This is in contrast to the approach in [12, 15] where the pressure–density relationship is inserted directly into the Reynolds equation, thereafter it is necessary to argue for that certain terms can be cancelled before they arrive at a complementarity formulation. Moreover, in addition to  the approach in Almqvist et al.  also covers compressible fluids. In the particular case when the pressure–density relationship is as in (2), a change of variables was introduced which transforms the problem such that the discrete formulation is an LCP.

This chapter extends the mathematical modelling presented in  to also include more general pressure–density relationships and allow the viscosity and film thickness to depend on the pressure. In particular this means that the model can be used to efficiently study elastohydrodynamic lubrication where cavitation is present. To demonstrate the applicability and evaluate the proposed model and the associated numerical solution method, four model problems are analysed.

## 2. Mathematical model

In this section, a model considering cavitation in thin film flow between two surfaces in relative motion is developed. The flow is regarded as compressible, the viscosity can be shear rate and pressure dependent (non-Newtonian and piezo-viscous) and elastic deformation of the contacting surfaces can be considered.

Let us start by defining the (three dimensional) fluid domain Ω between the two lubricated surfaces hl(lower) and hu(upper):

Ω={(x,z)|0xili,hlzhu}.E3

The fluid domain Ω is schematically illustrated in Figure 1 (left).

The hydrodynamic pressures that develop may be large enough to deform the surfaces. This implies that the film thickness h= huhldepends on the pressure, i.e., h= h(x,t,p(x,t)). The most frequently used approach for describing the pressure dependence of the film thickness, is based on standard contact mechanics results for line and point loading of elastic half-spaces. These are detailed in Chapters 2 and 3 in the book by Johnson  and in . In addition, when the surfaces move relative to each other Ω varies in time. Figure 1.Schematic illustration of the original (left) and the transformed (right) domain.

The classical thin film approximation, presented in e.g. [1, p. 147], can be obtained by employing the scaling

X1=x1/l1,X2=x2/l2,Z=z/h0,T=t/t0,μ¯=μ/μ0,ρ¯=ρ/ρ0,u¯=u/u0,v¯=v/v0,w¯=w/w0,P=p/μ0u0l1h02.E4

More precisely, using the fact that h0 << liand neglecting terms of the order h0/liand (h0/li)2 leads to that p= p(x,t) and the following relationship between the pressure and the velocity field in the lubricant

z(μ(x,z,t)uz(x,z,t))=p(x,t),E5

where u=(u,v)and =(/x1,/x2). Remark: an alternative, one-dimensional, thin film approximation was introduced in . Expressed in the terminology used here, the difference between the models is a consequence of Barus’ pressure–viscosity relationship; μ=μ0exp(αp)and the scaling of the vertical velocity component;

w¯=w/(u0h0l1).

In the following, the viscosity is assumed to be on the form μ(x,z,t)=μ0g(z,p(x,t)), where the reference parameter μ0 is typically taken as the viscosity at ambient pressure and gis a function which can be used to describe the variation of the viscosity across the film (z-dependence) and with the pressure. It should be noted that both μ, ρand hmay depend on the fluid pressure.

As pointed out above, Ω varies with time. In order to simplify the analysis, a transformation of the height coordinate zis introduced that leads to a fixed, time independent, domain Ω′. Indeed, let

z=zz1hE6

then the corresponding solution domain becomes

Ω={(x,z)|0xili,0z1},E7

irrespectively of the shape of the surfaces hland hu. The domain Ω and its corresponding dimensionless form Ω′ are depicted in Figure 1. In this case, the correspondence to Eq. (5) becomes

1h2(x,t)z(μ(x,z,t)uz(x,z,t))=p(x,t).E8

Integrating Eq. (8) twice with respect to z′ and by assuming no slip at the surfaces, i.e., u(x,1,t)=vu(x,t)at the upper surface and u(x,0,t)=vl(x,t)at the lower surface, yields the following closed form expression for the velocity field

u(x,z,t)=(f1(x,z,t)f1(x,1,t)f0(x,1,t)f0(x,z,t))h2(x)p(x)E9
+f0(x,z,t)f0(x,1,t)(vu(x,t)vl(x,t))+vl(x,t),

where

fi(x,z,t)=0zsiμ(x,s,t)dsE10

and vu=(uu(x,t),vu(x,t))and vl=(ul(x,t),vl(x,t)).

The flow is regarded as compressible according to the “arbitrary” pressure–density relation

ρ(p)=ρcf(p).E11

The function f(p) is strictly increasing for ppcand it satisfies f(p) = 1 when ppc.

With the aim set to derive an Reynolds type of equation for the pressure in the lubricant, the analysis continues by first formulating an expression for the mass flow and thereafter requiring continuity of the mass flow. By using Eq. (9) for the velocity field and Eq. (11) the mass flow, qbecomes

q(x,t)=01ρ(p(x,t))u(x,z,t)h(x,t)dzE12
=ρh2avρh312μ0bp,

where v=vu+vl=(uu+ul,vu+vl),

a(x,t)=(a100a2),E13
a1(x,t)=2(uu+ul)1(01f0(x,z,t)dzf0(x,1,t)(uuul)ul),a2(x,t)=2(vu+vl)1(01f0(x,z,t)dzf0(x,1,t)(vuvl)vl),1b(x,t)=0112μ0(f1(x,z,t)f1(x,1,t)f0(x,1,t)f0(x,z,t))dz.E14

With this expression for q, the requirement for conservation of mass reads

(ρh)t=q.E15

In order to incorporate the effect of cavitation, we assume that the following holds. In the full film zones, the pressure is larger than the cavitation pressure and the density is expressed by Eq. (11). In the cavitation zones, the pressure equals the cavitation pressure and the density is interpreted as degree of saturation δ. The saturation δis an unknown in the cavitation zones and satisfies that 0 ≤ δ≤ 1 while δ= 1 in the full film zones. Hence the pressure–density relationship in both the full film and the cavitated zones becomes

ρ(p)=ρc{f(p),p>pcδ,p=pc.E16

Note that both pand δare unknowns.

For computational purposes it is beneficial to introduce a change of variables p¯and ηthat satisfy the complementarity condition p¯η=0in the whole domain. This can be accomplished by defining p¯and ηas

p¯=ppc,p¯0E17

and

η=1δ={0,p¯>01δ,p¯=0.E18

In this notation, Eq. (16) becomes

ρ(p¯)=ρc{f(p¯+pc),p¯>01η,p¯=0.E19

Now, because of complementarity, this piecewise definition of ρcan be rewritten as a single expression valid throughout the whole domain, i.e.,

ρ=ρc((1η)+f(p¯+pc)1)=ρc(f(p¯+pc)η)E20

which then leads to the following expression for the mass flow

q=ρc(f(p¯+pc)η)h2avρc(f(p¯+pc)η)h312μ0bp
=ρc(fh2avfh312μ0bp¯ηh2av),p¯0.E21

Preservation of mass flow is ensured by inserting Eq. (21) and Eq. (20) into the continuity equation Eq. (15), which leads to the following mass preserving cavitation model:

(fhηh)t=(fh2avfh312μ0bp¯ηh2av)p¯0,0η1,p¯η=0.E22

The system in Eq. (22) can be solved numerically by the LCP-based solution procedure described in Section 3. When the solution (p¯,η)has been obtained, the fluid pressure pand the saturation δcan be found from (17) and (18).

## 3. Numerical solution procedure

In this section, we present a numerical solution procedure for the cavitation model Eq. (22) such that the standard theory for linear complementary problems (LCP) can be applied.

Let us start by introducing the notation

α=fh312μ0b   and   β=(β1,β2)=h2av.

while keeping in mind that f, h, aand bdepend on p. In this notation the system Eq. (22) reads

(fh)t(ηh)t=(fβαp¯)(ηβ)p¯0,0η1,p¯η=0.E23

A spatial finite difference discretization of the problem Eq. (23) can be obtained by dividing the domain into a uniform rectangular grid with N1×N2elements of size Δx1×Δx2. The following notation is adopted x1i=iΔx1, x2j=jΔx2,where i=0,,N1, j=0,,N2and

ui,j:=u(x1i,x2j,t).E251

Since the finite difference approximations of the partial derivatives w.r.t. x1 and x2 are obtained in a similar manner, the details (given below) are only provided for the partial derivatives w.r.t. x1. The problem at hand, is elliptic in the full film domain, where η= 0. In the caviatted regions, where p¯=0, we observe that the equation is hyperbolic in η. A central difference scheme is therefore used to approximate the derivatives of the first term in the right hand side of Eq. (23) while an upwind difference scheme is employed for the second term.

In order to present the finite differences we introduce the notation

ai±1/2,j=ai±1,j+ai,j2E250

and the approximation

ux1|i+1/2,jui+1,jui,jΔx1andux1|i1/2,jui,jui1,jΔx1.E252

In this notation, the first term in the right hand side of Eq. (23) becomes:

x1(fβ1αp¯x1)E253
(fβ1)i+1,j(fβ1)i1,j2Δx1αi+1/2,jp¯x1|i+1/2,jαi1/2,jp¯x1|i1/2,jΔx1
(fβ1)i+1,j(fβ1)i1,j2Δx1
(αi+1,j+αi,j2)(p¯i+1,jp¯i,jΔx1)(αi1,j+αi,j2)(p¯i,jp¯i1,jΔx1)Δx1
=(fβ1)i+1,j(fβ1)i1,j2Δx11Δx12(αi1/2,jp¯i1,j(αi1/2,j+αi+1/2,j)p¯i,j+αi+1/2,jp¯i+1,j).E24

The upwind discretization of the second term reads

x1(ηβ1)ηi,jβ1i,jηi1,jβ1i1,jΔx1.E25

By using Eq. (24), Eq. (25) and the corresponding finite difference approximations of the partial derivatives w.r.t. x2, the spatially discretized (which is still continuous in time) form of the system Eq. (23) can be written as

((fh)t)i,j((ηh)t)i,j=(fβ1)i+1,j(fβ1)i1,j2Δx1+(fβ2)i,j+1(fβ2)i,j12Δx21Δx12(αi1/2,jp¯i1,j(αi1/2,j+αi+1/2,j)p¯i,j+αi+1/2,jp¯i+1,j)1Δx22(αi,j1/2p¯i,j1(αi,j1/2+αi,j+1/2)p¯i,j+αi,j+1/2p¯i,j+1)(ηβ1)i+1,j(ηβ1)i1,j2Δx1(ηβ2)i,j+1(ηβ2)i,j12Δx2,E26

together with the conditions

p¯i,j0,0ηi,j1,p¯i,jηi,j=0.

Since hand fdepend on p, a numerical solution procedure can be posed as follows:

1. Guess a pressure distribution.

2. Compute the coefficients f, hαand β, in order to linearize Eq. (26) and pose it on the form

p¯=q+Mη,p¯0,0η1,p¯η=0,

where the vector qand the matrix Mare constants and where the unknowns p¯and ηare vectors with (N11)(N21)elements.

3. Solve the now obtained linear complementarity problem corresponding to Eq. (26), with the Lemke algorithm .

To fully discretize the problem at hand several approaches can be applied. For instance, first order forward (explicit) or backward (implicit) Euler, the second order implicit Crank-Nicolson method.

## 4. Numerical examples

In this section, the numerical solution procedure, for the cavitation model in Eq. (22), developed in Section 3, is examined by considering four different one-dimensional slider bearing examples. In all four examples, only the lower surface is moving, i.e. v=(ul,0), the lubricant is assumed to be Newtonian, the length of the bearing is Land the load the bearing supports is w0. In the first example, results obtained with the present approach are compared to results for the one-dimensional parabolic slider in . Two different pressure–density relationships are examined. The second example extends the first one in the sense that besides cavitation the bearing surfaces are also assumed to be linear elastic. In particular, this means that the film thickness in addition depends on the pressure. The film thickness, in this case, consists of three parts; a parameter related to the minimum film thickness h0, the bearing geometry g, and the elastic deformation of the bearing surfaces. The Boussinesq-Cerruti half space solution is used to model the elastic deformation . Indeed, his given by

h(x,p(x))=h0+g(x)4πEln|xs|p(s)ds,E27

where

2E=1νlEl+1νuEu,E28

where vland vuare the Poisson’s ratios for the lower and upper surfaces and Eland Euare the corresponding Young’s modulus.

The third example considers a double parabolic slider, which in addition to the single parabolic slider, exhibits reformation and highlight that mass is conserved.

The fourth and last example considers a quadruple parabolic slider bearing. The reason for choosing this configuration is to test the hypothesis that an elastically deformable bearing, in general, does not generate more film than the corresponding rigid one.

In all examples, the initial (undeformed) bearing geometry consists of 1, 2 or 4 parabolic parts. These bearing geometries, gn, can then be described via the L-periodic auxiliary function G

gn(x)=G(x/n),n=1,2or4,E29

where

G(x)=h1(2L)2(xL2)2,0xL.E30

In Table 1, parameters common to all four examples are listed. These are the bearing parameters Land h1, the speed of the lower surface ul, the applied load w0, the cavitation pressure pc, the dynamic viscosity μand the boundary conditions for the normalized density at x= 0 and x= L, denoted by ρ0/ρcand ρL/ρc, respectively.

### Example 1

In this example, a model problem with rigid surfaces and two different Newtonian lubricants is considered in order to compare with previous results presented in . Indeed, a single parabolic slider bearing of length Lwith rigid surfaces and film thickness of the form

h(x)=h0+g1(x),E31

is considered here. Two different Newtonian lubricants are studied, one which obeys the constant bulk modulus pressure–density relationship;

ρ=ρcexp((ppc)/β)E32

while the other one obeys the Dowson-Higginson pressure–density relationship;

ρ=ρcC1+C2(ppc)C1+(ppc).E33

The parameters for the two different lubricants are given in Table 2. As in , the case with constant bulk modulus, β, was solved for a fixed film thickness, hβ, where h0 = h1 = 25.4 μm. The pressure solution, pβ, obtained in this way corresponds to a load carrying capacity w0 ≈ 117 kN/m. In the case with the Dowson-Higginson compressibility, force-balance was included and the load carrying capacity obtained for the constant bulk modulus case was used as the applied load. This resulted in the film thickness hDHwith the minimum h0 = 25.8 μm and the pressure pDH. In Figure 2 the film thickness and pressure distributions corresponding to the two different cases are depicted. It can be seen that the pressure distribution pβis in good agreement with the results in .

Lh1ulw0pcμρ0/ρcρL/ρc
76.2 mm25.4 μm4.57 m/s117 kN/m100 kPa0.039 Pas1.00011

### Table 1.

Parameters common to all four examples.

βC1C2
0.069 GPa2.22 GPa1.66

### Table 2.

Parameters specific to Example 1. Figure 2.Film thickness (left) and pressure distributions (right) for two lubricants obeying constant bulk modulusβand Dowson-Higgison pressure-density relationships, respectively. The solutions correspond to the single parabolic slider defined in Example 1.

### Example 2

This example extends the previous one, by including elastic deformation of the bearing surfaces. To illustrate the effect of surface deformation, the combined elastic modulus, E′, was chosen to be 10% of the one for two steel surfaces with Poisson’s ratio vand Young’s modulus E, listed in Table 3. The lubricant was assumed to obey the constant bulk modulus, pressure –density relationship, with β= 3.34 as listed in Table 3. Both of the bearings carry the same load as in Example 1, i.e., w0 = 117 kN/m.

βvE
3.34 GPa0.3210 GPa

### Table 3.

Parameters specific to Examples 2, 3 and 4.

In Figure 3, the film thickness heand pressure distribution pefor the elastic bearing are depicted together with the solutions hrand prfor the rigid bearing. In particular, it can be seen that the elastic deformation leads to a thicker film, a smaller zone of cavitation and a lower maximum pressure. The minimum film thickness for the rigid bearing is 25.8 μm and for the bearing with elastic surfaces it is 30.4 μm. Figure 3.Film thickness (left) and pressure distributions (right) for the rigid and the elastic single parabolic slider bearings in Example 2.

### Example 3

The only difference between the problem studied in this example and the one studied in Example 2, is that the bearing geometry now corresponds to a double parabolic slider, g2, defined through Eq. (29). This example was chosen to illustrate that the proposed LCP-based numerical method renders a solution, which corresponds to a mass conserving flow.

The film thickness heand pressure distribution pefor the elastic bearing are depicted together with the solutions hrand prfor the rigid bearing in Figure 4. In a mass conserving flow in a rigid double parabolic slider (described with g2), the maximum compression occurs twice, once underneath each of the two parabolas. Since the compression is equal in these two points, the pressure, which is the maximum pressure, is also equal in these two points. It is clear from Figure 4 that the pressure solution for the rigid bearing fulfils these conditions. Figure 4.Film thickness (left) and pressure distributions (right) for the rigid and the elastic double parabolic slider bearings in Example 3.

Also, as in Example 1, the bearing with elastic surfaces generates a thicker lubricant film, a smaller zone of cavitation and a lower maximum pressure.

### Example 4

The only difference between this example and Examples 2 and 3, is that the bearing geometry now corresponds to a quadruple parabolic slider, g4. In both Examples 2 and 3, the elastic bearing generated thicker films, smaller zone of cavitation and lower maximum pressures than the corresponding rigid ones. This example was constructed to show that this is not true in general. Indeed, as shown in Figure 5, the film thickness hefor the elastic bearing is thinner, some of the zones of cavitation become longer and the pressure underneath the first and the last parabola is higher than for the rigid bearing. Figure 5.Film thickness (left) and pressure distributions (right) for the rigid and the elastic quadruple parabolic slider bearings in Example 4.

## 5. Concluding remarks

A new mathematical model for thin film lubrication has been derived. The model is quite general, e.g., it accounts for lubricant compressibility, cavitation, pressure dependent viscosity, non-Newtonian rheology and elastic deformation of the surfaces. The main novelty of the model is that cavitation of a compressible fluid is considered via a formulation of the mass flow, which ultimately results in a complementarity problem. Hence, standard methods developed for linear complementarity problems can be used in the numerical solution procedure. The model’s applicability has been demonstrated in several numerical examples.