Coupling multi-body dynamics and fluid dynamics to model lubricated spherical joints

Abstract

A new approach of coupling multibody dynamics and fluid dynamics is developed to model hydrodynamic lubrication of spherical clearance joints with thin fluid film and relative multidirectional motion. The model accounts for dynamics motion of articulating components as well as both squeeze- and wedge-film actions of the synovial fluid. Multibody dynamics methodology is employed to derive the motion equations and Reynolds equation governs the fluid dynamics. The finite difference method is utilized to discretize the governing equation of lubricant and the multi-grid method augments computational efficiency to acquire outcomes employing a Gauss–Seidel relaxation scheme. Fluid–structure interaction is incorporated into the methodology using a partitioned formulation embedded in a high-order Runge–Kutta time integrators for integrating the nonlinear equations of the coupled system over time of interest. A demonstrative example of total hip arthroplasty is considered and the developed model is assessed against outcomes available in the literature. The effect of initial conditions on the pressure, film thickness and dynamics of the lubricated spherical joint is analyzed and discussed. It is illustrated that maximum fluid pressure is undergone by the hip implant at the first walking cycle of movement due to an unstable state, which is strongly dependent upon the initial condition. Finally, the approach presented in this research work is a robust dynamic model to study hydrodynamic lubrication of spherical joints.

Appendices

Appendix A: Matrices in Eq. (6)

$$\begin{aligned} C_{1}^{i,j}= & {} \left( {\frac{\varepsilon _{i-0.5,j} }{\Delta \varphi ^{2}}} \right) \end{aligned}$$

(A1)

$$\begin{aligned} C_{2}^{i,j}= & {} -\frac{\varepsilon _{i-0.5,j} +\varepsilon _{i+0.5,j} }{\Delta \varphi ^{2}}\nonumber \\&\quad -\frac{\varepsilon _{i,j-0.5} sin\left( {\theta _{j-0.5} } \right) +\varepsilon _{i,j+0.5} sin\left( {\theta _{j+0.5} } \right) }{\Delta \theta ^{2}} \end{aligned}$$

(A2)

$$\begin{aligned} C_{3}^{i,j}= & {} \left( {\frac{\varepsilon _{i+0.5,j} }{\Delta \varphi ^{2}}} \right) \end{aligned}$$

(A3)

$$\begin{aligned} C_{4}^{i,j}= & {} sin\theta _{j} \left( {\frac{\varepsilon _{i,j-0.5} sin\left( {\theta _{j-0.5} } \right) }{\Delta \theta ^{2}}} \right) \end{aligned}$$

(A4)

$$\begin{aligned} C_{5}^{i,j}= & {} sin\theta _{j} \left( {\frac{\varepsilon _{i,j+0.5} sin\left( {\theta _{j+0.5} } \right) }{\Delta \theta ^{2}}} \right) \end{aligned}$$

(A5)

$$\begin{aligned} C_{6}^{i,j}= & {} -{\Omega }_{x} \left[ {\left( {sin\theta _{j} } \right) ^{2}sin\varphi _{i} \frac{H_{i,j} -H_{i,j-1} }{\Delta \theta }+sin\theta _{j} cos\theta _{j} cos\varphi _{i} \frac{H_{i,j} -H_{i-1,j} }{\Delta \varphi }} \right] \nonumber \\&\quad +{\Omega }_{y} \left[ {\left( {sin\theta _{j} } \right) ^{2}cos\varphi _{i} \frac{H_{i,j} -H_{i,j-1} }{\Delta \theta }-sin\theta _{j} cos\theta _{j} sin\varphi _{i} \frac{H_{i,j} -H_{i-1,j} }{\Delta \varphi }} \right] \nonumber \\&\quad +{\Omega }_{z} \left[ {\left( {sin\theta _{j} } \right) ^{2}\frac{H_{i,j} -H_{i-1,j} }{\Delta \varphi }} \right] \nonumber \\&\quad +2\left( {sin\theta _{j} } \right) ^{2}\frac{\partial H_{i,j} }{\partial T} \end{aligned}$$

(A6)

where

$$\begin{aligned} \varepsilon _{i,j+0.5} =\frac{\varepsilon _{i,j+1} +\varepsilon _{i,j} }{2} \end{aligned}$$

(A7)

$$\begin{aligned} \varepsilon _{i,j-0.5} =\frac{\varepsilon _{i,j-1} +\varepsilon _{i,j} }{2} \end{aligned}$$

(A8)

$$\begin{aligned} \varepsilon _{i+0.5,j} =\frac{\varepsilon _{i+1,j} +\varepsilon _{i,j} }{2} \end{aligned}$$

(A9)

$$\begin{aligned} \varepsilon _{i-0.5,j} =\frac{\varepsilon _{i-1,j} +\varepsilon _{i,j} }{2} \end{aligned}$$

(A10)

$$\begin{aligned} sin\left( {\theta _{j+0.5} } \right) =sin\left( {\frac{\theta _{j+1} +\theta _{j} }{2}} \right) \end{aligned}$$

(A11)

$$\begin{aligned} sin\left( {\theta _{j-0.5} } \right) =sin\left( {\frac{\theta _{j-1} +\theta _{j} }{2}} \right) \end{aligned}$$

(A12)

Appendix B: The multi-grid method

A fine mesh is commonly constructed to compute the lubricant pressure and film thickness using a relaxation method, e.g. Gauss-Seidel iterative method. However, such a solution method has two drawbacks: (1) a very slow convergence of the relaxation process especially for a large number of nodes; and (2) incapability to suppress error components in the solution that have wavelengths larger than the mesh size. The multigrid method allows us to accelerate the convergence of the relaxation process as well as suppress high wavelength errors by making use of coarse meshes within the original fine mesh. For further information about the multigrid method, the reader is referred to the following references [52,53,54]. The matrix form of Eq. (6) for a fine mesh with an even mesh size as $\Delta \theta = \Delta \varphi ={{\bar{h}}}$, is recast as follows:

$$\begin{aligned} L^{{{\bar{h}}}}\left( {p^{{{\bar{h}}}}} \right) =C_{6}^{{{\bar{h}}}} \end{aligned}$$

(B1)

where superscript (${{\bar{h}}})$ is to characterize its mesh size. Equation (B1) is commonly solved using a relaxation method, namely the Gauss-Seidel iterative methodology. Assume after a few sweeps, e.g. 2 or 3, an approximation, ${\tilde{p}}^{\bar{h}}$, to the solution $p^{{{\bar{h}}}}$ is achieved. The corresponding residual is calculated as follows:

$$\begin{aligned} r^{{{\bar{h}}}}=C_{6}^{{{\bar{h}}}} -L\left( {{\tilde{p}}}^{\bar{h}} \right) \end{aligned}$$

(B2)

Assuming an exact solution for Eq. (B1) is $p^{{{\bar{h}}}}$, the operator L is approximated in a proximity of ${\tilde{p}}$ as given by

$$\begin{aligned} L\left( {{\tilde{p}}} \right) =L\left( {{\tilde{p}}} \right) +\left. {\frac{\partial L}{\partial p}} \right| _{{\tilde{p}}} \left( {\bar{p}} -{{\tilde{p}}} \right) \end{aligned}$$

(B3)

As we are interested in the exact solution, ${\bar{p}}$, and know that the operator L is linear, the equation is recast as follows:

$$\begin{aligned} {{\bar{p}}}= & {} {\tilde{p}} + {\left( {{{\left. {\frac{{\partial L}}{{\partial p}}} \right| }_{{\tilde{p}}}}} \right) ^{ - 1}}\left( L\left( {{{\bar{p}}}} \right) \right. \left. - L\left( {{\tilde{p}}} \right) \right) = {\tilde{p}} + {\left( {{{\left. {\frac{{\partial L}}{{\partial p}}} \right| }_{{\tilde{p}}}}} \right) ^{ - 1}}L\left( {{{\bar{p}}} - {\tilde{p}}} \right) \nonumber \\= & {} {\tilde{p}} + {\left( {{{\left. {\frac{{\partial L}}{{\partial p}}} \right| }_{{\tilde{p}}}}} \right) ^{ - 1}}\left( {{C_6} - L\left( {{\tilde{p}}} \right) } \right) \end{aligned}$$

(B4)

in which the error is defined as $v={\bar{p}}-{\tilde{p}}$. After discretizing the cup surface with a uniform mesh of the size ${{\bar{h}}}$, the error and residual can be written on the grid as

$$\begin{aligned} v_{i,j}^{{{\bar{h}}}}= & {} {{\bar{p}}}_{i,j}^{{{\bar{h}}}} - {\tilde{p}}_{i,j}^{{{\bar{h}}}} \end{aligned}$$

(B5)

$$\begin{aligned} r_{i,j}^{{{\bar{h}}}}= & {} {\left( {C_6^{{{\bar{h}}}} - L{{\left( {{\tilde{p}}} \right) }^{{{\bar{h}}}}}} \right) _{i,j}} \end{aligned}$$

(B6)

Substituting Eqs. (B5) and (B6) into Eq. (B4), the following equation can be written for the error

$$\begin{aligned} {L^{{{\bar{h}}}}}\left( {v_{i,j}^{{{\bar{h}}}}} \right) = r_{i,j}^{{{\bar{h}}}} \end{aligned}$$

(B7)

The right side of Eq. (B7) is known, so it would be solved for the error $v_{i,j}^{{{\bar{h}}}}$ on the mesh grid. Using the Multigrid method enables us to solve Eq. (B7) for the error by switching the mesh from a fine mesh with the element size of ${{\bar{h}}}$ to a coarse mesh with the element size of 2${{\bar{h}}}$ that is called ${\bar{H}}$. Numbers of nodes are converted using the following relations, $i=2I-1$ and $j=2J-1$, in order to obtaining the $r_{I,J}^{{\bar{H}}}$, $(C_{6}^{{\bar{H}}} )_{I,J}$ and finally $(L\left( {{\tilde{p}}} \right) ^{{\bar{H}}})_{I,J}$ to rewrite Eq. (B7) in the coarse mesh grid, Eq. (B8).

$$\begin{aligned} L\left( v \right) _{I,J}^{{{\bar{H}}}} = r_{I,J}^{{{\bar{H}}}} \end{aligned}$$

(B8)

Equation (B8) is solved using the Gauss–Seidel iterative method until it converges and obtains the solution for the error, ${\bar{v}}^{{\bar{H}}}$. By this time, the error is written in the fine mesh grid using the following relations

$$\begin{aligned} {{{\bar{v}}}^{{{\bar{h}}}}}\left( {2I - 1,2J - 1} \right)= & {} {{{\bar{v}}}^{{{\bar{H}}}}}\left( {I,J} \right) \nonumber \\ {{{\bar{v}}}^{{{\bar{h}}}}}\left( {2I,2J - 1} \right)= & {} \frac{{{{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I,J} \right) + {{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I + 1,J} \right) }}{2} \nonumber \\ {{{\bar{v}}}^{{{\bar{h}}}}}\left( {2I - 1,2J} \right)= & {} \frac{{{{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I,J} \right) + {{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I,J + 1} \right) }}{2} \nonumber \\ {{{\bar{v}}}^{{{\bar{h}}}}}\left( {2I,2J} \right)= & {} \frac{{{{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I,J} \right) + {{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I + 1,J} \right) + {{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I,J + 1} \right) + {{{{\bar{v}}}}^{{{\bar{H}}}}}\left( {I + 1,J + 1} \right) }}{4} \end{aligned}$$

(B9)

Knowing ${\bar{v}}_{i,j}^{\bar{h}}$, the obtained error is added to the approximation for the solution, Eq. (B10). This process is repeated until the solution for Eq. (B1) converges.

$$\begin{aligned} {{\bar{p}}}_{i,j}^{new} = {{\bar{p}}}_{i,j}^{{{\bar{h}}}} + {{\bar{v}}}_{i,j}^{{{\bar{h}}}} \end{aligned}$$

(B10)