Numerical Simulation of the Cavitation in the Hydrodynamic Lubrication of Journal Bearings: a Parallel Algorithm

Numerical Simulation of the Cavitation in the Hydrodynamic Lubrication of Journal Bearings: a Parallel Algorithm

Mariangela Amendola

DPLPAG / FEAGRI

State University of Campinas, Campinas, SP. Brazil

[e-mail]

Carlos Antonio de Moura

IC / UFF

Fluminense Federal University, Niterói, RJ. Brazil

José Vitorio Zago

Petronio Pulino

Francisco de Assis M. Gomes Neto

DMA / IMECC

State University of Campinas, Campinas, SP. Brazil

In this paper we present an algorithm for the numerical simulation of the cavitation in the hydrodynamic lubrication of journal bearings. Despite the fact that this physical process is usually modelled as a free boundary problem, we adopted the equivalent variational inequality formulation.

We propose a two-level iterative algorithm, where the outer iteration is associated to the penalty method, used to transform the variational inequality into a variational equation, and the inner iteration is associated to the conjugate gradient method, used to solve the linear system generated by applying the finite element method to the variational equation. This inner part was implemented using the element by element strategy, which is easily parallelized.

We analyse the behavior of two physical parameters and discuss some numerical results. Also, we analyse some results related to the performance of a parallel implementation of the algorithm.

Keywords: Cavitation, parallel algorithm, variational inequalities, finite element method, penalty method

Introduction

In this work, we deal with the numerical simulation of cavitation in the hydrodynamic lubrication of journal bearings. This physical process is usually modelled as a free boundary problem.

Free boundary problems is a common designation for a class of mathematical problems used to represent physical processes whose numerical simulation requires not only the solution of a differential equation with boundary conditions, but also the determination of the boundary.

We present this work according to the following steps. First we describe the physical problem to be solved. Then, from the fluid mechanic classic laws, the Reynolds theory of hydrodynamic lubrication, the practical observations of the cavitation inside radial journal bearings and specifics considerations for the problem described, we show not only the usual differential formulation but also the equivalent variational inequality formulation, focusing their theoretical advantages. After that, we propose an algorithm to search the numerical solution of the variational inequality, keeping in mind its potential to be implemented in parallel. Then we discuss two ways of implementing the algorithm. Afther that, we show the numerical analysis of some physical parameters of the problem and the performance of the parallel algorithm. Finally, we present some conclusions.

Nomenclature

b = L/2r2

c = r₁ - r2

D = 2r2

e = excentricity

FD = 2c

h = radial lenght of Q

L = journal bearing lenght

p = lubricating fluid pressure

p₀ = athmospheric pressure

p_i = pressure p in each point i of the grid Q_h of Q

Q = region between S₁ and S2

Q_h = grid of Q

= Q È ¶Q

S₁ = external cylindrical surface

S₂ = internal cylindrical surface

r_i (i = 1,2) = cylinder's radius

Greek Letters

q = x/r₂ angle measured from a line that passes through the centers of both cylinders

e = e/c ratio of excentricity

m [Pa.s] = lubricating fluid viscosity

w [s^-1] = S₂ angular velocity

W = region of Q where p=0.

Physical Description of the Problem

The unknow of the problem is the pressure p[Pa] of a lubricating fluid with viscosity m [Pa.s] inside one spatial domain Q between two cylindrical surfaces S₁ and S₂, with length L[m] at atmospheric pressure at the ends. The axis of S₁ and S₂ are parallel but not coincident, and the internal cylinder S₂ rotates anticlockwise, with angular velocity w [s^-1]around its axis.

Before starting the movement, due to gravity, the cylinders are in contact (Fig. 1a).

At the begining of the movement, due to the adhesive properties of the lubricating fluid, a kind of dragging of the fluid starts, and S₂ starts to slip onto S₁, with a clockwise movement (Fig. 1b). As S₂ comes to more lubricated surfaces it starts sliding anticlockwise. After that, the movement of S₂ gradually generates a hydraulic pressure in the narrowest area between the cylinders that lifts S₂, eliminating that contact (Fig. 1c), and reducing the effort required to keep S₂ in movement. This process continues until an equilibrium is established and a certain clearance between the cylinders is obtained.

In normal conditions of operation, the pressure of the fluid reaches a minimum value associated to the vapor pressure. Below this pressure, cavitation occurs, and can lead to a damage to the journal bearings.

Mathematical Formulation of the Problem

The flow of a fluid in tridimensional movement is specified by its velocity U = (u, v, w) [m/s], pressure p and density p[kg.m^-3], all described in terms of x, y, z and t. To determine these five functions we have the same number of equations: the continuity equation, three movement equations and the thermodynamic state equation. For the special case of lubricating fluids some simplifications can be made on these equations since external forces need not to be considered and inertial forces can be disregarded when compared to viscous forces. As a result, the so called Reynolds equation is obtained (see Schlighting, 1968).

Further simplifications can be made to the model, due to the geometry of the journal bearings. Refering to S1, we can consider that the characteristic radial length h[m] of the clearance of region Q between the cylinders is very small if compared to the other radial lengths r1, r2 and L (h ~ (L/100)), where ri[m] is the radius of the cylinder Si. Consequently we can: 1) disregard curvature effects and use cartesian coordinates in the plane of the lubricating film; 2) consider that the variation of p through the film is small and can be disregarded; 3) suppose that the variation rate of any component of velocity in the plane of the film is smaller than that one through the film; 4) notice that the flow of the fluid in the film is mainly bidimensional, and the component v of the velocity of the fluid through the film can be disregarded.

After some mathematical manipulation that join and simplify all of these considerations, according to Amendola (1996), we obtain the equation

For the problem here described, it is convenient to rewrite this equation using cylindrical coordinates (q, z), taking the center of S₂ as the origin, its axis as direction z, and measuring the angle q anticlockwise from a line that passes through the centers of both cylinders, as shown in Fig.1.

So, establising the adimentional quantity

and introducing

where p₀ is the athmospheric pressure, e is the excentricity, c = r₁ – r₂, h = c + e cos q, e = e/c is the ratio of excentricity, D = 2r₂ and FD = 2c is the diametral clearance, we obtain the following Reynolds equation for the hydrodynamic pressure of the lubricating fluid of the journal bearings:

The associated boundary conditions are specified as

These boundary conditions, as well as the domain of the problem, are shown in Fig. 2.

Such conditions are not enough to solve Eq.(3), since we know that the fluid can not support subathmospheric pressure. So, we have to impose the condition of state transition on the boundary of the cavitation region. This is the so called Reynolds condition and it means that the pressure and its normal derivative are null on that boundary.

The differencial mathematical formulation of the process described above will be referred to as Problem 1.

Problem 1: Detrmine the function p = p (q, z) in and W Ì Q so that

where ¶ / ¶n is the normal derivative,

Q = {(q, z); -¥ < q < ¥ , ç z ç< b}Ì Â² ,

= Q ȶQ and

is the set of real and periodic continuous functions with continuous first and second derivative in .

However, even with this complete physical and mathematical formulation, there is not enough mathematical theory to assure the existence and unicity of the solution of the Problem 1. This is because the problem requires not only a mathematical methodology to find a solution for differential equations, but also an iterative procedure capable to locate the free boundary. As noticed by Orr and Scriven (1978), all iterative schemes employ a similar strategy: first a location for the free boundary is chosen, then the differential equation is solved, but only two of the three boundary conditions are satisfied, and finally the residual in the third boundary condition is used to decide how to alter the location of the free boundary. This process is repeated to convergence or frustation.

Understanding the general difficulty to assure convergence of this kind of scheme, we decided to choose another mathematical formulation to search the solution we are looking for.

Fortunately, this problem can be formulated by an equivalent variational inequality. In this case, Reynolds condition appears as a natural condition, it needs not to be imposed as in Eq.(6.3) (Rhode and McAllister, 1975).

This new formulation will be referred to as Problem 2.

Problem 2: Determine p Î K that

where

is the set of real and periodic continuous functions, and with continuous first derivatives in ,

is the set of functions of that are null in ¶Q,

is the completion of relative to the norm

, and

K = {v Î; v ³ 0 a.e in Q}convex and closed.

This problem has only one solution (Kinderleher and Stampacchia, 1980). Besides, it is not necessary to define an iterative method to fit the free boundary, since the cavitation region appears naturally at the solution.

Obviously, this choice is not cost free, as it demands a more complex mathematical way of writing and reeding the problem.

Numerical Solution of the Problem

In the process of choosing elementary numerical methods to compose our algorithm, one of the factors that we considered is the potential for their use in parallel.

Thus, focusing the use of the Finite Element Method (F.E.M.), we decided to transform the variational inequality into a variational equation in order to keep the problem ready to be solved by the F.E.M., suggested as ideal to be implemented in parallel (see Ciarlet and Lions (1991)). To perform this transformation we have chosen the Penalty Method (P.M.), as it was proposed in Moura and Amendola (1993).

The P.M. requires the introduction of an auxiliary functional that involves the definition of an Î parameter and a related function p_Î of p so that Î ® 0, p_Î ® p.

The F.E.M. gives us a systematic way of choosing basic functions for the definition of an subespace of finite dimension, denoted V_h, in order to give us the penalized solution p_hÎ as a linear combination of those functions over the grid Q_h of Q (see Becker, Carey and Oden (1981)).

Following a common F.E.M approach, here Q_h is taken as divided into triangles and the solution is sought as an expansion in terms of linear trial functions, which are defined within neighbouring subdomains (elements) and vanish elsewhere.

Using these ideas, our aim is to solve the equivalent penalized Galerkin problem, referred to as Problem 3.

Problem 3: Determine p_hÎ ÎV_hÌ where V_his a subespace of finite dimension, so that

where

Using the F.E.M. with f_j = f_j(q, z) as basic functions of V_h , Problem 3 is transformed into one system of algebric equations on p_i, where p_i represents p_hÎ computed at each point of the grid Q_h, i = 1,...N,. This system will be referred to as Problem 4.

Problem 4: Determine {p_i} Î Âⁿ so that

where

and the associated boundary conditions on the grid are supposed to be satisfied. The integrals are evaluated numerically by a Gaussian quadrature formula.

Computationally intensive components of the F.E.M. are the local calculations of the submatrices and right hand side vectors of the finite elements, the accumulation of these contributions to the global algebric system and the solution of this system. Thus, we have to be aware of these components when building an efficient algorithm.

In a quite different way from the usual implementation of the F.E.M., we followed the suggestion of Barragy and Carey (1988) and Carey et al. (1988), and strategically use the contribution to A of each element e, A_eindependently when solving linear systems, without explicitly assembling matrix A. Besides, the boundary conditions are also applied independently when computing the contributions to the global vector F. This will be referred to as the element by element (E.B.E.) strategy.

Using another notation, Problem 4 can be rewritten as

Problem 4': Determine p Î Â^N so that

where A is a symetric positive definite matrix.

The global system in Eq.(11) can be understood as Ãp = F, where à is formed adding the value 1/Î to some elements of the diagonal of A, according to the definition of .

Among the methods available to solve Eq.(11), we have chosen the Conjugate Gradient Method (C.G.M.) (see Hageman and Young (1981)), since its computationally intensive part, which consists in computing the product of matrix Ã. by a vector, can also be done using the E.B.E. strategy.

We also decided to use a diagonal preconditioner since it greatly improves the performance of the method and can be easily added according to E.B.E. strategy. Here, the C.G.M. with a diagonal preconditioner will be named P.C.G.M.

To clarify what we have discussed till now, we propose a two-level iterative algorithm, as we show below. The Algorithm 1, related to the P.M., is used as an external iterative processes and the Algorithm 2, related to the P.C.G.M. according to the E.B.E. strategy, appears as an internal iterative process.

Algorithm 1. (Penalty Method)

0. Set Î and l = 0

1. Form the right hand side vector F

2. p⁰¬ F

3. Repeat

3.1 l ¬ l + 1

3.2 Find p^l, an approximated solutuin for

3.3 Î ¬ Î /10

4. Until convergence is attained

Given NTEL, the number of finite elements of Q_h, the escalar a⁰, and vectors p⁰, r⁰ and u⁰, the following algorithm computes step 3.2 of Algorithm 1.

Algorithm 2. (P.C.G.M. according to the E.B.E strategy)

0. k ¬ 0

1. Repeat

1.1. For e = 1, ..., NTEL, do

1.1.1. Compute

1.1.2. Apply the boundary conditions to

1.1.3. Form the preconditioner diag

1.1.4.

1.2. end do

1.3. Assemble vector y and materiz M^k

1.4. g ¬ á y, u^kñ

1.5. step ¬ a ^k /g

1.6. p^k+1¬ p^k + step*u^k

1.7. r ^k+1¬ r k – step*y

1.8. w ¬ M^-1 * r^k+1

1.9. a ^{k +1}¬ ár ^k+1 , wñ

1.10. b ¬ a ^{k +1} / a ^k

1.11. u^{k +1}¬ w + b * u^k

1.12. k ¬ k+1

2. Until convergence is attained.

At the first iteration of the P.M., p⁰ ¬ F, r⁰¬ F - Ã * p⁰, u⁰¬ M^-1 * r⁰ and a ⁰¬ ár⁰ _, u⁰ñ . For an iteration l > 0, we use p, r, u and a generated at iteration l - 1 as initial values.

The P.C.G.M. stops when

or

For the P.C.G.M., that is the internal process, we define the convergence parameter, d_cg, as a function of the dimension of the grid, which is function of the number of finite elements we want in the direction q and z. This has to be done in terms of memory availability for the computer being used.

The P.M. stops when

For the P. M., that is the external iterative method, if in the beginning we take the penalty parameter,Î , as d_cg10², we can define the convergence parameter,d_p, as d_cg 10.

Two Ways of Implementing the Algorithm

To implement the algorithm described above, we must be aware of the availability of computer memory and processing time. Therefore, the dimension of the problem is decisive.

In case of small problems, we can keep the submatrices of all finite elements in memory, without forming and saving the global one. This will be called the E.B.E.-1 strategy.

For very large problems, we can not store the whole matrix A in memory so we need to use the E.B.E.-2 strategy, which means that we recompute the submatrices of the finite elements on each iteration of the P.C.G.M. This reduces the memory requirements but increases the processing time. This is the strategy we implemented in parallel.

To implement the E.B.E.-2 strategy, we have used the master-slave paradigm. In our parallel algorithm, one process, called master, works as a manager of the tasks that are to be actually done independent and simultaneously by the other processes, named slaves. The master and the slaves processes must execute specific instructions, in particular those related to the communication between the master and each slave.

The algorithm was validated for a simple problem specially formulated for this purpose.

Since the step 1.1 of Algorithm 2 takes about 98% of the total time spent by the whole program, and since the communication overhead is significant for the computer we are using, as it will become clear later, we decided to parallelize only this step.

Besides, we intend to run our algorithm on a heterogeneous parallel system, which means that we may have processors with different characteristics, so we should divide the work among the slaves according to their performance. Since we do not know the performance of each precessor in advance, we divide the finite elements into small batches that can to be processed independently and simultaneously. Initially, the master process sends one batch to each slave. After processing a batch, each slave returns the result to the master and ask for another one. This process is repeated until all of the finite element contributions are gathered by the master process. After that, the remaining steps of the algorithm are performed by the master as an usual sequential program.

Numerical Analysis of the Physical Parameters

For practical reasons, to analyse the physical parameters of the problem we have implemented the E.B.E.-1 strategy in FORTRAN on a 486DX 33 computer with 16Mb of RAM memory and 0.5 MFLOPS.

We have solved a problem with 8,481 variables and 256 x 32 finite elements in the directions q and z. For this case, starting with d_cg = 10^-5, d_p = 10^-4 and Î = 10^-3, we noticed that the number of iterations of the C.G.M. increases with the ratio of excentricity e/c while the performance of the P.C.G.M. is almost independent of this parameter, suggesting that the diagonal precoditioner is effective for this kind of problem. Some tests were made to determine the influence of the penalty parameter, Î, on the convergence of the conjugate gradient method. We tried some combinations of initial penalty parameters, Î₀, and reduction factors. Although one could expect that a small value of e would increse the condition number of matrix A and, therefore, reduce the rate of convergence of the C.G.M., we noticed that the algorithm was able to find the solution of this problem for a large variety of initial penalty parameters. Besides the rate of convergence was only marginally affected by the choice of Î_0..

This small problem was used to investigate two physical parameters, e and L/D. The remaining parameters were defined according to Duarte Jr. (1989) and are given below.

m = 0.001Pa . s, D = 0.01m, FD = 0.000112m, p_o = 102300Pa and w = 40 x 2ps^-1.

To analyse e a parameter that determines if the cylinders are more or less concentrical, for a fixed L/D = 0.5, we calculate the pressure distribution as function of the following e values: 0.1, 0.3, 0.5 and 0.9. Fig. 3 shows these results for z = 0. As we can see in this figure, the pressure p increases with e. Besides, the geometric shape of its distribution changes considerably, showing a compression and dragging towards the interface of the free boundary, confirming physically the great effort of the fluid to flow through a narrow region.

This fact also shows the numerical difficulty discussed before (the ill-conditioning of the Eq.(11) for some of these high values of e).

To analyse L/D, a parameter that determines if the journal bearing is long or short, for D previously established and a fixed e = 0.5, we calculate the pressure distribution as a function of the following L/D values: 0.5, 1.0 and 1.5. Fig. 4 shows these results for z = 0. As we can see in this figure, the extension of the cavitation region depends on the value of L/D. If we increase L, we obtain a small region, confirming the difficulty for the fluid to reach the cavitation pressure for a large L.

A numerical effect of this fact is the small influence of the boundary conditions in the solution in the case of a large L.

Analysis of the Performance of the Parallel Algorithm

To analyse the performance of the parallel algorithm we solved the same problem with 8,481 variables described in the previous section.

We have implemented the E.B.E.-2 strategy in FORTRAN, using the Parallel Virtual Machine (P.V.M.) as the communication library (see Geist et al.(1993)). The program was run at CENAPAD/SP-UNICAMP, using no more than 5 dedicated RISC/6000 processors, two with 12.7 MFLOPS and 32 Mbytes of RAM memory, and three with 20.72 MFLOPS and 64 Mbytes of RAM memory.

A key feature of our parallel algorithm is the definition of the batch variable, that defines the number of finite elements each process deals with at a time. Although, a good choice of batch can signicantly improve the performance of the algorithm, this parameter is hard to choose, since it depends on the size of the problem and the number of processors available.

To define the batch, some experiments were carried out using the best 3 of the 5 processors available. The results are shown in Table 1.

From Table 1, we noticed that the best value for the batch is 128. With this batch, the time spent using 3 processors (4,448.44 seconds) is significantly smaller than the time of the sequential program (9,273.53 seconds). However, even for the best batch size, about 30% of the total time is spent on master-slave communication.

The performance of the parallel algorithm can be measured by the speedup, here denoted by S(np) defined as

where T(np) is the time spent when np processors are used and T(1) is the time spent by the best sequential program.

The performance can also be measured by the efficiency, here denoted by E(np), defined as

To analyse the performance of our algorithm we carried out same experiments with batch equal to 128 and np varying from 2 to 5. The results are shown in Table 2.

The speedup and the efficiency obtained from these values are shown in Fig. 5. However, it is important to point out that, since the processors are of 2 different kinds, for the computation of the speedup, the time spent by the sequential program was calculated as a weighted average.

As we can see in Fig. 5 , the speedup increases almost linearly with np, while the efficiency is kept about 70% This results show the good scalability of the algorithm and suggest that a higher efficiency could by achieved for a computer with lower message passing overhead. Besides, Fig. 5 also confirms that our parallel approach based on the definition of batchs works well on computers with heterogeneous processors.

Finally, the parallel algorithm was also used to solve a problem with 33,924 variables. Although this problem is four times greater than the first one, the total number of conjugate gradient iterations was increased only by a factor of three, while the time spent per iteration was multiplied by four, so the total parellel processing time was about 12 times greater than the value obtained for the problem whith 8,481 variables. This result also confirms the good scalability of the algorithm and suggest that it can be successfully applied to very large problems.

Conclusions

We have built an algorithm to solve numerically the physical process of cavitation in the hydrodynamic lubrication of journal bearings.

The algorithm is basically composed of two iterative parts. The outer one is associated to a penalty method, and the inner to the conjugated gradient method, used to solve the linear systems generated by the finite element method. In the process of choosing the numerical methods for building this algorithm, we decided to use the element by element strategy, mainly because it can be easily parallelized. Also, as it is necessary due to the ill-conditioning, this strategy could be used to include the diagonal preconditioner. We notice that when the preconditioner is included, the algorithm is not sensible to an increase of the ratio of excentricity.

Two E.B.E. strategies were defined according to the problem dimension. For small problems, we keep the submatrices of all finite elements in memory; otherwise, we recompute the submatrices at each iteration of the P.C.G.M. This last alternative reduces the memory requirements but increases the processing time, and is well tailored for very large problems, so this is the strategy we implemented in parallel.

The E.B.E.-1 strategy was used to analyse the behavior of the physical parameters. In this analysis we could confirm the great effort of the fluid to flow through a narrow region and the difficulty for the fluid to reach the cavitation pressure for a large L.

The E.B.E.-2 strategy was implemented to analyse the scalability of the parallel algorithm. Some computational experiments show that the speedup increase according to number of processors in an almost linear way, so that the efficiency is almost constant.

Article received: June, 1999. Technical Editor: Angela O. Nieckele

• Amendola, M., 1996, "Numerical Solution of one Free Boundary Problem: Cavitation in the Hydrodynamic Lubrication of Journal Bearings"' (In Portuguese), Ph.D. Thesis, State University of Campinas, Campinas, S.P., Brazil, 170 p.
• Barragy, E. and Carey, G.F., 1988, "A Parallel Element - by - Element Solution Scheme'', Int. J. Numer. Meth. Engrg., Vol. 26, pp. 2367-2382.
• Becker, E.B., Carey, G.F. and Oden, J.T., 1981 "Finite Elements: An Introduction"', Prentice-Hall. New York, 258 p.
• Carey, G.F., Barragy, E., Mclay, R. and Sharma, M., 1988, "Element - by - Element Vector and Parallel Computations", Comm. Appl. Numer. Meth., Vol. 4, pp. 299-307.
• Ciarlet, P.G. and Lions, J.L. ed., 1991, "Handbook of Numerical Analysis", New York, North Holland (pt. 1) Vol. 2. Finite Element Methods, 928 p.
• Duarte Junior, D., 1989, "Mathematical Modeling and Simulation of Hydrodynamic Journal Bearings" (In Portuguese), Ph. D. Thesis, Polytechnics School -- USP, S\~{a}o Paulo, SP, Brazil, 179 p.
• Geist et. al., 1993, "P.V.M. User's Guide and Reference Manual". ORNL/TM - 12187.
• Hageman, L.A. and Young, D.M., 1981, "Applied Iterative Methods", Academic Press. New York, 386 p.
• Kinderlehrer, D. and Stampachia, G., 1980, "An Introduction to Variational Inequalities and their Applications", New York, Academic Press, 313 p.
• Moura, C.A. and Amendola, M., 1993, "A Parallel Algorithm for the Numerical Treatment of a Cavitation model in Lubrication". In: Chadan, J.M. and Rasmussen, H. (eds), Free Boundary Problems in Fluid Flow with Applications. New York. Longman, C. pp. 21-23 (Pitman Research Notes in Mathematics Series., Vol 282).
• Orr, F.M. and Scriven, L.E., 1978, "Rimming Flow: Numerical Simulation of Steady, Viscous, Free--Surface Flow with Surface Tension", J. Fluid Mech., Vol 84, part 1, pp. 145-165.
• Rhode, S.M. and Mc Allister, G.D., 1975, "A Variational Formulation for a Class of Free Boundary Problem Arising in Hydrodynamic Lubrication", Internat. J. Engrg. Sci. Vol. 13, pp. 841-850.
• Schlighting, H., 1968, "Boundary - Layer Theory", 7ed. New York, McGraw-Hill, 747 p.

Publication Dates

• Publication in this collection
  19 Aug 2002
• Date of issue
  2001

History

• Received
  June 1999