Implementation of a Finite Element Deformation Model Within an Elasto-Hydrodynamic Lubrication Numerical Solver for a Ball in Socket Tribopair

In the framework of the elasto-hydrodynamic lubrication simulation algorithms of lubricated tribopairs, a key role is played by the chosen deformation model, since it affects the surfaces’ separation, which guarantees the existence of a thin lubricant film thickness, even when the tribo-system is subjected to high loads. The aim of this article is to merge a finite element deformation model based on linear tetrahedra, previously developed by the same authors, within the Reynolds equation solver in the elasto-hydrodynamic mode, with reference to a generic ball in socket lubricated tribo-system. The main novelty of this research is the implementation of the finite element deformation model, allowing the authors to relate the deformation vector to the pressure one through an influence matrix which takes into account the spherical motion of the ball with respect to the socket. The computer code for the problem–solution was developed in a MATLAB environment and simulated a planar motion condition in terms of eccentricity and angular velocity vectors, in order to calculate the meatus fluid pressure field, surfaces’ separation, shear stress, deformation, and wear depth. The integration over time of the output fields led to the time evolution of the load vector, friction torque vector, and wear volume. Moreover, the lubrication algorithm takes into account the fluid non-Newtonian behavior and the surfaces’ progressive geometrical modification over time due to cumulated wear. The obtained results reproduced the classical elasto-hydrodynamic shapes of the involved quantities, following the meatus minimum thickness predicted by the Hamrock–Dowson model; furthermore, it provided information about the mechanical behavior of the whole bodies belonging to the spherical joint thanks to the finite element deformation model.

Introduction

The elasto-hydrodynamic lubrication is a lubrication mechanism established between two contacting deformable surfaces separated by a lubricant meatus. In general, the hydrodynamic lubrication related to two coupled surfaces provides their full separation thanks to the relative motion between them, in particular,

1) The sliding relative motion in the contact plane, when the surfaces are inclined along the motion direction forming a converging duct, leading the lubricant to generate a force acting on the surfaces orthogonally with respect to the contact plane in order to separate them;

2) The squeeze motion occurring when the surfaces approach each other along the contact normal direction, leading the lubricant to produce a separation force while trying to not be squeezed out of the surfaces’ gap.

When the loading closure force acting on the couple is too high to keep the separation between the surfaces, their deformability acts in order to permit full thin-film lubrication also in this case. The elasto-hydrodynamic lubrication mode produces a lubricant pressure shape along the relative motion direction very similar to the dry contact pressure evaluated by the well-established Hertz theory, but additionally, it is characterized by a minimum fluid film thickness located in correspondence of the contact zone outlet due to the relative motion, producing here the well-known pressure spike.

Actually, the elasto-hydrodynamic lubrication numerical simulation refers mostly to the Reynolds equation (Wang and Jin, 2004a; Wang and Jin, 2008; Gao, et al., 2009; Shettar, et al., 2018; Ruggiero and Sicilia, 2020a; Ruggiero and Sicilia, 2020b; Ruggiero, et al., 2020), which elaborates the fluid pressure field for a given separation field and entraining velocity: in this case, the deformation model contributes to the calculation of the surfaces’ separation, so its reliability affects the evaluated results.

In this work, the ball in socket kinematical couple is analyzed, due to its scientific relevance in the framework of the biomechanics (Mattei, et al., 2011; Ruggiero and Sicilia, 2020c). Several deformation models adapted to this spherical geometry were used by the researchers in the scientific literature: once the 2D lubrication grid on the socket internal surface is defined, the common purpose is to evaluate the radial deformation in correspondence of each grid point due to the fluid pressure acting on all the nodes belonging to the grid domain. Three usual approaches can be identified:

1) The proportional models, which assume a proportionality between the local deformation and the local pressure, such that each nodal deformation on the analyzed lubrication grid does not depend on other nodes [e.g., the constraint column model or independent spring model (Jalali-Vahid, 2001; Jalali-Vahid, et al., 2001; Jalali-Vahid, et al., 2003; Ruggiero and Sicilia, 2020b), the elastic foundation model (Fregly, et al., 2003; Halloran, et al., 2005; Pèrez-Gonzàlez, et al., 2008; Mukras, et al., 2010; Askari and Andersen, 2018; Srivastava, et al., 2021), etc.]; it can be used when the two materials composing the couple are characterized by very different stiffness, such as to consider only the deformation of the softer body;

2) The convolution method, which consists in solving the mechanical problem associated with the evaluation of the surface deformation of an elastic half-space due to a distributed surface load (i.e., the pressure field) by applying the discrete Fourier transform or the boundary element method (Evans and Hughes, 2000; Wang, et al., 2003; Wang and Jin, 2004b; Wang and Jin, 2004a; Wang and Jin, 2008; Gao, et al., 2009; Wang, et al., 2009; Shettar, et al., 2018); in this case, each lubrication grid point communicates with the other nodes through the definition of an influence matrix which relates the radial surface deformation field to the pressure one;

3) The finite element method, aims to elaborate displacements and constraint reactions of the whole analyzed bodies, by discretizing them in several finite elements and solving the weak formulation of the mechanical equilibrium problem (Jagatia and Jin, 2001; Halloran, et al., 2005; Ruggiero, et al., 2018; Ruggiero, et al., 2019; Ruggiero and Sicilia, 2021); this approach ensures accurate deformation results, however, it can be very expensive in terms of computational time and could require an additional dedicated software with respect to the one used to solve the lubrication problem.

The purpose of this work is to implement a deformation model based on the discretization of the ball in socket joint in linear tetrahedra finite elements (Ruggiero and Sicilia, 2021), in order to assemble an influence matrix which relates the surfaces’ radial deformation vector to the pressure one generated by the fluid meatus within the joint gap. The objective is to fully merge the finite element model with the lubrication algorithm into a comprehensive MATLAB code which works to solve the lubrication problem at the interface between the ball and socket surfaces, keeping information about the whole bodies composing the kinematical couple, in order to obtain a fluid–solid interaction system in which both the computational fluid dynamics and computational solid model are implemented.

Materials and Methods

The Finite Elements Deformation Model

The spherical joint comprises the ball, i.e., a solid sphere of radius $r$, and the socket, i.e., a solid hollow hemisphere of internal radius $R$ and thickness $H$, respectively, described by the nodal positions ${\mathit{x}}_b$ and ${\mathit{x}}_s$ in Eq. 1 through the spherical coordinates $\rho$, $\theta$, and $\phi$.

$\begin{array}{c}\begin{array}{cc}\begin{array}{c}{\mathit{x}}_b=\left[\begin{array}{c}\rho \sin \theta \cos \phi \\ \rho \sin \theta \sin \phi \\ \rho \cos \theta \end{array}\right]\\ \begin{array}{c}0\le \rho \le r\\ 0\le \theta \le \pi \\ 0\le \phi \le 2\pi \end{array}\end{array}& \begin{array}{c}{\mathit{x}}_s=\left[\begin{array}{c}\rho \sin \theta \cos \phi \\ \rho \sin \theta \sin \phi \\ \rho \cos \theta \end{array}\right]\\ \begin{array}{c}R\le \rho \le R+H\\ 0\le \theta \le \pi \\ 0\le \phi \le \pi \end{array}\end{array}\end{array}\end{array}(1)$

With reference to a previous work of the same authors (Ruggiero and Sicilia, 2021), the linear tetrahedra mesh was created for both the bodies and is shown in Figure 1, where the fixed nodes are highlighted by black triangles and the pressure loaded nodes are identified by red arrows oriented toward the radial loading direction. The pressure vector acts simultaneously on the internal socket surface and external ball surface, while the socket is fixed on the external surface and the ball is characterized by fixed nodes with the $y$-coordinate below a given value.

FIGURE 1

FIGURE 1. Socket (A) and ball (B) linear tetrahedra discretization.

With the developed approach (Ruggiero and Sicilia, 2021), the equilibrium equation associated with the linear elastic mechanical behavior of the two structures is converted in Eq. 2 into the relationship which links the surface radial deformation vector $\mathit{\delta }$ due to the acting surface pressure vector $\mathit{p}$ by the definition of the influence matrix $\mathit{C}$, coming from the original stiffness matrix $\mathit{K}$, the nodal displacements $\mathit{q}$, and the nodal forces $\mathbf{\Phi }$: the procedure is feasible and is achieved by matrix manipulation since the surface radial deformation $\mathit{\delta }$ is hidden in the free part of the nodal displacement $\mathit{q}$, while the pressure nodal forces ${\mathbf{\Phi }}_p$ are directly calculated through the surface pressure vector $\mathit{p}$.

$\mathit{Kq}=\mathbf{\Phi }+{\mathbf{\Phi }}_{\mathit{p}}\to \mathit{\delta }=\mathit{Cp}(2)$

As described next, the Reynolds lubrication equation refers to a surface grid domain fixed with the socket's internal surface, so, while the constitutive Eq. 2 referring to the socket does not need further transformations, the one referring to the ball has to be rotated. The rotation is obtained through the cubic interpolation matrices ${\mathit{J}}_{bs}$ and ${\mathit{J}}_{sb}$ [described in (Ruggiero and Sicilia, 2021)] which, respectively, convert a surface quantity vector defined in the ball surface grid ${\mathit{u}}_b$ into its version defined in the socket surface grid ${\mathit{u}}_s$ and vice versa, as written in Eq. 3. In addition to the possibility of translating from a grid to another one, the two matrices ${\mathit{J}}_{\mathit{bs}}$ and ${\mathit{J}}_{sb}$ contain information about the relative angular motion of the ball with respect to the socket; however, since in the next paragraph the model will be tested for a planar motion with high angular velocity, in this work, the grid on the ball's external surface is considered fixed in the space.

${\mathit{u}}_s={\mathit{J}}_{bs}{\mathit{u}}_b{\mathit{u}}_b={\mathit{J}}_{sb}{\mathit{u}}_b(3)$

Then, writing Eq. 2 for both the spherical joint bodies and rotating all the involved quantities in the socket’s internal surface grid through Eq. 3, the final constitutive deformation model which relates to the total surface deformation vector $\mathit{\delta }$ generated by the lubricant pressure vector $\mathit{p}$ within the gap through the global influence matrix $\mathit{C}$ is obtained in the Eq. 4.

$\{\begin{array}{l}{\mathit{\delta }}_s={\mathit{C}}_s{\mathit{p}}_s\\ {\mathit{\delta }}_b={\mathit{C}}_b{\mathit{p}}_b\\ \mathit{\delta }={\mathit{\delta }}_s+{\mathit{J}}_{bs}{\mathit{\delta }}_b\to \mathit{\delta }=\left({\mathit{C}}_s+{\mathit{J}}_{bs}{\mathit{C}}_b{J}_{sb}\right)\mathit{p}=\mathit{Cp}\\ {\mathit{p}}_s=\mathit{p}\\ {\mathit{p}}_b={\mathit{J}}_{sb}{\mathit{p}}_s\end{array}(4)$

The Elasto-Hydrodynamic Lubrication Model

Under the classical hypotheses, the unsteady Reynolds equation is reported in its general form in Eq. 5: it is a partial differential equation in which the lubricant pressure $p$ field over time $t$ is determined by knowing the surfaces’ separation $h$, the entrainment velocity vector $\mathit{v}$, and the lubricant rheological properties, i.e., the density $\rho$ and viscosity $\mu$. Due to the elasto-hydrodynamic mode, the surfaces’ separation $h$ is directly dependent on pressure $p$ through the deformation model, therefore the Reynolds equation is highly nonlinear.

$\mathbf{\nabla }\cdot \left(\frac{\rho h^3}{12\mu }\mathbf{\nabla }p\right)=\mathbf{\nabla }\cdot \left(\rho h\mathit{v}\right)+\frac{\partial }{\partial t}\left(\rho h\right)(5)$

The quantities involved in the Reynolds equation are characterized by the following (Ruggiero and Sicilia, 2020b), considering that the spatial coordinates related to the analyzed spherical contact are referred to the spherical angles $\theta$ and $\phi$ belonging to the socket's internal surface grid, coming from the scheme shown in Figure 2.

FIGURE 2

FIGURE 2. Ball in socket lubrication system scheme.

The rheological properties are considered as functions of the pressure through:

1) The Dowson–Higginson density model in Eq. 6, in which it is assumed that the lubricant density $\rho$ increases with pressure $p$ in dependence on the coefficients $a_{\rho }$ and $b_{\rho }$ starting from the nominal density ${\rho }_0$ (van Leeuwen, 2009);

$\begin{array}{c}\rho ={\rho }_0\frac{a_{\rho }+b_{\rho }p}{a_{\rho }+p}\end{array}(6)$

2) The product between the cross non-Newtonian model and the Barus one in Eq. 7, which consider, respectively, that the lubricant viscosity $\mu$ varies in a finite range delimited by the given viscosity values ${\mu }_0$ and ${\mu }_{\infty }$ due to the shear-thinning or shear-thickening effects with respect to the shear rate (approximated with the ratio between the sliding velocity $v_{sl}$ and the surfaces’ separation $h$) in dependence on the power-law coefficients $k_{\mu }$ and $n_{\mu }$ (Gao, et al., 2016), and that, moreover, the viscosity experiences the elasto-hydrodynamic exponential increase with pressure $p$ governed by the coefficient ${\alpha }_{\mu }$ (van Leeuwen, 2009).

$\begin{array}{c}\mu =\left[{\mu }_{\infty }+\frac{{\mu }_0-{\mu }_{\infty }}{1+k_{\mu }{\left(\frac{v_{sl}}{h}\right)}^{n_{\mu }}}\right]e^{{\alpha }_{\mu }p}\end{array}(7)$

The sub-models described in the Eqs 6, 7 are very useful in the context of the biomechanical tribological pairs represented by the artificial hip joint, in which the high pressure level reached within the gap could be responsible for an appreciable increase of the synovial fluid density, while the low relative angular speed causes the aggregation of the macromolecules composing the fluid, generating a considerable increase in its viscosity.

The surfaces’ separation $h$ field is obtained in Eq. 8 by the composition of:

1) the geometrical approach $h_g$ due to the eccentricity vector $\mathit{e}$ of the ball center with respect to the socket one, which is projected on the radial unit vector $\widehat{\mathit{r}}$ and subtracted from the radial clearance $c$ (Ruggiero and Sicilia, 2020b);

2) the total surfaces’ radial deformation field $\delta$, which is directly obtained from its vector version $\mathit{\delta }$ coming from the finite element deformation model described above (Ruggiero and Sicilia, 2021);

3) the penetration wear depth field $u_w$ cumulated until the analyzed time instant $t$, which is obtained through the Archard wear model by integrating over time the linear wear rate, i.e., the product between the wear factor function $k_w$ (scaling the nominal wear factor $k_{w_0}$ by a power ${\alpha }_w$ of the ratio between the separation $h$ and the surface arithmetic roughness $R_a$, in order to take advantage of the Archard wear model over the whole analyzed domain, even in that lubricated zones in which the film thickness is much larger than the roughness), the lubricant pressure $p$, and the sliding velocity $v_{sl}$ (Gao, et al., 2017; Ruggiero and Sicilia, 2020b).

$\{\begin{array}{l}c=R-r\to h_g=c-{\mathit{e}}^T\widehat{\mathit{r}}\\ \mathit{\delta }=\mathit{Cp}\to \delta \\ k_w=k_{w_0}{\left(\frac{h}{R_a}\right)}^{-{\alpha }_w}\to u_w={\displaystyle \underset{0}{\overset{t}{\int }}{k}_{w}p{v}_{sl}dt}\end{array}\to h=h_g+\delta +u_w(8)$

The entraining velocity vector $\mathit{v}$ in spherical coordinates is represented by the arithmetic average between the velocities of two points belonging to the contacting bodies (Ruggiero and Sicilia, 2020b), and it is written in Eq. 9 that it comes from the composition of the translational term due to the time derivative of the eccentricity vector $\dot{\mathit{e}}$ and the rotational term being governed by the angular velocity vector $\mathit{\omega }$ of the ball with respect to the socket; finally, the spherical rotation matrix ${\mathit{R}}_s$ is applied in order to move from the Cartesian space to the spherical one. The sliding velocity $v_{sl}$ can be obtained directly from the components of $\mathit{v}$ lying in the contact plane, i.e., $U_{\theta }$ and $U_{\phi }$.

$\begin{array}{c}\mathit{v}=\frac{1}{2}{\mathit{R}}_s^T\left\{\dot{\mathit{e}}+\mathit{\omega }\times \left[\left(R-h\right)\widehat{\mathit{r}}-\mathit{e}\right]\right\}=\left[\begin{array}{c}{U}_{\theta }\\ U_{\phi }\\ U_{\rho }\end{array}\right]\to v_{sl}=2\sqrt{U_{\theta }^2+U_{\phi }^2}\end{array}(9)$

Once all the quantities involved in the Reynolds equation have been identified, it is reported in Eq. 10, turned in to spherical coordinates, and coupled with the boundary pressure $p_0$ in correspondence with the grid domain limits (Ruggiero and Sicilia, 2020b).

$\{\begin{array}{l}\sin \theta \frac{\partial }{\partial \theta }\left(\sin \theta \frac{\rho h^3}{12\mu }\frac{\partial p}{\partial \theta }\right)+\frac{\partial }{\partial \phi }\left(\frac{\rho h^3}{12\mu }\frac{\partial p}{\partial \phi }\right)=R\sin \theta \left[\frac{\partial }{\partial \theta }\left(\sin \theta \rho h{U}_{\theta }\right)+\frac{\partial }{\partial \phi }\left(\rho h{U}_{\phi }\right)+R\sin \theta \frac{\partial }{\partial t}\left(\rho h\right) \right]\\ p\left(0,\phi \right)=p\left(\pi ,\phi \right)=p_0\\ p\left(\theta ,0\right)=p\left(\theta ,\pi \right)=p_0\end{array}(10)$

Eq. 10 is discretized with the finite differences on the socket’s internal surface grid by the definition of some groups, vectors and matrices explained in the authors’ previous work. The objective is to write the Reynolds equation $R_I$ in correspondence to each grid point identified by index $I$, resulting in a nonlinear closed system $\mathit{R}$ which is characterized by the five-diagonal band square matrix $\mathit{A}$ and the vector $\mathit{b}$, both dependent on the pressure vector $\mathit{p}$.

$\begin{array}{c}{\mathit{f}}_s^T{\mathit{D}}_{2_s}{\mathit{p}}_s-\left({\mathit{D}}_{{\theta }_s}^T{\mathit{u}}_{{\theta }_s}+{\mathit{D}}_{{\phi }_s}^T{\mathit{u}}_{{\phi }_s}+{\mathit{D}}_{t_s}^T{\mathit{u}}_{t_s}\right)=R_I\left(\mathit{p}\right)=0\to \mathit{R}\left(\mathit{p}\right)=\mathit{A}\left(\mathit{p}\right)\mathit{p}-\mathit{b}\left(\mathit{p}\right)\end{array}(11)$

Even though the analytical forms of all the quantities involved in the Reynolds equation are known, from a numerical point of view, it is not convenient to use the Newton iterative method to solve the system Eq. 11 as done in Ruggiero and Sicilia (2020b), because the related Jacobian matrix in this case is a full matrix due to the presence of all the non-zero elements belonging to the influence matrix $\mathit{C}$ assembled by the finite element model, and since matrix $\mathit{A}$ is a band matrix, the more convenient under-relaxed iterative scheme reported in Eq. 12 is used, where $\lambda$ is the relaxation coefficient.

$\begin{array}{c}{\mathit{p}}^{\left(k+1\right)}={\mathit{p}}^{\left(k\right)}+\lambda \left({\mathit{A}}^{-1}\left({\mathit{p}}^{\left(k\right)}\right)\mathit{b}\left({\mathit{p}}^{\left(k\right)}\right)-{\mathit{p}}^{\left(k\right)}\right)\end{array}(12)$

The working algorithm logic is shown in Figure 3 and is supplied by the input data of eccentricity $\mathit{e}\left(t\right)$ and angular velocity $\mathit{\omega }\left(t\right)$ vectors' time trends and elaborates in correspondence with each time instant the pressure field $p\left(\theta ,\phi ,t\right)$; then, the new linear wear rate is integrated over time $t$ and is added to the next surfaces’ separation, in order to consider the progressive geometry modification being updated due to material loss.

FIGURE 3

FIGURE 3. Elasto-hydrodynamic lubrication model flowchart.

Once the lubrication algorithm has reached convergence over all the analyzed time periods, the interesting integral quantities of load, friction torque, and wear volume can be evaluated. The load vector over time $\mathit{N}\left(t\right)$ is obtained by integrating the lubricant pressure over the socket's internal surface area, as written in Eq. 13.

$\begin{array}{c}\mathit{N}\left(t\right)={\displaystyle \underset{0}{\overset{\pi }{\int }}{\displaystyle \underset{0}{\overset{\pi }{\int }}p\widehat{\mathit{r}}{R}^2\sin \theta d\theta d\phi }}\end{array}(13)$

The friction torque vector calculation requires the knowledge of the shear stress fields due to the lubricant acting on the socket internal surface; they are collected in vector $\mathit{\tau }$ depending on the surface pressure gradient $\mathbf{\nabla }p$ and on the sliding velocity vector ${\mathit{v}}_{sl}$. Turning the involved quantities in to spherical coordinates, the two surface shear stress fields ${\tau }_{\theta }$ and ${\tau }_{\phi }$ are given in Eq. 14.

$\begin{array}{c}\left[\begin{array}{c}\begin{array}{c}{\tau }_{\theta }\\ {\tau }_{\phi }\end{array}\\ 0\end{array}\right]=\mathit{\tau }=\frac{h}{2}\mathbf{\nabla }p+\mu \frac{{\mathit{v}}_{sl}}{h}=\left[\begin{array}{c}\begin{array}{c}\frac{h}{2R}\frac{\partial p}{\partial \theta }+\mu \frac{\left(-2U_{\theta }\right)}{h}\\ \frac{h}{2R\sin \theta }\frac{\partial p}{\partial \phi }+\mu \frac{\left(-2U_{\phi }\right)}{h}\end{array}\\ 0\end{array}\right]\end{array}(14)$

Then, the friction torque vector $\mathit{T}\left(t\right)$ related to the Cartesian reference frame is obtained in Eq. 15 by integrating the torque generated by the shear stress fields over the socket's internal surface, considering the transformation produced by the application of the spherical rotation matrix ${\mathit{R}}_s$.

$\begin{array}{c}\mathit{T}\left(t\right)={\displaystyle \underset{0}{\overset{\pi }{\int }}{\displaystyle \underset{0}{\overset{\pi }{\int }}R\widehat{\mathit{r}}\times {\mathit{R}}_s\mathit{\tau }{R}^2\sin \theta d\theta d\phi }}\end{array}(15)$

Finally, the wear volume loss over time $V_w\left(t\right)$ is given in Eq. 16 by the integration over the socket's internal surface of the cumulated penetration wear depth field $u_w$.

$\begin{array}{c}{V}_w\left(t\right)={\displaystyle \underset{0}{\overset{\pi }{\int }}{\displaystyle \underset{0}{\overset{\pi }{\int }}{u}_w{R}^2\sin \theta d\theta d\phi }}\end{array}(16)$

As stated during the finite element deformation model description, additional interesting information about the entire tribo-system body is achieved by returning to the mechanical equilibrium reported in Eq. 2, that is, the solution of the elasto-hydrodynamic lubrication problem provides the time evolution of the lubricant pressure vector $\mathit{p}$, which is used to evaluate the pressure nodal forces vector ${\mathbf{\Phi }}_p$ acting on the nodes belonging to the ball or socket structures; therefore, solving Eq. 2 for the two bodies in the unknown nodal displacement vector $\mathit{q}$, the stress vector $\mathit{\sigma }$ and strain vector $\mathit{\epsilon }$ for each tetrahedron element are evaluated by introducing the elemental deformation matrix $\mathit{B}$ and elasticity matrix $\mathit{E}$. Considering the nodal displacement vector associated with the $i$ tetrahedral element ${\mathit{q}}_i$, the associated stress and strain vectors are given in Eq. 17.

$\begin{array}{c}\begin{array}{cc}\left[\begin{array}{c}\begin{array}{c}{\epsilon }_{x_i}\\ {\epsilon }_{y_i}\\ {\epsilon }_{z_i}\end{array}\\ {\gamma }_{x{y}_i}\\ {\gamma }_{y{z}_i}\\ {\gamma }_{x{z}_i}\end{array}\right]={\mathit{\epsilon }}_i={\mathit{B}}_i{\mathit{q}}_i& \left[\begin{array}{c}\begin{array}{c}{\sigma }_{x_i}\\ {\sigma }_{y_i}\\ {\sigma }_{z_i}\end{array}\\ {\tau }_{x{y}_i}\\ {\tau }_{y{z}_i}\\ {\tau }_{x{z}_i}\end{array}\right]={\mathit{\sigma }}_i={\mathit{E}}_i{\mathit{\epsilon }}_i={\mathit{E}}_i{\mathit{B}}_i{\mathit{q}}_i\end{array}\end{array}(17)$

The two fields over time evaluated in this work are the Von Mises equivalent stress ${\sigma }_{VM}$ and strain energy $U$, with reference to Eq. 18, where the first is calculated by the principal stress coming from the diagonalization of the rearranged stress tensor, while the second is given as the scalar product between the stress and strain vectors and by the tetrahedron volume $V$.

$\begin{array}{c}\begin{array}{c}\left[\begin{array}{ccc}{\sigma }_{x_i}& {\tau }_{x{y}_i}& {\tau }_{x{z}_i}\\ {\tau }_{x{y}_i}& {\sigma }_{y_i}& {\tau }_{y{z}_i}\\ {\tau }_{x{z}_i}& {\tau }_{y{z}_i}& {\sigma }_{z_i}\end{array}\right]\to \left[\begin{array}{ccc}{\sigma }_{I_i}& 0& 0\\ 0& {\sigma }_{I{I}_i}& 0\\ 0& 0& {\sigma }_{II{I}_i}\end{array}\right]\to {\sigma }_{V{M}_i}=\frac{\sqrt{{\left({\sigma }_{I_i}-{\sigma }_{I{I}_i}\right)}^2+{\left({\sigma }_{I{I}_i}-{\sigma }_{II{I}_i}\right)}^2+{\left({\sigma }_{I_i}-{\sigma }_{II{I}_i}\right)}^2}}{\sqrt{2}}\\ U_i=\frac{1}{2}\left({\mathit{\sigma }}_i^T{\mathit{\epsilon }}_i\right)V_i\end{array}\end{array}(18)$

The whole merged model flowchart is shown in Figure 4, in which all the input and output quantities are highlighted.

FIGURE 4

FIGURE 4. Input/output scheme of the developed model.

Results and Discussions

The developed elasto-hydrodynamic lubrication simulator was tested for a typical biotribological system considered in the human total joint replacements: a Ceramic ball in an ultra-high-molecular-weight-polyethylene (UHMWPE) socket. The joint kinematics was selected, without loss of generality, as a planar motion (Figure 5) in the $xy$ plane, in which only the $z$-component of the relative angular velocity vector $\mathit{\omega }$ is different from zero and equal to 300 rad/s, while only the $y$-component of the eccentricity vector $\mathit{e}$ goes from zero to 1.2 times the radial clearance $c$ during 1 s ($\mathit{n}=\mathit{e}/c$ is the eccentricity vector in dimensionless form). The other input data are reported in Table 1. For this type of simulation, a processor Intel(R) Xeon(R) Gold 5218 CPU @ 2.30 GHz was used to analyze a 100 × 100 lubrication grid mesh in addition to the structures discretization in 8,005 elements for the socket and 13,937 elements for the ball along a time vector made of 40 steps.

FIGURE 5

FIGURE 5. Spherical joint kinematics input: eccentricity (A); angular velocity (B).

TABLE 1

TABLE 1. Parameters input data.

The tribological quantities' output are shown in their following dimensionless form, by using the maximum Hertzian pressure and shear stress $p_{Hz}$, ${\tau }_{Hz}$ (defined in Eq. 19 through the vertical load $N_y$, the equivalent curvature radius $R_x$, and the equivalent Young modulus $E^{\ast }$) and the radial clearance $c$ as reference parameters.

$\begin{array}{c}\{\begin{array}{c}\frac{1}{R_x}=\frac{1}{r}-\frac{1}{R}\\ \frac{1}{E^{\ast }}=\frac{1-{\nu }_b^2}{E_b}+\frac{1-{\nu }_s^2}{E_s}\end{array}\to a=0.908\sqrt[3]{\frac{N_y{R}_x}{E^{\ast }}}\to \begin{array}{cc}{p}_{Hz}=\frac{3}{2}\frac{N_y}{\pi a^2}& {\tau }_{Hz}=\frac{1}{3}{p}_{Hz}\end{array}\end{array}(19)$

The calculated fluid pressure $p$ and the surface separation $h$ fields are shown in correspondence of the last time instant in Figure 6: the squeeze action, due to the $y$-component of the eccentricity vector $\mathit{e}$, causes the pressure bell-shape within the socket surface and the deformation needed to keep the surface separated. The $z$-component of the angular velocity vector $\mathit{\omega }$ is responsible for the sharp pressure growth toward the outlet of the contact zone, characterized by minimum separation. The time evolution of the two described fields is available in Supplementary Video S1.

FIGURE 6

FIGURE 6. Lubricant pressure (A) and surface separation (B) fields.

The classical elasto-hydrodynamic profiles of the pressure and separation are better visualized in Figure 7, showing the complete time evolution of the same quantities depicted in Figure 6; here, they are reported in correspondence with the pressure peak position, in which each line is characterized by proportional transparency going from 100% to 0% along with time. Due to the only relative rotational motion around the $z$-axis, the pressure spike and minimum film thickness are established along the angle $\phi$-direction, while the angle $\theta$-direction is characterized by symmetrical shapes. The smooth time evolution of the related graph can be visualized in Supplementary Video S2.

FIGURE 7

FIGURE 7. Lubricant pressure and surface separation profiles evolution along $\theta$ (A) and $\phi$ (B).

During the analyzed time instant, the fields of the surface radial deformation $\delta$ and the wear depth $u_w$ are reported in Figure 8. Since the deformation is the response of the spherical joint to the pressure action, in this case, the surface deformation field absorbs the pressure spike, providing an almost symmetrical bell shape, while the wear depth field is affected by the intensity of the pressure spike which produces material loss cumulating over time in the outlet zone and resulting in a nonsymmetrical shape. All the related time steps can be viewed in Supplementary Video S3.

FIGURE 8

FIGURE 8. Surface radial deformation (A) and penetration wear depth (B) fields.

The shear stress fields are shown in Figure 9 and in Supplementary Video S4 over the socket's internal surface domain, in order to analyze their action in terms of friction; the shear stress acting along the $\theta$-direction, i.e., ${\tau }_{\theta }$, is characterized by an anti-symmetrical behavior with respect to the $xy$ plane, and it does not produce appreciable friction resistance to the imposed rotational motion around the $z$-axis; the largest contribution to the resistance against the rotational motion is given by the high negative values reached by the ${\tau }_{\phi }$ field.

FIGURE 9

FIGURE 9. Shear stress fields over the socket spherical domain along $\theta$ (A) and $\phi$ (B) directions.

The integration over the socket surface of the pressure and of the shear stress led to the time evolution of load $\mathit{N}$ and of the friction torque $\mathit{T}$ vectors shown in Figure 10; as expected, the combined action of the eccentricity $e_y$ and of the angular velocity ${\omega }_z$ produces a loading force lying in the $xy$ contact plane, with the $y$-component rapidly growing, while the only appreciable friction action is generated around the $z$-axis with a torque initially constant and subsequently quickly increasing as the contact area increases its dimensions.

FIGURE 10

FIGURE 10. Load (A) and friction torque (B) time evolution.

The surface integration over the analyzed domain of the penetration wear depth shown in the Figure 8 causes the wear volume loss time evolution reported in Figure 11.

FIGURE 11

FIGURE 11. Wear volume time evolution.

The quantities describing the mechanical behavior of the spherical joint with respect to the analyzed kinematics can be viewed both for the entire body and for its internal structure. In particular, the Von Mises stress and the strain energy fields are depicted in Figure 12 in correspondence to the last time instant for both the ball and socket, by viewing the contact interface (the complete time evolution of the visualized quantities is available in Supplementary Video S5); on the contacting surfaces, the socket is characterized by higher values of both Von Mises stress and strain energy than the ball because of its lower stiffness; in general, the ball is more loaded on its back (negative $y$-coordinate) because of the imposed fixed nodes.

FIGURE 12

FIGURE 12. Von Mises stress [socket (A); ball (B)] and strain energy fields [socket (C); ball (D)].

The same quantities referring to the internal structure of the spherical joint are shown, respectively, in Figures 13, 14, highlighting their propagation inside the bodies in the $xy$ plane at the $z$-coordinate equal to zero (the time evolution is available in Supplementary Videos S6, 7). In both cases, the action of the lubricant pressure is clearly observable: the two evident features are the ball’s capacity to better distribute contact stress inside the body and the higher load in correspondence to its fixation.

FIGURE 13

FIGURE 13. Von Mises stress field inside of the spherical joint.

FIGURE 14

FIGURE 14. Strain energy field inside of the spherical joint.

The reliability of the model was discussed by comparing the minimum fluid meatus thickness $h_{\min }$ reached within the joint surfaces with the one calculated through the well-known Hamrock–Dowson formula (Hamrock and Dowson, 1976; Lubrecht, et al., 2009) reported in Eq. 20.

$\begin{array}{c}\frac{h_{\min }}{R_x}=3.63{\left(\frac{\overline{\mu }U}{2E^{\ast }{R}_x}\right)}^{0.68}{\left(2{\alpha }_{\mu }{E}^{\ast }\right)}^{0.49}{\left(\frac{N_y}{2R_x^2{E}^{\ast }}\right)}^{-0.073}\left(1-e^{-0.68k}\right)\end{array}(20)$

where the viscosity surface integral mean $\overline{\mu }$ is used and the other parameters are defined in Eq. 21.

$\begin{array}{c}\begin{array}{ccc}U=\frac{{\omega }_{z}R}{2}& \frac{1}{R_x}=\frac{1}{R_y}& k=\mathrm{1,0339}{\left(\frac{R_y}{R_x}\right)}^{0.636}\end{array}\end{array}(21)$

The comparison is depicted in Figure 15, where the minimum fluid meatus thickness is plotted against the vertical load $N_y$ in their dimensionless form. Obviously, in the first part of the graph, the two approaches do not match because the load is too low, which does not produce enough amount of deformation; then, in correspondence to higher loads, the pressure spike takes place within the surfaces, such that the minimum thickness simulated by the proposed model overlaps on the one predicted by the Hamrock–Dowson formula, demonstrating the reliability of the developed algorithm.

FIGURE 15

FIGURE 15. Minimum film thickness comparison with the Hamrock–Dowson model.

Conclusion

The objective of this work was to integrate a finite element deformation model already developed by the authors into an elasto-hydrodynamic lubrication algorithm, in order to characterize the deformation of the surfaces belonging to the analyzed tribo-system through a complex and accurate fluid–solid interaction tool.

The mechanical and tribological behavior of the ball in socket joint was discussed and described through the analytical theories coming from the Reynolds equation and the finite element method in the case of linear elastic materials.

In particular, the surface deformation model referred to the spherical geometry of the ball and socket bodies and, starting from the mechanical equilibrium which links the nodal displacements to the nodal forces through the stiffness matrix, led to the definition of the influence matrix which directly related the total surfaces’ radial deformation to the pressure established within the interface between the ball and the socket; following the proposed approach, the interfacial influence matrix contains information about the mechanical and geometrical properties of the whole body comprising the spherical joint and, furthermore, the finite element model allows analyses of the tensional state of the global structure once the lubrication algorithm has elaborated the interfacial lubricant pressure.

The evaluated deformation directly communicates to the summation which composes the surfaces’ separation that is needed to solve the Reynolds equation referred to the elasto-hydrodynamic mode. The lubricant rheology properties of density and dynamical viscosity are defined as pressure functions and the viscosity, moreover, is considered as a function of the lubricant shear rate in order to simulate its non-Newtonian shear-thinning or -thickening behavior. The surface separation geometry varies over time even because its composition accounts for a loss contribution equal to the penetration wear depth, which is updated over each time instant. The resulting highly nonlinear partial differential equation is discretized over the socket’s internal surface through the finite differences and solved by an iterative technique. The lubricant pressure solution due to the relative motion input in terms of eccentricity and angular velocity of the ball with respect to the socket is used to evaluate the shear stress acting on the socket surface, then the integral quantities of load, friction torque, and wear volume are elaborated. Thanks to the used finite element deformation model, the Von Mises stress and the strain energy related to the spherical joint bulk are computed by the knowledge of the nodal displacements.

The developed algorithm related to the integration of the two described models is written in a MATLAB environment, and it is used to analyze the elaborated output quantities referred to as a ceramic-in-UHMWPE spherical joint during a basic planar motion. The results computed by the model reproduced the classical elasto-hydrodynamic shapes referring to the lubricant pressure and to the surface separation, which presented the pressure spike in correspondence to the minimum film thickness and provided the time evolution of the load lying in the relative motion plane and of the friction torque consistent with the allowed relative rotation, showing a good computational behavior, demonstrated by the comparison with the Hamrock–Dowson model. Furthermore, the propagation of the lubricant pressure action along with the ball and socket both superficial and sub-superficial structures was visualized in terms of Von Mises stress and strain energy, providing information about the modality to absorb the stress of the materials involved in the construction of the spherical joint.

Due to the high information content provided by the developed model, it could be exploited in the framework of the tribological characterization or design of lubricated joints for both mechanical and biomechanical purposes (e.g., the total hip replacements). From a purely mechanical point of view, a Stribeck curve describing the friction coefficient in correspondence with a chosen kinematics can be depicted, in order to, for example, define a map that can compute the friction torque due to the ball and socket relative motion and loading configurations, which could be very useful within the multibody simulation algorithms. In order to reach more accurate results, the main improvements in the authors’ perspective could regard

1) the adaptation of the entire model to the mixed lubrication regime, dividing the socket's internal surface into lubricated areas and direct contact zones;

2) the introduction of a damping matrix within the finite element model, in order to consider the possibility of having a viscoelastic mechanical behavior;

3) the coupling of the Reynolds equation to the energy conservation equation, such that the computed interfacial temperature field could be used to define a further contribution to the nodal forces due to the consequent thermal expansion of the finite elements.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

AR: conceptualization, formal analysis, methodology, supervision, reviewing and editing, and funding. AS: formal analysis, programming, calculations, writing—original draft, and writing—reviewing, and editing.

Funding

This research was funded by MIUR, PRIN 2017 BIONIC.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmech.2022.909156/full#supplementary-material

References

Askari, E., and Andersen, M. S. (2018). A Closed-form Formulation for the Conformal Articulation of Metal-On-Polyethylene Hip Prostheses: Contact Mechanics and Sliding Distance. Proc. Inst. Mech. Eng. H. 232 (12), 1196–1208. doi:10.1177/0954411918810044

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, H. P., and Hughes, T. G. (2000). Evaluation of Deflection in Semi-infinite Bodies by a Differential Method. Proc. Institution Mech. Eng. Part C J. Mech. Eng. Sci. 214 (4), 563–584. doi:10.1243/0954406001523911

CrossRef Full Text | Google Scholar

Fregly, B. J., Bei, Y., and Sylvester, M. E. (2003). Experimental Evaluation of an Elastic Foundation Model to Predict Contact Pressures in Knee Replacements. J. biomechanics 2003 (36), 1659–1668. doi:10.1016/s0021-9290(03)00176-3

CrossRef Full Text | Google Scholar

Gao, L., Dowson, D., and Hewson, R. W. (2016). A Numerical Study of Non-newtonian Transient Elastohydrodynamic Lubrication of Metal-On-Metal Hip Prostheses. Tribol. Int. 93, 486–494. doi:10.1016/j.triboint.2015.03.003

CrossRef Full Text | Google Scholar

Gao, L., Dowson, D., and Hewson, R. W. (2017). Predictive Wear Modeling of the Articulating Metal-On-Metal Hip Replacements. J. Biomed. Mat. Res. 105 (3), 497–506. doi:10.1002/jbm.b.33568

CrossRef Full Text | Google Scholar

Gao, L., Wang, F., Yang, P., and Jin, Z. (2009). Effect of 3D Physiological Loading and Motion on Elastohydrodynamic Lubrication of Metal-On-Metal Total Hip Replacements. Med. Eng. Phys. 31 (6), 720–729. doi:10.1016/j.medengphy.2009.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Halloran, J. P., Easley, S. K., Petrella, A. J., and Rullkoetter, P. J. (2005). Comparison of Deformable and Elastic Foundation Finite Element Simulations for Predicting Knee Replacement Mechanics. J. Biomech. Eng. 127, 813–818. doi:10.1115/1.1992522

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamrock, B. J., and Dowson, D. (1976). Isothermal Elastohydrodynamic Lubrication of Point Contacts: Part III – Fully Flooded Results. Cleveland, OH: National Aeronautics and Space Administration Technical Note.

Google Scholar

Jagatia, M., and Jin, Z. M. (2001). Elastohydrodynamic Lubrication Analysis of Metal-On-Metal Hip Prostheses under Steady State Entraining Motion. Proc. Inst. Mech. Eng. H. 215 (6), 531–541. doi:10.1243/0954411011536136

PubMed Abstract | CrossRef Full Text | Google Scholar

Jalali-Vahid, D., Jagatia, M., Jin, Z. M., and Dowson, D. (2001). Prediction of Lubricating Film Thickness in a Ball-In-Socket Model with a Soft Lining Representing Human Natural and Artificial Hip Joints. Proc. Institution Mech. Eng. Part J J. Eng. Tribol. 215 (4), 363–372. doi:10.1243/1350650011543600

CrossRef Full Text | Google Scholar

Jalali-Vahid, D., Jin, Z. M., and Dowson, D. (2003). Elastohydrodynamic Lubrication Analysis of Hip Implants with Ultra High Molecular Weight Polyethylene Cups under Transient Conditions. Proc. Institution Mech. Eng. Part C J. Mech. Eng. Sci. 217 (7), 767–777. doi:10.1243/095440603767764417

CrossRef Full Text | Google Scholar

Jalali-Vahid, D., and Jin, Z. M. (2001). Transient Elastohydrodynamic Lubrication Analysis of Ultra-high Molecular Weight Polyethylene Hip Joint Replacements. Proc. Institution Mech. Eng. Part C J. Mech. Eng. Sci. 216 (4), 409–420. doi:10.1243/0954406021525205

CrossRef Full Text | Google Scholar

Lubrecht, A. A., Venner, C. H., and Colin, F. (2009). Film Thickness Calculation in Elasto-Hydrodynamic Lubricated Line and Elliptical Contacts: The Dowson, Higginson, Hamrock Contribution. Proc. Institution Mech. Eng. Part J J. Eng. Tribol. 223 (3), 511–515. doi:10.1243/13506501jet508

CrossRef Full Text | Google Scholar

Mattei, L., Di Puccio, F., Piccigallo, B., and Ciulli, E. (2011). Lubrication and Wear Modelling of Artificial Hip Joints: A Review. Tribol. Int. 44 (5), 532–549. doi:10.1016/j.triboint.2010.06.010

CrossRef Full Text | Google Scholar

Mukras, S. M., Kim, N. H., Schmitz, T. L., and Sawyer, W. G. (2010). Evaluation of Contact Force and Elastic Foundation Models for Wear Analysis of Multibody Systems. Int. Des. Eng. Tech. Conf. Comput. Inf. Eng. Conf. 44106, 1565–1577. doi:10.1115/detc2010-28750

CrossRef Full Text | Google Scholar

Pèrez-Gonzàlez, A., Fenollosa-Esteve, C., Sancho-Bru, J. L., Sànchez-Marìn, F. T., and Vergara, M. a. R.-C. P. J. (2008). A Modified Elastic Foundation Contact Model for Application in 3D Models of the Prosthetic Knee. Med. Eng. Phys. 30 (3), 387–398.

Google Scholar

Ruggiero, A., D’Amato, R., and Affatato, S. (2019). Comparison of Meshing Strategies in THR Finite Element Modelling. Materials 12 (14), 2332. doi:10.3390/ma12142332

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruggiero, A., Merola, M., and Affatato, S. (2018). Finite Element Simulations of Hard-On-Soft Hip Joint Prosthesis Accounting for Dynamic Loads Calculated from a Musculoskeletal Model during Walking. Materials 11 (4), 574. doi:10.3390/ma11040574

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruggiero, A., and Sicilia, A. (2020b). A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints. Lubricants 8 (7), 72. doi:10.3390/lubricants8070072

CrossRef Full Text | Google Scholar

Ruggiero, A., and Sicilia, A. (2020c). A Novel Explicit Analytical Multibody Approach for the Analysis of Upper Limb Dynamics and Joint Reactions Calculation Considering Muscle Wrapping. Appl. Sci. 10 (21), 7760. doi:10.3390/app10217760

CrossRef Full Text | Google Scholar

Ruggiero, A., Sicilia, A., and Affatato, S. (2020). In Silico total Hip Replacement Wear Testing in the Framework of ISO 14242-3 Accounting for Mixed Elasto-Hydrodynamic Lubrication Effects. Wear 460, 203420. doi:10.1016/j.wear.2020.203420

CrossRef Full Text | Google Scholar

Ruggiero, A., and Sicilia, A. (2020a). Lubrication Modeling and Wear Calculation in Artificial Hip Joint during the Gait. Tribol. Int. 142, 105993. doi:10.1016/j.triboint.2019.105993

CrossRef Full Text | Google Scholar

Ruggiero, A., and Sicilia, A. (2021). Mathematical Development of a Novel Discrete Hip Deformation Algorithm for the In Silico Elasto-Hydrodynamic Lubrication Modelling of Total Hip Replacements. Lubricants 9 (4), 41. doi:10.3390/lubricants9040041

CrossRef Full Text | Google Scholar

Shettar, B. M., Hiremath, P., and Bujurke, N. (2018). A Novel Numerical Scheme for the Analysis of Effects of Surface Roughness on EHL Line Contact with Couple Stress Fluid as Lubricant. Sadhana 43 (8), 1–12. doi:10.1007/s12046-018-0906-y

CrossRef Full Text | Google Scholar

Srivastava, G., Christian, N., and Fred Higgs, C. (2021). A Predictive Framework of the Tribological Impact of Physical Activities on Metal-On-Plastic Hip Implants. Biotribology 25, 100156. doi:10.1016/j.biotri.2020.100156

CrossRef Full Text | Google Scholar

van Leeuwen, H. (2009). The Determination of the Pressure-Viscosity Coefficient of a Lubricant through an Accurate Film Thickness Formula and Accurate Film Thickness Measurements. Proc. Institution Mech. Eng. Part J J. Eng. Tribol. 223 (8), 1143–1163. doi:10.1243/13506501jet504

CrossRef Full Text | Google Scholar

Wang, F. C., Jin, Z. M., McEwen, H. M. J., and Fisher, J. (2003). Microscopic Asperity Contact and Deformation of Ultrahigh Molecular Weight Polyethylene Bearing Surfaces. Proc. Inst. Mech. Eng. H. 217 (6), 477–490. doi:10.1243/09544110360729117

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F. C., and Jin, Z. M. (2004a). Prediction of Elastic Deformation of Acetabular Cups and Femoral Heads for Lubrication Analysis of Artificial Hip Joints. Proc. Institution Mech. Eng. Part J J. Eng. Tribol. 218 (3), 201–209. doi:10.1243/1350650041323331

CrossRef Full Text | Google Scholar

Wang, F., and Jin, Z. (2004b). Lubrication Modelling of Artificial Hip Joints: From Fluid Film to Boundary Lubrication Regimes. Eng. Syst. Des. Analysis 417731, 605–611. doi:10.1115/esda2004-58077

CrossRef Full Text | Google Scholar

Wang, F., and Jin, Z. (2008). Transient Elastohydrodynamic Lubrication of Hip Joint Implants. J. Tribol. 130, 6200. doi:10.1115/1.2806200

CrossRef Full Text | Google Scholar

Wang, W.-Z., Wang, F.-C., Jin, Z.-M., Dowson, D., and Hu, Y.-Z. (2009). Numerical Lubrication Simulation of Metal-On-Metal Artificial Hip Joint Replacements: Ball-In-Socket Model and Ball-On-Plane Model. Proc. Institution Mech. Eng. Part J J. Eng. Tribol. 223 (7), 1073–1082. doi:10.1243/13506501jet581

CrossRef Full Text | Google Scholar

Nomenclature

Abbreviations

$r$ Ball radius

$\mathrm{R}$ Socket radius, Reynolds equation function vector

$H$ Socket width

$x$ Cartesian position vector

$\rho ,\theta ,\phi$ Spherical coordinates

$\delta$ Deformation vector

$p$ Pressure vector Lubricant pressure

$K$ Finite Element stiffness matrix

$q$ Nodal displacement vector

$\Phi$ Nodal force vector

$C$ Influence matrix

$J_{bs},J_{sb}$ Cubic interpolation matrices

$\rho$ Lubricant density

$\mu$ Lubricant dynamical viscosity

$h$ Film thickness

$v$ Entraining velocity vector

$t$ Time

$p$ Pressure vector Lubricant pressure

$v_{sl}$ Sliding velocity

$h_g$ Surfaces’ geometrical approach

$c$ Radial clearance

$e$ Eccentricity vector

$\widehat{r}$ Radial unit vector

$k_w$ Wear factor

$u_w$ Penetration wear depth

$\omega$ Angular velocity vector

$R_s$ Spherical rotation matrix

$N$ Load vector

$\tau$ Shear stress vector

$T$ Friction torque vector

$V_w$ Wear volume

$\sigma$ Stress vector

$\epsilon$ Strain vector

${\sigma }_{VM},U$ Von Mises stress and strain energy

$p_{Hz},{\tau }_{Hz}$ Hertzian pressure and shear stress