Modelling Cavitation in (Elasto)Hydrodynamic Lubrication

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.

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 v_u and the velocity of the lower surface is v_l, then under the assumptions outlined in e.g. [1], the following Reynolds equation is obtained

$\frac{\partial }{\partial t}\left(\rho h\right)=\nabla \cdot \left(\frac{\rho h^3}{12\mu }\nabla p-\frac{1}{2}\rho h\left(v_u+v_l\right)\right)\mathrm{,}$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 [2] 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 p_c.

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 [3–5], Elrod’s work [6] 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 [7] 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

$\rho ={\rho }_{c}exp\left(\left(p-p_c\right)/\beta \right)\mathrm{,}$E2

where ρ_c is the density at cavitation pressure p_c, and β is the bulk modulus of the lubricant. One example of a generalization of the results in [6] and [7] is the paper [8] 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 [6] 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. [12]. 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. [15]. 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. [14], 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 [12] the approach in Almqvist et al. [14] 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 [14] 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 h_l (lower) and h_u (upper):

$\Omega =\left\{\left(x\mathrm{,}z\right)|0\le x_i\le l_i\mathrm{,}{h}_l\le z\le h_u\right\}\mathrm{.}$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 = h_u − h_l depends 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 [16] and in [2]. 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

$\begin{array}{c}{X}_1=x_1/l_1\mathrm{,}{X}_2=x_2/l_2\mathrm{,}Z=z/h_0\mathrm{,}T=t/t_0\mathrm{,}\overline{\mu }=\mu /{\mu }_0\mathrm{,}\overline{\rho }=\rho /{\rho }_0\mathrm{,}\\ \\ \overline{u}=u/u_0\mathrm{,}\overline{v}=v/v_0\mathrm{,}\overline{w}=w/w_0\mathrm{,}P=p/\frac{{\mu }_0{u}_0{l}_1}{h_0^2}\mathrm{.}\end{array}$E4

More precisely, using the fact that h₀ << l_i and neglecting terms of the order h₀/l_i and (h₀/l_i)² leads to that p = p(x,t) and the following relationship between the pressure and the velocity field in the lubricant

$\frac{\partial }{\partial z}\left(\mu \left(x\mathrm{,}z\mathrm{,}t\right)\frac{\partial u}{\partial z}\left(x\mathrm{,}z\mathrm{,}t\right)\right)=\nabla p\left(x\mathrm{,}t\right)\mathrm{,}$E5

where $u=(u\mathrm{,}v)$ and $\nabla =(\partial /\partial x_1\mathrm{,}\partial /\partial x_2)$. Remark: an alternative, one-dimensional, thin film approximation was introduced in [17]. Expressed in the terminology used here, the difference between the models is a consequence of Barus’ pressure–viscosity relationship; $\mu ={\mu }_{0}exp\left(\alpha p\right)$ and the scaling of the vertical velocity component;

$\begin{array}{c}\overline{w}=w/\left(u_0\frac{h_0}{l_1}\right)\mathrm{.}\end{array}$

In the following, the viscosity is assumed to be on the form $\mu \left(x\mathrm{,}z\mathrm{,}t\right)={\mu }_{0}g\left(z\mathrm{,}p\left(x\mathrm{,}t\right)\right)$, where the reference parameter μ₀ is typically taken as the viscosity at ambient pressure and g is 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 h may depend on the fluid pressure.

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

$\begin{array}{l}{z}^{\prime }=\frac{z-z_1}{h}\hfill \end{array}$E6

then the corresponding solution domain becomes

${\Omega }^{\prime }=\left\{\left(x\mathrm{,}z\right)|0\le x_i\le l_i\mathrm{,0}\le z^{\prime }\le 1\right\}\mathrm{,}$E7

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

$\frac{1}{h^2\left(x\mathrm{,}t\right)}\frac{\partial }{\partial z^{\prime }}\left(\mu \left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)\frac{\partial u}{\partial z^{\prime }}\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)\right)=\nabla p\left(x\mathrm{,}t\right)\mathrm{.}$E8

Integrating Eq. (8) twice with respect to z′ and by assuming no slip at the surfaces, i.e., $u\left(x\mathrm{,1,}t\right)=v_u\left(x\mathrm{,}t\right)$ at the upper surface and $u\left(x\mathrm{,0,}t\right)=v_l\left(x\mathrm{,}t\right)$ at the lower surface, yields the following closed form expression for the velocity field

$u\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)=\left(f_1\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)-\frac{f_1\left(x\mathrm{,1,}t\right)}{f_0\left(x\mathrm{,1,}t\right)}{f}_0\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)\right)h^2\left(x\right)\nabla p\left(x\right)$E9

$+\frac{f_0\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)}{f_0\left(x\mathrm{,1,}t\right)}\left(v_u\left(x\mathrm{,}t\right)-v_l\left(x\mathrm{,}t\right)\right)+v_l\left(x\mathrm{,}t\right)\mathrm{,}$

where

$f_i\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)={\displaystyle {\int }_0^{z^{\prime }}}\frac{s^i}{\mu \left(x\mathrm{,}s\mathrm{,}t\right)}ds$E10

and $v_u=(u_u\left(x\mathrm{,}t\right)\mathrm{,}{v}_u\left(x\mathrm{,}t\right))$ and $v_l=(u_l\left(x\mathrm{,}t\right)\mathrm{,}{v}_l\left(x\mathrm{,}t\right))$.

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

$\begin{array}{l}\rho \left(p\right)={\rho }_{c}f\left(p\right)\mathrm{.}\hfill \end{array}$E11

The function f(p) is strictly increasing for p ≥ p_c and it satisfies f(p) = 1 when p ≤ p_c.

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, q becomes

$q\left(x\mathrm{,}t\right)={\displaystyle {\int }_0^1}\rho \left(p\left(x\mathrm{,}t\right)\right)u\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)h\left(x\mathrm{,}t\right)d{z}^{\prime }$E12

$\begin{array}{l}=\frac{\rho h}{2}av-\frac{\rho h^3}{12{\mu }_{0}b}\nabla p\mathrm{,}\hfill \end{array}$

where $v=v_u+v_l=(u_u+u_l\mathrm{,}{v}_u+v_l)$,

$a\left(x\mathrm{,}t\right)=\left(\begin{array}{ll}{a}_1\hfill & 0\hfill \\ 0\hfill & a_2\hfill \end{array}\right)\mathrm{,}$E13

$\begin{array}{ll}{a}_1\left(x\mathrm{,}t\right)\hfill & =2{\left(u_u+u_l\right)}^{-1}\left(\frac{{\displaystyle {\int }_0^1}{f}_0\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)d{z}^{\prime }}{f_0\left(x\mathrm{,1,}t\right)}\left(u_u-u_l\right)-u_l\right)\mathrm{,}\hfill \\ \hfill & \hfill \\ a_2\left(x\mathrm{,}t\right)\hfill & =2{\left(v_u+v_l\right)}^{-1}\left(\frac{{\displaystyle {\int }_0^1}{f}_0\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)d{z}^{\prime }}{f_0\left(x\mathrm{,1,}t\right)}\left(v_u-v_l\right)-v_l\right)\mathrm{,}\hfill \\ \hfill & \hfill \\ \frac{1}{b\left(x\mathrm{,}t\right)}\hfill & ={\displaystyle {\int }_0^1}12{\mu }_0\left(f_1\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)-\frac{f_1\left(x\mathrm{,1,}t\right)}{f_0\left(x\mathrm{,1,}t\right)}{f}_0\left(x\mathrm{,}{z}^{\prime }\mathrm{,}t\right)\right)d{z}^{\prime }\mathrm{.}\hfill \end{array}$E14

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

$\frac{\partial \left(\rho h\right)}{\partial t}=\nabla \cdot q\mathrm{.}$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

$\rho \left(p\right)={\rho }_c\{\begin{array}{ll}f\left(p\right)\hfill & \mathrm{,}p>p_c\hfill \\ \delta \hfill & \mathrm{,}p=p_c\hfill \end{array}\mathrm{.}$E16

Note that both p and δ are unknowns.

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

$\overline{p}=p-p_c\mathrm{,}\overline{p}\ge 0$E17

and

$\eta =1-\delta =\{\begin{array}{ll}0\hfill & \mathrm{,}\overline{p}>0\hfill \\ 1-\delta \hfill & \mathrm{,}\overline{p}=0\hfill \end{array}\mathrm{.}$E18

In this notation, Eq. (16) becomes

$\rho \left(\overline{p}\right)={\rho }_c\{\begin{array}{ll}f\left(\overline{p}+p_c\right)\hfill & \mathrm{,}\overline{p}>0\hfill \\ 1-\eta \hfill & \mathrm{,}\overline{p}=0\hfill \end{array}\mathrm{.}$E19

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

$\rho ={\rho }_c\left(\left(1-\eta \right)+f\left(\overline{p}+p_c\right)-1\right)={\rho }_c\left(f\left(\overline{p}+p_c\right)-\eta \right)$E20

which then leads to the following expression for the mass flow

$\begin{array}{l}q=\frac{{\rho }_c\left(f\left(\overline{p}+p_c\right)-\eta \right)h}{2}av-\frac{{\rho }_c\left(f\left(\overline{p}+p_c\right)-\eta \right)h^3}{12{\mu }_{0}b}\nabla p\hfill \end{array}$

$={\rho }_c\left(\frac{fh}{2}av-\frac{f{h}^3}{12{\mu }_{0}b}\nabla \overline{p}-\frac{\eta h}{2}av\right)\mathrm{,}\overline{p}\ge 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:

$\begin{array}{l}\frac{\partial \left(fh-\eta h\right)}{\partial t}=\nabla \cdot \left(\frac{fh}{2}av-\frac{f{h}^3}{12{\mu }_{0}b}\nabla \overline{p}-\frac{\eta h}{2}av\right)\hfill \\ \hfill \\ \overline{p}\ge \mathrm{0,}0\le \eta \le \mathrm{1,}\overline{p}\eta =0.\hfill \end{array}$E22

The system in Eq. (22) can be solved numerically by the LCP-based solution procedure described in Section 3. When the solution $\left(\overline{p}\mathrm{,}\eta \right)$ has been obtained, the fluid pressure p and 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

$\alpha =\frac{f{h}^3}{12{\mu }_{0}b}\text{ and }\beta =({\beta }_1,{\beta }_2)=\frac{h}{2}av.$

while keeping in mind that f, h, a and b depend on p. In this notation the system Eq. (22) reads

$\begin{array}{l}\frac{\partial \left(fh\right)}{\partial t}-\frac{\partial \left(\eta h\right)}{\partial t}=\nabla \cdot \left(f\beta -\alpha \nabla \overline{p}\right)-\nabla \cdot \left(\eta \beta \right)\hfill \\ \hfill \\ \overline{p}\ge \mathrm{0,}0\le \eta \le \mathrm{1,}\overline{p}\eta =0.\hfill \end{array}$E23

A spatial finite difference discretization of the problem Eq. (23) can be obtained by dividing the domain into a uniform rectangular grid with $N_1\times N_2$ elements of size $\Delta x_1\times \Delta x_2$. The following notation is adopted $x_1^i=i\Delta x_1$, $x_2^j=j\Delta x_2\mathrm{,}$ where $i=\mathrm{0,}\dots \mathrm{,}{N}_1$, $j=\mathrm{0,}\dots \mathrm{,}{N}_2$ and

$u_{i\mathrm{,}j}:=u\left(x_1^i\mathrm{,}{x}_2^j\mathrm{,}t\right)\mathrm{.}$E251

Since the finite difference approximations of the partial derivatives w.r.t. x₁ and x₂ are obtained in a similar manner, the details (given below) are only provided for the partial derivatives w.r.t. x₁. The problem at hand, is elliptic in the full film domain, where η = 0. In the caviatted regions, where $\overline{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

$a_{i\pm 1/\mathrm{2,}j}=\frac{a_{i\pm \mathrm{1,}j}+a_{i\mathrm{,}j}}{2}$E250

and the approximation

${\frac{\partial u}{\partial x_1}|}_{i+1/\mathrm{2,}j}\approx \frac{u_{i+\mathrm{1,}j}-u_{i\mathrm{,}j}}{\Delta x_1}\text{and}{\frac{\partial u}{\partial x_1}|}_{i-1/\mathrm{2,}j}\approx \frac{u_{i\mathrm{,}j}-u_{i-\mathrm{1,}j}}{\Delta x_1}\mathrm{.}$E252

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

$\frac{\partial }{\partial x_1}\left(f{\beta }_1-\alpha \frac{\partial \overline{p}}{\partial x_1}\right)$E253

$\approx \frac{{\left(f{\beta }_1\right)}_{i+\mathrm{1,}j}-{\left(f{\beta }_1\right)}_{i-\mathrm{1,}j}}{2\Delta x_1}-\frac{{{\alpha }_{i+1/\mathrm{2,}j}\frac{\partial \overline{p}}{\partial x_1}|}_{i+1/\mathrm{2,}j}-{{\alpha }_{i-1/\mathrm{2,}j}\frac{\partial \overline{p}}{\partial x_1}|}_{i-1/\mathrm{2,}j}}{\Delta x_1}$

$\approx \frac{{\left(f{\beta }_1\right)}_{i+\mathrm{1,}j}-{\left(f{\beta }_1\right)}_{i-\mathrm{1,}j}}{2\Delta x_1}$

$-\frac{\left(\frac{{\alpha }_{i+\mathrm{1,}j}+{\alpha }_{i\mathrm{,}j}}{2}\right)\left(\frac{{\overline{p}}_{i+\mathrm{1,}j}-{\overline{p}}_{i\mathrm{,}j}}{\Delta x_1}\right)-\left(\frac{{\alpha }_{i-\mathrm{1,}j}+{\alpha }_{i\mathrm{,}j}}{2}\right)\left(\frac{{\overline{p}}_{i\mathrm{,}j}-{\overline{p}}_{i-\mathrm{1,}j}}{\Delta x_1}\right)}{\Delta x_1}$

$\begin{array}{c}=\frac{{\left(f{\beta }_1\right)}_{i+\mathrm{1,}j}-{\left(f{\beta }_1\right)}_{i-\mathrm{1,}j}}{2\Delta x_1}\\ \\ -\frac{1}{\Delta x_1{}^2}\left({\alpha }_{i-1/\mathrm{2,}j}{\overline{p}}_{i-\mathrm{1,}j}-\left({\alpha }_{i-1/\mathrm{2,}j}+{\alpha }_{i+1/\mathrm{2,}j}\right){\overline{p}}_{i\mathrm{,}j}+{\alpha }_{i+1/\mathrm{2,}j}{\overline{p}}_{i+\mathrm{1,}j}\right)\mathrm{.}\end{array}$E24

The upwind discretization of the second term reads

$\frac{\partial }{\partial x_1}\left(\eta {\beta }_1\right)\approx \frac{{\eta }_{i\mathrm{,}j}{\beta }_1{}_{i\mathrm{,}j}-{\eta }_{i-\mathrm{1,}j}{\beta }_1{}_{i-\mathrm{1,}j}}{\Delta x_1}\mathrm{.}$E25

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

$\begin{array}{c}{\left(\frac{\partial \left(fh\right)}{\partial t}\right)}_{i\mathrm{,}j}-{\left(\frac{\partial \left(\eta h\right)}{\partial t}\right)}_{i\mathrm{,}j}=\frac{{\left(f{\beta }_1\right)}_{i+\mathrm{1,}j}-{\left(f{\beta }_1\right)}_{i-\mathrm{1,}j}}{2\Delta x_1}+\frac{{\left(f{\beta }_2\right)}_{i\mathrm{,}j+1}-{\left(f{\beta }_2\right)}_{i\mathrm{,}j-1}}{2\Delta x_2}\\ \\ -\frac{1}{\Delta x_1{}^2}\left({\alpha }_{i-1/\mathrm{2,}j}{\overline{p}}_{i-\mathrm{1,}j}-\left({\alpha }_{i-1/\mathrm{2,}j}+{\alpha }_{i+1/\mathrm{2,}j}\right){\overline{p}}_{i\mathrm{,}j}+{\alpha }_{i+1/\mathrm{2,}j}{\overline{p}}_{i+\mathrm{1,}j}\right)\\ \\ -\frac{1}{\Delta x_2{}^2}\left({\alpha }_{i\mathrm{,}j-1/2}{\overline{p}}_{i\mathrm{,}j-1}-\left({\alpha }_{i\mathrm{,}j-1/2}+{\alpha }_{i\mathrm{,}j+1/2}\right){\overline{p}}_{i\mathrm{,}j}+{\alpha }_{i\mathrm{,}j+1/2}{\overline{p}}_{i\mathrm{,}j+1}\right)\\ \\ -\frac{{\left(\eta {\beta }_1\right)}_{i+\mathrm{1,}j}-{\left(\eta {\beta }_1\right)}_{i-\mathrm{1,}j}}{2\Delta x_1}-\frac{{\left(\eta {\beta }_2\right)}_{i\mathrm{,}j+1}-{\left(\eta {\beta }_2\right)}_{i\mathrm{,}j-1}}{2\Delta x_2}\mathrm{,}\end{array}$E26

together with the conditions

${\overline{p}}_{i\mathrm{,}j}\ge \mathrm{0,}0\le {\eta }_{i\mathrm{,}j}\le \mathrm{1,}{\overline{p}}_{i\mathrm{,}j}{\eta }_{i\mathrm{,}j}=0.$

Since h and f depend 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

$\overline{p}=q+M\eta \mathrm{,}\overline{p}\ge \mathrm{0,}0\le \eta \le \mathrm{1,}\overline{p}\cdot \eta =\mathrm{0,}$

where the vector q and the matrix M are constants and where the unknowns $\overline{p}$ and η are vectors with $\left(N_1-1\right)\cdot \left(N_2-1\right)$ elements.

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

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=(u_l\mathrm{,0})$, the lubricant is assumed to be Newtonian, the length of the bearing is L and the load the bearing supports is w₀. In the first example, results obtained with the present approach are compared to results for the one-dimensional parabolic slider in [8]. 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 h₀, 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 [2]. Indeed, h is given by

$h\left(x\mathrm{,}p\left(x\right)\right)=h_0+g\left(x\right)-\frac{4}{\pi E^{\prime }}{\displaystyle {\int }_{-\infty }^{\infty }}ln\left|x-s\right|p\left(s\right)ds\mathrm{,}$E27

where

$\frac{2}{E^{\prime }}=\frac{1-{\nu }_l}{E_l}+\frac{1-{\nu }_u}{E_u}\mathrm{,}$E28

where v_l and v_u are the Poisson’s ratios for the lower and upper surfaces and E_l and E_u are 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, g_n, can then be described via the L-periodic auxiliary function G

$g_n\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}G\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo>/</mo>}n\mathrm{<mo\; stretchy="false">)</mo>,}n\mathrm{<mo>=</mo>1,2}\text{or}\mathrm{4,}$E29

where

$G(x)=h_1{\left(\frac{2}{L}\right)}^2{\left(x-\frac{L}{2}\right)}^2\mathrm{,}0\le x\le L\mathrm{.}$E30

In Table 1, parameters common to all four examples are listed. These are the bearing parameters L and h₁, the speed of the lower surface u_l, the applied load w₀, the cavitation pressure p_c, the dynamic viscosity μ and the boundary conditions for the normalized density at x = 0 and x = L, denoted by ρ₀/ρ_c and ρ_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 [8]. Indeed, a single parabolic slider bearing of length L with rigid surfaces and film thickness of the form

$h\left(x\right)=h_0+g_1\left(x\right)\mathrm{,}$E31

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

$\rho ={\rho }_c\text{exp}(\left(p-p_c\right)/\beta )$E32

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

$\rho ={\rho }_c\frac{C_1+C_2\left(p-p_c\right)}{C_1+\left(p-p_c\right)}\mathrm{.}$E33

The parameters for the two different lubricants are given in Table 2. As in [8], the case with constant bulk modulus, β, was solved for a fixed film thickness, h_β, where h₀ = h₁ = 25.4 μm. The pressure solution, p_β, obtained in this way corresponds to a load carrying capacity w₀ ≈ 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 h_DH with the minimum h₀ = 25.8 μm and the pressure p_DH. 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 [8].

L	h1	ul	w0	pc	μ	ρ0/ρc	ρL/ρc
76.2 mm	25.4 μm	4.57 m/s	117 kN/m	100 kPa	0.039 Pas	1.0001	1

Table 1.

Parameters common to all four examples.

β	C1	C2
0.069 GPa	2.22 GPa	1.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 v and 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., w₀ = 117 kN/m.

β	v	E
3.34 GPa	0.3	210 GPa

Table 3.

Parameters specific to Examples 2, 3 and 4.

In Figure 3, the film thickness h_e and pressure distribution p_e for the elastic bearing are depicted together with the solutions h_r and p_r for 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, g₂, 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 h_e and pressure distribution p_e for the elastic bearing are depicted together with the solutions h_r and p_r for the rigid bearing in Figure 4. In a mass conserving flow in a rigid double parabolic slider (described with g₂), 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, g₄. 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 h_e for 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.

Acknowledgments

The authors want to acknowledge The Swedish Research Council (VR) for the financial support through the grants no. 2013-4978 and 2014-4894.