Review on the Bubble Dynamics Based Cavitation Dynamics for the Negative Squeeze Motion in Lubricated Contacts

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. The approach with bubble dynamics has been introduced to the lubrication community about 20 years ago by Someya, and it is based on the coupling of Reynolds equation and Rayleigh–Plesset equation. It has been used for journal bearings, squeeze film dampers, and 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. More than a decade ago, in 2009, the first paper on the negative squeeze motion with bubble dynamics—allowing numerical calculations of tensile stresses in the lubricant—had been published. The application in mind is the simulation of mixed lubrication for rough surfaces. The negative squeeze motion is then understood as the motion of asperities (on smaller length scales). The paper at hand summarizes some of the research on the dynamics of cavitation in lubricated contacts from different research groups from the last 10–15 years and sketches key topics for further research on the topic. The roadmap is centered around the three key issues that remained from the previous research of the author: (a) numerical stability of the calculations for curved plates, (b) characteristic time scale for separation of plates, and (c) experimental evidence for validating the calculation results.


INTRODUCTION
Cavitation is a typical phenomenon in lubricated contacts and appears in different forms. For journal bearings, film rupture and reformation is typically connected to air ingestion respectively feeding new oil. In squeeze film dampers, oil stiction problems or mixed lubrication on the other hand, pressure induced growth and collapse of bubbles is most relevant.
In most simulations of lubricated contacts (e.g., journal bearings) tensile stresses (negative pressure) are neglected. With respect to pressure p cavitation is then characterized by p = p v or even simpler p = 0 (p v being the vapor pressure at the given temperature).
Nonetheless, negative pressure appears even in journal bearings. To the authors knowledge, Someya (2003) was the first to do numerical calculations for journal bearings that included negative pressure 1 . He used a bubble dynamics approach including a dilatational viscosity. For discretization, he used the Finite Difference method. Someya reports a good agreement between experiment and calculations based on bubble dynamics.
Tensile stresses cannot be neglected when simulating the negative squeeze motion of lubricated contacts (i.e., the separation of two plates, see Figure 1). At first sight, the negative squeeze motion of a macroscopic plate does not seem relevant for practical purposes. Far from it! An important practical application is the oil stiction, which is highly relevant for the performance of reed valves in compressors or of switching valves in hydraulic applications (e.g., Resch and Scheidl, 2014;Yoshizumi et al., 2018). Switching times as low as 1 ms require sophisticated models for the lubricant's behavior in valves.
Certainly the negative squeeze motion is also relevant when looking at smaller length scales. Then we look at asperities moving relative to each other and thereby again and again separating from each other. Tensile stresses on the microscopic level resulting from negative squeeze motion can contribute to fatigue and wear in mixed lubrication. Indeed, the simulation of rough surfaces including negative pressure was the application in mind when the author started his research work on lubricated contacts.
There is experimental evidence for tensile stresses in lubrication films for the setup under study for more than five decades. But it was only in 2009 that calculations with tensile stresses have been reported for the negative squeeze motion Geike and Popov (2009a). To reproduce tensile stresses in numerical simulations bubble dynamics needs to be taken into account. Traditional methods based on Reynolds equation and static cavitation conditions do not yield tensile stresses. Bubble dynamics is also required for getting the relevant time scale for the separation of plates.
The key question is which system of equations shall constitute the numerical simulation model. The answer was and still is the coupling of Reynolds equation for compressible fluid flow and Rayleigh-Plesset equation for bubble dynamics.
In order to identify the key topics for further research, it is worthwhile to look at the three questions that remained when the author's research came to a temporary end in early 2008.
1. How can the numerical stability problems be solved for the case of spherical plates (finite radius of curvature)? 2. What is the correct time scale for the plate separation? 3. What experimental evidence is available for a quantitative validation of numerical calculations?
The paper at hand looks at relevant contributions to the topic, particularly in the last 12 years 2 , and suggests potential corner stones for future research on the topic. The papers is organized along the three questions from above.

The Authors Approach to the Cavitation Dynamics for the Negative Squeeze Motion
To the author's knowledge, the papers from January 2009 Geike and Popov (2009a,b) are the first to report on numerical calculations of the negative squeeze motion based on Reynolds equation and bubble dynamics, thereby allowing tensile stresses to develop in the lubricant. The author used Reynolds equation and Rayleigh-Plesset-equation as the basis of the simulation model. p v is the vapor pressure, γ is the surface tension constant, η liq is the viscosity, and ρ liq the density of the lubricant (liquid). Reynolds equation is the partial differential equation for the pressure in the lubricant and can be obtained from Navier-Stokes equation and continuity equation on assumptions such as negligible inertia, Newtonian fluid of constant viscosity and thin film geometry. Rayleigh-Plesset equation describes the growth in radius of a spherical bubble in an incompressible fluid of infinite extent. It is used to approximately calculate the vapor fraction and thus the density, which is needed in the Reynolds equation. It is important to understand that the simulation does not consider single bubbles in the lubricant. Instead, the simulation is based on a characteristic bubble radius at each location r that is connected to the vapor fraction. Rayleigh-Plesset equation is used to derive a partial differential equation for vapor fraction, which is then used in the simulation model.
For the published results the first term with the second derivative of R had been neglected. Further simulations with the full Equation (2) yielded almost the same results.
The partial differential Equation (1) had been transformed into a system of ordinary differential equations based on the Differential Quadrature method (Shu, 2000). The final system of equations is a differential-algebraic system of differentiation index 1. One of the challenges in solving the problem numerically is the strong transient behavior. Within a short time interval significant negative pressure appears and disappears again, while vapor fraction increases quickly.
The numerical solution converged for flat plates (infinite radius of curvature) for a variation of grid types. The results are-at least qualitatively-correct. Husen et al. (2016) came in their paper from 2015 to the same results using a Finite Element method, thereby validating our discretization method and computer implementation of the numerical calculations.
For spherical plates (finite radius of curvature) numerical simulations did not converge for no obvious reason. The numerical integration of the initial value problem (that results from the spatial discretization) fails to converge. However, for the mixed lubrication of rough surfaces it is key to master the spherical plate problem.
From first-hand experience, the author knows that in the Differential Quadrature method stability of numerical calculations sometimes depend on the grid type used (e.g., uniform grid vs. Chebyshev-Gauss-Lobatto grid), but the stability problems remained no matter which grid had been used. The stability issue's root cause and appropriate countermeasures have not been identified yet.

Calculations for Squeeze Film Dampers With Dilatational Viscosity Term
Few months later, in April 2009, a paper by Gehannin et al. (2009) was published on a quite similar topic-cavitation modeling based on the coupling of Reynolds equation and Rayleigh-Plesset equation. The setup and the application in mind is however quite different: they study squeeze film dampers, i.e., they look at a tangential motion instead of the vertical motion (separation). Bubbles contain vapor and undissolved gas thus the pressure inside the bubble is the sum of the partial pressure of the gas and the vapor pressure. Consequently, they write p B instead of p v in equation (2), with In addition, following Someya (2003), Gehannin et al. used an additional term on the right-hand side of the Rayleigh-Plesset Equation (2) with the dilatational viscosity, They finally solve the system of equations using a Finite Volume approach and conclude that for the squeeze film damper the bubble pressure term and the dilatational viscosity term are most relevant. The agreement between their calculations and the experimental results from Adiletta and Pietra (2006) will be discussed below.
Gehannin et al. do not report on stability issues in their numerical approach. This might be caused by one of the three major differences.
• The setup is different (tangential instead of vertical motion).
• Additional terms have been used in the Rayleigh-Plesset equation. • A different method for discretization has been used (Finite volume instead of Differential Quadrature).
In future research, the relevance of the additional terms for the negative squeeze motion may be studied first. A later paper from the same group (Braun et al., 2017) reports on a simulation model for journal bearings that is extended by an energy equation to model temperature effects and heat transfer. Jaramillo et al. (2018) study the coupled equations from the mathematicians perspective, i.e., they study the existence and stability of stationary solutions. They study the case with zero vertical relative velocity (e.g., the standard procedure for journal bearing simulation) and take the dilatational viscosity into account. For journal bearings they perform numerical experiments based on the Finite Volume method and a backward Euler scheme. They conclude that for the eccentricity ǫ < 0.41 the transient solution converges toward the stationary solution; for higher values of ǫ time-convergence is no longer obtained. One of the hypothesis for stability identified by Jaramillo et al. is indeed violated for ǫ > 0.41. This hypothesis is centered around the term

Mathematical Studies on the Well-Posedness of Reynolds-Rayleigh-Plesset Coupling for Journal Bearings
Even though the setup under study is different, this result is understood as an additional hint to look at the additional terms in the Rayleigh-Plesset equation when doing further research on the the negative squeeze motion.

The Way Forward: Extend the Model Equations
For further numerical simulations, it seems most promising to include the dilatational viscosity term and extend the pressure term to allow for other pressure than the vapor pressure inside the bubble. The stability issue can be revisited based on the extended model. In case that stability issues remain, the method of discretization can be changed, from Differential Quadrature method to Finite Volume method or Finite Element method. In addition to numerical experiments it would be beneficial to partner up with mathematicians to study the negative squeeze motion from the mathematical perspective. This could result in statements on the parameter values for which stable numerical calculations can be expected.

The Authors Results Based on Bubble Dynamics
The numerical results published in 2009 are qualitatively reasonable. The pressure distribution starts at the onset of motion with the distribution that is expected without cavitation. Thus, significant tensile stresses can be observed for a very short time. Tensile stresses quickly drop while vapor fraction increases. However, the characteristic time scale for plate separation seems rather short and needs attention.

Calculations of the Oil Stiction Force Without Bubble Dynamics
As mentioned above, an important practical application for the plate separation is the oil stiction in technical devices as valves in hydraulic applications or reed valves in compressors. Oil stiction influences the opening time and therefore the performance of valves in hydraulic or pneumatic applications. Roemer et al. (2015) state that for small initial distances the assumption of negligible tensile strength of the liquid is not applicable. They propose a model with a fluid tensile strength as additional material property to deal with the tensile stresses that develop when two plates quickly separate. Also, a distinct cavitation zone is modeled. The cavitation zone appears in the event of the tension exceeding the threshold. Roemer et al. report of one example where a time of about 12 ms is required before a cavitation zone is formed. Scheidl and Gradl (2016) study the problem of two separating plates based on Reynolds equation and cavitation. Their calculation is based on static cavitation conditions, i.e., is without bubble dynamics, and disregards tensile stresses 3 . Another paper from the research group (Resch and Scheidl, 2014) reports tensile stresses for short periods of time-in the range of tenth of milliseconds.

Calculations With Bubble Dynamics and Experiments for the Oil Stiction of Reed Valves
The paper by Yoshizumi et al. (2018) from 2018 reports on experiments and numerical calculations for the opening of a reed valve. As reed and valve seat are flat and the motion is vertical this situation is relatively close to the negative squeeze motion. However, reed and valve seat are parallel only at the beginning.
As the reed bends while opening, the motion of the two parts is not parallel.
The authors include the elastic deformation of the reed as well as the dynamics of the oil film. For the latter, they use a coupling of Reynolds equation and Rayleigh-Plesset equation. Again, the dilatational viscosity is included and the inertia terms are neglected. The Finite Volume method is used for discretization. For comparison, not only the dynamic cavitation model with bubble dynamics is used for calculation but also two other models (one with static cavitation, one without cavitation). They do not report any stability issues in their calculations.
For the reed valve delay time (a very relevant performance indicator) Yoshizumi et al. got similar results from experiments and from numerical calculations with the dynamic cavitation model 4 . As one would expect, the delay time is far too high for the model without cavitation and far too small for the static model. This is a wonderful hint that the dynamic cavitation model with bubble dynamics is just the right choice-it yields correct results, and simpler choices for modeling are insufficient for a full understanding.

The Way Forward: Extend the Model Equations
The results for reed valves from Yoshizumi et al. indicate a dependency between delay time and the cavitation model used for calculations (see above). Even more important and not as obvious: the delay time significantly depends on the dilatational viscosity. In particular, the lower the value of the dilatational viscosity the lower the delay time. Therefore, a model without dilatational viscosity might be insufficient for a correct modeling of the time scale.
Using the above mentioned extensions to the simulation model, first of all including the dilatational viscosity, the time scale needs to be looked at closely again. Based on the results from Yoshizumi et al. it is quite save to assume an increase in the characteristic time for the negative squeeze motion.
Having the simulation of mixed lubrication of rough surfaces in mind, Roemer et al.'s (2015) approach with a tensile stress should also be considered. The advantage of this approach is the lower computational effort. Geike and Popov (2009a) give an overview on experimental evidence for tensile stresses in the negative squeeze film motion and discuss why the published data are not sufficient for model validation. In particular, they point to the work of Hays and Feiten (1964), Parkins and May-Miller (1984), Chen et al. (2004), and Wang et al. (2005), who had studied the timedependent cavitation experimentally in a simple parallel-plate squeeze film configuration.

GETTING FURTHER EXPERIMENTAL EVIDENCE
It seems that no experimental results for the negative squeeze motion of spherical plates-sufficient for a quantitative validation of simulation models-have been published yet. The work of Sun et al. (2008) seems promising for the model validation for flat plates in oscillatory motion. The experiments from Yoshizumi et al. (2018) are also relevant as their setup is somewhat close to the negative squeeze motion and a characteristic time scale is studied. It is also worthwhile to study other contributions on the oil stiction topic to find further experimental data.
Most experimental studies, however, focus on other setups 5 . The focus is mostly on journal bearings or squeeze film dampers. As an example, the paper by Adiletta and Pietra (2006) from 2006 studies the squeeze film damper. Their results show that negative values for the absolute pressure appear in squeeze film dampers too. There seems to be still a difference between numerical results from Gehannin et al. (2009) and the experimental results from Adiletta and Della Pietra with respect to the absolute pressure inside the cavitation zone. Gehannin et al. conclude that it is essential to include the dilatational viscosity to get negative pressure in the studied squeeze film damper. They also emphasize that a model based on Reynolds and Rayleigh-Plesset equations relies on many parameters, temperature dependent lubricant properties and also initial conditions on bubble size and vapor fraction.
For the future, a joint research between partners who on the one side undertake the necessary experimental studies and on the other side work on the numerical simulation model seems beneficial to answer the open questions.

WHAT IS NEXT?
In the author's opinion the research on the cavitation dynamics for mixed lubrication should be centered around the above mentioned key topics-stability of numerical calculation, clarifying time scale and finding further experimental evidence.
As mentioned above, it seems most promising to set the focus on extending the Rayleigh-Plesset equation (bubble dynamics). Hopefully, the questions around the stability issues and the characteristic time scale can then be answered. Involving partners for (i) the mathematical study on the stability issue and (ii) for further experiments would be beneficial for future research.
From todays perspective, the next leap forward would be the simulation of the elasto-hydrodynamic problem of rough surfaces including the cavitation dynamics. For this the boundary element method (BEM) seems 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 methods that require the discretization of the entire volume. The tool and experience of Popov's research group Pohrt and Li (2014) could be used here.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.