Squeeze-Film Lubrication Using Biphasic Fluid-Structure Interactions

Squeeze-Film Lubrication Using Biphasic Fluid-Structure Interactions

You are here:
< Back


Squeeze-film lubrication examines how a thin lubricant film pressurizes to support load across two bearing surfaces. If one or both of these bearing surfaces is porous, permeable and deformable, allowing the lubricant to flow in or out of it, it is possible to use FEBio’s biphasic fluid-structure interaction material to analyze this type of lubrication. This problem was first solved semi-analytically in the study of Hou et al. Hou-1992Hou JS, Mow VC, Lai WM, Holmes MH. An analysis of the squeeze-film lubrication mechanism for articular cartilage. J Biomech. 1992 Mar;25(3):247-59. doi: 10.1016/0021-9290(92)90024-u. PMID: 1564060. and reprised in the study of Hlavácek Hlavacek-1993aHlavácek M. The role of synovial fluid filtration by cartilage in lubrication of synovial joints–I. Mixture model of synovial fluid. J Biomech. 1993 Oct;26(10):1145-50. doi: 10.1016/0021-9290(93)90062-j. PMID: 8253819.. Here, we perform a similar analysis using the finite deformation framework of FEBio. This type of lubrication analysis may be of interest for understanding the mode of lubrication between articular surfaces of diarthrodial joints.

Finite Element Model


The geometry employed here is the same as the one used in the case study on elastohydrodynamic squeeze-film lubrication. We examine a fluid film being squeezed between a rigid cylinder of radius $R$ and a biphasic layer bonded to a rigid substrate. The rigid cylinder is subjected to a load intensity $W$, squeezing out the fluid sideways, with potential fluid exchanges with the biphasic layer. 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$).

Figure 1. Finite element mesh for biphasic-FSI domain representing a tissue layer (green) and fluid-FSI domain representing the lubricant (purple). The top surface of the fluid domain is cylindrical with radius $R=500$ mm. The width of the biphasic domain is $20$ mm and its height is $1$ mm. The initial minimum film thickness $h$ is 0.1 mm.

Material Models

The tissue layer was modeled with a biphasic-FSI material with a solid volume fraction set to $0.2$, a Newtonian viscous fluid with bulk modulus equal to $2200 \,\text{MPa}$ (which is the bulk modulus of water at room temperature and atmospheric pressure), and shear viscosity equal to $10^{-9}\,\text{MPa}\cdot\text{s}$. The solid matrix was neo-Hookean with Young’s modulus set to $0.5\,\text{MPa}$ and Poisson’s ratio set to $0$. A constant isotropic hydraulic permeability model was adopted, with $k=0.001\,\text{mm}^4/\text{N}\cdot\text{s}$. Since biphasic-FSI analyses are inherently dynamic analyses, applying a step load on a biphasic-FSI material may produce an oscillatory dynamic response similar to a spring-mass-damper system. To prevent this type of behavior and hew more closely to traditional squeeze-film lubrication analyses, the mass density of the solid and fluid materials were set to zero.

The lubricant was modeled using a fluid-FSI material, using a Newtonian viscous fluid. The fluid mass density was set to zero (for the same reason as given above) and the bulk modulus was set to $2200\,\text{MPa}$. The shear viscosity was set to $\mu=10^{-9}\,\text{MPa}\cdot\text{s}$ (comparable to the viscosity of water).


The 8-node hexahedral mesh for this geometry was generated in Cubit (https://cubit.sandia.gov). It consists of two domains, one for the biphasic-FSI material representing the tissue layer, 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). In addition, 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 tissue and lubricant domains.

A rigid body was interfaced with the top surface of the fluid-FSI lubricant 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 compressive load equal to $5 \, \text{N}/\text{mm}$ along the $Y-$axis.

A FSI interface traction was set between the lubricant domain and the rigid body, to transfer the prescribed load $W$ to the fluid.

Analysis Settings

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, with 5 maximum reformations and 50 maximum updates. 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 20. The optimal iterations parameter was set to 50.


The analysis ran until time $t=0.79\,\text{s}$, when it finally failed to converge after the requested maximum reformations and updates (total solution time was 58 minutes on a high-end desktop computer with 18 cores).

The central film thickness ($h_{0}$) 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, then another rapid decrease (Figure 2). The analysis ended with $h_{0}=6\times 10^{-7} \, \text{mm}$ at $t=0.79\,\text{s}$.

Figure 2. Central film thickness $h_{0}$ versus time, for biphasic squeeze-film lubrication analysis.

The profiles of the rigid cylinder and the biphasic bearing surface are presented at three time points, corresponding to the times when $h_{0}=50 \, \mu\text{m}$, $h_{0}=10 \, \mu\text{m}$ and $h_{0}=1 \, \mu\text{m}$ (Figure 3).

Figure 3. Profiles of the rigid cylinder and biphasic bearing surface for three different values of $h_{0}$ ($50 \, \mu\text{m}$ in green, $10 \, \mu\text{m}$ in orange, and $1 \, \mu\text{m}$ in blue). Note the uneven scales along $x$ and $y$.

The fluid film pressure, which is evaluated from the dilatation $e^f$ using the fluid’s bulk modulus, is displayed at the same three time points in Figure 4. It is notable that the fluid pressure becomes slightly negative at the edge of its footprint, during intermediate time points as the fluid film shrinks in thickness, as seen in the case of $h_0=10\,\mu\text{m}$ in this figure. However, as the film thickness reduces further, the fluid pressure becomes positive or zero across the entire fluid film domain. Oscillations in the fluid pressure at $h_0=1\,\mu\text{m}$ are indicative of the numerical stiffness of the governing equations as the film thickness decreases. Though not shown here, the amplitude of oscillations does not increase further as the film thickness shrinks to its lowest value shown in Figure 2.

Figure 4. Fluid film pressure $p$ at three time points corresponding to $h_{0}=50$, $10$ and $1 \, \mu\text{m}$.

Since fluid may be exchanged at the interface between the fluid lubricant film and the biphasic layer, we may also examine the fluid flux $w_y$ in the direction normal to the interface. The normal flux is shown in Figure 5 for the time points corresponding to a central film thickness of $h_{0}=50 \, \mu\text{m}$, $h_{0}=10 \, \mu\text{m}$ and $h_{0}=1 \, \mu\text{m}$. These results show that fluid flows from the film into the biphasic layer in the central region of the pressure footprint, while it flows from the biphasic layer into the film in the outer region.

Figure 5. Fluid flux $w_y$ normal to the interface between the fluid film and biphasic layer, at time points corresponding to $h_{0}=50 \, \mu\text{m}$, $h_{0}=10 \, \mu\text{m}$ and $h_{0}=1 \, \mu\text{m}$. A positive value of $w_y$ indicates fluid flowing from the biphasic layer into the film.


This case study illustrates that biphasic hydrodynamic lubrication problems may be solved using the fluid-FSI module in FEBio, with a biphasic-FSI material representing the porous, permeable and deformable tissue, and a fluid-FSI material representing the lubricant. The temporal evolution of the central film thickness (Figure 2) demonstrates three phases in the evolution of the lubricant film: First, there is a rapid loss of film thickness due to lateral squeezing of the lubricant out of the film; as the film thickness reduces considerably during that phase, the $x-y$ shear rate correspondingly rises, leading to a significant increase in the resistance to flow along the lateral direction; in the third phase, as the film thickness reduces much further, the path of least resistance for the lubricant film is to flow into the biphasic bearing layer, causing a rapid decrease in film thickness until the film collapses, leading to boundary contact. These results are qualitatively consistent with the literature reports in the Overview section.

The material properties employed in this analysis are relevant to articular carttilage. The surface roughness of healthy cartilage is on the order of $0.5\,\mu\text{m}$. Therefore, when the film thickness has reduced to $1\,\mu\text{m}$ we expect that the fluid-film (hydrodynamic) lubrication regime has transitioned to mixed lubrication, where asperities on the bearing surfaces may start coming into contact, possibly trapping lubricant pools in small pockets. According to Figure 2, $h_0$ reduces to $1\,\mu\text{m}$ at approximately $0.16\,\text{s}$ after load application, long before the third phase of rapid film depletion takes place. Indeed, the bearing surface profiles become highly congruent at that juncture, as shown in Figure 3. Therefore, any further decrease in film thickness that occurs during the last phase evidently leads to boundary contact of the bearing surfaces, thus transitioning to the boundary lubrication regime.

The fluid film pressure reported in Figure 4 also represents the interstitial fluid pressure in the biphasic layer, since the fluid pressure must be continuous across porous interfaces in biphasic theory. As the central film thickness decreases to $h_0=1\,\mu\text{m}$, the fluid pressure footprint loses its “tail” at the outer edge of the pressure footprint. This result suggests investigating the limiting case when the bearing surfaces come into direct contact, without an intervening fluid film. This type of boundary contact analysis may be simulated in FEBio using the biphasic module, employing the exact same mesh and material properties for the biphasic layer. A comparison of the fluid-film pressure in the squeeze-film analysis near the point of film depletion ($h_{0}=1\,\mu\text{m}$), against the contact pressure of the biphasic contact analysis shows excellent agreement (Figure 7a), as well as for the deformation of the biphasic bearing surface (Figure 7b). This agreement confirms that the fluid-film pressure in the squeeze-film analysis approaches the contact pressure as the film thickness reduces, consistent with classical findings in elastohydrodynamic lubrication. It is also interesting to note that the interstitial fluid pressure in the biphasic contact analysis (green curve in Figure 7a) is slighly smaller in magnitude than the contact pressure at the center of the contact footprint, though it slightly exceeds the contact pressure at the outer edge. Integrating the contact pressure over the contact surface recovers the load intensity $W=5\,\text{N}/\text{mm}$, whereas integrating the biphasic interstitial fluid pressure produces a fluid load intensity $W^p=4.92\,\text{N}/\text{mm}$. This result indicates that the solid matrix of the biphasic layer contributes to supporting a small fraction of the applied load when boundary contact is taking place, in contrast to fluid-film lubrication where all of the load is supported by the fluid lubricant.


Figure 7. Comparison of biphasic fluid-film lubrication and biphasic contact analyses for (a) fluid-film versus contact pressure and (b) biphasic surface deformation. Fluid-film lubrication results correspond to the time point when $h_{0}=1\,\mu\text{m}$, approaching the limit of depletion of the fluid film.


  • Hou-1992
    Hou JS, Mow VC, Lai WM, Holmes MH. An analysis of the squeeze-film lubrication mechanism for articular cartilage. J Biomech. 1992 Mar;25(3):247-59. doi: 10.1016/0021-9290(92)90024-u. PMID: 1564060.
  • Hlavacek-1993a
    Hlavácek M. The role of synovial fluid filtration by cartilage in lubrication of synovial joints–I. Mixture model of synovial fluid. J Biomech. 1993 Oct;26(10):1145-50. doi: 10.1016/0021-9290(93)90062-j. PMID: 8253819.
Was this article helpful?
0 out of 5 stars
5 Stars 0%
4 Stars 0%
3 Stars 0%
2 Stars 0%
1 Stars 0%
How can we improve this article?
Please submit the reason for your vote so that we can improve the article.
Table of Contents
Go to Top