Elastohydrodynamic Squeeze-Film Lubrication Using Fluid-Structure Interactions
Overview
Elastohydrodynamic lubrication (EHL) examines how a thin lubricant film pressurizes to support load across two deformable elastic bearing surfaces. In classical tribology, hydrodynamic (fluid film) lubrication is generally studied by solving the Reynolds equation for the fluid pressure $p$ in a film whose variable thickness $h$ is given. The Reynolds equation is reduced from the general Navier-Stokes equation for incompressible Newtonian fluids, in the limit when the film thickness is much smaller than the characteristic size of the footprint of the pressurized fluid. If the bearing surfaces deform significantly under the action of the fluid pressure $p$ (i.e., if the deformation is comparable to or greater than the film thickness $h$), it becomes necessary to solve for the bearing surface deformations, typically using linear isotropic elasticity theory for the bearing material. This is the essence of EHL, which combines the solution of the Reynolds equation for fluid flow and pressurization with the elasticity solution for the bearing deformation. For metal bearings, the film thickness may be on the order of one micron or less, whereas the metal dimensions are on the order of millimeters or greater. Therefore, the deformations of the bearing material remain in the range of infinitesimal strains.
In biomechanics, EHL may be of interest for understanding the mode of lubrication between articular surfaces of diarthrodial joints, or between a tendon and its sheath, and other similar applications. Since biological soft tissues typically undergo large deformations, it may become necessary to use finite deformation elasticity for such applications. Since biological lubricants are typically non-Newtonian, exhibiting shear-thinning properties, the Reynolds equation is not strictly applicable for such problems.
Fortunately, FEBio provides the necessary tools for solving EHL problems under finite deformation of the bearing surfaces, using non-Newtonian fluids when desired. These capabilities exist in the fluid-FSI module, which implements fluid-structure interactions. In this case study, we used FEBio to examine the elastohydrodynamic squeeze-film lubrication between a rigid cylindrical surface and a deformable elastic solid layer under the action of a step load. Unlike traditional EHL analyses that rely on solving the Reynolds equation while using infinitesimal strain elasticity for the bearing surface deformations, FEBio solves the full set of momentum and mass balance equations for the fluid and solid domains, for any desired constitutive relation.
Finite Element Model
Geometry
In this case study we examine a fluid film being squeezed between a rigid cylinder of radius $R$ and an elastic layer bonded to a rigid substrate. The rigid cylinder is subjected to a load intensity $W$, squeezing out the fluid. The fluid pressure in the film is transiently supporting the load $W$ until it gets depleted. How long it takes for the film to get depleted represents an essential aspect of this problem.
For computational efficiency, we chose to perform a two-dimensional plane strain analysis in the $XY-$plane (Figure 1). We also took advantage of the symmetry about the $YZ-$plane to cut the model in half. In this case, the prescribed load represents a load intensity $W$ (load per unit depth along the $Z-$axis). We modeled the rigid cylinder with a radius $R=500\,\text{mm}$. The elastic layer had a half-width of $20\,\text{mm}$ along the $X-$axis, a thickness of $1\,\text{mm}$ along the $Y-$axis, and a depth of $1 \,\text{mm}$ along the $Z-$axis. The initial minimum separation between the cylinder and elastic solid bearing surface was set to $0.1 \,\text{mm}$ (representing the minimum film thickness at the initial time $t=0$).
Material Models
The solid was modeled with a Mooney-Rivlin uncoupled elastic material with $c_1=0.25 \,\text{MPa}$, $c_2=0$ and $k=250\,\text{MPa}$. Since fluid-FSI analyses are inherently dynamic analyses, applying a step load on an elastic solid may produce an oscillatory dynamic response similar to a spring-mass system. To prevent this type of behavior and hew more closely to traditional EHL, the mass density $\rho$ of the solid material was set to zero.
Two models were used for the fluid-FSI material representing the lubricant, one where the fluid was modeled as Newtonian viscous, and the other where it was modeled as a Carreau (shear-thinning) fluid. For both models, the fluid mass density $\rho$ was set to zero (for the same reason as given above) and the bulk modulus was set to $2.2\,\text{GPa}$ (which is the bulk modulus of water at room temperature and atmospheric pressure). For the Newtonian model, the shear viscosity was set to $\mu=0.001 \,\text{Pa}\cdot\text{s}$ (comparable to the viscosity of water). For the Carreau fluid, properties of synovial fluid were obtained from the literature Tirtaatmadja-1984Tirtaatmadja, V., Boger, D.V. and Fraser, J.R.E., 1984. The dynamic and steady shear properties of synovial fluid and of the components making up synovial fluid. Rheologica Acta, 23(3), pp.311-321.: $\mu_{0}=5\,\text{Pa}\cdot\text{s}$, $\mu_{\infty}=0.001\,\text{Pa}\cdot\text{s}$, $\lambda=55 \,\text{s}$ and $n=0.41$ (Figure 2). These properties show that the viscosity of synovial fluid may vary over three orders of magnitude with increasing shear rate $\dot{\gamma}$.
Mesh
The 8-node hexahedral mesh for this geometry was generated in Cubit (https://cubit.sandia.gov). It consists of two domains, one for the solid material and the other for the fluid-FSI material representing the lubricant domain (Figure 1). Each domain had 200 uniformly distributed elements along $X$, 12 elements along $Y$ using a double bias with a factor of 1.25, and 1 element along $Z$.
Boundary Conditions
Boundary conditions consisted of fixing the displacement along $X$ on the left faces, along $X$ and $Y$ on the bottom face, and along $Z$ on the front and back faces (Figure 1). On the fluid domain only, the fluid (relative) velocity was set to zero along $X$ at the left face, along $X$ and $Y$ at the top and bottom faces (to model the no-slip condition), and along $Z$ at the front and back faces. The fluid dilatation was set to zero on the right face of the fluid domain.
A rigid body was interfaced with the top surface of the fluid domain. The rigid body was constrained from translating along $X$ and $Z$, and rotating about $X$, $Y$ and $Z$. The load intensity $W$ was prescribed on the rigid body as a step load equal to $-5 \, \text{N}/\text{mm}$ along the $Y-$axis.
A fluid-FSI interface was set between the fluid and solid domains, to transfer the fluid traction components to the elastic layer, and between the fluid domain and rigid body, to transfer the prescribed load $W$ to the fluid.
Analysis Settings
For both models, the analysis settings requested $10^6$ time steps with step size of $10^{-6} \,\text{s}$, for a total analysis time of $1\,\text{s}$. The step size was allowed to increase by one decade for every increasing decade of time starting from $10^{-4}\,\text{s}$. Broyden’s quasi-Newton method was used in both cases, with 5 maximum reformations and 50 maximum updates. For the Newtonian model, the parameters $\mathtt{diverge\_reform}$ and $\mathtt{reform\_each\_time\_step}$ were set to $\mathtt{false}$ and the maximum number of retries was set to 5. For the Carreau case, due to its highly nonlinear behavior, the parameter $\mathtt{diverge\_reform}$ was set to $\mathtt{true}$ and the maximum number of retries was set to 20. Both models set the optimal iterations parameter to 50.
Results
The Newtonian fluid case ran until time $t=0.41\,\text{s}$, when it finally failed to converge after the requested maximum reformations and updates (total solution time was 38 minutes on a high-end desktop computer with 18 cores). The Carreau fluid case ran to completion until time $t=1\,\text{s}$ (total solution time of nearly nine hours).
The minimum film thickness ($h_{\min}$) was evaluated by subtracting the positions of nodes at the top and bottom of the left face of the fluid domain. When plotted on a semi-log scale, results showed a rapid initial decrease in film thickness, followed by a slower decay. The film thickness decreased much more rapidly with the Newtonian lubricant, when compared to the Carreau lubricant (Figure 3). The Newtonian analysis ended with $h_{\min}=0.3 \, \mu\text{m}$ at $t=0.41\,\text{s}$, whereas the Carreau analysis ended with $h_{\min}=0.9 \, \mu\text{m}$ at $t=1\,\text{s}$.
The profiles of the rigid cylinder and the elastic solid bearing surface are presented at three time points, in the case of the Newtonian lubricant, corresponding to the times when $h_{\min}=50 \, \mu\text{m}$, $h_{\min}=10 \, \mu\text{m}$ and $h_{\min}=1 \, \mu\text{m}$ (Figure 4).
The fluid film pressure at same three time points is shown in Figure 5 for the Newtonian and Carreau lubricants.
Figure 5. Fluid film pressure $p$ for (a) the Newtonian lubricant and (b) the Carreau lubricant, at three time points corresponding to $h_{\min}=50$, $10$ and $1 \, \mu\text{m}$.
In the Carreau fluid analysis, the value of the shear viscosity $\mu$ can be examined throughout the fluid domain, at various time points. Since the shear viscosity is dependent on the shear rate $\dot{\gamma}$, and since the shear rate varies along $X$ and $Y$, contour maps of the viscosity are presented in Figure 6 for the time points corresponding to a minimum film thickness of $h_{\min}=50 \, \mu\text{m}$, $h_{\min}=10 \, \mu\text{m}$ and $h_{\min}=1 \, \mu\text{m}$. These maps show that the shear viscosity was lowest at the earliest time points, increasing by two orders of magnitude as the film thickness $h$ decreased.
Figure 6. Contour maps of the shear viscosity $\mu$ of the Carreau lubricant, over the (undeformed) fluid domain, at three different time points corresponding to (a) $h_{\min}=50 \, \mu\text{m}$, (b) $h_{\min}=10 \, \mu\text{m}$ and (c) $h_{\min}=1 \, \mu\text{m}$. The viscosity increases by two orders of magnitude over time.
Discussion
This case study illustrates that EHL problems may be solved using the fluid-FSI module in FEBio. The results for the Newtonian lubricant converge well and replicate the expected temporal evolution of the minimum film thickness (Figure 3), the bearing surface deformation (Figure 4) and the fluid film pressure (Figure 5a). The Newtonian lubricant analysis terminated before reaching $t=1 \,\text{s}$ because the fluid mesh became highly compressed, reducing from $h_{\min}=100 \, \mu\text{m}$ down to $h_{\min}=0.3 \, \mu\text{m}$, which is greater than a 300-fold reduction in height. Of course, squeeze-film lubrication analyses of this type start from an arbitrarily chosen initial film thickness. Had the analysis been initiated at $h_{\min}=10 \, \mu\text{m}$, it is possible that the final film thickness achieved during the analysis could have been reduced by a similar factor.
Though the elastic bearing surface deformation was on the order of $30 \,\mu\text{m}$ (Figure 4), it would be erroneous to assume that these small deformations remained in the range of infinitesimal strains. In these analyses the peak value of the maximum principal Lagrangian strain in the elastic solid layer was greater than 0.16 (corresponding to a stretch ratio of 1.15 or engineering strain of 15%).
In principle, higher load intensities $W$ may be applied for these analyses. However, the application of a step load in squeeze-film lubrication problem may become problematic; computationally, the resulting instantaneous reduction in film thickness may cause element inversions; furthermore, from a physics-based argument, inertial effects would prevent any load application to be instantaneous. Therefore, for higher magnitudes of $W$ it may become necessary to prescribe a time-dependent increase of its value from 0. For example, $W$ may be ramped up linearly from 0 to the desired value over a $10 \,\mu\text{s}$ duration.
The analysis of squeeze-film EHL using a Carreau model with synovial fluid properties represents a significant innovation, both computationally and from the perspective of understanding articular cartilage lubrication. Computationally, it is remarkable that a solution could be found, despite the highly nonlinear variation of the shear viscosity $\mu$ with varying shear rate $\dot{\gamma}$ (Figure 2). However, this solution came at great computational expense, requiring a much longer solution time than that of the Newtonian lubricant. Nevertheless, the fact that a computational solution could be obtained remains a significant achievement.
From a physics perspective, it is interesting that the higher viscosity of SF fluid appears to be more influential as the squeeze-film response progresses in response to a step load, with less of an effect immediately upon loading (Figure 3). Despite the fact that the synovial fluid viscosity is relatively very high under the initial rest conditions ($\mu_0 = 5 \text{Pa}\cdot\text{s}$), the early time response exhibits a much lower viscosity, closer to that of water (Figure 6a). This occurs because the peak fluid film velocity is on the order of $15 \, \text{m}/\text{s}$ and the shear rate reaches $\dot{\gamma}=8.8\times10^5\,\text{s}^{-1}$ when $h_{\min}=50 \,\mu\text{m}$. Thus, despite the fact that the fluid film is thicker in the early time response, the shear rate is very large. Conversely, as the fluid film becomes thinner, e.g., when $h_{\min}=1 \,\mu\text{m}$, the peak velocity drops to less than $10\,\text{mm}/\text{s}$ while the shear rate drops to $\dot{\gamma}=3.3\times10^4\,\text{s}^{-1}$, producing a higher shear viscosity $\mu$ (Figure 6c). Of course, the shear rate varies considerably over the fluid film domain, implying that the calculations presented here are relevant to the minimum values of the shear viscosity $\mu$ over the fluid film domain. The contour maps presented in Figure 6 provide a more comprehensive description of $\mu$ over the entire domain.
A film thickness of $1\,\mu\text{m}$ is comparable or smaller to the characteristic surface roughness of articular cartilage. This film thickness is achieved in slightly less than $1\,\text{s}$ with the synovial fluid lubricant (Figure 3). For lower extremity joints, a duration of $1\,\text{s}$ is comparable to the duration of a gait cycle, whereas upper extremity joints may sustain longer durations of loading under activities of daily living, leading to further reduction of the fluid film thickness. By the conventional standards of EHL theory, this result implies that the lubrication mode of articular cartilage is mixed (i.e., a mix of fluid-film and boundary lubrication) under short-term loading durations. Nevertheless, it is very interesting to see that the behavior of a Carreau lubricant deviates quite significantly from that of a Newtonian lubricant as time increases (Figure 3), leading to conditions that are more favorable to EHL lubrication.
For the types of nonlinear analyses reported in this case study, there is no available closed-form solution that may serve as a verification of the computational results. However, an important limiting case can be examined, representing the condition when the fluid film has been completely depleted. In that case, the rigid cylindrical bearing surface indents the elastic solid layer under “dry” contact conditions, which can also be simulated in FEBio, using the structural mechanics module. A comparison of the fluid-film pressure in the EHL analysis, near the point of film depletion ($h_{\min}=1\,\mu\text{m}$), against the dry contact pressure shows excellent agreement (Figure 7a), as well as for the deformation of the elastic bearing surface (Figure 7b). This agreement provides confidence in the validity of the EHL analysis in FEBio.
Figure 7. Comparison of Newtonian EHL and dry contact analyses for (a) fluid-film versus contact pressure and (b) elastic solid surface deformation. EHL results correspond to the time point when $h_{\min}=1\,\mu\text{m}$, approaching the limit of depletion of the fluid film.
In summary, this case study demonstrates that FEBio can solve challenging problems in the field of elastohydrodynamic lubrication in biological structures, where tissues may undergo large deformations and lubricants may exhibit strong non-Newtonian behaviors.
Footnotes
- Tirtaatmadja-1984Tirtaatmadja, V., Boger, D.V. and Fraser, J.R.E., 1984. The dynamic and steady shear properties of synovial fluid and of the components making up synovial fluid. Rheologica Acta, 23(3), pp.311-321.