FUME 7343 FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 19, No 1, 2021, pp. 115 - 124 https://doi.org/10.22190/FUME210112027G © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper BUBBLE DYNAMICS-BASED MODELING OF THE CAVITATION DYNAMICS IN LUBRICATED CONTACTS Thomas Geike Beuth Hochschule für Technik Berlin, Berlin, Germany Abstract. Cavitation is a common phenomenon in fluid machinery and lubricated contacts. In lubricated contacts, there is a presumption that the short-term tensile stresses at the onset of bubble formation have an influence on material wear. To investigate the duration and magnitude of tensile stresses in lubricating films using numerical simulation, a suitable simulation model must be developed. The chosen simulation approach with bubble dynamics is based on the coupling of the Reynolds equation and Rayleigh-Plesset equation (introduced about 20 years ago by Someya).Following the basic approach from the author’s earlier papers on the negative squeeze motion with bubble dynamics for the simulation of mixed lubrication of rough surfaces, the paper at hand shows modifications to the Rayleigh-Plesset equation that are required to get the time scale for the dynamic processes right. This additional term is called the dilatational viscosity term, and it significantly influences the behavior of the numerical model. Key Words: Cavitation, Mixed Lubrication, Oil Stiction, Negative Squeeze Motion, Bubble Dynamics, Negative Pressure, Wear 1. INTRODUCTION Cavitation is a common phenomenon in fluid machinery and lubricated contacts. In the context of fluid machinery, the focus is on cavitation erosion, in which pressure shocks lead to material failure [1]. In lubricated contacts, on the contrary, there is a presumption that the short-term tensile stresses at the onset of bubble formation have an influence on material wear. To investigate the duration and magnitude of tensile stresses in lubricating films using numerical simulation, a suitable simulation model must be developed. The author proposed such a cavitation model for numerical simulations of mixed lubrication that includes tensile stresses in the lubrication film at the onset of bubble growth [2]. Received January 12, 2021 / Accepted March 07, 2021 Corresponding author: Thomas Geike Beuth Hochschule für Technik Berlin, Luxemburger Str. 10, 13353 Berlin, Germany E-mail: Thomas.Geike@beuth-hochschule.de 116 T. GEIKE Generally speaking, simulation models for the cavitation dynamics in lubricated contacts can be roughly clustered into two groups: either without or with bubble dynamics, the first one being the standard case for most fluid film bearing calculations. Cavitation is a crucial phenomenon in the calculation of the load carrying capacity of plain bearings. The literature available on this subject is extensive. For most bearing calculations, the tensile stresses present for a short time are not relevant; in this respect, “pressure equals zero” is used for calculations when there is no positive pressure. The presented simulation model belongs to the second group – having the bubble dynamics as one key feature included in the model. The aim here is to calculate, besides macroscopic forces between moving rough surfaces, also the tensile stresses that may affect the wear in mixed lubrication applications. Also, for adhesion problems in lubricated contacts (the so-called oil stiction problem), the tensile stresses in the lubricating film are necessarily taken into account when two surfaces separate. The original model of the author yielded qualitatively reasonable results, but three key issues remained [3]. 1. The characteristic time for the decline of tensile stress seems rather short. There is no reliable source of information on the characteristic time scale for the processes studied but results from the oil stiction problem and experiments hint to larger characteristic times than those observed with the original model. This will be discussed in detail in Sections 2.2 and 4.3. 2. Numerical simulations work well with circular plates, but spherical plates cause numerical stability problems in the original model. 3. Experimental evidence with quantitative results for time evolution of pressure on very short time scales are missing for the squeeze motion studied here. The paper at hand deals with the first two open issues from this list. 2. CAVITATION IN LUBRICATION FILMS – SOME BASICS 2.1. Cavitation means formation and destruction of cavities In lubricated contacts, conditions can occur in which the lubricant becomes temporarily discontinuous. In the lubricating film, a coexistence of liquid and cavities is created and perishes. The phenomenon called cavitation can be divided into gas and vapor cavitation [4,5]. Cavitation occurs when the pressure locally falls below a certain limit. Gas cavitation can be observed in aerated lubricants at pressures below the ambient pressure. The cavities then contain dissolved gases. Vapor cavitation occurs at dynamic load conditions when the pressure falls below the vapor pressure. Vapor of the lubricating fluid is then present in the cavities. Vapor cavitation is a dynamic process; vaporization does not occur abruptly. Therefore, lubricants can also transmit tensile stresses for a short time. 2.2. Tensile stresses are observed in negative squeeze motion The classical solution for normal force F and pressure p when squeezing a viscous fluid (dynamic viscosity ) between two parallel circular plates (Fig. 1, velocity V, plate radius L, distance h) is Bubble Dynamics-Based Modeling of the Cavitation Dynamics in Lubricated Contacts 117 𝐹 = 3𝜂𝜋𝑉𝐿4 2ℎ3 (1) and 𝑝 = 𝑝amb − 3𝜂𝑉 ℎ3 (𝐿2 − 𝑟2) (2) with ambient pressure 𝑝amb [6] and radial coordinate r. When the upper plate approaches the base plate, the pressure is always positive (squeeze motion). When the plates are separated (negative squeeze motion), the pressure at the center of the plate falls below the vapor pressure 𝑝v for 𝑉 > ℎ3 3𝜂𝐿2 (𝑝amb − 𝑝v) (3) and cavitation starts. Fig. 1 Negative squeeze motion for the setup with circular plate (pull-off experiment) Time-dependent cavitation in a simple arrangement of parallel plates has been experimentally studied by Hays and Feiten [7], Parkins and May-Miller [8], and Chen et al. [9]. Hays and Feiten considered the case of constant velocity, while the other two groups studied the case of periodic motion. The experiments confirm that the lubricating film can transmit a tensile stress. Due to the tensile stress, bubble growth occurs and finally the macroscopically visible disruption of continuity. The tensile stress approaches zero within a short time. In the experiments, the cavitation region arises in the center and disappears there. All experiments focus on the characterization of cavitation patterns; in particular, no data are provided that can be used to check numerical investigations sufficiently well. The experiments under practical conditions have shown that cavitation can occur in bearings, e.g. in connecting rod bearings in internal combustion engines, and that this is sometimes associated with considerable damage [10]. 118 T. GEIKE 2.3. Tensile stresses are reproduced in simulation models with bubble dynamics An overview of research results on the simulation of cavitation dynamics can be found in the publications [2,3]. There is one common feature of all simulations known to the author: they do not reproduce the fine spatial structures as observed in experiments. As mentioned before, simulation models for the cavitation dynamics in lubricated contacts can be roughly clustered into two groups: either without or with bubble dynamics. Approaches without bubble dynamics usually use the Reynolds equation complemented with some additional conditions to define the cavitated areas, and do not reproduce tensile stresses in the fluid. The approach with bubble dynamics, on the other hand, reproduces tensile stresses in the fluid film in negative squeeze motion. The building blocks of bubble dynamics-based simulation models are the Reynolds equation and Rayleigh-Plesset equation. The idea of combining the two equations for lubrication problems goes back to Someya about 20 years ago [11]. It has been used for journal bearings and squeeze film dampers; it is essentially required for correct numerical calculations of the negative squeeze motion (i.e. the separation of two plates) or the oil stiction problem. 3. MODELING THE CAVITATION DYNAMICS IN THE SQUEEZE MOTION 3.1. The idea in a nutshell: Reynolds equation and bubble growth As said before, a simulation model that reproduces tensile stresses must account for the dynamics of bubble growth and decay. Starting from the Rayleigh-Plesset equation [12,13] for the bubble radius R 𝑅 𝑑2𝑅 𝑑𝑡2 + 3 2 ( 𝑑𝑅 𝑑𝑡 ) 2 = 𝑝v − 𝑝 𝜌liq − 4𝜂liq 𝜌liq𝑅 𝑑𝑅 𝑑𝑡 − 2𝛾 𝜌liq𝑅 (4) and the Reynolds equation for compressible fluids 1 𝑟 𝜕 𝜕𝑟 ( 𝑟𝜌ℎ3 𝜂 𝜕𝑝 𝜕𝑟 ) = 12ℎ 𝜕𝜌 𝜕𝑡 + 12𝜌𝑉 (5) partial differential equations can be derived for density  and pressure p in the pull-off experiment. Both quantities depend on time t and radial coordinate r. The subscript liq indicates the liquid state.  characterizes the surface tension. This approach for the cavitation problem in lubricating films in the pull-off experiment was first presented in [2]. The work of Someya [11] was not known to the author at that time. Instead, the approach was inspired by the work of Sauer on cavitation in fluid machinery [14]. The spatial discretization of the partial differential equations with the differential quadrature method [15] leads to a system of differential and algebraic equations solved with a solver for differential-algebraic systems of equations (DAE). As shown by numerical experiments, a simulation model based on the above equations yields extremely short times for the pressure drop [2] and numerical stability problems for the spherical cap (unpublished research results from the author). It will be shown that adding an additional term (dilatational viscosity term) solves both problems. The stability problem will be discussed in a subsequent publication, dealing with the spherical cap in more detail. The dilatational viscosity term has been used by different authors in the context Bubble Dynamics-Based Modeling of the Cavitation Dynamics in Lubricated Contacts 119 of cavitation in lubricants, beginning with Someya [11]. According to Someya, the dilatational viscosity term is related to the Marangoni effect. When the bubble expands under negative pressure, the dilatational viscosity increases the resistance against bubble growth and rupture. In the paper by Gehannin et al. [16], two terms change compared to the above-mentioned Rayleigh-Plesset equation. First, an additional term is added with the dilatational viscosity (following Someya’s approach), and second, the constant vapor pressure is replaced by a more complex expression. With respect to the dilatational viscosity, Gehannin et al. follow the explanations of Someya [11] and add the term − 4𝜅 𝜌liq𝑅 2 𝑑𝑅 𝑑𝑡 to the right-side of the Rayleigh-Plesset equation, with  characterizing the dilatational viscosity term. Adding this term to the simulation model has a significant impact on the time scale on which tensile stresses exist during the negative squeeze motion. Numerical experiments with a single bubble and the additional term, recently undertaken by the author, give three main insights: (1) the characteristic time scale for the bubble dynamics is about 1000 times larger than in earlier numerical experiments. (2) The inertia term with the second derivative of the bubble radius and the surface tension term can be neglected. (3) When solving the Rayleigh-Plesset equation for the first derivative of the bubble radius, catastrophic cancellation may occur if not dealt with appropriately. 3.2. Continuum with microstructure - coupling bubble dynamics with the Reynolds equation The Rayleigh-Plesset equation describes the behavior of a single, spherical bubble in an incompressible fluid of infinite extension. It applies to the growth of the bubble in the first phase, where mechanical effects dominate: inertia, pressure difference, viscosity, and surface tension. Modeling cavitation by bubble growth assumes that the bubbles in the lubricating film are far enough apart for the Rayleigh-Plesset equation to be applicable in good approximation. Individual bubbles, characterized by position and radius, are not described. Instead, a density  or a vapor fraction  is assigned to a location, more precisely to a radial coordinate r, in the framework of a continuous description method. Between density and vapor fraction exists the relation 𝜌 = 𝛼𝜌vap + (1 − 𝛼)𝜌liq (6) where index vap indicates to vapor state. Correspondingly, for the time derivative 𝜕𝜌 𝜕𝑡 = 𝜕𝛼 𝜕𝑡 (𝜌vap − 𝜌liq) (7) Formally, via 𝑅 = √ 3 4𝜋𝑛0 𝛼 1 − 𝛼 3 (8) each location is assigned a bubble radius, with 𝑛0 characterizing the concentration of bubble nuclei. From the Rayleigh-Plesset equation follows the time evolution of the bubble radius and hence the time evolution of the density. 120 T. GEIKE 4. NEGATIVE SQUEEZE MOTION FOR THE CIRCULAR PLATE Now the pull-off experiment (negative squeeze motion) shown in Fig. 1 will be studied in detail. A circular plate is pulled upward with constant velocity in the simplest case. At the boundary, the pressure is always the ambient pressure. This simple arrangement is studied by Boedo and Booker [17] using the finite element method and a very simple cavitation model. The arrangement has practical significance for the study of the adhesion behavior of lubricated contacts (oil stiction problem), as it is relevant in the simulation of rapidly opening valves. 4.1. Dimensionless model equations: Coupling of the Rayleigh-Plesset equation and Reynolds equation The equations constituting the simulation model are made dimensionless as follows. All quantities with a bar on top are corresponding dimensionless quantities. For the lengths, the reference quantity is the plate radius. There is one exception: for the bubble radius is 𝑅 = 𝛬 𝑅 𝐿 (9) although in the practical implementation =1 has always been used so far. For pressure, the ambient pressure is taken as reference value. For the reference time 𝑇ref = 𝜂liq 𝑝amb (10) is set. For the dimensionless bubble radius and the dimensionless concentration of bubble nuclei, we have 𝑅 = √ 3 4𝜋𝑛0   𝛼 1 − 𝛼 3 (11) and 𝑛0  =   𝑛0 𝐿  3 𝛬3 (12) The dimensionless velocity is �̅� = 𝑉 𝑇ref 𝐿 (13) The differential equation for the bubble radius, neglecting the term with the second derivative of the bubble radius and the surface tension term, can be reduced to the form 𝑅′ = − 4 3𝑅 (𝜁1 + 𝜁2 𝑅 ) + √ 16 9𝑅 2 (𝜁1 + 𝜁2 𝑅 ) 2 2 3 𝜁1 + (𝑝v − 𝑝) (14) Eq. (14) is the dimensionless form of Eq. (4) considering all mentioned simplifications to the Rayleigh-Plesset equation. Reynolds equation (5) is not repeated here, as its form does not change by making the equation dimensionless. The dimensionless constants are defined as 𝜁1 = 𝜂liq 2 𝛬2 𝑝amb𝐿 2𝜌liq (15) and Bubble Dynamics-Based Modeling of the Cavitation Dynamics in Lubricated Contacts 121 𝜁2 = 𝜅𝜂liq𝛬 3 𝑝amb𝐿 3𝜌liq (16) For practical implementation, again note that for small bubble radius, the difference is formed from two numbers of nearly equal size. To avoid the catastrophic cancellation effect, the program code needs to be written appropriately. For the time derivative of the vapor volume fraction the relation 𝜕𝛼 𝜕𝑡 = 3𝛼(1 − 𝛼) 1 𝑅 𝑑𝑅 𝑑𝑡 + ℎ2 12𝜂 𝜕𝑝 𝜕𝑟 𝜕𝛼 𝜕𝑟 (17) is used according to Eq. (12) from [2]. This gives the change in density over time needed for the Reynolds equation. 4.2. Dynamics of tensile stresses are plausibly reproduced Selected simulation results for the circular plate are discussed below. For the selected dynamic viscosity, 𝑇ref = 0.35 μs. Initial distance is ℎ = 10 −2at 𝑡 = 0. The lubricant film thickness is thus 1% of the plate radius. The velocity is constant 𝑉 = 10−4. For the start value of vapor fractionα0 two different values, 10 -8 and 10-3 are chosen. In diagrams with the spatial coordinate as horizontal axis (i.e. Fig. 2 and Fig. 4), the curves show different time moments between 0 to 200 as indicated in the diagrams. For 𝑡 = 0 the system starts with the theoretical pressure distribution without cavitation (Fig. 2 for α0 = 10 −8). The pressure distribution is not given as an initial condition but results from the given α distribution as part of the solution of the DAE system. With increasing time, the tensile stress decreases more and more. The massive decrease of the tensile stress happens in a time of the order of 0.1 millisecond. In the context of this work, the characteristic time for pressure drop is defined as the time required for the decay of force to 1% of the initial value. The initial value α0 = 10 −3 gives a similar pressure distribution (quite contrary to previously published calculations without dilatational viscosity). Fig. 2 Pressure distribution for the negative squeeze motion 122 T. GEIKE For 𝑡 = 400, the dimensionless pressure in the center is about -1.4. The dimensionless force drops from -463 at 𝑡 = 0 to -3.7 for 𝑡 = 400 (0.14ms); the drop in force to 1% of the initial value occurs within 0.13ms (Fig. 3).Thus the characteristic time for pressure drop is here 0.13 ms. Fig. 4 shows the corresponding time evolution of vapor fraction and bubble radius as function of radial coordinate and time. Vapor ratio, bubble radius and thus density do not change significantly during the phase of rapid decrease of tensile stress. Essentially, there are two phases. In the first phase the tensile stress drops rapidly, then in the second phase the density changes appreciably, while the pressure is nearly zero. Fig. 3 Macroscopic force as function of time for the negative squeeze motion Fig. 4 Spatial distribution of vapor fraction (left) and bubble radius (right) Bubble Dynamics-Based Modeling of the Cavitation Dynamics in Lubricated Contacts 123 4.3. Additional term in the Rayleigh-Plesset equation leads to significantly larger characteristic times In the earlier simulations, the pressure drop occurs in a time interval of order 𝑡 ≈ 10−1 or 𝑡 ≈ 0.035  μs [2]. In the new simulations considering dilatational viscosity, a comparable pressure drop occurs about 1000 times slower! The now observed characteristic times of the tensile stress drop in the order of 0.1ms matches the times reported by Resch and Scheidl for the oil stiction problem [18]. A characteristic time of the order of 0.1 ms is also more plausible when thinking of the experimental observations [7-9]. Thus, including dilatational viscosity in the simulation model is key. 5. CONCLUSIONS An existing model for the numerical simulation of cavitation in mixed lubrication contacts that reproduces tensile stresses on short time scales has been developed further. As assumed earlier adding the dilatational viscosity term to the Rayleigh-Plesset equation significantly increases the characteristic time of the existence of tensile stresses in the negative squeeze motion. The larger time scale seems more realistic when looking at the related problem of oil stiction and early experimental evidence. Furthermore, initial simulations for the spherical cap indicate that the previously stated numerical stability problems can be overcome by the additional term with the dilatational viscosity. Further work needs to be done to study the spherical cap in negative squeeze motion in detail, including a comparison of density evolution on larger time scales. From today’s perspective, the two major directions of further research are (1) the simulation of the elasto-hydrodynamic problem of rough surfaces including the cavitation dynamics and (2) studies on the effect of short-time tensile stresses in lubricants on the material wear. For the first topic, the boundary element method (BEM) seems to be the right choice for modeling the elastic part. Only the surface is discretized in the BEM. Hence the method allows to model surface roughness with very fine meshes and it is often more efficient for contact problems than the methods requiring the discretization of the entire volume. To make BEM simulations with cavitation dynamics efficient, it is also advisable to undertake further studies of the negative squeeze motion regarding the relevance of terms (inertia, viscosity, surface tension). REFERENCES 1. Prandtl, L., 1957, Strömungslehre, Friedr. Vieweg & Sohn. 2. Geike, T., Popov, V.L., 2009, A bubble dynamics based approach to the simulation of cavitation in lubricated contacts, ASME J. Tribol., 131(1), 011704. 3. Geike, T., 2020, Review on the bubble dynamics based cavitation dynamics for the negative squeeze motion in lubricated contacts, Front. Mech. Eng., 6, 33. 4. Dowson, D., Taylor, C.M., 1979, Cavitation in bearings, Annu. Rev. Fluid Mech., 11, pp. 35-66. 5. Szeri, A.Z., 2005, Fluid Film Lubrication, Cambridge University Press. 6. Popov, V.L., 2015, Kontaktmechanik und Reibung, Springer Vieweg. 7. Hays, D.F., Feiten, J.B., 1964, Cavities between moving parallel plates, in: Davies, R. (Ed.), Cavitation in Real Liquids, Elsevier, New York, pp. 122-127. 124 T. GEIKE 8. Parkins, D.W., May-Miller, R., 1984, Cavitation in an oscillatory oil squeeze film, ASME J. Tribol., 106, pp. 360-367. 9. Chen, X., Sun, M., Wang, W., Sun, D.C., Zhang, Z., Wang, X., 2004, Experimental investigation of time dependent cavitation in an oscillatory squeeze film, Sci. China Ser. G, 47, pp. 107-112. 10. Optasanu, V., Bonneau, D., 2000, Finite element mass-conserving cavitation algorithm in pure squeeze motion - validation/application to a connecting-rod small end bearing, ASME J. Tribol., 122, pp. 162-169. 11. Someya, T., 2003, Negative pressure in the oil-film of journal bearing, Rotrib 03 National Tribology Conference, University of Galati, Galati, pp. 215-220. 12. Plesset, M.S., Prosperetti, A., 1977, Bubble dynamics and cavitation, Annu. Rev. Fluid Mech., 9, pp. 145-185. 13. Feng, Z.C., Leal, L.G., 1997, Nonlinear bubble dynamics, Annu. Rev. Fluid Mech., 29, pp. 201-243. 14. Sauer, J., 2000, Instationär kavitierende Strömungen - ein neues Modell basierend auf Front Capturing und Blasendynamik, PhD Thesis, Universität Karlsruhe, Karlsruhe, Germany. 15. Shu, C., 2000, Differential Quadrature and Its Applications in Engineering, Springer, Berlin. 16. Gehannin, J., Arghir, M., Bonneau, O., 2009, Evaluation of Rayleigh-Plesset equation based cavitation models for squeeze film dampers, ASME J. Tribol., 131, 024501. 17. Boedo, S., Booker, J.F., 1995, Cavitation in normal separation of square and circular plates, ASME J. Tribol., 117, pp. 403-410. 18. Resch, M., Scheidl, R., 2014,A model for fluid stiction of quickly separating circular plates, Proc. IMechE C J. Mech. Eng. Sci., 228, pp. 1540–1556.