A Model of Synovial Fluid with a Hyaluronic Acid Source: A Numerical Challenge

Abstract

Initially motivated by the analysis of the flow dynamics of the synovial fluid, taken as non-Newtonian, this paper also reports on a numerical challenge which occurred unexpectedly while solving the momentum equation of the model. The configuration consists of two infinitely long horizontal parallel flat plates where the top plate is sheared at constant speed and the bottom plate is fixed. The synovial fluid shows a shear-thinning rheology, and furthermore it thickens with the hyaluronic acid (HA) concentration, i.e., it is also chemically-thickening. Accordingly, a modified Cross model is employed to express the shear rate and concentration-dependent viscosity, whose parameter values are determined from experimental data. Another significance of the study is the investigation of the effect of an external stimulus on the flow dynamics via a HA source term. The resulting flow exhibits peculiar features resulting from extremely large and small, but positive, numerical quantities, such as the viscosity and the shear rates. This requires constructing a parametrized zero-machine level solver, up to 300 accurate digits or so, for capturing the correct length scales of the flow physics. As a conclusion, the physical model, although simple, but original, leads to interesting results whose numerical determination turns out to be successful only once the real cause of the numerical trap is identified.

1. Introduction

The synovial joints and, as a result, the synovial fluid (SF) are responsible for the mobility of the skeletal system. The SF resides in the joints, and it is enclosed by the tissues surrounding the joint space, namely, the synovium [1]. The SF is responsible for the viscous, elastic, and lubricant properties in the joints [2], and its viscoelastic properties are measured [3]. King [4] determined the normal stresses for SF. Distortions in these mentioned SF properties may result in various joint dysfunctions such as osteoarthritis and rheumatoid arthritis. There are some tribological modeling studies in the presence of synovial fluid with the aim of total hip replacement [5]. The SF can be basically approximated as blood plasma, except that it is free from large proteins and enriched with hyaluronic acid (HA) molecules, which are synthesized locally in the joints [6]. The response of the SF to instantaneous external stimuli is non-Newtonian when it is healthy, whereas the response is Newtonian when it is inflamed. The presence of HA is the main reason behind that behavior, and the response is affected by the concentration and the molecular weight of the HA in the joints [6,7]. The normal HA concentration of a healthy joint is approximately between 2 and 4 mg/mL [1]. The compositions of the healthy and of the unhealthy SF are in Table 1 in [1]. A review of the rheological properties of the SF in terms of the HA concentration, shear rate dependency, and elasticity was done in [8]. More recent reviews of the SF rheology focusing on various pathologies can be found in [1,9]. Tandon et al. [10] modeled the SF as a Bingham fluid. However, most studies used a power law model. Rudraiah et al. [11] modeled the SF as a power law shear-thinning fluid with constant consistency and flow behavior indices. Hron et al. [6] employed two different power law equations for modeling SF: one where a function of the HA concentration multiplies the consistency index, and a second one where the flow index is concentration dependent. In this second model, the HA concentration influences the shear-thinning behavior of the SF, which is also observed in the experiments [7]. Fam et al. [7] modeled the viscosity of the SF with a modified Cross model [12]. A summary of the theoretical approaches on the lubrication modeling of the synovial joints is given in [13].

The Newtonian SF within an artificial hip joint is modeled between two parallel plates where one plate is oscillated [14]. An analytical solution for the unsteady flow was obtained via Duhamel’s theorem. Pekkan et al. [15] also considered oscillatory flow of a Newtonian SF first between infinite parallel plates, then in two-dimensional joint geometries, and including a geometry with cartilage curvature. The effects of magnetic or acoustic fields on the SF are also studied for practical reasons. Khan et al. [16] considered the SF either shear-thinning or shear-thickening and subject to magnetic field while it is on peristaltic motion. Romanishina et al. [17] assumed the SF as a two-phase fluid with solid and liquid phases. VijayaKumar and Ratchagar [18] modeled lubricating SF between parallel plates in the presence of a reaction.

Modeling biofluids is very complicated as there are many changes that take place. Consequently, one method of modeling a biofluid is to consider it as a chemically reacting fluid and focus on the concentration of the compound-of-interest [6,19,20,21,22], e.g., hyaluronic acid. Bridges et al. [20] studied the pulsating synovial flow in an annular region between two cylinders. They assumed that the viscosity depends on the concentration of HA to represent the shear-thinning and chemically-thickening behavior of the SF, adopting the model presented in [6].

In this work, the SF is placed between two parallel plates, with one of the plates being driven at constant speed. The SF is assumed to exhibit both shear-thinning and chemically-thickening behaviors. Consequently, a modified version of the Cross model is used to show the rheology of the SF. A source term, which represents the response of the SF to a stimulus, is employed in the model. The source, i.e., the increase in the acid concentration, is given with a linear function. The resulting nonlinear differential equations are solved using a Chebyshev spectral collocation (CSC) [23,24,25] method. Although simple, this model leads to very peculiar physical behaviors. In particular, as the HA concentration increases, the viscosity reaches extremely high values leading to almost constant velocity in the core and very large gradients near the plates. The associated strain rate must anyhow be positive across the flow domain, however small it is. Capturing this behavior with numerics turns out to be very challenging. It will then be seen in this paper that the physical reliability of the numerical results is controlled by employing the necessary number of accurate digits, which cannot be anticipated before facing the resulting issue. In other words, for some parameters, any double-precision solver fails to provide acceptable results, and leads to nonphysical behaviors.

There are examples of calculations where there is a need for higher accuracy than double precision [26,27,28,29,30]. Khanna [31] adjusted the number of digits for solving a hyperbolic differential equation. Historically, the use of accurate digits has increased from single to double precision with 64-bit computers. After the birth of 128-bit computers, it would not be surprising that the quadruple precision will be the standard use. To detect the physical phenomena using a numerical tool, sometimes there is a need for more than double precision. Consequently, there are some software like Advanpix and Multiprecision Computing Toolbox for Matlab (The MathWorks Inc., Massachusetts, USA) that are currently in use. Here, we present a practical problem, synovial fluid flow, where a parametrized machine-zero level is employed to numerically capture the main features of the flow.

In Section 2, the model, solution method, and procedure are presented. In Section 3, the physical results are presented as functions of the strength of the source r. Next, in Section 4, the necessity to use controlled number of accurate digits is shown. Furthermore, the required number of digits is given as a function of r. Section 5 presents the conclusion. Appendix A presents the numerical treatment of the problem.

2. Model

The physical system is depicted in Figure 1. It consists of two infinitely long parallel plates, with the top plate moving at constant speed, $U_w$, in its plane and the bottom one being at rest. First, the governing equations are presented, then the viscosity model is introduced.

2.1. Equations

The governing equations express the balances of momentum and concentration as

$$\rho \frac{\partial u}{\partial t}=\frac{\partial }{\partial y}\left[\eta \left(\frac{\partial u}{\partial y},c\right)\frac{\partial u}{\partial y}\right],$$

(1)

and

$$\frac{\partial c}{\partial t}=D\frac{{\partial }^{2}c}{\partial y^2}+rg\left(c\right),$$

(2)

where u and c represent the velocity in the x-direction ($\mathrm{m}/\mathrm{s}$) and the mass concentration ($\mathrm{kg}/{\mathrm{m}}^3$), respectively. The diffusion coefficient D (${\mathrm{m}}^2/\mathrm{s}$) and the density $\rho$ ($\mathrm{kg}/{\mathrm{m}}^3$) are assumed to be constant. Moreover, $g\left(c\right)$ represents a source term and r is the strength of that term, a positive (zero if there is no source) control variable (see in [25], p. 173, for a similar treatment in a different problem, namely, ignition in a solid). The unit of r depends on the choice of the function $g\left(c\right)$. In addition, $\eta$ is the viscosity ($\mathrm{kg}/\mathrm{m}/\mathrm{s}$), which depends on the strain rate as well as the concentration.

The problem is rendered dimensionless using the following scales for the length, velocity, time, and concentration: h, $U_w$, $h/U_w$, and $c_{jel}$ = 20 mg/mL, which is the jellification concentration of the SF [6]. The dimensionless equation for the momentum balance is

$$Re\frac{\partial \tilde{u}}{\partial \tilde{t}}=\frac{\partial }{\partial \tilde{y}}\left(\tilde{\eta }\frac{\partial \tilde{u}}{\partial \tilde{y}}\right).$$

(3)

The tilde represents the scaled variable. Here, the Reynolds number $Re$ is defined as

$$Re=\frac{\rho U_{w}h}{{\beta }_1},$$

(4)

where ${\beta }_1$, the zero-shear viscosity of the SF for $\tilde{c}$ = 0, is given in Section 2.2. The dimensionless concentration equation reads

$$Pe\frac{\partial \tilde{c}}{\partial \tilde{t}}=\frac{{\partial }^2\tilde{c}}{\partial {\tilde{y}}^2}+\tilde{r}\tilde{g}\left(\tilde{c}\right).$$

(5)

Here, $\tilde{g}$ is the dimensionless source term and the Péclet number $Pe$ is defined as

$$Pe=\frac{U_{w}h}{D}.$$

(6)

The dimensionless numbers $Re$ and $Pe$ are 10 and ${10}^3$, respectively [6]. The above equations are subject to the following boundary conditions:

$$\tilde{u}\left(-1,\tilde{t}\right)=0\tilde{u}\left(1,\tilde{t}\right)=1$$

(7)

and

$$\tilde{c}\left(-1,\tilde{t}\right)=0.15\tilde{c}\left(1,\tilde{t}\right)=0.16.$$

(8)

The system of equations, i.e., Equations (3) and (5), will be considered with the source term $\tilde{g}$ taken as linear, i.e., with $\tilde{g}=\tilde{c}$.

2.2. Viscosity Model

The viscosity of the SF is assumed to follow the modified Cross model [7], given as

$$\eta =\frac{{\eta }_0}{1+{\left(\lambda \frac{\partial u}{\partial y}\right)}^{\left(1-n\right)}}$$

(9)

where $\frac{\partial u}{\partial y}$ represents the strain rate, ${\eta }_0$ (kg/m/s) the zero-shear viscosity, $\lambda$ (s) the time constant, and n (dimensionless) the flow behavior index. It is known from experimental data that the viscosity of the SF depends on the HA concentration. Consequently, this work incorporates, as it was done in [6] in their viscosity model, this dependency into Equation (9) through the above-mentioned coefficients.

The viscosity is of course a positive quantity, which is satisfied by the value of ${\eta }_0$ given in Equation (9), provided that the denominator is positive (see Section 2.4). To fix the dependency of the parameters, Equation (9) is fitted to the experimental data wherein the viscosity is measured as a function of the shear rate at various constant concentration values [1]. Therefore, for each constant concentration dataset, a set of ${\eta }_0$, $\lambda$, and n is determined. This led to

$${\eta }_0\left(\tilde{c}\right)={\beta }_{1}exp\left(\widetilde{{\beta }_2}\tilde{c}\right),$$

(10)

$$\tilde{\lambda }\left(\tilde{c}\right)=\widetilde{{\beta }_3}{\tilde{c}}^{{\beta }_4},$$

(11)

$$n\left(\tilde{c}\right)=-\widetilde{{\beta }_5}\tilde{c}+{\beta }_6.$$

(12)

Here, $\beta$′s are constant coefficients. Note that all parameters and dependencies in Equations (10)–(12), except for ${\beta }_1$ and ${\eta }_0$, are dimensionless. The dimensionless form of the viscosity model reads as

$$\tilde{\eta }=\frac{\eta }{{\beta }_1}=\frac{exp\left(\widetilde{{\beta }_2}\tilde{c}\right)}{1+{\left(\tilde{\lambda }\frac{\partial \tilde{u}}{\partial \tilde{y}}\right)}^{(1-n)}}.$$

(13)

The coefficients are determined using Matlab’s lsqcurvefit function by fitting Equation (9) to the concentration data as

$$({\beta }_1,\widetilde{{\beta }_2},\widetilde{{\beta }_3},{\beta }_4,\widetilde{{\beta }_5},{\beta }_6)=(0.0114\frac{\mathrm{kg}}{\mathrm{m}\mathrm{s}},25.9131,0.0105,3.16,0.9724,0.5216).$$

The fitted curves together with the experimental data are shown for four values of the concentration in Figure 2.

2.3. Analysis of the Model: Steady State

The system of Equations (3) and (5) is coupled by the $\tilde{c}$ field. It is not a strong coupling as the concentration evolution is not affected by the momentum dynamics. Based on the comments made in Section 2.4, it can be concluded that this system can only converge to a steady solution, at any value of $\tilde{r}$. We can then decide that $\tilde{c}$ is fixed to its steady distribution with $\tilde{y}$, i.e., it satisfies the steady version of Equation (5),

$$\frac{{\partial }^2\tilde{c}}{\partial {\tilde{y}}^2}+\tilde{r}\tilde{g}\left(\tilde{c}\right)=0.$$

(14)

The analytical solution is given by

$$\tilde{c}\left(\tilde{y}\right)=\left\{\begin{array}{cc}0.005\tilde{y}+0.155,\hfill & \tilde{r}=0\hfill \\ \frac{0.155}{\cos \left(\sqrt{\tilde{r}}\right)}\cos \left(\sqrt{\tilde{r}}\tilde{y}\right)+\frac{0.005}{\sin \left(\sqrt{\tilde{r}}\right)}\sin \left(\sqrt{\tilde{r}}\tilde{y}\right),\hfill & \tilde{r}>0.\hfill \end{array}\right.$$

(15)

The concentration $\tilde{c}\left(\tilde{y}\right)$ has to be positive in the flow domain. This imposes that the parameter $\tilde{r}$ is bounded by $\frac{{\pi }^2}{4}\approx 2.46$. Beyond this value, the solution for the concentration is no longer positive over the entire domain. We have thus limited the $\tilde{r}$ value to $0\le \tilde{r}<\frac{{\pi }^2}{4}$. The steady-state concentration field is henceforth inserted into the momentum equation, Equation (3).

The steady-state solution of the momentum balance equation, Equation (3), has to verify the following relationship:

$$\tilde{\eta }\left(\frac{\partial \tilde{u}}{\partial \tilde{y}},\tilde{c}\right)\frac{\partial \tilde{u}}{\partial \tilde{y}}=\tilde{K},$$

(16)

where $\tilde{K}$, the dimensionless shear stress, is constant for a given $\tilde{r}$. As the upper wall is sheared in the positive x-direction and the viscosity is positive, $\tilde{K}$ is positive everywhere. With the constant $\tilde{K}$ being positive, the strain rate $\frac{\partial \tilde{u}}{\partial \tilde{y}}$ must also be strictly positive in the entire flow domain.

2.4. Stability Consideration

The momentum balance Equation (3) is nonlinear in $\tilde{u}$, but only through the viscosity which has to be always positive. Therefore, the right-hand-side of Equation (3) has only negative eigenvalues. In other words, it is a pure diffusion problem whose solution has to be steady as the boundary conditions are time-independent. No physical instability nor sensitivity to the initial condition can then be expected from this configuration, provided the initial condition satisfies the boundary conditions together with $\frac{\partial \tilde{u}\left(\tilde{y},\tilde{t}=0\right)}{\partial \tilde{y}}>0$. As will be seen in Section 3, this strain rate for large values of $\tilde{r}$ may become vanishingly small, e.g., reaching ${10}^{-300}$ for $\tilde{r}=2.45$. However, it is not zero, and it should not be zero; otherwise, the physical quantities like velocity or stress are not computed correctly.

3. Results

The unsteady momentum balance Equation (3) is solved by using the CSC spatial discretization and by marching in time, see Appendix A, until the criterion given in Equation (A4) is satisfied. It is proceeded with increasing values of $\tilde{r}$, taking the same initial condition, i.e., $\tilde{u}(\tilde{y},\tilde{t}=0)=(\tilde{y}+1)/2$, and the same time-step size, $\delta \tilde{t}=1$.

The physical results are presented in terms of three dimensionless variables: the concentration $\tilde{c}$ in Figure 3, the velocity $\tilde{u}$ in Figure 4, and the viscosity $\tilde{\eta }$ in Figure 5. All variables are plotted as functions of position $\tilde{y}$ for five values of $\tilde{r}$, the strength of the hyaluronic acid (HA) source term, viz., $\tilde{r}=[0,0.5,1.5,2.2,2.45]$. The concentrations and the viscosities are given as scaled by their values at $\tilde{y}=0$, which are $\tilde{c}(\tilde{y}=0)=[0.155,0.204,0.457,2.86,27.9]$, and $\tilde{\eta }(\tilde{y}=0)=[55.5,1.96\times {10}^2,1.39\times {10}^5,1.54\times {10}^{32},2.35\times {10}^{314}]$, respectively.

The steady concentration profiles are obtained analytically via Equation (15) and are plotted in Figure 3. The linear concentration profile, obtained for $\tilde{r}=0$, originates from the different values of the concentration on the walls as given in Equation (8). For $\tilde{r}>0$, the concentrations present an almost parabolic shape whose maximum is slightly shifted from $\tilde{y}=0$. As $\tilde{r}$ increases, the concentration amplitude also increases, but only moderately.

Increasing $\tilde{r}$ results in a dramatic change in the velocity profile near the walls, as seen in Figure 4. This behavior of the velocity profile is due to the shear-thinning and chemically-thickening viscosity. Note that the velocity profile for $\tilde{r}=0$ is close to a linear (Newtonian) one, but it is still a nonlinear profile. As $\tilde{r}$ increases, the velocity becomes nearly constant in the core of the domain with very sharp gradients close to the plates. For $\tilde{r}=2.45$, which is just shy of the critical value of ${\pi }^2/4$ (see Equation (15)), the velocity magnitude shows almost a step change going from the boundary to the core. In other words, the flow evolves from an almost constant-viscosity creeping flow for $\tilde{r}=0$ to a plug flow for $\tilde{r}\ge 1.5$.

The increase in the HA concentration affects the center of the domain the most, and as a result the viscosity around the center becomes very large, whereas near the walls it is relatively small (Figure 5). Increasing the concentration has a major effect on the viscosity, mainly due to the numerator of the Cross equation, Equation (9), i.e., ${\eta }_0$, which has an exponential dependency. This leads to a significant contrast in viscosity, of almost Gaussian type, not only in amplitude, but also in its localization around $\tilde{y}=0$.

The final quantity of interest is the dimensionless stress, $\tilde{K}$. First of all, it is checked that $\tilde{K}$ is constant with $\tilde{y}$ for all values of $\tilde{r}$. An example for $\tilde{r}=2.2$ is given in Figure 6. Then, Figure 7 shows the strain rate (blue dots) and the viscosity (orange dots) at $\tilde{y}=0$ as functions of $\tilde{r}$. Their product leads to $\tilde{K}$, defined in Equation (16), and it is superimposed in Figure 7b (green dots). Figure 7a offers a zoom of the stress $\tilde{K}$ plot as a function of $\tilde{r}$. It shows an exponential increase with $\tilde{r}$, while not exhibiting as drastic of an evolution as the strain rate and the viscosity do with $\tilde{r}$. The viscosity contrast leads this pure momentum diffusion configuration to exhibit a boundary layer behavior on the plates, with very thin layers when $\tilde{r}$ exceeds 2, together with an almost constant velocity in the flow core. This “almost constant” is a peculiarity and it has to be captured numerically.

4. Discussion about Numerics

The physical results presented in the previous section are in full agreement with the mathematical key points mentioned in the model analysis such as that the strain rate is always positive and the shear stress is always positive and constant.

The numerical determination of these results requires to adjust the grid refinement with the parameter $\tilde{r}$. Figure 8, obtained with double precision, shows the evolution of the minimum spectral cut-off N, i.e., the necessary number of grid points with $\tilde{r}$ (see Equation (A1)). However, satisfying this requirement is not sufficient and in fact not even the main criterion for obtaining the aforementioned results.

The stress $\tilde{K}$ is theoretically constant at any position (see Equation (16)), and it is obtained accordingly in Figure 6. It results from the product of two numerically sensitive quantities: the strain rate and the viscosity. As $\tilde{r}$ becomes larger than 2, the viscosity starts to increase exponentially at the centerline, i.e., $\tilde{y}=0$, whereas the velocity slope decreases exponentially, reaching amplitudes on the order of ${10}^{-300}$ (Figure 7). The solver must then be able to handle numerical values with an arbitrary number of accurate digits in order to capture this wide range of dynamics. This is the numerical challenge. Otherwise, the numerical scheme will get off from the mathematical key points of the model to be solved and would lead to nonphysical results identified by a non-constant $\tilde{K}$ or even a $\tilde{K}$ that blows up.

Going back to Figure 4, one notices that the $\frac{\partial \tilde{u}}{\partial \tilde{y}}$ data cover a wide range of numerical dynamics when $\tilde{r}$ becomes large. Actually, the ratio between the smallest and the largest slopes crosses the double-precision limit for $\tilde{r}$ of approximately $2.1$. Beyond this value of $\tilde{r}$, any double-precision solver loses numerical control of the system. The machine double precision has to be extended to a parametrized machine-zero level to have a correct digital representation of the slopes. Let us introduce the ratio, $N_d$,

$$N_d\left(\tilde{r}\right)={\log }_{10}\left(\frac{\max \left(\frac{\partial \tilde{u}}{\partial \tilde{y}}\right)}{{\left(\frac{\partial \tilde{u}}{\partial \tilde{y}}\right)}_{\tilde{y}=0}}\right).$$

(17)

It defines the minimum number of accurate digits required for numerically capturing the slopes for any value of $\tilde{r}$. It is plotted in Figure 9.

To this end, the solver of the SF model is coded in the framework of the software Mathematica 12 (from Wolfram Research), which allows the use of the number of accurate digits at any precision. All the physical results come from the time-marching approach with $\delta \tilde{t}=1$, and a time convergence is achieved after 15–20 time steps.

As cited in Section 1, there are many papers in the literature which need more than double precision for an accurate calculation. Even though the present physical system is simple, it shows a difficult trap, which cannot be solved with a refinement of the grids. This is a new and original configuration where an arbitrary number of accurate digits is required to get a solution for a diffusion problem.

5. Conclusions

The synovial fluid (SF) flow between two parallel plates is investigated in the presence of hyaluronic acid (HA) added as a source term. This source, which represents the response of the SF to a stimulus, is given as a linear function of the concentration. Experimental viscosity data, given for different values of the strain rate and HA concentration, are fitted to a modified Cross model whose parameters are taken as concentration dependent.

The evolution of the flow dynamics with the strength of the source term r becomes very striking, especially for large values of r. The viscosity becomes extremely large, even though the concentration magnitude increases only moderately. This results in distinct velocity profiles with sharp gradients near the plates and an almost constant level in the core. The dimensionless strain rates get as low as ${10}^{-300}$, and the dimensionless viscosities are as large as ${10}^{300}$, whose product leads, however, to a constant value of the stress all throughout the fluid, as it is dictated by the model. Correctly capturing these very small and very large numbers turns out to be essential, and in fact it is the only way to obtain physically meaningful results, including the stress. We show that this can only be achieved by controlling the machine-zero level, i.e., employing the necessary number of digits in the computation.