An arbitrary Lagrangian–Eulerian formulation for the Reynolds equation considering the JFO boundary condition

In the realm of lubrication theory, the hydrodynamics of the fluid lubricant are commonly elucidated through the Reynolds equation. Its establishment stems from using the Eulerian description method, restricting the fluid domain to a fixed spatial configuration. However, certain scenarios necessitate coupling the Reynolds equation with solid deformation, wherein the fluid boundary undergoes alterations concurrent with the displacement or deformation of the solid. In such cases, the employment of a movable fluid domain becomes imperative. To address this requirement, this contribution introduces the arbitrary Lagrangian–Eulerian method, which effectively transforms the reference system of the Reynolds equation into a fluid domain amenable to arbitrary movement. Subsequently, the weak formulation associated with the movable fluid domain is derived and synergistically combined with a preexisting cavitation formulation. This integration facilitates an efficient solution of the Reynolds equation utilizing the finite element method while duly accounting for the mass‐conserving cavitation effects. Ultimately, a transient hydrodynamic lubrication example is presented to demonstrate the feasibility of the newly developed formulations.


INTRODUCTION
The investigation of lubricated contacts typically entails the consideration of two solid surfaces that may come into contact, with the interstitial region filled by a fluid lubricant.The lubricant flow within the narrow gap between the solid surfaces can be characterized by the classical Reynolds equation [1].It is imperative to highlight that the phenomenon known as cavitation may arise under specific circumstances.When the fluid pressure drops below a critical threshold, frequently referred to as the cavitation pressure, the fluid can release dissolved gas bubbles or draw in ambient air into the interstice.This leads to a significant alteration in the pressure distribution across the fluid film.To account for this, diverse boundary conditions have been introduced in the lubrication theory, including the Gümbel, Swift-Stieber, and Jakobsson-Floberg-Olsson (JFO) boundary conditions.Among these, the JFO boundary condition aligns most closely with physical reality as it considers the conservation of mass for the mixture flow containing bubbles.On the other hand, it also poses greater numerical challenges in implementation regarding its complementary constraints.Further comprehensive discussions on the Reynolds equation and cavitation boundary conditions can be found in-depth within [1].Recently, the authors of this paper have developed a cavitation formulation [2] based on the augmented-Lagrangian (AL) method for the Reynolds equation with the JFO boundary condition.This formulation transforms the original inequalityconstrained problem into an unconstrained one, which can be efficiently solved using the Newton method.Nevertheless, this formulation employs the Eulerian description method, where the reference system consists of fixed spatial points.This approach remains adequate when the motion of the solid surfaces can be disregarded, as typically observed in the stationary hydrodynamic lubrication regime.In contrast, in the transient elastohydrodynamic lubrication (EHL) regime, the motion or substantial deformation of the solid surface profoundly impacts the fluid's geometric boundaries.To capture the evolution of physical quantities influenced by the changing boundaries, a movable fluid domain proves to be more effective.
The methodology employed to transfer the reference system to a moving fluid domain is referred to as the arbitrary Lagrangian-Eulerian (ALE) method.Initially proposed within the framework of finite difference and finite volume methods, pertinent contributions can be found in [3] and [4].Subsequently, the ALE method has been extended to encompass the field of finite element analysis [5,6], and an exhaustive review of the ALE method is presented in [7].
However, a model of the Reynolds equation incorporating the JFO boundary condition on a moving fluid domain remains currently absent.Therefore, this work will employ the ALE method for the Reynolds equation, enabling the transfer of its reference system to a fluid domain that can be moved or deformed arbitrarily.The weak formulation will also be derived, allowing it to be seamlessly integrated with the preexisting cavitation formulation and subsequently solved through the finite element method.Moreover, this novel formulation can be coupled with solid deformations and contact mechanics [8] in future studies to systematically investigate the contact problem involving partially filled lubricants [9,10] and establish comparisons with experimental data [11,12].

REYNOLDS EQUATION CONSIDERING THE JFO BOUNDARY CONDITION
This section provides a concise overview of the Reynolds equation and elucidates the mass-conserving constraints for cavitated bubbles, known as the JFO boundary condition.A more comprehensive understanding can be referred to [2].For illustrative purposes, Figure 1(A) portrays a narrow gap between two solid surfaces (with velocity   = [  ,   ,   ]  ), accommodating both the fluid lubricant and cavitated bubbles.In the lubrication theory, the primary focus lies on the rupture of the pressure distribution within the fluid film caused by the presence of cavitated bubbles rather than on the intricate dynamics of the bubbles themselves.Consequently, the two-phase flow is regarded as a mixture flow.
Given that the film thickness ℎ (in the -direction) within this fluid domain is considerably smaller compared to its dimensions in the and -directions, and taking into account the inherent viscous flow of the lubricant, it becomes feasible to disregard the flow velocity in the -direction [13].Therefore, the pressure  can be approximated as constant along the -direction.Additionally, the flow velocity  = [, ]  (in the and -directions) can be effectively averaged across the film thickness, thus defining an average flow velocity vector ū, derived from the Navier-Stokes equation: Herein, the lubricant is assumed to be incompressible with a viscosity denoted as , and Ū = [ 1 +  2 ,  1 +  2 ]  represents the entraining velocity vector.
Concurrently, the continuity equation guarantees the mass conservation of the mixture, whereby an averaging process along the film thickness yields in which the density ρ corresponds to the average mixture density in the -direction.Substituting the derived average flow velocity (1) into the averaged continuity equation ( 2) results in the famous Reynolds equation (for more details, please refer to Dowson [13]): It effectively reduces the original three-dimensional fluid domain to Ω  defined within the -plane, further partitioned into two distinct regions, the full-film region Ω +  and the cavitated region Ω −  , based on the presence of bubbles (refer to Figure 1B).
It is essential to acknowledge that the application of the Reynolds equation ( 3) is invalid at the boundary Γ  between Ω +  and Ω −  , as the pressure gradient may become indefinite due to film rupture caused by the presence of bubbles.Hence, it becomes necessary to introduce a new equation on Γ  to ensure the conservation of mass: In this context,  = ρ ūℎ represents the mass flow rate, while  signifies the unit normal vector on Γ  .Additionally, the superscripts + and − indicate the limiting values of the physical quantities adjacent to the Ω +  and Ω −  sides of Γ  , respectively.  denotes the moving velocity (distinct from the flow velocity) of Γ  in the normal direction.
Moreover, by assuming a constant pressure   within the cavitated region and a density of  0 for the incompressible lubricant, the following constraints can be derived: These constraints, along with the mass conservation equation ( 4), collectively constitute the JFO boundary conditions for the Reynolds equation [1].

AN ALE FRAMEWORK FOR THE REYNOLDS EQUATION
It is worth emphasizing that the Reynolds equation, as introduced in Section 2, employs the conventional Eulerian description, whereby fixed spatial points serve as the reference system.While this approach effectively captures the evolution of physical field quantities within a stationary fluid domain, it is inadequate in EHL problems that involve fluid flow and solid motion or deformation.
In such scenarios, the solid surface's motion or deformation alters the fluid domain's boundaries.A varying control volume or fluid domain (see Figure 2) is often employed to accurately track the fluid's physical field quantities near specific characteristic points (e.g., maximum deformation) on the solid surface.In this context, the Eulerian description proves inadequate as it necessitates the construction of a fluid domain large enough to encompass the entirety of the solid's trajectory.Alternatively, the ALE description method, which regards the fluid domain in motion as the reference system, offers a more suitable approach.

Reynolds equation using ALE description
Under the assumption that at the initial time  0 , every point within the reduced fluid domain Ω 0 (depicted in Figure 2) corresponding to the Reynolds equation can be represented by a two-dimensional parameter vector  ∈ Ξ, it follows that at any subsequent time  ∈  = [ 0 ,   ], each parametrical point  can be mapped onto the two-dimensional Euclidean space (-plane) through a function , which characterizes the motion or deformation of the fluid domain: The parameter vector  of each point remains fixed over time, akin to a material point in continuum mechanics.However, it represents not a fluid particle but a geometric point within the fluid domain.Therefore, the moving velocity v of a parametric point  can be defined as with the symbol |  signifying that the parameter vector  remains constant.
In order to shift the reference frame of the Reynolds equation from the spatial points to the parameter points residing on the moving fluid domain, a scalar field (, ) defined on a two-dimensional Euclidean space can initially be examined.This scalar field corresponds to a distribution function  ( , ) constructed over the parameter vector  , from which the following relation can be readily derived: Performing partial time derivatives on both sides leads to By assigning ρℎ as the value of  in Equation ( 9) and substituting it into the Reynolds Equation (3) utilizing the Eulerian description method, the ALE form of the Reynolds equation can be derived as The reference frame of this equation can be arbitrary depending on the selected motion velocity v for the fluid domain.

Weak formulation for the moving fluid domain
It is necessary to establish appropriate finite element spaces to obtain a weak form of the strong formulation (10) compatible with the finite element method.Assuming that Dirichlet boundary conditions are imposed on both the unknowns  and ρ, the solution spaces   and  ρ , as well as a weight space   , can be defined for any given time  ∈  as where ℍ 1 denotes the Hilbert space over the moving fluid domain Ω  .According to the Lax-Milgram theorem, the solution to Equation (10) must adhere to the following weak formulation: Find  ∈   and ρ ∈  ρ , such that ∀  ∈   , Furthermore, applying the Reynolds transport theorem to the moving fluid domain leads to As the weight function   remains invariant for a specific point  over time, Equation ( 15) can be reformulated in the following manner: By substituting it into Equation ( 14), the weak formulation of the Reynolds equation within the ALE framework can be derived as Regarding the cavitation constraints (5), it is permissible to set the pressure   within the cavitated region to zero, as it is solely the pressure gradient that influences Equation (17).The void fraction  = 1 − ρ∕ 0 can be introduced to describe the volume fraction occupied by the bubbles in the film thickness direction.Consequently, the cavitation constraints ( 5) can be expressed in a complementary manner as To address these constraints, a solution space and a weight space can be defined as A weak formulation based on the AL method was developed in Tong et al. [2] to replace the complementary constraints (5): where  is a positive predefined penalty parameter.This formulation effectively transforms the original constrained problem into an unconstrained one.Additionally, substituting ρ = (1 − ) 0 into Equation ( 17) and performing partial integrations yields It is crucial to note that the second integral in the weak formulation (22) corresponds to a convective term, which, when discretized using the continuous Galerkin method, can lead to nonphysical oscillations in the obtained results.Therefore, an appropriate stabilization technique, such as the streamline-upwind/Petrov-Galerkin (SUPG) or the Galerkin-leastsquares (GLS) method, becomes imperative.In this study, the SUPG method is employed as an illustrative example.It introduces a perturbation term   to modify the weight function associated with the convective terms, defined as with  representing an element-size dependent parameter (for a more comprehensive exposition, please refer to [2,14]).
The resulting stabilized weak formulation, Π , , can be expressed as Furthermore, it has been demonstrated in [2] that by discretizing the weak formulation (17), the mass conservation condition (4) on Γ  can be automatically satisfied using the continuous Galerkin method, eliminating the need for redundant procedures.Finite-dimensional subspaces  ℎ  ⊂   ,  ℎ  ⊂   ,  ℎ  ⊂   , and  ℎ  ⊂   can be established to achieve this, employing continuous basis functions to estimate   ,   ,   , and   , respectively.Assuming that all subspaces are estimated using the same set of basis functions  = [ 1 ,  2 , … ,   ]  , which follows the concept of isoparametric formulation, all unknowns and corresponding weight functions can be expressed in discretized form: By substituting them into Π , and Π  , the discretized weak formulations Π ℎ , and Π ℎ  can be obtained.Their solution can be accomplished by employing the classical Newton's method.The detailed derivation is extensively discussed in [2], thus avoiding redundancy in this current discussion.
The formulation presented in this paper does not impose specific requirements on the selection of basis functions, which can encompass both Lagrange polynomials commonly employed in traditional finite element analysis and nonuniform rational B-spline (NURBS) functions utilized in isogeometric analysis.Nevertheless, in terms of the robustness of higher order basis functions and the combination with the mortar method, NURBS functions possess a notable advantage over Lagrange polynomials, as shown in a detailed comparison provided in [2].

Numerical example
This section aims to demonstrate the practical implementation of the newly developed formulation in the context of transient hydrodynamic lubrication.As depicted in Figure 3, a small spherical section is appended to the upper flat plate, inducing a constant velocity of   = 10 m/s in the -direction.The lower flat plate remains fixed, resulting in a timevarying gap height distribution that squeezes the fluid lubricant toward the rear.Given the focus of this study on exploring the applicability of the ALE formulation for the Reynolds equation while considering the JFO boundary condition, elastic deformations are initially disregarded, and both bodies are modeled as rigid entities.
To accurately capture the pressure distribution of the fluid and the volume fraction of cavitation bubbles near the spherical section, a rectangular fluid domain is constructed for the Reynolds equation, employing NURBS basis functions.This NURBS surface utilizes a denser arrangement of control points near the initial position of the spherical section.As the upper solid surface initiates its motion, the control points in its vicinity move synchronously so that the ALE velocity in Equation ( 22) is computed as Ū∕2 − v = [−  ∕2, 0]  , ensuring that the points within the fluid domain track the motion of the spherical section.The outcomes are recorded in Figure 4, where Figures 4(A)-(C) illustrate the temporal evolution of the pressure distribution within the moving fluid domain.Notably, the peak pressure in the region ahead of the spherical section is rapidly developed and remains relatively constant at the selected time samples.The forward propulsion of the spherical section instigates the formation of a substantial cavitated region in its wake, ultimately leading to the development of the trailing flow displayed in Figure 4(F).Upon examination, these findings align with those using the Eulerian description.

CONCLUSION
The present work employs the ALE method to modify the Reynolds equation, enabling the selection of an arbitrarily moving fluid domain as the reference system.Furthermore, the weak formulation for the fluid domain in motion is derived and integrated with a well-established cavitation formulation, resulting in a novel ALE-based cavitation formulation.This formulation can be efficiently solved using the finite element method and effectively addresses two-phase flow problems associated with moving boundaries in the field of lubrication.In future investigations, the formulation can be extended to accommodate large solid deformations, allowing for the resolution of complex transient EHL phenomena.

A C K N O W L E D G M E N T S
The authors gratefully acknowledge the financial support provided by the German Research Foundation (DFG) -project number 390252106.
Open access funding enabled and organized by Projekt DEAL.

F I G U R E 1
The two-phase flow in a narrow gap between two solid bodies.(A) The cross section along the -plane, the fluid domain can be reduced to Ω  for the Reynolds equation.(B) Top view of Ω  , divided in full-film region Ω +  and cavitated region Ω −  .

F I G U R E 2
Moving control volume following the relative motion of the solid surface leads to the temporally varying reduced fluid domain for the Reynolds equation: Ω 0 at time  0 and Ω  at time .

F I G U R E 3
Moving mesh adaption for the fluid domain tracking around a moving sphere section, in order to capture the nearby field quantities (pressure and void fraction).F I G U R E 4 Results evaluated on the moving fluid domain for the Reynolds equation.(A)-(C) Evolution of the pressure distributions at  = 0.2 ms,  = 1 ms, and  = 2 ms, respectively.(D)-(F) Evolution of the void fraction at  = 0.2 ms,  = 1 ms, and  = 2 ms, respectively.