On the Modeling of the Synovial Fluid

Abstract

Synovial fluid is a polymeric liquid which generally behaves as a viscoelastic fluid due to the presence of hyaluronan molecules. We restrict ourselves to the regime in which the fluid responds as a viscous fluid. A novel generalized power-law fluid model is developed wherein the power-law exponent depends on the concentration of the hyaluronan. Such a model will be adequate to describe the flows of such fluids as long as they are not subjected to instantaneous stimuli. Assuming two different structures for the form of the power law exponent, both in keeping with physical expectations, we numerically solve for the flow of the synovial fluid (described by the constraint of incompressibility, the balance of linear momentum, and a convection-diffusion equation for the concentration of hyaluronan) in a rectangular cavity. The solutions obtained with our models are compared with the predictions of those based on a model that has been used in the past to describe synovial fluids. While all the three models seem to agree well with available experimental results, one of the models proposed by us seems to fit the data the best; it would, however, be hasty to pass judgment based on this one particular experimental correlation.

1. Introduction

Synovial fluid (SF) is a polymeric solution that is viscoelastic. That this is indeed so was documented by Ogston and Stanier [1], and confirmed later by Gibbs et al. [2] who examined the effects of the external stimuli on the viscoelastic properties of SF and by Thurston and Greiling [3] who studied the change in the viscoelastic properties in relation to pathological conditions of SF. More recent work concerning the viscoelastic properties of SF under various experimental conditions has been carried out by Gomez and Thurston [4] and Rewi et al. [5]. Despite clear evidence pertaining to the viscoelastic nature of SF, one can yet find many publications describing the SF as a viscous fluid. One reason for this is the fact the SF behaves as a generalized viscous fluid as long as it is not subject to instantaneous external stimuli and thus in processes that are steady or do not have instantaneous components its response is that of a generalized viscous fluid. In this short study, we try to delineate the conditions under which the viscous-like approximation is reasonable and to formulate a suitable mathematical model of such behavior, namely, the ability to capture the shear-thinning of the fluid (see, e.g., Ogston and Stanier [1] and King [6]).

Most biological fluids undergo a plethora of chemical reactions that leads to either degrading or enhancing their mechanical properties. To ignore such biochemical reactions and to merely model the fluid as an inert body essentially makes the studies irrelevant to most problems of physiological interest as it is the pathological situations that need to be understood and modeled appropriately. For instance, if one were to think of a biological material such as blood, it is the biochemical reactions that lead to a variety of blood disorders that present challenges with regard to the modeling; not that the fact that the flow takes place in a complex geometry is unimportant, merely that it concerns itself more with the solution methodology than the development of an appropriate mathematical model. Interestingly, in the case of SF, the response is non-Newtonian when it is healthy and normal and close to the classical Newtonian fluid when it is inflamed (see page 4). Our interest here is with regard to the mathematical modeling rather than the study of a specific equation in a complex geometry. The approach that we use, namely, that of coupling the balance of linear momentum for the body of interest with a convection-diffusion equation for a constituent which affects the properties of the body, has been adopted to study the degradation and enhancement of the flow characteristics of fluids (see Bridges and Rajagopal [7]) and the load carrying capacity of solids by several authors (the relevant references and other ad hoc approaches to the problem are discussed in detail in Shah et al. [8] and Darbha and Rajagopal [9]).

From the biochemical point of view the synovial fluid is the ultrafiltrate of blood plasma (namely, plasma free from large proteins) enriched with the locally synthesized polysaccharide molecules, called hyaluronan (HA). For closer biochemical description of SF and HA see Chapter 11 in Voet and Voet [10]. The presence of HA in the synovial fluid is important since the HA is the constituent responsible for the non-Newtonian behavior of the synovial fluid, which determines the functioning of the whole synovial joint system. The molecule of hyaluronan in SF comprises of a very long unbranched and strongly anionic polymer of as many as 25000 repeated disaccharide units, which consist of N-acetylglucosamine and D-glucuronic acid linked alternately by 𝛽(1→4) and 𝛽(1→3) bonds. The anionicity of the HA chain plays a significant role in its hydrophilic properties. Roughly speaking, the outwardfacing hydrophilic carboxylate groups and the carboxamide and hydroxyl groups of adjacent sugar residues provide high number of free potential bound places for the water molecules. Since in aqueous solutions the HA chains are heavily hydrated, they occupy a large volume, more than 1000-fold larger than in a dry state and they are able to bind water up to 10000 times their own volume. This allows the chains to create in the SF expanded coil configurations which come into contact with each other even at relatively low mass concentration of 0.1%. At higher concentration the chains create flexible molecular networks throughout the solution which is fastened by noncovalent (and thus relatively weak) bonds with the same character as the water-chain bonds. In this point of view, the fastening bonds are actually the concurrents of the water-chain bonds. Which of the bonds is preferred depends mainly on the flow conditions. During loading which is relatively fast, the chains of HA remain linked to each other randomly organized in the solution, and thus the SF response is elastic-like and consequently the synovial joint can bear the intensive loading. On the other hand, during slow and prolonged flow conditions, the chains are forced to align and thus the bonds between them are being replaced by the water molecules. This leads to the partial separation of the HA chains which results in their ability to slip past each other more freely and the SF responds like a viscous fluid. For more details, see Bothner and Wik [11], Scott et al. [12], or Voet and Voet [10]. We shall be concerned with this viscous fluid-like response.

2. Model

2.1. Model Proposal

Since the response of the fluid, whether elastic effects or viscous effects are predominant, depends on the nature of the flow, the model for SF must depend on the “dynamics” of the flow. Higher shear rates imply higher alignment of the chains and thus a decrease in the viscosity. On the other hand, the influence of concentration works contrariwise because higher concentration of HA implies higher enlacement of the chains, which increases the viscosity. In previous studies, the viscous behavior of SF has been mathematically modeled by a shear-thinning fluid with constant concentration (Rudraiah et al. [13], Lai et al. [14]) or the effects of concentration and the shear rate on the viscosity were separated (Morris et al. [15]). This is contrary to the results of the experiments which show that the concentration influences the shear rate response, which could imply that the shear rate index could depend on the concentration for the power-law model. The restriction of constant concentration is not appropriate for modeling the SF behavior under physiological conditions. In real joints, the concentration of HA varies. For example, it has been shown (see Coleman et al. [16]) that HA creates some kind of a boundary layer near the synovium with concentration five times higher than in the central parts of the synovial joint cavity (this is the consequence of the varying HA production in the synovium combined with the flow conditions). Thus we develop a model for the generalized viscosity for the SF that depends on shear rate and concentration of HA.

We consider models for the viscosity 𝜇 that belong to the power-law class. We compare the model with those introduced in the literature (e.g., Laurent et al. [17], Lai et al. [14])(Model 1) 𝜇=𝜇₀𝑒^𝛼𝑐(1+𝛾|𝐃|²)^𝑛, (1)

with our new model that takes into account the varying concentration influencing the shear-thinning effects in the following form:(Model 2) 𝜇=𝜇₀(𝛽+𝛾|𝐃|²)^𝑛(𝑐). (2)

In the both models the parameters 𝛼,𝛽,𝛾, and 𝑛 are unknown and they have to be experimentally determined. Since the SF with zero concentration of HA is basically blood plasma, the parameter 𝜇₀ should represent the plasma viscosity. The tensor 𝐃 stands for the symmetric part of the velocity gradient. Because the real physiological values of concentration of HA occurs normally between 0 mg/mL and 15 mg/mL and for concentrations higher than 20 mg/mL the solution starts to turn into a gel (for which our approach is inappropriate), we make the parametrization of the concentration in the following fashion: 𝑐=𝑐_real/𝑐_lim, where 𝑐_lim=20 mg/mL and 𝑐 is the parametrized concentration of HA with values between 0 and 1. For model 2 we need to also specify the dependence of the shear-thinning index on the concentration 𝑐.

To find a suitable form for 𝑛(𝑐) in model 2, we have to satisfy the natural conditions for the shear-thinning index, which are (1) the values of 𝑛(𝑐) have to remain in the range between (−0.5,0) to ensure the shear-thinning behavior, (2) the dependence on the concentration has to be monotonic, and (3) 𝑛(0)=0 as the SF becomes Newtonian for zero concentration. In order not to make the formula too complicated, we decided to use the exponential behavior with one free parameter(Model 2a) 1𝑛(𝑐)=2(𝑒^−𝛼𝑐−1), (3)

and a simple rational function with two free parameters(Model 2b) 1𝑛(𝑐)=𝜔(𝛼𝑐²+1−1), (4)

which satisfy these conditions and mainly, the fitting procedures, which are the topic of the next section, and lead to better results than for any other simple function with only one or two free parameters.

2.2. Unknown Parameters Identification

Each of the models introduced contains some unknown parameters. For model 1 they are 𝛼,𝛾,𝑛, in model 2a we have 𝛼,𝛽,𝛾, and in model 2b we have 𝛼,𝛽,𝛾,𝜔. We find these parameters by a fitting technique, more precisely by using the least squares method to fit the available experimental data taken from Ogston and Stanier [1], which represent the range of the mean concentration of HA in the normal SF. (Pathological SF or SF where inflammation has been introduced usually exhibit decline in concentration of HA, or the HA chain length is significantly shorter and thus these fluids do not exhibit such dominant non-Newtonian characteristics any longer. By normal SF we mean SF with rheological responses and biochemical composition as the SF of a healthy young individual.) The fits were computed using the MATLAB software. One can see the plotted results in Figure 1 and fitted values of the model parameters in Table 1.

Table 1 

Fitted values of the parameters for the models; 𝜇₀=0.01𝑃𝑎⋅𝑠.

(a) Model 1

(b) Model 2a

(c) Model 2b

Figure 1 

Relative viscosity against shear rate for different concentrations. Graphs of the relative viscosity for all the models show the fitted curves and the experimental data (points) which were taken for the fitting procedure. Here we use the notation from Ogston and Stanier [1] for 𝜇rel=𝜇/𝜇ref, where 𝜇ref is defined by the viscosity of the glycerol solution.

Even though the models of class 2 fit the experimental data more accurately, the model 1 can fit the data reasonably well for some specific applications in the range of the concentrations that is being used, this means in the range of 0.14–0.25. To compare all models for higher concentration we use the results portrayed in Figure 2. From these graphs, it is easy to observe that the model of class 2 predicts realistic values of viscosities even for higher concentrations, while the model 1 can not be used outside the range of concentrations in which it was fitted. This does not mean that the model of class 2 extrapolations give accurate viscosity values, for these can be validated only by experiments for extended range of concentrations.

(a) Model 1

(b) Model 1; logarithmic scale

(c) Model 2a

(d) Model 2a; logarithmic scale

(e) Model 2b

(f) Model 2b; logarithmic scale

Figure 2 

Relative viscosity against shear rate for higher concentrations, plots (a), (c), (e), and the same graphs in the logarithmic scale, plots (b), (d), (f). Here we use notation from Ogston and Stanier [1] for 𝜇rel=𝜇/𝜇ref.

For comparison, we plotted the shear-thinning index behavior of all the models that have been introduced, see Figure 3. We found that when the concentration tended towards unity, when the fluid behaves close to gel, the power-law index for the viscosity tended to −0.5, the parametric value that determines the limit of shear thinning.

Figure 3 

Shear-thinning index for all models. The highlighted area represents the concentration range of the experimental data, thus the mean concentrations in SF under normal physiological conditions.

3. Numerical Simulations

3.1. Governing Equations

Even though SF is complex biological material, under normal conditions, it can be approximated as an incompressible homogeneous single constituent fluid. This homogenization can be done because the physiological mass concentration of HA is very low (usually less than 1%) and we do not assume any presence of high-molecular proteins, aggressors, or inflammation. In general we would have to use a mixture theory approach and record the balance equations for each of the constituents. We however ignore the balance laws for HA and merely write down a convection-diffusion equation for it. The tacit assumption is that the inertial effects due to the motion of HA and the interaction terms are negligible. The only effect of HA is the influence its concentration has on the properties of the plasma. Such an approach has been adopted by Bridges and Rajagopal [7] to study of pulsatile flow of chemically reacting fluids. Further, we describe the flow of the SF fluid in the terms of velocity field 𝐯 and pressure field 𝑝 which are governed by the generalized Navier-Stokes equations. The concentration (𝑐) distribution which influences the flow only through the generalized viscosity is described by a convection-diffusion equation. One could even extend the convection-diffusion equation by the volumetric production/destruction term which would stand for eventual enzymatic reactions in the inflammatory stage of pathological diseases of the SF. Since we do not allow any kind of pathological behavior of the SF in our model, and the HA chains are produced only by the synovial membrane, we assume this volumetric term to be zero. The system of the governing equations takes the form of𝜚div𝐯=0,𝜕𝐯𝜕𝑡+𝜚[grad𝐯]𝐯=−grad𝑝+2div(𝜇(𝑐,|𝐃|)𝐃),𝜕𝑐𝜕𝑡+(grad𝑐)⋅𝐯=div𝐅(𝑐), (5) where 𝜚 represents the SF density, 𝐃=1/2(∇𝐯+(∇𝐯)^𝑇) and 𝐅 stands for the diffusive flux of the concentration 𝑐. Since the concentrations of HA are very low, we assume that the density 𝜚 is not influenced by the varying concentration 𝑐 and thus it remains constant. We specify the diffusive flux by the constitutive equation 𝐅=𝐷_𝑐grad𝑐, (6) with constant diffusivity 𝐷_𝑐, whose physiological value was taken from Coleman et al. [16]. Even though experimentally it has been found that the diffusivity can change with the concentration and shear-rate (see Rudraiah et al. [13], Gribbon et al. [18]), we consider it to be constant since the variations and the value itself are small compared to dominant convection.

It will be useful to recast (5) in terms of nondimensional variables, defined by 𝐱𝐗=𝐿^∗𝐯,𝐕=𝑉^∗𝑐,𝐶=𝐶^∗,𝑝𝑃=𝑃^∗𝜇,𝑀=𝑀^∗, (7) where 𝐿^∗ and 𝑉^∗ are the characteristic length and velocity, respectively. We use the parametrization of concentration which has been introduced in Section 3, and for consistency, we choose 𝑃^∗=𝜚(𝑉^∗)²,𝑀^∗=𝜇₀. (8)

From now, for simplicity of the notation, we use instead of capitals the small letters for the nondimensional variables. Thus, the system of governing equations is transformed into div𝐯=0,𝜕𝐯1𝜕𝑡+[grad𝐯]𝐯=−grad𝑝+Re2div(𝜇(𝑐,|𝐃|)𝐃),𝜕𝑐1𝜕𝑡+(grad𝑐)⋅𝐯=Pediv(grad𝑐), (9) with standard notation for reduced Reynolds number Re=𝜚𝐿^∗𝑉^∗/𝜇₀ and Péclet number Pe=𝐿^∗𝑉^∗/𝐷_𝑐.

The mathematical analysis of such kinds of system can be found in Bulíček et al. [19] or in Antontsev and Rodrigues [20]. While Bulíček et al. [19] consider the flow of an incompressible fluid with concentration dependent viscosity which corresponds to our model 1, the paper of Antontsev and Rodrigues [20] treats the flow of incompressible power-law model for viscosity with power-law index dependent on temperature which could be adapted for our model of class 2.

3.2. Geometry, Initial, and Boundary Conditions

For the needs of numerical simulations we consider a simple two-dimensional rectangular domain denoted by Ω, whose geometry is shown in Figure 4. On the walls we prescribe the no-slip boundary conditions for the velocity 𝐯 and Neumann or Dirichlet boundary conditions for concentration 𝑐 in the following way: 𝑐|_Γ1=0.1+0.1𝑡𝑡∈⟨0,2⟩,0.3𝑡>2,𝜕𝑐|Γ𝑣𝜕𝐧=0,_𝜏|_Γ1=𝑣₀𝑣𝑡𝑡∈⟨0,1⟩,₀𝑣𝑡>1,^𝑛|_Γ1Γ=0,𝐯|=𝟎, (10) where Γ=Γ₂⋃Γ₃⋃Γ₄ (see Figure 4) and 𝑣₀ is the time constant parabolic profile of the bottom wall velocity of the form 𝑣₀=(1+𝑥)(1−𝑥);𝑥∈Γ₁. By 𝑣_𝜏 and 𝑣_𝑛 we denote the tangential and the normal component of the velocity vector 𝐯, respectively. The nonzero boundary conditions are prescribed only on the bottom wall. The initial condition is the rest state with uniform concentration 𝐯(0,𝐱)=0,𝑐(0,𝐱)=0.1. (11)

Figure 4 

Geometry of the domain Ω=⟨−1,1⟩×⟨−0.1,0.1⟩. Here, we indicate the boundary parts of Ω by symbols Γ1,Γ2,Γ3, and Γ4.

We decided to choose such geometry and boundary and initial conditions to demonstrate and compare the basic flow behavior of the models that are being considered. The boundary conditions are made to be continuous along the whole boundary and to exhibit continuous evolution in time from the rest state, which is advantageous in the numerical implementation.

3.3. Numerical Method

In order to use the finite element method we introduce the governing equations system (9) in weak form: we seek for triplet (𝐯,𝑐,𝑝) such that 𝐯∈(𝑊^1,2(Ω))², 𝑐∈𝑊^1,2(Ω), and 𝑝∈𝐿²(Ω), 𝐯 and 𝑐 satisfy the Dirichlet boundary conditions from (10) and the triple (𝐯,𝑐,𝑝) satisfies _Ωdiv𝐯𝜁d𝑥=0∀𝜁∈𝑃,_Ω𝜕𝐯𝜕𝑡⋅𝝋d𝑥+_Ω=([∇𝐯]𝐯)⋅𝝋d𝑥_Ω2𝑝div𝝋d𝑥−Re_Ω𝜇(𝑐,|𝐃|)𝐃∶∇𝝋d𝑥∀𝝋∈𝐕,_Ω𝜕𝑐𝜕𝑡𝜓d𝑥+_Ω1∇𝑐⋅𝐯𝜓d𝑥=−Pe_Ω∇𝑐⋅∇𝜓d𝑥∀𝜓∈𝐶, (12) where the spaces 𝐕,𝐶, and 𝑃 are defined as 𝐕=(𝑊₀^1,2(Ω))²,𝐶={𝑐∈𝑊^1,2(Ω);𝑐|_Γ1=0}, and 𝑃=𝐿²(Ω).

The discretization of the problem then follows standard finite element procedure. First we use the quadrilateral mesh 𝒯_ℎ to approximate our domain Ω by domain Ω_ℎ and then we define finite dimensional approximation to the spaces 𝐕,𝐶, and 𝑃𝐕_ℎ={𝐯∈(𝐶(Ω_ℎ))²∶𝐯|_𝐾∈𝑄₂(𝐾)∀𝐾∈𝒯_ℎ𝐶},_ℎ={𝑐∈𝐶(Ω_ℎ)∶𝑐|_𝐾∈𝑄₁(𝐾)∀𝐾∈𝒯_ℎ𝑃},_ℎ={𝑝∈𝐿²(Ω_ℎ)∶𝑝|_𝐾∈𝑃^dis₁(𝐾)∀𝐾∈𝒯_ℎ}, (13) where 𝑄₂(𝐾) and 𝑄₁(𝐾) denote the space of biquadratic, bilinear functions on the quadrilateral element 𝐾, respectively, and 𝑃^dis₁(𝐾) denotes the space of linear functions on 𝐾, without the requirement of continuity between adjacent elements. This combination of velocity and pressure spaces satisfies the Babuška-Brezzi stability condition and hence guaranties the solvability of the coupled discrete system.

Since the reduced Reynolds number is of order 10, the Navier-Stokes system discretization need not to be stabilized. On the other hand, since the Péclet number is of order 10³ the discretization of convection-diffusion equation for the concentration has to be stabilized in order to prevent spurious oscillations in the numerical solution. We have used the classical streamline diffusion stabilization which decreases the presence of the spurious oscillation by introducing additional numerical diffusion only in the direction of the flow.

For the time discretization we use the Crank-Nicolson scheme and implicit treatment of the incompressibility constraint and pressure with time step 𝜏. The resulting nonlinear algebraic system can be written in the following form: 𝐁𝐯^𝑛+1=0,(𝐌_𝐯+𝜏2𝐀_𝐯(𝐯^𝑛+1,𝐜^𝑛+1))𝐯^𝑛+1+𝜏𝐁^𝑇𝐩^𝑛+1=𝐹(𝐯^𝑛,𝐜^𝑛),(𝐌_𝐜+𝜏2𝐀_𝐜(𝐯^𝑛+1))𝐜^𝑛+1=𝐺(𝐯^𝑛,𝐜^𝑛), (14) where 𝐌_𝐯 and 𝐌_𝐜 are the finite element mass matrices, 𝐁 is the discrete gradient operator, 𝐁^𝑇 the discrete divergence operator, 𝐀_𝐯 and 𝐀_𝐜 are the corresponding discrete convection-diffusion operators. The vector (𝐯^𝑛,𝐜^𝑛,𝐩^𝑛) refers to the solution of the finite element approximation at the time level 𝑛.

This nonlinear algebraic system is solved by the damped Newton method and linear subproblems by using direct linear solver.

3.4. Computational Results

We computed numerical solutions corresponding to all the models for the boundary conditions (10) and initial conditions (11). To show the influence of shear-thinning on the flow, we also solved the problem described by the Navier-Stokes equations with constant viscosity (𝜇₀) for the same settings (in the text denoted as the N-S model). All plotted graphs show the numerical solution at time 𝑡=100. We introduce both area distribution and profiles on central lines of the computed solutions.

First, we present the distribution of the velocity magnitude in the whole domain (Figure 5), and for more detailed comparison we include the 𝑣_𝑥 and 𝑣_𝑦 velocity profiles on 𝑦- and 𝑥-cuts of the domain (Figure 6). Next, we show the concentration (Figure 7) and dimensional viscosity distributions (Figures 9 and 11 in logarithmic scale) and their profiles on the same domain cuts as for the velocity (Figures 8, 10, and 12). We do not plot the viscosity distribution for the case of the N-S model since it remains constant for all time. As the last set of figures we present the change of shear-thinning index for models of class 2, (Figure 13).

(a) N-S: 𝜇=𝜇0

(b) Model 1: 𝜇=𝜇0𝑒𝛼𝑐(1+𝛾|𝐃|2)𝑛

(c) Model 2a: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)(1/2)(𝑒−𝛼𝑐−1)

(d) Model 2b: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)𝜔(1/(𝛼𝑐2+1)−1)

Figure 5 

Velocity magnitude for all models.

(a) Velocity profiles on line 𝑥=0

(b) Velocity profiles on line 𝑦=0

Figure 6 

Comparison of velocity profiles on 𝑦- and 𝑥-cuts.

(a) N-S: 𝜇=𝜇0

(b) Model 1: 𝜇=𝜇0𝑒𝛼𝑐(1+𝛾|𝐃|2)𝑛

(c) Model 2a: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)(1/2)(𝑒−𝛼𝑐−1)

(d) Model 2b: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)𝜔(1/(𝛼𝑐2+1)−1)

Figure 7 

Concentration distribution.

(a) Concentration profiles on line 𝑥=0

(b) Concentration profiles on line 𝑦=0

Figure 8 

Comparison of concentration profiles on 𝑦- and 𝑥-cuts.

(a) Model 1: 𝜇=𝜇0𝑒𝛼𝑐(1+𝛾|𝐃|2)𝑛

(b) Model 2a: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)(1/2)(𝑒−𝛼𝑐−1)

(c) Model 2b: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)𝜔(1/(𝛼𝑐2+1)−1)

Figure 9 

Viscosity distribution.

(a) Viscosity profiles on line 𝑥=0, zoom

(b) Viscosity profiles on line 𝑦=0, zoom

Figure 10 

Comparison of viscosity profiles on 𝑦- and 𝑥-cuts.

(a) Model 1: 𝜇=𝜇0𝑒𝛼𝑐(1+𝛾|𝐃|2)𝑛

(b) Model 2a: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)(1/2)(𝑒−𝛼𝑐−1)

(c) Model 2b: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)𝜔(1/(𝛼𝑐2+1)−1)

Figure 11 

Viscosity distribution in logarithmic scale.

(a) Viscosity profiles on line 𝑥=0, logarithmic scale

(b) Viscosity profiles on line 𝑦=0, logarithmic scale

Figure 12 

Comparison of viscosity profiles on 𝑦- and 𝑥-cuts.

(a) Model 2a: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)(1/2)(𝑒−𝛼𝑐−1)

(b) Model 2b: 𝜇=𝜇0(𝛽+𝛾|𝐃|2)𝜔(1/(𝛼𝑐2+1)−1)

Figure 13 

Shear-thinning index.

Comparing the flow fields, one can see significant difference between the solution of the N-S model, model 1, and models 2a, 2b. The computation for N-S model gives nearly symmetric solution along 𝑦-axis while the solutions for the other models are asymmetric (see Figure 6(b)). As one can see from the velocity profiles in Figure 6(a), while the maximum of the back flow in the N-S case is located above the 𝑥-axis, for the model of class 2 it is slightly smaller and located around the 𝑥-axis, whereas for model 1 the maximum is the smallest, flattened and located below the 𝑥-axis.

As for the velocity profiles, the concentration profiles for model 2a and 2b are similar. The transport of concentration is fastest for the N-S model due to the higher speed of the flow compared to the other models. In Figure 8(a) one can see the presence of numerical oscillations in the solution through very small diffusivity and high gradient of the concentration, which causes the decrease of concentration below the value of 0.1.

Inspecting Figure 9 one can observe, that the computed viscosity for model 1 attains higher values compared to model of class 2, mainly in the areas of higher concentration. On the other hand, Figure 11 shows the apparent difference in the viscosity distribution.

The last set of graphs in Figure 13 presenting the shear-thinning index distributions shows that the model 2b exhibits slightly larger spread of values, it extends to a higher maximum and lower minimum than model 2a, in this example.

4. Summary

We have modeled the synovial fluid as a generalized viscous fluid as such a model is a good approximation for a class of flows of the synovial fluid provided the processes do not involve instantaneous response. We have compared three different models, with that of the Newtonian model. Model 1 combines in an artificial way two characteristics introduced in the literature, namely, exponential dependence of viscosity on the concentration and the shear-thinning ability, and the newly introduced models 2a, 2b which describe the SF as a power-law type fluid with concentration dependent exponent. While all three models fit the experimental data reasonably well it is quite clear that the model 2b provides the best fit. Of course, this might just be happenstance for this particular boundary value problem. Unfortunately additional experimental data are not available in order to corroborate the model further. In light of the results that we have established, it would be reasonable to suppose that the model 2b predicts the response of the synovial fluid the best amongst models that have been used thus far.

We numerically solved the flow in a simple configuration for all three models. There are certain differences between solutions for model 2a, 2b even though the fitted results were very close to each other. On the other hand, the solution for model 1 exhibits distinct patterns with regard to the flow field and viscosity distribution. This means that even though the viscosities for the models are close, the small difference is sufficient to produce significant differences in the flow pattern.

Since in the system of governing equations there is a convection-diffusion equation with very small diffusivity, on using standard finite element method the solutions exhibit large spurious oscillations, primarily in the solution for the concentration and consequently in the viscosity. These oscillations were reduced by the streamline diffusion method, however, it will be useful to implement better numerical stabilization methods (internal penalty method of Turek and Ouazzi [21], TVD method by Kuzmin and Turek [22]) to obtain solution without any artificial oscillations.

We have solved the problem in a geometry that does not reflect the complexity of the actual physiological problem. Our aim was to merely obtain an understanding of the flow characteristics of the fluid. In the future, we plan to study the flow in a geometry which is more physiologically relevant.

Acknowledgments

J. Hron was supported by the project LC06052 (Jindřich Nečas Center for Mathematical Modeling) financed by MŠMT. J. Málek’s contribution is a part of the research project MSM 0021620839 financed by MŠMT; the support of GAČR 201/09/0917 is also acknowledged. P. Pustějovská was supported by the project LC06052 (Jindřich Nečas Center for Mathematical Modeling), financed by MŠMT and by GAUK Grant no. 2509/2007. K. R. Rajagopal thanks the National Science Foundation for its support.