Keywords

13.1 Introduction

Wavelets are recognized as a powerful mathematical tool in signal and image processing, time-series analysis, geophysics, computer graphics, etc. The last few years have witnessed tremendous amount of interest and activity in the application of the wavelet theory and its associated multi-resolution analysis to different areas of science and technology. Lately, an area in which wavelets are gaining currency is numerical simulations of ordinary and partial differential equations. What makes them particularly appealing to the solution of partial differential equations is their hierarchical structure and localization property, coupled with their time and scale invariance. Wavelet-multigrid method is applied to solve modified Reynolds equation arising in biomedical engineering, an illustration of the advantages of wavelet methods over traditional methods. The unique feature of this method is combining classical multigrid scheme with the modern wavelet theory in such a way that each benefits from the other.

We turn our attention to synovial joints [1] of the human body. Recently, the research on synovial joints has focused on two aspects: To investigate the fundamental lubrication processes occurring in the natural joints of the human body and the development of artificial, replacement joints based on theoretical results. Mathematical models of human joints serve to predict quantities which are difficult to measure experimentally and to simulate changes to the physiological conditions. With particular reference to the human knee joint, many researchers are involved in modeling of the bio-mechanism of the joint lubrication, proposing sophisticated models that involve complex numerical computation in order to obtain the desired solutions.

The main aim of this chapter is to propose an analytical approximate squeeze film lubrication model of the human knee joint for the quick assessment of the synovial fluid pressure and the associated load carrying capacity. The proposed model is based on a modified Reynolds equation; the solution of which gives the fluid pressure and consequently the load-carrying capacity. Normal synovial fluid is viscous, and its functions are nutrition and lubrication of cartilage, load bearing and shock absorption, ensuring efficient lubrication and proper functioning of the synovial joints.

The modified Reynolds equation which incorporates randomized surface roughness structure as well as the elastic nature of the articular cartilage with viscous and non-Newtonian synovial fluid as lubricant is derived. A simplified mathematical model has been developed for understanding the combined effects of surface roughness, poroelasticity, and lubrication aspects of viscous and non-Newtonian fluid of human knee joint.

Retrospectively, the simplest and the most successful linear biphasic model for articular cartilage has been developed by Mow and others [2]. This model includes small deformation of poroelastic material which corresponds to Biot [3] theory for soil consolidation. Using this (Biot’s) theory, the governing equations for cartilage deformation and motion of interstitial fluid were formulated. Monsour and others [4] modeled the joint as porous permeable, deformable elastic material (articular cartilage) filled with a single layer of homogenous fluid. The tissue secretes viscous and highly non-Newtonian fluid called synovial fluid which mainly consists of hyaluronic acid, nutrients, glycoprotein, etc. This synovial fluid bathes and supplies nutrients to both surfaces of the cartilage. Hou and others [5] analyzed the squeeze film lubrication of articular cartilage by assuming synovial fluid to be linearly viscous fluid. Detailed analyses about articular cartilage and non-Newtonian characteristics of synovial fluid are given in [6].

Poroelasticty is a continuum theory for the analysis of porous media with elastic matrix consisting of interconnected fluid filled pores. When fluid permeates into a poroelastic material, the drag force between the fluid and the porous medium may cause deformation in the porous matrix. This leads to volumetric changes in the pores. Since the pores are filled with fluid, the presence of the fluid not only acts as a stiffener of the material but also results in the flow of pore fluid between regions of higher and lower pore pressure. A successful model for cartilage with interstitial fluid has been developed by Mow and his coworkers [6]. This simplest linear version of biphasic mixture includes the small deformation of the porous elastic matrix, which corresponds to Biot’s model for soil consolidation. The above biphasic model for a homogeneous and isotropic articular cartilage was used in a series of papers to model the fluid flow across the articular surface in geometrically idealized step-loaded synovial joints [5–7]. Various aspects of articular cartilage and non-Newtonian characteristics of the synovial fluid are presented by Torzilli and Mow [6]. Collins [8] considered Biot’s theory to model poroelastic matrix for cartilage which is assumed to satisfy generalized form of Darcy’s law for unsteady flow in an elastic porous medium. Later, modified and corrected forms of the constitutive equations were given in [9–12]. Most of the studies on synovial joint mechanics and lubrication have used elastic single-phase models of cartilage and a Newtonian single-phase model for synovial fluid. Recently, Mercer and Barry [13] gave a numerical method on finite difference approximations for the calculation of deformation, pressure, and flow within a finite two-dimensional poroelastic medium by considering slip and no-slip boundaries.

Sayles and others [14] experimentally revealed that the cartilage surface is rough and that roughness-height distribution is Gaussian in nature. This motivates us to investigate the effect of roughness in the cartilage surface. Christensen [15] developed the stochastic theory to investigate the effect of roughness in hydrodynamic lubrication on the assumption that roughness can be represented as a randomly varying quantity. It is assumed that in classical hydrodynamic lubrication theory, the roughness heights are small compared to the film thickness. This theory consists of two types of roughness structures, namely, longitudinal and transverse roughness. The former one has the roughness striations in the form of long ridges and narrows in horizontal direction whereas the latter one in the vertical direction.

The present chapter is organized as follows. In order to make it as self-contained as possible, the anatomy of the human knee joint and its bio-mechanism is given in Sect. 13.2. The simplified modified constitutive relations of poroelastic material are considered in the Sect. 13.3. In Sect. 13.4, the modified Reynolds equation which holds in the film region is derived for two roughness patterns, namely, longitudinal and transverse roughness structures. In Sect. 13.5, we describe briefly wavelet-multigrid method, recently developed by the author, for the solution of elliptic partial differential equations. Section 13.6 is devoted to the discussion of the results obtained in the previous sections for various parameters. Section 13.7 summarizes the major findings of the present study.

13.2 Anatomy and Bio-Mechanism of a Human Knee Joint

The knee joint is one of the most important joints of our body. It plays an essential role in movement related to carrying the body weight in horizontal (running and walking) and vertical (jumps) directions. The knee joint is one of the largest and most complex joints in the body. The knee joins the thighbone (femur) to the shinbone (tibia). The smaller bone that runs alongside the tibia (fibula) and the kneecap (patella) are the other bones that make the knee joint. Tendons connect the knee bones to the leg muscles that move the knee joint. Ligaments join the knee bones and provide stability to the knee in preventive and self-corrective ways. The anterior cruciate ligament prevents the femur from sliding backward on the tibia (or the tibia sliding forward on the femur).The posterior cruciate ligament prevents the femur from sliding forward on the tibia (or the tibia from sliding backward on the femur). The medial and lateral collateral ligaments prevent the femur from sliding side to side.

Two C-shaped pieces of cartilage called the medial and lateral menisci act as shock absorbers between the femur and tibia. Numerous synovial fluid-filled sacs help the knee move smoothly. When a person walks or stands still, the force L acting through pushes the femoral condyles and the tibial plateau together when L  >  0 and to create a squeeze film effect when L  <  0. The synovial fluid is sucked into or sucked out of the cartilage. Cartilage is avascular (free from blood supply) and the flow of synovial fluid into the cartilage is one of the ways it receives nutrients. However, when a person stands still for an extended period of time, it is believed that the synovial fluid will eventually be squeezed out of the gap between tibial plateau and femoral condyles, causing direct contact between cartilage-coated surfaces. This is an undesirable condition, as it results in cartilage degradation and promotes in the long run osteoarthritis.

The physical configuration of the problem is shown in [16] to predict the performance of the knee joint. The bone ends are covered by articular cartilage to prevent natural abrasion, which is in a sac containing fluid for lubricating the two surfaces. A tough fibrous capsule together with the muscles, ligaments, intra-articular structures, etc. encloses the normal joint cavity. The inner lining of this capsule, the synovial membrane, secretes viscous and highly lubricating fluid called synovial fluid. This fluid bathes both articular surfaces and intra-articular structures. Following Walker and Erkman [17], as the load-bearing area of the synovial knee joint is small, two articular surfaces may be considered to be parallel under high loading conditions. For mathematical simplicity, the average of the three layers of the cartilages is modeled as a single poroelastic layer. So, the problem considered is that of squeeze film lubrication between two rectangular surfaces with finite dimensions.

13.3 Mathematical Modeling and Computer Simulation

We present a simple mathematical model for the human knee joint. The governing equations are the principles of mass and momentum balances applied to the simplified geometry under the following hypotheses.

  1. 1.

    The lubricant in the film region is modeled as Newtonian, i.e., linearly viscous and incompressible fluid.

  2. 2.

    For simplicity, both the tibial plateau and the femoral condyles are modeled as rigid impermeable flat surfaces and the effects of menisci are neglected.

  3. 3.

    Classical lubrication theory is used to describe the flow of synovial fluid in the thin gap.

  4. 4.

    The flow is assumed to be steady state, laminar, incompressible, and three-dimensional.

  5. 5.

    Viscosity is kinematical and depends on the hyaluronic acid concentration.

The application of hypotheses from 1 to 5 to the conservative laws leads to a modified Reynolds lubrication equation giving fluid pressure in the joint which supports the load avoiding the direct contact between the solid surfaces.

The upper rigid rough impervious spherical indenter approaches the lower poroelastic matrix normally with a constant velocity \( dH/dt. \) The film thickness between two plates is characterized by:

$$ H=h(x,y)+{{h}_{s}}(x,y,\xi ), $$
(13.1)

where \( h(x,y)={{h}_{0}}+( {{x}^{2}}+{{y}^{2}} )/2R \) is the nominal smooth part of the geometry, \( {{h}_{0}} \) is the minimum thickness, R is the radius of the indenter in x–y plane, \( {{h}_{s}} \) is the height of the surface asperities measured from the nominal level which is a randomly varying quantity of zero mean and \( \xi\) is the index parameter determining a definite roughness structure. All the articulations of knee joint under fluid film lubrication involve cartilage–viscous fluid–cartilage interactions. The lubricant in the film region is assumed to be Newtonian fluid, i.e., linearly viscous. So, the problem considered would be that of three-dimensional squeeze film lubrication between the upper rigid spherical indenter and the lower smooth poroelastic matrix. With usual assumptions of lubrication, the governing equations are:

$$ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0, $$
(13.2)
$$ \frac{{{\partial }^{2}}u}{\partial {{z}^{2}}}=\frac{1}{\mu }\frac{\partial p}{\partial x}, $$
(13.3)
$$ \frac{{{\partial }^{2}}v}{\partial {{z}^{2}}}=\frac{1}{\mu }\frac{\partial p}{\partial y} $$
(13.4)
$$ 0=\frac{\partial p}{\partial z}, $$
(13.5)

where \(u,v,\text{ and }w\) are the velocity components in \(x,y,\ \text{and}\ \text{z}\) directions, respectively, \(p\) is the pressure, and \(\mu \) is the viscosity of the fluid.

13.3.1 Poroelastic Region

Following Mow and others [2], the poroelastic material is assumed to be made of a linearly elastic solid phase (deformable cartilage matrix) and an inviscid fluid phase in which these two phases are intrinsically incompressible. For the cartilage matrix deformation and viscous fluid contained in its pores, we write modified coupled equations of motion as in [9, 12, 17]:

$$ \mathop{Matrix:\ \rho }_{m}\frac{\mathop{\partial }^{2}\mathbf{U}}{\partial \mathop{t}^{2}}=div\mathop{\sigma }_{m}-{{( {{\phi }^{f}} )}^{2}}\frac{\mu }{k}\left( \frac{\partial \mathbf{U}}{\partial t}-\mathbf{V} \right) $$
(13.6)
$$ \text{Fluid:}\ {{\rho }_{f}}\frac{\partial \mathbf{V}}{\partial t}=div\mathop{\sigma }_{f}+{{( {{\phi }^{f}} )}^{2}}\frac{\mu }{k}\left( \frac{\partial \mathbf{U}}{\partial t}-\mathbf{V} \right)\ \text{and} $$
(13.7)
$$ \nabla \cdot \left( {{\phi }^{s}}\frac{\partial \mathbf{U}}{\partial t}+{{\phi }^{f}}\mathbf{V} \right)=0, $$
(13.8)

where \({{\phi }^{f}}\) is the porosity and \({{\phi }^{s}}=1-{{\phi }^{f}}\) is the solidity of the poroelastic material. \({{\rho }_{m}}\ \text{and}\ {{\rho }_{f}}\) denote the densities of solid matrix and fluid, respectively, \(\mathbf{U}\) is the corresponding displacement vector, \(k\) is the permeability of the cartilaginous matrix to the fluid. The left hand terms denote the local forces (mass times acceleration), which are counterbalanced by the right hand terms namely the surface forces,\(div\sigma ,\) and the porous medium driving forces (Darcy’s law), respectively. These two component equations may be simply viewed as generalized forms of Darcy’s law for unsteady flow in a deformable porous medium in terms of relative velocity \(( \partial \mathbf{U}/\partial t-\mathbf{V} )\) between the moving cartilage and the fluid contained in its pores. Also, Eqs. 6 and 7 denote force balances for the linear elastic solid and viscous fluid components of the cartilage, respectively. The classical stress tensor \(\sigma \) for a continuous homogeneous medium may be expressed for the matrix (cartilage) and fluid (synovial). Thus, the constitutive relations for the solid and fluid phases, respectively are

$$ {{\sigma }_{m}}=-{{\phi }^{s}}P\mathbf{I}+2Ne+{{A}^{'}}e\mathbf{I}, $$
(13.9)
$$ {{\sigma }_{f}}=-{{\phi }^{f}}P\mathbf{I}+{{E}^{'}}e\mathbf{I}, $$
(13.10)

in terms of the elastic parameters \(N,{{E}^{'}},\ \text{and}\ {{A}^{'}}\) of the cartilage and the hydrostatic pressure \(P\) and \(\mathbf{I}\) the identity tensor, e the cartilage dilation. The inertial terms in Eqs. 13.6 and 13.7 are neglected, because in the balance of the momentum equation, the fluid–fluid viscous stress is negligible compared with the drag between the fluid and the solid matrix [10]. After neglecting the inertia terms in Eqs. 13.6 and 13.7 and using this result into Eq. 13.8, we get:

$$ {{\nabla }^{2}}e-\frac{\mu }{k(A+{E}'+2N)}\frac{\partial e}{\partial t}=0, $$
(13.11)

where \(e=\nabla \cdot \mathbf{U}\). The typical stress-strain curve obtained from confined compression tests an articular cartilage under loading stresses is essentially a linear relationship beyond the initial stress level [18]. Following Collins [8], their results may be characterized by a simple linear equation in terms of the corresponding average bulk modulus K and pressure Pas:

$$ e={{e}_{0}}+\frac{P}{K}. $$
(13.12)

Using this empirical relation (Eq. 13.12) in Eq. 13.11, we get the equation describing pressure in the poroelastic region:

$$ {{\nabla }^{2}}P=0. $$
(13.13)

13.3.2 Boundary Conditions

The relevant boundary conditions for the velocity field (0  <  z  <  H) are:

$$ u(x,y,0)=u(x,y,H)=v(x,y,0)=v(x,y,H)=0, $$
(13.14)
$$ w(x,y,0)=-{{w}_{n}},\text{}w(x,y,H)=-\frac{dH}{dt}, $$
(13.15)

where \({{w}_{n}}\) represents the normal component of the relative velocity of the fluid at the cartilage surface. Conditions (Eq. 13.14) are no-slip velocity conditions.

13.4 Solution Procedure

Integrating Eq. 13.13 with respect to z over the porous layer thickness (\(-\delta \) , 0) and using the solid backing boundary condition \(\partial P/\partial z=0\ \text{at}\ z=-\delta \) we get:

$$ {{\left. \frac{\partial P}{\partial z} \right|}_{z=0}}=-\int\limits_{-\delta }^{0}{\left( \frac{{{\partial }^{2}}P}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}P}{\partial {{y}^{2}}} \right)}dz. $$
(13.16)

Here \(\delta \)is the thickness of the poroelastic layer. Using the Morgan–Cameron approximation [19] valid for the poroelastic layer to be small and incorporating the pressure continuity condition (p = \({{p}^{*}}\)) at the porous interface z = 0, we get:

$$ {{\left. \frac{\partial P}{\partial z} \right|}_{z=0}}=-\delta \left( \frac{{{\partial }^{2}}p}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}p}{\partial {{y}^{2}}} \right). $$
(13.17)

After neglecting inertia terms, Eq. 13.7 may be arranged in terms of relative velocity in the form:

$$ \left( \mathbf{V}-\frac{\partial \mathbf{U}}{\partial t} \right)=-\frac{k}{\mu {{( {{\phi }^{f}} )}^{2}}}( {{\phi }^{f}}\nabla P-{{E}^{'}}\nabla e ). $$
(13.18)

Elimination of \(e\) through Eqs. 13.12 and 13.18 gives:

$$ \left( \mathbf{V}-\frac{\partial \mathbf{U}}{\partial t} \right)=-\frac{k}{\mu {{( {{\phi }^{f}} )}^{2}}}\nabla P\left( {{\phi }^{f}}-\frac{{{E}^{'}}}{K} \right). $$
(13.19)

The normal component of the relative fluid velocity at the cartilage surface is

$$ {{w}_{n}}=-{{\left( \mathbf{V}-\frac{\partial \mathbf{U}}{\partial t} \right)}_{n}}=-\frac{k}{\mu {{( {{\phi }^{f}} )}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right){{\left. \frac{\partial P}{\partial z} \right|}_{z=0}}. $$
(13.20)

Using Eq. 13.17 in Eq. 13.20, we get:

$$ {{w}_{n}}=\frac{k\delta }{\mu {{( {{\phi }^{f}} )}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right)\left( \frac{{{\partial }^{2}}p}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}p}{\partial {{y}^{2}}} \right). $$
(13.21)

Equations 13.3 and 13.4 can be integrated for u and v with respect to zusing boundary conditions (Eq. 13.14). Substituting u and v in Eq. 13.2 and integrating across the film thickness from \(z=0\ \text{to}\ z=H\) with respect to z using boundary conditions (Eq. 13.15), we obtain modified Reynolds equation

$$ \begin{aligned}& \frac{\partial }{\partial x}\left\{ \left[ {{H}^{3}}-\frac{24\delta k}{{{\left( {{\phi }^{f}} \right)}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial p}{\partial x} \right\}+\frac{\partial }{\partial y}\left\{ \left[ {{H}^{3}}-\frac{24\delta k}{{{\left( {{\phi }^{f}} \right)}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial p}{\partial y} \right\}=-12\mu \frac{dH}{dt}, \\& \quad \quad \quad\\ \end{aligned} $$
(13.22)

for the pressure distribution in the joint cavity. For including roughness features, taking the stochastic average of Eq. 13.22 with respect to \({{h}_{s}},\) we get:

$$ \begin{array}{*{35}{l}} \frac{\partial }{\partial x}\left\{ E\left[ \left[ {{H}^{3}}-\frac{24\delta k}{{{\left( {{\phi }^{f}} \right)}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial p}{\partial x} \right] \right\}+\\ {}\\\end{array}\frac{\partial }{\partial y}\left\{ E\left[ \left[ {{H}^{3}}-\frac{24\delta k}{{{\left( {{\phi }^{f}} \right)}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial p}{\partial y} \right] \right\}=-12\mu \frac{dE( H )}{dt}, $$
(13.23)

where expectancy operator \(E(\bullet )\)is defined by:

$$ E(\bullet )=\int\limits_{-\infty }^{\infty }{(\bullet )f({{h}_{s}})d{{h}_{s}}}, $$
(13.24)

\(f\) is the probability density function of the stochastic film thickness \({{h}_{s}}.\) In most of the problems, surfaces show a roughness height distribution which is Gaussian in nature. Therefore, polynomial form which approximates the Gaussian is chosen in the analysis. Such a probability density function is [12]:

$$ f({{h}_{s}})=\left\{ \begin{aligned}& \frac{35}{32{{c}^{7}}}{{\left( {{c}^{2}}-h_{s}^{2} \right)}^{3}},\text{}-c<{{h}_{s}}<c \\& \text{ 0,elsewhere} \\ \end{aligned} \right. $$

where \(c=3{{\sigma }_{1}}\ \text{and}\ {{\sigma }_{1}}\) is the standard deviation.

13.4.1 Longitudinal Roughness

In the context of rough surfaces, there are two types of roughness patterns, which are of special interest. The one-dimensional longitudinal structure where the roughness has the form of long narrows ridges and valleys running in the x-direction and the one-dimensional transverse structure where roughness striations are running in the y-direction in the form of narrow and valleys. However, the present study is restricted to one-dimensional longitudinal roughness since the one roughness structure can be obtained from other by just rotation of the coordinate axes. For the one-dimensional longitudinal roughness pattern, the film thickness assumes the form:

$$ H={{h}_{0}}(t)+{{h}_{s}}(x,\xi ), $$

then, modified Stochastic Reynolds equation (Eq. 13.23) takes the form:

$$ \begin{aligned} &\frac{\partial }{\partial x}\left\{ [ E\left( {{H}^{3}} )-\frac{24\delta k}{{{\left( {{\phi }^{f}}\right)}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right)\right]\frac{\partial E(p)}{\partial x} \right\}+\\& \frac{\partial }{\partial y}\left\{ \left[ \frac{1}{E\left( {1}/{{{H}^{3}}}\; \right)}-\frac{24\delta k}{{{\left( {{\phi }^{f}} \right)}^{2}}}\left( \frac{{{E}^{'}}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial E(p)}{\partial y} \right\}=-12\mu \frac{dE( H )}{dt}. \end{aligned}$$
(13.25)

13.4.2 Transverse Roughness

In this structure, the roughness striations are running in the y-direction in the form of narrow ridges and valleys. In this case, the film thickness takes the form:

$$ H=h+{{h}_{s}}(y,\xi ) $$

and modified Reynolds equation becomes:

$$ \frac{\partial }{\partial x}\left\{ \left[ \frac{1}{E\left( {{H}^{-3}} \right)}-\frac{24k\delta }{{{({{\phi }^{f}})}^{2}}} \right]\frac{\partial E(p)}{\partial x} \right\}+\frac{\partial }{\partial y}\left\{ \left[ E({{H}^{3}})-\frac{24k\delta }{{{({{\phi }^{f}})}^{2}}} \right]\frac{\partial E(p)}{\partial y} \right\}=-12\mu \frac{dE(H)}{dt} $$
(13.26)

For the distribution function given by Eq. 13.24, we have:

$$ E(H)=h. $$
(13.27)
$$ E({{H}^{3}})={{h}^{3}}+\frac{h{{c}^{2}}}{3} $$
(13.28)
$$ E( {\scriptstyle{}^{1}\!\!\diagup\!\!{}_{{{H}^{3}}}\;} )=\frac{35}{32{{c}^{7}}}[ 3( 5{{h}^{2}}-{{c}^{2}} )\left( {{c}^{2}}-{{h}^{2}} )\log \left( \frac{h+c}{h-c} \right)+2ch( 15{{h}^{2}}-13{{c}^{2}} \right) ] $$
(13.29)

Therefore, Reynolds equation (Eq. 27) for longitudinal roughness takes the form:

$$ \begin{aligned}& \frac{\partial }{\partial x}\left\{ \left[ E({{H}^{3}})-\frac{24k\delta }{{{({{\phi }^{f}})}^{2}}}\left( \frac{{{E}'}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial E(p)}{\partial x} \right\}+\frac{\partial }{\partial y}\left\{ \left[ \frac{1}{E({{H}^{-3}})}-\frac{24k\delta }{{{({{\phi }^{f}})}^{2}}}\left( \frac{{{E}'}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial E(p)}{\partial y} \right\}=-12\mu \frac{dh}{dt}. \\&\\ \end{aligned} $$
(13.30)

Introducing the following nondimensional parameters and variables, in Eq. 13.30:

$$\begin{aligned}& \bar{x}=\frac{x}{a},\text{ }\bar{y}\text{ = }\frac{y}{a},\quad \bar{H}=\frac{H}{{{h}_{0}}}=1+\frac{{{{\bar{x}}}^{2}}+{{{\bar{y}}}^{2}}}{2{{R}_{1}}{{h}_{1}}}+{{{\bar{h}}}_{s}},\text{}{{R}_{\text{1}}}=\frac{R}{a},\text{ }{{h}_{\text{1}}}=\frac{{{h}_{0}}}{a},{{{\bar{h}}}_{s}}=\frac{{{h}_{s}}}{{{h}_{0}}},\\&\quad \bar{k}=\frac{k\delta }{h_{0}^{3}},\ \bar{p}=\frac{E(p)h_{0}^{3}}{\mu _{0}^{2}{dh}/{dt}\;},\text{ }c=\frac{C}{{{h}_{0}}},\end{aligned}$$

we get

$$ \frac{\partial }{\partial \bar{x}}\left\{ \left[ E({{{\bar{H}}}^{3}})-\frac{24\bar{k}}{{{({{\phi }^{f}})}^{2}}}\left( \frac{{{E}'}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial \bar{p}}{\partial \bar{x}} \right\}+\frac{\partial }{\partial \bar{y}}\left\{ \left[ E({{{\bar{H}}}^{-3}})-\frac{24\bar{k}}{{{({{\phi }^{f}})}^{2}}}\left( \frac{{{E}'}}{K}-{{\phi }^{f}} \right) \right]\frac{\partial \bar{p}}{\partial \bar{y}} \right\}=-12. $$
(13.31)

Here:

$$ \begin{aligned}& E({{{\bar{H}}}^{3}})={{{\bar{h}}}^{3}}_{0}+\frac{\bar{h}{{C}^{2}}}{3} \\& E({{{\bar{H}}}^{-3}})=\frac{35}{32{{C}^{7}}}\left( 3(5{{{\bar{h}}}^{2}}-{{C}^{2}})({{C}^{2}}-{{{\bar{h}}}^{2}})log\left( \frac{\bar{h}+C}{\bar{h}-C} \right)+2C\bar{h}(15{{{\bar{h}}}^{2}}-13{{C}^{2}}) \right), \\ \end{aligned} $$

and for pressure distribution, the boundary conditions are:

$$ \begin{aligned}& \bar{p}=0\ \text{at}\ \bar{x}=\pm 1, \\& \bar{p}=0\ \text{at}\ \bar{y}=\pm 1. \\ \end{aligned} $$
(13.32)

13.5 Wavelet-Multigrid Method

Since Eq. 13.32 is an elliptic equation with variable coefficients, it is difficult to solve analytically; we propose to solve it numerically using recently developed wavelet-multigrid method [20]. Using the first order finite difference scheme for the derivatives in Eq. 13.32, we write:

$$ {{A}_{0}}{{\bar{p}}_{i+1,j}}+{{A}_{1}}{{\bar{p}}_{i-1,j}}+{{A}_{2}}{{\bar{p}}_{i,j+1}}+{{A}_{3}}{{\bar{p}}_{i,j-1}}+{{A}_{4}}{{\bar{p}}_{i,j}}={{A}_{5}}, $$
(13.33)

where

$$ \begin{array}{*{35}{l}} {{A}_{0}}=\frac{{{\alpha }_{i+1/2,j}}}{{{(\Delta x)}^{2}}},\text{ }{{\text{A}}_{\text{1}}}=\frac{{{\alpha }_{i-1/2,j}}}{{{(\Delta x)}^{2}}},\text{ }{{\text{A}}_{\text{2}}}=\frac{{{\beta }_{i,j+1/2}}}{{{(\Delta y)}^{2}}},\text{ }{{\text{A}}_{\text{3}}}=\frac{{{\beta }_{i,j-1/2}}}{{{(\Delta y)}^{2}}},\text{ }\\ {{A}_{4}}=-({{A}_{0}}+{{A}_{1}}+{{A}_{2}}+{{A}_{3}}),\text{ }{{\text{A}}_{\text{5}}}=-12,\\ {{\alpha }_{i,j}}=E(\bar{H}_{i,j}^{3})-\frac{24\bar{k}}{{{({{\phi }^{f}})}^{2}}}\left( \frac{{{E}'}}{K}-{{\phi }^{f}} \right),\text{ }{{\beta }_{i,j}}=E(\bar{H}_{i,j}^{-3})-\frac{24\bar{k}}{{{({{\phi }^{f}})}^{2}}}\left( \frac{{{E}'}}{K}-{{\phi }^{f}} \right).\\\end{array} $$

Let the matrix formulation of Eq. 13.33 be:

$$ A\mathbf{x}=\mathbf{y}. $$
(13.34)

By considering D-4 family of wavelets, the corresponding transform matrix (just for illustrative purpose) is:

$$ W=\frac{1}{\sqrt{2}}\left[ \begin{matrix} {{c}_{0}} & {{c}_{1}} & {{c}_{2}} & {{c}_{3}} & 0 & 0 & 0 & 0\\ {{d}_{0}} & {{d}_{1}} & {{d}_{2}} & {{d}_{3}} & 0 & 0 & 0 & 0\\ 0 & 0 & {{c}_{0}} & {{c}_{1}} & {{c}_{2}} & {{c}_{3}} & 0 & 0\\ 0 & 0 & {{d}_{0}} & {{d}_{1}} & {{d}_{2}} & {{d}_{3}} & 0 & 0\\ 0 & 0 & 0 & 0 & {{c}_{0}} & {{c}_{1}} & {{c}_{2}} & {{c}_{3}}\\ 0 & 0 & 0 & 0 & {{d}_{0}} & {{d}_{1}} & {{d}_{2}} & {{d}_{3}}\\ {{c}_{2}} & {{c}_{3}} & 0 & 0 & 0 & 0 & {{c}_{0}} & {{c}_{1}}\\ {{d}_{2}} & {{d}_{3}} & 0 & 0 & 0 & 0 & {{d}_{0}} & {{d}_{1}}\\\end{matrix} \right], $$
(13.35)

where \({{c}_{i}}\ \text{and}\ {{d}_{i}},\)’s are the filter coefficients [21].

Now the fast wavelet transform (FWT) is performed on A and y of Eq. 13.35 recursively till the coarsest level is reached at level − J. Then matrix equation \({{\hat{A}}_{l}}{{\mathbf{\hat{x}}}_{l}}={{\mathbf{\hat{y}}}_{l}}\) is used to obtain \({{\mathbf{\hat{x}}}_{l}}\) at the coarsest level using the Gauss–Jordan method. Finally, integer wavelet transform (IWT) is performed on \(\mathbf{\hat{x}}\) to obtain \({{\mathbf{\hat{x}}}_{0}}\) level − 1 which gives the distribution of fluid film pressure \(\bar{p}.\) In the calculation, for all numerical experiments D-4 wavelets are employed. However, one has the freedom and flexibility to choose any wavelet basis. To test the accuracy, the problem is solved at resolutions \({{2}^{4}}\text{ and }{{\text{2}}^{\text{8}}}.\) It is observed that there is marginal increase in the accuracy of the solution; better accuracy can be achieved by increasing the resolution and/or the order of the wavelet family. It is also observed that the amount of computational effort is considerably less than that of classical multigrid method.

Once we have obtained the fluid film pressure by using the wavelet-multigrid method, the load carrying capacity can be evaluated straightforwardly. We can find the nondimensional load carrying capacity \(\bar{W}\) per unit area of the joint surface using the following formula.

$$ \bar{W}=\int\limits_{-1}^{1}{\int\limits_{-1}^{1}{\bar{p}(\bar{x},\bar{y})d\bar{x}d\bar{y}}}. $$
(13.36)

13.6 Results and Discussion

A simplified mathematical model has been developed for the study of effect of surface roughness on the poroelastic bearing by modeling articulate cartilage as poroelastic matrix and synovial fluid as a linearly viscous fluid. The problem considered is that of the squeeze film lubrication between the rough spherical indenter and smooth poroelastic material. The modified Reynolds equation incorporates the constitutive equations of poroelastic matrix. This Reynolds equation is solved using wavelet-multigrid method. The values of \({E}',K,\ \text{and}\ \bar{k}\) are taken from Torzilli [6], which are associated with healthy human cartilage with normal functioning. The parameters \(\bar{k}\) and \({{\phi }^{f}}\) are chosen in such way that the load-carrying capacity that can be sustained by the joints should be more. The pressure \(\bar{p}\ \text{and}\ \bar{w}\) are functions of nondimensional parameters \(C,\text{ }\bar{k},\ \text{and}\ \text{ }\bar{H}\) and constant parameters \({{\phi }^{f}}\ \text{and}\ {E}'/K.\) As the radius of the upper rigid indenter increases, it becomes relatively flat and uniform, which leads to increase in the area of the film region. This wide film area avoids exit of lateral fluid and is responsible for retaining the large amount of fluid which enhances the pressure and the load-carrying capacity. For large radius, the spherical indenter becomes relatively flat and thus the study reduces to the squeeze film lubrication between two parallel plates.

The distribution of pressure \(\bar{p}\) with rectangular coordinates \(\bar{x}\;{\text{and}}\;\bar{y}\) is given in Figs. 13.1a and 13.1b. For C = 0.3, the roughness effects are more pronounced when compared with C = 0.1. In Fig. 13.2, it is observed that the effect of roughness is to increase the pressure distribution increasing values of C. The roughness increases the pressure in the film region compared with the smooth case. This is because, the presence of surface asperities reduces the fluid flow and therefore large fluid accumulates in the film region, which enhances the pressure built up.

Fig. 13.1
figure 1

a Fluid film pressure distribution for C = 0.3 with \(\bar{k}=7.65\times {{10}^{-5}},{{\phi }^{f}}=0.8\ \text{and}\ \mathbf{{E}'}/\mathbf{K=1}\mathbf{.0}\). b Fluid film pressure distribution for C = 0.1 with \(\bar{k}=7.65\times {{10}^{-5}}{{\phi }^{f}}=0.8\ \text{and}\ \mathbf{{E}'}/\mathbf{K=1}\mathbf{.0}\)

Full size image
Fig. 13.2
figure 2

Variation of pressure \(\bar{p}\) with \(\bar{x}\) for different values of C with \(\mathbf{\bar{k}=7}\mathbf{.65\times 1}{{\mathbf{0}}^{\mathbf{-5}}},{{\mathbf{f}}^{\mathbf{f}}}\mathbf{=0}\mathbf{.8}\ \text{and}\ \mathbf{{E}'}/\mathbf{K=1}\mathbf{.0}\)

Full size image

The variation of load \(\bar{w}\) that can be sustained by a joint, with roughness parameter C for different elastic parameter \({E}'/K\) is shown in Fig. 13.3. The load-carrying capacity increases with increasing \({E}'/K\) for different roughness parameters C. Presence of hyaluronic acid, molecular, and other molecular weight substances in the synovial fluid results in the formation of a thick dense substance on the cartilage surface during squeezing action. This leads to increase in the pressure distribution and in turn increase in the load-carrying capacity. This trend is preserved for all values of \({E}'/K\).

Fig. 13.3
figure 3

Variation of load \(\mathbf{\bar{W}}\) with C for different values of \(\mathbf{{E}'}/\mathbf{K}\) with \(\mathbf{\bar{k}=7}\mathbf{.65\times 1}{{\mathbf{0}}^{\mathbf{-5}}},{{\mathbf{f}}^{\mathbf{f}}}\mathbf{=0}\mathbf{.8}\)

Full size image

13.7 Summary

Wavelet-multigrid method is found to be accurate for the solution of modified Reynolds equation, which incorporates surface roughness and poroelastic nature of articular cartilage. The method is proven to be more efficient than the classical multigrid method. Comparing the rates of convergence of the two methods, it is found that Wavelet-multigrid method converges faster compared with the multigrid scheme. It provides the ability to carry out calculations focusing only in regions where it is needed. It can as well be generalized to nonuniform grids, which is a unique feature of this approach. In classical scheme we find that as the grid size decreases the condition number increases; a serious disadvantage whereas in the present Wavelet-multigrid scheme the conditioning of the matrix is incorporated in the procedure itself without requiring separate analysis for it. It is observed that the effect of roughness is to increase the pressure built up (and hence the load-carrying capacity) compared with the classical case. The governing equations describing complex structure of cartilage and synovial fluid are complicated because of nonlinearity and joints having a wide range of articulating features. However, the proposed simplified model does predict some of the salient features of bearing characteristics, which would enable the biomedical engineers in selecting suitable design parameters. In addition to being of scientific and medical interest in its own right, the present study has achieved the objective of gaining deeper understanding of the lubrication in the knee joint which is of crucial importance in the total knee-replacement technologies. The results obtained could guide the experimentation with new material for knee replacement with mechanical characteristics analyzed here. Currently, the most promising compliant materials are hydro gels, which have similar properties to those of natural cartilage.