Synovial Joint Study

Nawal Odah Al-Atawi, Daoud Suleiman Mashat

Department of Mathematics, Numerical Analysis, King Abdulaziz University, Jeddah, KSA.

DOI: 10.4236/ajcm.2022.121002                

Abstract

The impact of certain separate characteristics, including the porosity parameter, reaction rate parameter, and viscoelastic parameters of steady convective diffusion across a rectangular channel, has been investigated in this article. The model’s momentum and concentration equations were developed using the similarities technique, and the numerically finite volume method was combined with the Beavers and Joseph slip conditions. Various graphs have been used to get insight into various parameters of the problem on velocity and concentration. The cartilage surfaces are assumed to be porous, and the viscosity of synovial fluid varies with hyaluronate (HA) content.

Keywords

Lubrication Theory, Rheology of Synovial Fluid, Velocity & Concentra-tion, Similarity Parameters, Finite Volume Method

1. Introduction

The synovial joints are crucial in both human and animal movement. Under normal physiological conditions (see Figure 1), these joints can withstand very high loads and have very low friction [1]. Articular cartilage is a porous gel of proteoglycan aggregates embedded in a water-swollen network of collagen fibrils [1]. When cartilage is compressed, the interstitial fluid is forced to flow relative to the solid organic matrix and exude from it [2]. The vessels that transport nutrients are present in immature articular cartilage [2] [3]. The extracellular fluid transports nutrients from the synovial fluid via diffusion and convection [4]. The process of dispersion is critical in both chemical and biological systems [5] [6]. Hyaluronic acid [HA], glycoprotein, and other macromolecular components move from synovial fluid to cartilage in synovial joints [5] [6]. In healthy synovial fluid, the diffusion coefficient of hyaluronan was on average 30% slower than expected by sample viscosity [5] [6]. While dispersing from synovial fluid

Figure 1. Shows physical information of synovial joint (Human Knee joint) [9].

to cartilage, HA and other components may also perform microcirculation [7]. Endothelium in arteries and cartilage in synovial joints are porous cell layers that can or cannot deform [7]. Synovial fluid is essential for joint lubrication, nutrient delivery, and metabolite removal from avascular articular cartilage [7] [8]. Increased lubrication, imbibition, and exudation increase the concentration of hyaluronic acid molecules in synovial fluid [7] [8]. According to Wegamir’s findings, an increase in hyaluronic acid concentration causes an increase in synovial fluid viscosity.

1.1. Models of Lubrication

The most prominent theory of lubrication regarding human joints was regulated by hydrodynamic. Well-built investigations were carried out to find the similarities between physical wedges and articulating surfaces in human joints shown by Reynolds which is necessary for hydrodynamic or fluid film lubrication [10]. The intraarticular cartilage (synovial fluid in knee joint) behaves as a slight inclination of the coupled surfaces to cater wedge shaped lubricant films, which is suggested by MacConail [10] [11]. After his research, there was conclusion about joints to view as hydrodynamically lubricated systems by performing the first detailed experiments on friction in naturally acting joints, observing a horse stifle mounted as the fulcrum in a pendulum machine. A very low coefficient of friction ( $f=0.02$ ) in case of sliding was observed when the joint was lubricated by synovial fluid whereas an unlubricated joint revealed creaking noise after 8000 cycles when applied a constant load of 445 N [12] [13]. After 4 hours, observations were the joint heating, debris being thrown, bone granting, stream rising, etc. [12] [13]. After such experiment he highlighted the essential role of lubrication in natural joints. So hydrodynamic theory at that time played a fundamental role in the study of natural human joint lubrication [12] [13] [14]. Even if this topic remains under investigation, to summarize and highlight the scientific progress coupled with mathematical equations made through the years in the understanding and modelling of lubrication mechanisms in natural human synovial joints, covering more popular and relevant theories from earlier research in this field up to the latest ones, can be viewed by introducing briefly some milestones from the 1930’s to today [12] [13] [14].

1.2. Boundary Lubrication

The conclusion was drown which was based on several pendulum experiments (in which a linear decay in amplitude was observed) that hydrodynamic action could not exist due to low sliding velocities under the heavy loads acting in the human joints. An alternative to MacConail theory, a boundary lubrication model was proposed by Charnley in 1959 by his experiment which recorded friction coefficient values [13] [14]. Such lubrication exists between unlubricated sliding surfaces and fluid-film lubrication. It is also characterized as condition in which the friction between surfaces is determined by features of the surfaces and properties of the lubricant other than viscosity [13] [14]. Barnett and Cobbold made critical observations on Charnley’s theory, in which linear decay in amplitude was attributable to the dissection of the joint [15]. Such decay showed an essential linear relationship with time when replaced the joint in the pendulum fulcrum with a hydrostatic bearing. In 1968, Linn and Radin conducted experiments on animal joints at a constant load attracted a lot of interests in which extraneous dynamic force component give rise to friction force [16]. The most important conclusion was the animal joints operate within mixed film region of lubrication.

1.3. Weeping Lubrication

McCutchen considered the porosity and elasticity of articular cartilage of joint lubrication problem for the first time in 1966 [17]. According to his views, the pressured synovial fluid travels through the porous cartilage which acts as a sponge like material [17]. Weeping lubrication is a form of fluid lubrication in which the load bearing surfaces are held apart by a film of lubricant that is maintained under pressure [16] [17] [18]. Weeping implied that the lubricant film was sweated into the high-pressure zone between opposing cartilages, while the boundary lubrication effect between the contact surfaces remained [16] [17] [18]. McCutchen came up with the term ‘‘weeping lubrication” because of bearing materials that accomplish it weep fluids when pressurized. An external pump is frequently used to supply pressure in engineering aspects while the pump action in the human body is provided by muscle contractions around the joint or by weight bearing compression of articular cartilage, which causes the cartilage to distort and “weep” fluid, forming a fluid film across the articular surfaces [16] [17] [18]. Because the impervious layer of calcified cartilage prevents the fluid from being driven into the subchondral bone, it can only move into the joint. When the stress is removed, osmotic pressure causes the fluid to flow back into the articular cartilage [16] [17] [18]. When cartilage weeps under force, the weight is absorbed by the hydrostatic pressure of the fluid, even though the liquid content of the porous cartilage would be evacuated and the friction coefficient would grow if the squeeze is sustained indefinitely. It also appears that wringing water out of the porous matrix takes longer and is more difficult, whereas resoaking is simple [16] [17] [18]. This sort of lubrication is most effective when the load is high, but it can be used in a variety of situations.

1.4. Elastohydrodynamic Lubrication

The protective fluid film is maintained at an adequate thickness in the elastohydrodynamic lubrication model by the elastic deformation of the articular surfaces [19]. Elastohydrodynamic lubrication is a hydrodynamic process in which the fluid film pressure produces elastic deformation of the bearing surfaces, altering the pressure in the film region. To put it another way, the elastic cartilage deforms slightly to keep a layer of fluid between the opposing joint surfaces. Under strong loading, the elastohydrodynamic motion can keep a fluid film intact [19]. The fluid film thickness varies between 10⁻⁵ and 10⁻⁴ cm [19]. In general, elastic distortion offers greater geometrical conformity in the contact zone than is ordinarily present, resulting in a substantially thicker lubricating film for a given load [19]. Furthermore, the high pressure may produce a significant rise in the viscosity of the lubricant. The thickness of the lubricating film is further increased by this viscosity effect [19]. With the addition of the effective elastic modulus, the minimum film thickness is a function of the same factors as in hydrodynamic lubrication [19]. In 1963, Dintenfass firstly, showed failure of the hydrodynamic lubrication mode in which he completely neglected deformations while he took into account the deformability of articular cartilage and led to Elastrohydrodynamic lubricayion theory (EHD) [20]. The main finding was that the operating film thickness in highly loaded lubricated contacts could be up to 100 times greater than predicted by conventional lubrication theory, because the human joint requires the correct operating conditions for the cartilage to be separated by a synovial film and thus boundary lubrication does not exist [20]. In 1966, Tanner theoretical values suggested the possibility of lubricating film in the human joints when surfaces are subject to relative motion [21]. Later on, in 1986 confirmed EHD lubrication as prominent mode of lubrication of soft bearing or highly deformable systems, such as synovial joints [21].

1.5. Squeeze-Film Lubrication

When the bearing surfaces are moving perpendicularly towards one other, squeeze-film lubrication occurs. The movement of articular surfaces that are perpendicular to one another creates pressure in the fluid film. The fluid film is squeezed out of the area of approaching contact as the opposing surfaces draw closer together. The viscosity of the fluid in the gap between the surfaces creates pressure, which causes the lubricant to be forced out. The surface is kept separated by the pressure caused by the fluid’s viscosity [22]. This sort of lubrication is best for heavy loads that are only used for a short time. This device has the ability to carry heavy loads for brief periods of time [22]. The layer of fluid lubricant thins when the fluid is driven out, and the joint surfaces come into contact. Under particular conditions, the film’s height may be measured, and the time it takes for the film to shrink from a given height to zero can be calculated [22]. A squeeze film mechanism of lubrication was suggested by Dowson and Popov in 1966 which based on assumption that the analysis concerning a loaded rigid cylinder approaching a rigid plate [23]. A rigid plane opposes a rigid bone covered by a porous, elastic layer of articular cartilage lubricated by synovial fluid with relative velocities are both sliding and squeezing [23]. Dowson used the primary results of hydrodynamic and EHD theories in a simple model for human joints for the first time in tribology history [22] [23]. Theoretical and experimental values of film thickness were found to be in excellent agreement. Higginson and Norman (1974) and Higginson (1977) looked at a pair of compliant surfaces where the lubricant was squeezed out. However, the non-Newtonian behaviour of the synovial fluid was not taken into account in the computations [24]. An asymptotic analysis of a lubrication problem for a model of articular cartilage and synovial fluid under squeeze-film conditions, was proposed an original analytical approximate model for the synovial pressure field determination in the ankle joint in a pure squeeze motion, accounting for the non-Newtonian behaviour of synovial fluid and cartilage porosity [23] [24].

2. The Rheology of Synovial Fluid

● Power Law Model:

Ostwald-de Waele relationship which is called power law model can be related to viscosity and shear rate in a steady shear flow,

$\begin{array}{l}\mu =\mathbb{K}{\left(\stackrel{\dot{}}{\gamma }\right)}^{n-1}\\ \mu =\mathbb{K}{\left(\frac{\partial u}{\partial y}\right)}^{n-1}\mathrm{,}\end{array}$ (1)

where $\mu$ is viscosity, $\stackrel{\dot{}}{\gamma }$ shear rate, $\mathbb{K}$ and n are cure fitting coefficients [25].

● Newtonian Fluid in Generalized Form:

A general type of Newtonian fluid for shear thinning can be described by Ostwald-de Waele relationship which satisfies the Rheological equation,

$\begin{array}{l}\tau =\left(\stackrel{\dot{}}{\gamma }\right)\mu \left(\stackrel{\dot{}}{\gamma }\right),\\ \tau =\mathbb{K}{\left(\stackrel{\dot{}}{\gamma }\right)}^n.\end{array}$ (2)

An idealized fluid is generalized Newtonian fluid in which shear stress is a function of shear rate at a particular time whereas independent of deformation history [25]. With the help of power law for fluid, equation of Rheology becomes,

$\tau =\mathbb{K}{\left(\frac{\partial u}{\partial y}\right)}^n\mathrm{,}$ (3)

● Stokes Coupled Stress Fluid:

Stokes theory (1966) is used to model synovial fluid to build up link with the momentum and continuity equations [25]. According to this theory, defined equations are,

$\begin{array}{l}\rho \frac{D\mathbb{V}}{Dt}=-\nabla \mathbb{P}+\mu \nabla ℂ,\\ \mathbb{Q}=\rho \mathbb{F}+\frac{1}{2}\rho \nabla \times ℂ-\eta {\nabla }^4\mathbb{V},\\ \nabla \cdot \mathbb{V}=0.\end{array}$ (4)

where $\rho$ is density of synovial fluid, vectors $\mathbb{V}$, $\mathbb{F}$ and $ℂ$ represent velocity, body force and body couple per unit mass. $\mu$ is viscosity and $\eta$ represents coupled stress constant [25].

3. Kinematics

In term of mathematical equations for fluid flow is the conservation laws of physics.

● The mass of a fluid is conserved (continuity equation).

● The rate of change of momentum equals the sum of the forces on a fluid particle (Newton’s 2nd law of motion) [26].

Steady, two-dimensional mass conservation or continuity equation at a point in a compressible fluid is,

$\frac{\partial \rho }{\partial t}+\text{div}\left(\rho \mathbb{V}\right)=\mathrm{0,}$

where the first term is the rate of change of density with respect to time while second term represents the net flow of mass out of the element across its boundaries and is called the convective term. Our assumption in this article is that the density is constant (incompressible flow), so above equation can be written as,

$\begin{array}{l}\text{div}\left(\mathbb{V}\right)=\mathrm{0,}\\ \Rightarrow \frac{\partial e_I}{\partial x}+\frac{\partial e_J}{\partial y}=0.\end{array}$ (5)

Conservation of momentum equation in x-direction is defined as,

$\frac{\partial e_I}{\partial t}+e_I\frac{\partial e_I}{\partial x}+e_J\frac{\partial e_I}{\partial y}=\frac{\partial }{\partial x}\left(-\mathbb{P}+{\tau }_{xx}\right)+\frac{\partial }{\partial y}\left({\tau }_{yx}\right)+S_{M_x},$

Conservation of momentum equation in y-direction is defined as,

$\frac{\partial e_J}{\partial t}+e_I\frac{\partial e_J}{\partial x}+e_J\frac{\partial e_J}{\partial y}=\frac{\partial }{\partial x}\left({\tau }_{xy}\right)+\frac{\partial }{\partial y}\left(-\mathbb{P}+{\tau }_{yy}\right)+S_{M_y},$

where ${\tau }_{xx}$ and ${\tau }_{yy}$ are viscous stress components which is composed of the linear deformation rate and the volumetric deformation rate. Liquids are incompressible so the mass conservation equation is continuity and viscous stresses are twice the local rate of linear deformation times the dynamic viscosity [26] [27].

3.1. Assumptions on Momentum Equation

● Steady, two dimensional and incompressible.

● $e_I$ is x component of the seepage velocity in fluid film region.

● ${\bar{e}}_I$ is average velocity.

● $S_{M_x}$ and $S_{M_y}$ are source terms include contribution due to body forces only and $\mathbb{P}$ is pressure. We consider $S_{M_x}=S_{M_y}=0$ and linear deformation rate.

● Neglect the variation of pressure normal to very thin film of lubrication.

● n is power law index with $n=1$.

Above assumptions led to the following expression,

$0=-\frac{\partial \mathbb{P}}{\partial x}+\mu \frac{\partial }{\partial y}\left[{\left|\frac{\partial e_I}{\partial y}\right|}^{n-1}\frac{\partial e_I}{\partial y}\right].$ (6)

$0=-\frac{\partial \mathbb{P}}{\partial y}.$ (7)

Let us consider a special slip boundary conditions at lower and upper permeable surfaces.

$e_I=-\frac{\sqrt{\mathbb{K}}}{\alpha }\frac{\partial e_I}{\partial y},\text{at}y=h,$ (8)

where $\mathbb{K}$ and $\alpha$ are permeability of porous medium and slip parameter respectively [26] [27]. The boundary of porous medium with a free viscous fluid is the continuity of pressure and the normal mass flux which give rise to wipe up the tangential velocity of the synovial fluid [28] [29]. Such boundary condition was proposed by Beavers and Joseph (BJ) in 1967 as a empirical formula. The analysis showed that BJ condition is macroscopic law which describes a phenomena at a macroscopic scale whereas characteristic length is much larger than the pore characteristic size [28] [29]. Similarity transformation is useful to scale out the problem.

$e_I=0,\text{at}y=0.$ (9)

Convection diffusion phenomena has great potential for applications in mathematical biology as well as pattern formation. The hyaluronic acid (HA) molecules in synovial fluid linked by convection diffusion phenomena is simple mass balance between changes in concentration $ℂ$ of solutes [28] [29].

$\frac{\partial ℂ}{\partial t}+\left(e_I-{\bar{e}}_I\right)\frac{\partial ℂ}{\partial x}+\left(e_J-{\bar{e}}_J\right)\frac{\partial ℂ}{\partial y}={\mathbb{D}}_{Diff}\left(\frac{{\partial }^2ℂ}{\partial x^2}+\frac{{\partial }^2ℂ}{\partial y^2}\right)-{ℝ}_1ℂ\mathrm{,}$ (10)

where $ℂ$, ${ℝ}_1$ and ${\mathbb{D}}_{Diff}$ are concentration of solutes, chemical reaction with first order rate (such reaction occurs in fully developed flow after diffusion) and diffusion constant respectively.

3.2. Assumptions on Concentration Equation

● Transverse diffusion is very much higher than longitudinal diffusion.

$\Rightarrow \frac{{\partial }^2ℂ}{\partial x^2}\ll \frac{{\partial }^2ℂ}{\partial y^2}\mathrm{,}$

● Steady convection diffusion phenomena.

We apply the following boundary conditions on Equation (10),

$ℂ=c_0,\text{at}y=h,$ (11)

$ℂ=c_1,\text{at}y=-h.$ (12)

4. Governing System for Synovial Fluid

The physical configuration to proposed problem as incompressible synovial fluid flowing two parallel surfaces which is apart by 2h separated distance [28]. The synovial fluid has been represented by viscoelastic fluid and its elasticity is pertinent in joint lubrication and bounded by porous layer [28].

$\begin{array}{l}0=-\frac{\partial \mathbb{P}}{\partial x}+\mu \frac{{\partial }^2{e}_1}{\partial y^2},\\ 0=-\frac{\partial \mathbb{P}}{\partial y}\end{array}\}{e}_I=-\frac{\sqrt{\mathbb{K}}}{\alpha }\frac{\partial e_I}{\partial y},\text{at}y=h\&e_I=0,\text{at}y=0.$ (13)

$\left(e_I-{\bar{e}}_I\right)\frac{\partial ℂ}{\partial x}={\mathbb{D}}_{Diff}\left(\frac{{\partial }^2ℂ}{\partial x^2}\right)-{ℝ}_1ℂ\}ℂ=c_0,\text{at}y=h\&ℂ=c_1,\text{at}y=-h.$ (14)

5. Methodology

Partial differential equations occur in almost every branch of science and engineering. Exact analytical solutions may be obtained only for very few cases among them. In most cases one has to adopt some numerical method to solve the problems. As we know that the momentum and continuity equations were first time solved by German engineer H. Blasius, which is done by transforming the partial differential equations into ordinary differential equations [30]. He introduced a new independent variable which is called similarity variable (nondimensionalization) [30].

5.1. Nondimensionalization

5.1.1. Velocity Distribution Profile

For non-dimensionalization technique, let us consider $\eta =\frac{y}{h}$ which is accommodated in Equation (13) [30] [31] [32].

$\frac{1}{\mu }\frac{\partial \mathbb{P}}{\partial x}=\frac{1}{h^2}\frac{{\partial }^2{e}_1}{\partial {\eta }^2}\}{e}_I=-\frac{\sqrt{\mathbb{K}}}{\alpha h}\frac{\partial e_I}{\partial \eta },\text{at}\eta =1\&e_I=0,\text{at}\eta =0.$ (15)

Solution to problem (15) is,

$\begin{array}{l}{e}_1=\frac{h^2}{2\mu }{\eta }^2\frac{\text{d}\mathbb{P}}{\text{d}x}+A\eta +B,\\ B=0,\\ A=-\frac{\text{d}\mathbb{P}}{\text{d}x}\frac{\frac{h^2}{2\mu }+\frac{\sqrt{\mathbb{K}}h}{\alpha \mu }}{1+\frac{\sqrt{\mathbb{K}}}{\alpha h}}\end{array}\}\text{where}A\text{and}B\text{are}\text{constants}\text{of}\text{integration}.$ (16)

$e_1=\frac{h^2}{2\mu }{\eta }^2\frac{\text{d}\mathbb{P}}{\text{d}x}-\eta \frac{\text{d}\mathbb{P}}{\text{d}x}\frac{\frac{h^2}{2\mu }+\frac{\sqrt{\mathbb{K}}h}{\alpha \mu }}{1+\frac{\sqrt{\mathbb{K}}}{\alpha h}}$ (17)

As synovial fluid velocity is distributed randomly in the cartilage, so we can find the average velocity by introcing the following expression,

${\bar{e}}_1=\frac{1}{2h}{\displaystyle {\int }_{-h}^h}{e}_1\left(y\right)\text{d}y$ (18)

Now we try to find $e_1^{\star }$ which can be used in concentration equation.

$\begin{array}{l}{e}_1^{\star }=\frac{e_1-{\bar{e}}_1}{{\bar{e}}_1},\\ e_1^{\star }=6\mu \left[\frac{{\eta }^2}{2\mu }-\eta \frac{\frac{1}{2\mu }+\frac{\sqrt{\mathbb{K}}}{\alpha \mu h}}{1+\frac{\sqrt{\mathbb{K}}}{\alpha h}}\right]-1,\\ e_1^{\star }=3{\eta }^2-6\eta \beta -1.\end{array}\}\text{where},\beta =\frac{\frac{1}{2}+\sigma }{1+\sigma }\text{and}\sigma =\frac{\sqrt{\mathbb{K}}}{\alpha h}.$ (19)

5.1.2. Concentration Distribution Profile

Apply non-dimensionalization technique to Equation (14) [30] [31] [32], such as

${ℂ}^{\star }=\frac{c}{c_0},y^{\star }=\frac{y}{h},t^{\star }=\frac{t}{\bar{t}},\xi =\frac{x-{\bar{e}}_{1}t}{L},{\beta }_1^2=\frac{h^2{ℝ}_1}{D_{Diff}}$ (20)

which led to the following equation,

$\begin{array}{l}\begin{array}{l}\frac{h^2}{D_{Diff}}\left(\frac{D_{Diff}{c}_0}{h^2}\right)\frac{{\partial }^2{ℂ}^{\star }}{\partial y^{\star }{}^2}-\frac{h^2}{D_{Diff}}{c}_0{ℝ}_1{ℂ}^{\star }=\frac{h^2}{L{D}_{Diff}}{c}_0{e}_I^{\star }\frac{{ℂ}^{\star }}{\partial \xi },\\ \frac{{\partial }^2ℂ}{\partial y^2}-{\beta }_1^2ℂ=\mathbb{Q}{e}_I^{\star },\\ \mathbb{Q}=\frac{h^2}{L{D}_{Diff}}\frac{ℂ}{\partial \xi }\end{array}\}\\ ℂ=c_0\text{at}y=1\text{and}ℂ=c_1,\text{at}y=-1.\end{array}$ (21)

Solution to problem (21) is,

$\begin{array}{l}\begin{array}{l}ℂ=-\frac{1}{{\beta }_1^2}\mathbb{Q}{e}_I^{\star }+A_1\cosh {\beta }_{1}y+B_1\sinh {\beta }_{1}y,\\ A_1=\frac{1}{\cosh {\beta }_1}\left[\frac{c_0+c_1}{2}+\frac{\mathbb{Q}{e}_I^{\star }}{{\beta }_1^2}\right],\\ B_1=\frac{1}{\sinh {\beta }_1}\left[\frac{c_0-c_1}{2}\right],\end{array}\}\\ \text{where}{A}_1\text{and}{B}_1\text{are}\text{constants}\text{of}\text{integration}\mathrm{.}\end{array}$ (22)

$ℂ=\frac{1}{\cosh {\beta }_1}\left[\frac{c_0+c_1}{2}+\frac{\mathbb{Q}{e}_I^{\star }}{{\beta }_1^2}\right]\cosh {\beta }_{1}y+\frac{1}{\sinh {\beta }_1}\left[\frac{c_0-c_1}{2}\right]\sinh {\beta }_{1}y-\frac{1}{{\beta }_1^2}\mathbb{Q}{e}_I^{\star },$ (23)

5.2. Finite Volume Method

Computational methods like finite difference, finite volume, finite element or their variants, that are most frequently used in computational fluid dynamics. The present work emphasises mainly the finite volume method for computational solution to aforementioned problems (13) & (14). Implementation steps are as follows:

● Divide the domain into the finite sized subdomains (finite control volumes) and each subdomain is represented by a finite number of grid points (like Nodes) [33],

● integrate the governing differential Equation (GDE) over each subdomain [33],

● consider a profile assumption for the dependent variable (like, interpolation function) to evaluate the above integral which express the result as an algebraic quantity at the grid points [33].

${\displaystyle {\int }_{\Delta V}}{u}_{\eta \eta }\text{d}V+{\displaystyle {\int }_{\Delta V}}Sourc{e}_{\eta }\text{d}V=\mathrm{0,}$ (24)

above Equation (24) can be written as,

${\displaystyle {\int }_{\Delta V}}{d}_{\eta }\left(u_{\eta }\right)\text{d}V+{\displaystyle {\int }_{\Delta V}}Sourc{e}_{\eta }\text{d}V=0.$ (25)

which implies,

${\left[n{u}_{\eta }\right]}_n-{\left[A_s{u}_{\eta }\right]}_s+Sourc{e}_{\eta }A\delta \eta =\mathrm{0,}$ (26)

profile assumptions for $u_{\eta }$ at control volume faces n and s are as follows,

$SouthFaceFlu{x}_{\eta }=\frac{A_s\left(u_P-u_S\right)}{\delta {\eta }_{SP}}\mathrm{,}$ (27)

$NorthFaceFlu{x}_{\eta }=\frac{A_n\left(u_N-u_P\right)}{\delta {\eta }_{PN}}\mathrm{,}$ (28)

Combining Equations (26) $\to$ (28) to get the following system,

$\left(\frac{A_s}{\delta {\eta }_{SP}}+\frac{A_n}{\delta {\eta }_{PN}}-S_P\right)u_P=\frac{A_s}{\delta {\eta }_{SP}}{u}_S+\frac{A_n}{\delta {\eta }_{PN}}{u}_N+S_u$ (29)

5.2.1. Algorithm Implementation for Velocity Profile

Scheme procedure can be deduced from Equation (29), which is as follows:

Interior Boundary Nodes:

$\begin{array}{l}{a}_P{u}_P=a_S{u}_S+a_N{u}_N+S_u,\\ a_S=\frac{A_s}{\delta {\eta }_{SP}},\\ a_N=\frac{A_n}{\delta {\eta }_{PN}},\\ a_P=a_S+a_N-S_P.\end{array}\}\text{where}\delta {\eta }_{SP}=\delta {\eta }_{PN}=\delta \eta ,\&A_s=A_n=A$ (30)

Left Boundary Node: For first boundary node, we use $u=0$ at $\eta =0$.

$\frac{A}{\delta \eta }\left(u_N-u_P\right)-\frac{2A}{\delta \eta }\left(u_P-u_a\right)+Sourc{e}_{\eta }A\delta \eta =0$ (31)

$\begin{array}{l}{a}_P{u}_P=a_S{u}_S+a_N{u}_N+S_u,\\ a_S=0,\\ a_N=\frac{A}{\delta \eta },\\ a_P=a_S+a_N-S_P,\\ S_P=-\frac{2A}{\delta \eta },\\ S_u=\frac{2A}{\delta \eta }{u}_a+Sourc{e}_{\eta }A\delta \eta .\end{array}\}$ (32)

Right Boundary Node: For last (right) boundary node, we use Equation (15) at $\eta =1$.

$A\left(-\frac{\alpha h}{\sqrt{\mathbb{K}}}{u}_P\right)-\frac{A}{\delta \eta }\left(u_P-u_S\right)+Sourc{e}_{\eta }A\delta \eta =0$ (33)

$\begin{array}{l}{a}_P{u}_P=a_S{u}_S+a_N{u}_N+S_u,\\ a_S=\frac{A}{\delta \eta },\\ a_N=0,\\ a_P=a_S+a_N-S_P,\\ S_P=A\left(-\frac{\alpha h}{\sqrt{\mathbb{K}}}\right),\\ S_u=Sourc{e}_{\eta }A\delta \eta .\end{array}\}$ (34)

5.2.2. Algorithm Implementation for Concentration Profile

Finite volume method to Concentration Equation (14) is defined as;

${\displaystyle {\int }_{\Delta V}}{C}_{yy}\text{d}V+{\displaystyle {\int }_{\Delta V}}{\beta }_1^{2}C\text{d}V+{\displaystyle {\int }_{\Delta V}}Sourc{e}_C\text{d}V=\mathrm{0,}$ (35)

above Equation (35) can be written as,

${\displaystyle {\int }_{\Delta V}}{d}_y\left(C_y\right)\text{d}V-{\displaystyle {\int }_{\Delta V}}{\beta }_1^{2}C\text{d}V+{\displaystyle {\int }_{\Delta V}}Sourc{e}_C\text{d}V=\mathrm{0,}$ (36)

which implies,

${\left[n{C}_y\right]}_n-{\left[A_s{C}_y\right]}_s-{\beta }_1^{2}A\delta y{C}_P+Sourc{e}_{C}A\delta y=0.$ (37)

Profile assumptions for $C_y$ at control volume faces n and s are as follows,

$SouthFaceFlu{x}_y=\frac{A_s\left(C_P-C_S\right)}{\delta y_{SP}}\mathrm{,}$ (38)

$NorthFaceFlu{x}_y=\frac{A_n\left(C_N-C_P\right)}{\delta y_{PN}}\mathrm{,}$ (39)

Combining Equations (37) $\to$ (39) to get the following system,

$\left(\frac{A}{\delta y}+\frac{A}{\delta y}-S_P\right)C_P=\frac{A}{\delta y}{C}_S+\frac{A}{\delta y}{C}_N+S_c$ (40)

Scheme procedure can be deduced from Equation (40), which is as follows:

Interior Boundary Nodes:

$\begin{array}{l}{a}_P{C}_P=a_S{C}_S+a_N{C}_N+S_c,\\ a_S=\frac{A}{\delta y},\\ a_N=\frac{A}{\delta y},\\ a_P=a_S+a_N-S_P,\\ S_P=-{\beta }_1^{2}A\delta y,\\ S_c=ASourc{e}_C\delta y\end{array}\}$ (41)

Left Boundary Node: For first boundary node, we use $C=c_o$ at $y=1$.

$\frac{A}{\delta y}\left(C_N-C_P\right)-\frac{2A}{\delta y}\left(C_P-C_a\right)+Sourc{e}_{C}A\delta y-{\beta }_1^{2}A\delta y{C}_P=0$ (42)

which implies,

$\begin{array}{l}{a}_P{C}_P=a_S{C}_S+a_N{C}_N+S_c,\\ a_S=0,\\ a_N=\frac{A}{\delta y},\\ a_P=a_S+a_N-S_P,\\ S_P=-\frac{2A}{\delta y}-{\beta }_1^{2}A\delta y,\\ S_u=\frac{2A}{\delta y}{C}_a+Sourc{e}_{C}A\delta y.\end{array}\}$ (43)

Right Boundary Node: For last (right) boundary node, we use $C=c_1$ at $y=-1$.

$\frac{2A}{\delta y}\left(C_b-C_P\right)-\frac{A}{\delta y}\left(C_P-C_S\right)+Sourc{e}_{C}A\delta y-{\beta }_1^{2}A\delta y{C}_P=0$ (44)

which implies,

$\begin{array}{l}{a}_P{C}_P=a_S{C}_S+a_N{C}_N+S_c,\\ a_S=\frac{A}{\delta y},\\ a_N=0,\\ a_P=a_S+a_N-S_P,\\ S_P=-\frac{2A}{\delta y}-{\beta }_1^{2}A\delta y,\\ S_u=\frac{2A}{\delta y}{C}_b+Sourc{e}_{C}A\delta y.\end{array}\}$ (45)

6. Results

In this section, we discussed the part of numerical results of the weakly coupled problem through graphs and expressions which acquired from non-dimensionalization of velocity and concentration equations. These results are calculated numerically by utilizing MATLAB software. Figure 2 shows results for velocity profile which is deduced from Equation (29) by using Finite volume scheme. Figure 3 shows results for velocity profile error which deduced from Equation (29) by using Finite volume scheme at different grid settings. Figure 4 shows results for concentration profile which deduced from Equation (40) by using Finite volume scheme. Figure 5 shows results for concentration profile error which deduced from Equation (40) by using Finite volume scheme. Figure 6 shows results for velocity profile which deduced from Equation (29) by using Finite volume

Figure 2. Shows numerical and analytical results comparison for velocity Equations (28) & (30).

Figure 3. Shows error results for velocity Equation (30).

Figure 4. Shows numerical and analytical results comparison for concentration Equation (41).

Figure 5. Shows error results for concentration Equation (41).

Figure 6. Shows result for velocity Equation (13) with different values of Be = Bejian number.

Figure 7. Shows result for velocity Equation (13) with different values of $\alpha$.

Figure 8. Shows result for concentration Equation (14) with different values of ${\beta }_1$.

scheme at different $Be=\text{Bejan}\text{Number}=\frac{h^2}{\mu }\frac{\partial \mathbb{P}}{\partial x}$ with $\mu =1$, while Figure 7 shows results for velocity profile by using Finite volume scheme at different $\alpha$ values [34] [35]. Figure 8 shows results for concentration profile which deduced from Equation (29) by using Finite volume scheme at different ${\beta }_1$ values.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Walker, P.S., Dowson, D., Longfield, M.D. and Wright, V. (1968) “Boosted Lubrication” in Synovial Joints by Fluid Entrapment and Enrichment. Annuls of the Rheumatic Diseases, 27, 512-520. https://doi.org/10.1136/ard.27.6.512
[2]	Bali, R. and Shukla, A.K. (2000) Rheological Effects of Synovial Fluid on Nutritional Transport. Tribology Letters, 9, 233-239. https://doi.org/10.1023/A:1018825425258
[3]	Jin, Z.M., Dawson, D. and Fisher, J. (1992) The Effect of Porosity of Articular Cartilage on the Lubrication of a Normal Human Hip Joint. Proceeding of Institution of Mechanical Engineers H: Journal of Engineering and Medicine, 206, 117-124. https://doi.org/10.1243/PIME_PROC_1992_206_279_02
[4]	Shiroky, A., Alexander, V., Volkov, V. and Novochadov, V. (2014) Crucial Processes Interaction during the Renewal of Articular Cartilage: the Mathematical Modeling. Alexander European Journal of Molecular Biotechnology, 4, 86-94.
[5]	Videcoq, P., Steenkeste, K., Bonnina, E. and Garnier, C. (2013) A Multi Scale Study of Enzyme Diffusion in Macromolecular Solutions and Physical Gels of Pectin Polysaccharides. Soft Matter, 9, 5110-5118. https://doi.org/10.1039/c3sm00058c
[6]	Chandra, P. and Agarwal, R.P. (1983) Dispersion in Simple Microfluid Flows. International Journal of Engineering Science, 21, 431-442. https://doi.org/10.1016/0020-7225(83)90093-9
[7]	Taylor, G.I. (1953) Dispersion of Soluble Matter in Solvent Flowing Slowly through tube. Proceedings of the Royal Society of London. Series A, 219, 186-203. https://doi.org/10.1098/rspa.1953.0139
[8]	Aris, R. (1956) On the Dispersion of a Solute in Fluid Flow through a Tube. Proceedings of the Royal Society of London. Series A, 235, 67-77. https://doi.org/10.1098/rspa.1956.0065
[9]	Tandon, P.N. and Chaurasia, A. (1991) Application of Magnetic Fields to Synovial Joints. Computers & Mathematics with Applications, 22, 33-45. https://doi.org/10.1016/0898-1221(91)90145-T
[10]	MacConaill, M.A. (1932) The Function of Intra Articular Fibrocartilages, with Special Reference to the Knee and Inferior Radio-Ulnar Joints. Journal of Anatomy, 66, 210-227.
[11]	MacConaill, M.A. (1956) Studies in the Mechanics of Synovial Joints. Irish Journal of Medical Science, 31, 353-364. https://doi.org/10.1007/BF02951170
[12]	MacConaill, M.A. (1960) Lubrication of Mammalian Joints. Nature, 185, 920. https://doi.org/10.1038/185920a0
[13]	Jones, E.S. (1936) Joint Lubrication. Lancet, 227, 1043-1045. https://doi.org/10.1016/S0140-6736(01)37157-X
[14]	Charnley, J. (1959) The Lubrication of Animal Joints. Proceedings of IMechE Conference on. Biomechanics, London, 12-19. https://doi.org/10.1136/ard.19.1.10
[15]	Barnett, C.H. and Cobbold, A.F. (1962) Lubrication within Living Joints. The Journal of Bone and Joint Surgery, 44B, 662-674. https://doi.org/10.1302/0301-620X.44B3.662
[16]	Linn, F.C. and Radin, E.L. (1968) Lubrication of Animal Joints. III. The Effect of certain Chemical Alterations of the Cartilage and Lubricant. Arthritis & Rheumatism, 11, 674-682. https://doi.org/10.1002/art.1780110510
[17]	McCutchen, C.W. (1966) Paper 1: Physiological Lubrication. Proceedings of the Institution of Mechanical Engineers, Conference Proceedings, 181, 55-62. https://doi.org/10.1243/PIME_CONF_1966_181_207_02
[18]	McCutchen, C.W. (1973) A Note on Weeping Lubrication. In: Kenedi, R.M., Ed., Perspectives in Biomedical Engineering, Palgrave Macmillan, London, 109-110. https://doi.org/10.1007/978-1-349-01604-4_18
[19]	McCutchen, C.W. (1983) Lubrication of and by Articular Cartilage. In: Hall, B.K., Ed., Cartilage, Academic Press, Cambridge, 87-107. https://doi.org/10.1016/B978-0-12-319503-6.50009-8
[20]	Dintenfass, L. (1963) Lubrication in Synovial Joints. Nature, 197, 496-497. https://doi.org/10.1038/197496b0
[21]	Tanner, R.I. (1966) An Alternative Mechanism for the Lubrication of Synovial Joints. Physics in Medicine & Biology, 11, 119-127. https://doi.org/10.1088/0031-9155/11/1/312
[22]	Dowson, D. (1966) Paper 12: Modes of Lubrication in Human Joints. Proceedings of the Institution of Mechanical Engineers, Conference Proceedings, 181, 45-54. https://doi.org/10.1243/PIME_CONF_1966_181_206_02
[23]	Popova, E. and Popov, V.L. (2014) On the History of Elastohydrodynamics: The Dramatic Destiny of Alexander Mohrenstein-Ertel and His Contribution to the Theory and Practice of Lubrication. ZAMM: Journal of Applied Mathematics and Mechanics, 95, 652-663. https://doi.org/10.1002/zamm.201400050
[24]	Higginson, G.R. and Norman, R. (1974) The Lubrication of Porous Elastic Solids with Reference to the Functioning of Human Joints. Journal of Mechanical Engineering Science, 16, 250-257. https://doi.org/10.1243/JMES_JOUR_1974_016_045_02
[25]	Ruggiero, A. (2020) Milestones in Natural Lubrication of Synovial Joints. Frontiers in Mechanical Engineering, 6, Article No. 52. https://doi.org/10.3389/fmech.2020.00052
[26]	Versteeg, H.K. and Malalasekera, W. (2007) An Introduction to Computational Fluid Dynamics. 2nd Edition, Pearson, London, Chapter 03.
[27]	Löhner, R. (2008) Applied Computational Fluid Dynamics Techniques. John Wiley & Sons, Ltd., Hoboken, Chapter 02-06. https://doi.org/10.1002/9780470989746
[28]	Meenapriya, P. and Uma Maheswari, K. (2021) An Analytical Study of Aerosol Dispersion in the Atmosphere Bounded by Porous Layers. In: Rushi Kumar, B., Sivaraj, R. and Prakash, J., Eds., Advances in Fluid Dynamics, Springer, Singapore, 67-79. https://doi.org/10.1007/978-981-15-4308-1_5
[29]	Auriault, J.L. (2010) About the Beavers and Joseph Boundary Condition. Transport in Porous Media, 83, 257-266. https://doi.org/10.1007/s11242-009-9435-9
[30]	Ejaz, A., Abbas, I., Nawaz, Y., Shoaib Arif, M., Shatanawi, W. and Abbasi, J. (2020) Thermal Analysis of MHD Non-Newtonian Nanofluids over a Porous Media. Computer Modeling in Engineering & Sciences, 125, 1119-1134. https://doi.org/10.32604/cmes.2020.012091
[31]	Narsu, S., Kumar, R.B. (2018) Comparative Study of Chemically Reacting Blasius and Sakiadis Unsteady MHD Radiated Flow with Variable Conductivity. Journal of Physics: Conference Series, 1000, Article ID: 012149. https://doi.org/10.1088/1742-6596/1000/1/012149
[32]	Vijaya Kumar, R. and Ratchagar, N.P. (2021) Mathematical Modeling of Synovial Joints with Chemical Reaction. Journal of Physics: Conference Series, 1724, Article ID: 012051. https://doi.org/10.1088/1742-6596/1724/1/012051
[33]	Versteeg, H.K. and Malalasekera, W. (2007) An Introduction to Computational Fluid Dynamics. 2nd Edition, Pearson, London, Chapter 04 & 05.
[34]	Ascher, U.M. and Ruuth, S.J. (1995) Implicit Explicit Methods for Time Dependent Partial Differential Equations. SIAM Journal on Numerical Analysis, 32, 797-823. https://doi.org/10.1137/0732037
[35]	Li, J. (2008) Computational Partial Differential Equations Using MATLAB. Chapman and Hall/CRC, New York. https://doi.org/10.1201/9781420089059