An analytical study of the squeezing ﬂow of synovial ﬂuid

– The problem of the squeezing of a ﬂuid ﬁlm with a complex rheology is studied. The ﬂuid is represented by the Phan-Tien and Tanner (PTT) constitutive equation. The goal of the paper is to develop an analytical approach which allows one to describe and quantify this process, which is often present in classical lubrication problems, but particularly in human joints. The present approach can potentially supplement time-consuming and costly purely numerical methods which are also susceptible to misleading artifacts. A parametric study reveals the eﬀect of ﬂuid elasticity compared to viscous eﬀects and some explanation is suggested as to their possible links with joint diseases


Introduction
The tribological effectiveness of a joint is a major factor in its longevity.There are a wide variety of proposed theories as to the mode of lubrication in human joints.Experimental studies have led researchers to associate the four phases of the gait cycle to regimes of lubrication of the knee joint, Dumbleton [1].The squeezing effect is dominant in the phase of the cycle, called the "heel strike", with significant hydrodynamic load carrying.The porosity of the cartilage can boost the squeezing effect, by leaking of the aqueous phase of synovial fluid through the matrix.The effective viscosity is thereby increased and helps maintain a film thickness in the micrometer range [2][3][4].
According to the theory of elastohydrodynamic lubrication, the squeezing effect plays an important role in maintaining the unsteady film thickness during the walking cycle [5].During this unsteady phase, the linear viscoelastic properties of synovial fluid in prostheses cannot be ignored, according to Mazzucco et al. [6].
Synovial fluid, clearly a non-Newtonian fluid, has been the subject of a multitude of research studies, see Fung [7].From early on, its viscoelastic properties have been identified, Ogston and Stainer [8]; and its shear thinning behavior described, King [9].Similarly the effects of normal stress have been characterized, e.g., by McCutchen [10], and Mow and Ateshian [11].Measurements of elastic modulus and loss modulus have shown that the viscoelasticity a Corresponding author: Benyebka.Bou-Said@insa-lyon.fr of synovial fluid is due to its high content of hyaluronic acid, a complex macromolecule.Studies on bovine synovial fluid were performed to demonstrate the role of rheology on its behavior, Oates et al. [12].Hyaluronic acid has been well-studied in both natural and synthetic settings and its viscoelastic properties are clearly demonstrated, Szwajczak et al. [13].Oda and Sugishita [14] and many other authors [11] have shown that important lubricating properties of synovial fluid in natural joints can be attributed to its viscoelastic properties.Synovial fluid is an excellent shock absorber during rapid movements, but also has lubricating properties in relatively slow motion, Ghosh and Guidolin [15].For an articulated joint almost at rest, the period of deformation is long during which the cross linked molecular chains are broken, and the rheological properties are predominantly viscous.At high deformation rate the adjustment of molecular chains of hyaluronic acid in short periods of oscillatory deformation cannot take place, which forces the fluid to behave in a partially elastic manner (Safari et al. [16] and Kobayashi et al. [17]).Modeling the behavior of synovial fluid in articulated joints continues as a subject of ongoing research [18][19][20].
An additional understanding of the rheology of synovial fluid can be useful for better understanding of the tribological behavior of artificial joints, and their improved design.While the viscoelastic behavior of synovial fluid is widely accepted, the choice of a constitutive model is still a cause for debate.The model used in this study, introduced by Phan-Thien and Tanner [21] (PTT), can be derived from the theory of polymer networks at the molecular level.The law is considered one of the most realistic First convected time derivative (= Dimensionless quantity, see Equation ( 31 in describing the behavior of polymer solutions because it describes the main characteristics of non-Newtonian fluids in shear and elongation.In the field of lubrication one can cite references [22][23][24] in which the research described is based on this model.

Analysis
This work is intended as a contribution to the study of lubrication of the knee joint by using the PTT model.The effect of squeezing of the synovial fluid is considered to be the predominant mechanism governing the joint behavior.During the squeezing phase, the synovial fluid is not instantly expelled from between the two joint surfaces and since the fluid resists this extensional motion, a finite time is needed for the surfaces to come in contact.In this time interval, a pressure builds up and the load is carried by the squeezed film; see, for example, Fein [25], Torzilli [26], and Timothy and David [27].
The contact between articular surfaces as well as the behavior of synovial fluid have been considered in different ways, depending on the particular study.Sinha et al. [28] used a model of the rigid contact of a sphere with a porous surface, neglecting the sliding movement.Bujurke [29] considered a slider bearing configuration.Nabhani et al. [30] have modeled the synovial fluid as a couple stress fluid and considered the flow in the porous matrix with a generalized approach including viscous and inertial effects.
The knee joint is schematically represented, as shown in Figure 1, by two horizontal surfaces, separated by a lubricant fluid film.The effect of squeezing between two disks is interesting for several reasons.This configuration captures much of the basic physics of the knee joint, and is among the most fundamental and widely studied unsteady flow.
The PTT constitutive equation, describing the stress in terms of the fluid velocity field, is expressed in the notation of Bird et al. [31] as: where τ = τ i with The extra stress τ is described as the superposition of partial stresses τ i .The rate-of-strain is denoted γ, and the subscript (1) refers to the first convected time derivative meaning the upper convective Maxwell derivative.The relaxation time λ, the viscosity η, and the function Z are described by the spectra {λ 1 , λ 2 , . ..}, etc.
The originality of the PTT model comes from the introduction of an effective velocity gradient tensor (L = v − ( 1 2 )ξ γ) and the function Z i (tr(τ i )), (tr denoting the trace of the tensor), characterized respectively by the parameters ξ and ε.These parameters are usually determined by experimental data.The effective velocity gradient tensor term corrects for disturbance to the assumed background motion of the solvent by the solute macromolecules.The hydrodynamic interaction of a macromolecule with the surrounding macromolecules modifies the flow around the macromolecule, which causes a motion different from what would otherwise be imposed by the solvent.The function Z represents the ratio of the rate of destruction of junctions at particular time to that of equilibrium conditions.
Phan-Thien and Tanner designate ξ as a slip parameter, representing a measure of the difference between the imposed flow and the flow sensed by the macromolecule.A non-zero value of ξ is linked to the existence of a second normal stress difference, as well as the amplitude of stress over shoot.When ε and ξ are zero, the upper convective Maxwell model is recovered.When ξ = 1, the PTT model reduces to the co-rotational Maxwell model.For ξ = 2, the PTT model becomes the lower convective Maxwell model.In fact, experimental data have shown that not exceed the value of 0.4 [32].
In the theory of lubrication, scale of the coordinate direction across the film thickness z, is small compared to the other two directions.Thus, the order of magnitude of the various terms are, v u, w; and ∂ ∂z ∂ ∂r , ( 1 r ) ∂ ∂θ ; with u, v, and w being the flow velocity components in the radial, axial and azimuthal directions, respectively.
Subject to this scaling, and omitting the numerical i index, for clarity, Equation (1) can be written in cylindrical coordinates: The flow is considered sufficiently slow that angular distortion does not occur and thus τ rθ = 0.
By analogy to the corresponding Newtonian problem, the pressure and normal stresses are functions of r, θ and t; while the remaining shear stresses τ rz and τ θz are functions of r, z, θ and t.The momentum equation is now presented in the three cylindrical coordinate directions r, z and θ.Neglecting the effect of gravitational forces and fluid inertia, gives: In viscoelastic flow pressure p does not have an obvious physical meaning.If the flow is incompressible, it is the total stress π which exerts a physical force on a bounding surface.The extra stress τ is uniquely determined by the kinematics of the flow and the rheological model, and is an additive portion of the total stress π.Isotropic pressure p is simply their difference, π = τ + pδ.
In principle, the initial movement is impulsive in nature.At t = 0 − , the fluid is at rest, and the stress and velocities are zero at any point of the fluid.At t = 0 + , the upper disk is set in motion.The movement of the upper disk is transmitted to adjoining particles through viscous effects and eventually all the fluid is set in motion.During this time, the dynamics equations are not defined, the hydrostatic pressure and the acceleration of the upper disk are indeterminate.Several studies have addressed this point [22,25,31].In our case, we assume that at the initial moment, we have: The boundary conditions are such that: Thus, the flow is described by six equations which yield the six components of the extra stress and two equations of momentum.These shear flows can be classified as "slow" or "creeping" meaning fluid particles possess negligible inertia.The microstructure of the total solventpolymer mixture does not undergo significant deformation.Phan-Thien and Tanner showed that in this case it is not necessary to introduce variation of the function Z(tr(τ )) of Equations ( 1)- (7).The ε parameter has no special significance in a flow where the shear is dominant, and in a simple shear it is zero.By the same token, the material parameter ξ does not have great significance in a stretching flow.If, in addition, it is assumed that the flow is quasisteady, which is a reasonable assumption except for the brief interval just after the impulsive start, the tangential stresses can be found algebraically as: with ) in which η r and η θ are the apparent viscosities in the r and θ directions, respectively, in the plane normal to the direction z.

Modified Reynolds equation
We seek an approximate solution based on the Newtonian velocity field and following Tichy and Bou-Said [33] we can approximate the square of velocity gradients (shear rates) according to the average value of the Newtonian velocities across the film thickness: Where u N and w N represent the radial and axial Newtonian velocity components, respectively.The normal stresses thus become: Consequently, the assumption that the normal stresses are independent of the axial coordinate z is respected.In addition, the apparent viscosities are no longer a function of the axial coordinate z.From these relationships, the momentum Equation (8) becomes: (24) and the momentum Equation ( 9) can be expressed as: Following the classical method for deriving the Reynolds equation, we obtain: This is a new modified Reynolds equation of a PTT fluid for squeezing between two surfaces in cylindrical coordinates.

Axisymmetric squeezing
Both articular surfaces are considered as two impermeable parallel rigid disks of radius R separated by a horizontal distance h(t).The upper disk is moving at a velocity dh dt while the lower disk is stationary.The flow of the lubricating film between the two disks is axisymmetric ( ∂ ∂θ = 0) and the squeezing begins from an initial height h 0 .If the boundary conditions are set to zero except for V 2 (t) = dh dt , the modified Reynolds Equation ( 27) is simplified considerably to: and the average shear rates become: The radial component u N of the Newtonian velocity, is The boundary conditions on the normal stress at a point on the edge of the disk (at z = h and r = R) are a source of debate.Several authors use the condition p = p atm when there is no effect of extra stress, and similarly others use π zz = p atm .Lee et al. [34] recommend the integral condition h 0 π rr dz = p atm which has intuitive physical merit.We have used π zz = p atm which results in the simplest implementation of the equations: at r = 0: ∂π zz ∂r = 0, at z = h and r = R : π zz = p atm (30) The integration of the modified Reynolds equation with the axisymmetric boundary conditions (30), gives the dimensionless axial normal stress: The dimensionless load capacity of the squeeze film is determined by integrating the axial normal stress over the disk area: The load capacity is determined by two parameters, namely the Weissenberg number W e, which represents the elastic forces versus the viscous forces W e = λv0 h0 and the aspect ratio h0 h0 = h0 R , which characterizes the initial fluid film thickness h 0 versus the radius R of the disks.
Likewise we obtain the dimensionless shear stress: τrz (r, t) = − 3 h3 The corresponding dimensionless Newtonian load capacity and dimensionless Newtonian shear stress are W N ( t) , respectively.The various dimensionless quantities are defined as follows: We recall that v 0 is the initial axial squeeze velocity.

Results and discussion
The values of the rheological parameters, geometry and mass of the modeled physical system refer to those of healthy synovial fluid in a typical knee joint: η is the viscosity of 0.086 Pa.s, the slip parameter ξ is 0.4, the radius R of the upper and lower disks is 1 cm.Two values of the initial film thickness are used, h 0 = 10 and 100 μm, and the initial squeeze velocity v 0 is 0.1 m.s −1 .
The study of squeezing at a constant speed is used to evaluate the load capacity as a function of the Weissenberg number, reflecting the potential importance of the elasticity of a healthy synovial fluid compared to the diseased case of osteoarthritis.The load is explicitly applied and the interval is divided in 101 equidistant points.
In the case of a purely Newtonian fluid (W e = 0), the load carrying capacity of the fluid film increases hyperbolically.For a Weissenberg number W e = 5 × 10 −3 (Fig. 2), the load capacity is very low.The fluid is too elastic to be able to generate significant pressure.At W e = 5× 10 −4 (Fig. 2), the lift created increases to reach a maximum and then drops to zero.This can be explained by a competition between the viscous forces and elastic forces.For the largest value of Weissenberg number, the shearing of the lubricating film does not cause a significant load bearing capacity.At the second highest value of the Weissenberg number we observe a slight increase in load capacity, but the elastic forces apparently dominate and reduce the carrying capacity to zero.For the lowest values of Weissenberg number, the viscous forces prevail causing a high load capacity, with the appearance of a hyperbolic increase (also Fig. 2).However, towards the end of the motion the elastic forces tend to reduce the bearing capacity to zero.Be aware that in this figure increasing time is from right to left.
For the lowest values of Weissenberg number, the load capacity approaches that of the Newtonian case and we can observe that the Newtonian behavior is more or less followed by the values of the Weissenberg number (Fig. 3).
At the micro level, the macromolecule at rest is a dense chain -tangled and folded on itself.The shear deformation tends to disrupt this entanglement in the direction of flow.If the shear rate is relatively lower the elasticity (relaxation time) high, the macromolecules remain entangled as the forces of interaction between the solvent and polymers are weaker than the cohesive forces between macromolecules.In this case, the apparent viscosity drops rapidly causing a lower load capacity (Fig. 4).At a high shear or low elasticity, the action of solvent tends to unravel the macromolecule in the direction of flow, and the viscosity drops much less rapidly creating a larger load.This orientation effect predominates as the elasticity effect of the macromolecules cannot counteract the direction of these macromolecules.
The evolution of the relative friction coefficient f f N as a function of the relative squeeze height is represented in Figure 5 for three values of the Weissenberg number and two values of the film thickness aspect ratio.There is a higher relative friction coefficient at higher Weissenberg numbers mainly because of the decrease of the carrying capacity with the increase of the Weissenberg number.
To study the squeezing process in greater detail is necessary to determine the evolution of the film thickness h(t) from the mechanical system dynamics.The following is the equation of motion of the upper disk in dimensionless form: where m is the mass of the upper disk and g is the gravity acceleration.
To calculate the instantaneous dimensionless height of the upper disk, a method of finite difference discretization is used.The initial rate of dimensionless squeezing is equal to minus one (the negative sign indicating squeezing).The instantaneous film thickness at succeeding time increments h( ti+1 ) is obtained by the bisection method, (a) (b) 0,0 0, 2 0 ,4 0 ,6 0,8 Dimensionless film thickness Newtonian case 0,0 0,2 0,4 0,6 0,8 The two graphs of Figure 6 represent the variation of dimensionless film thickness h( t), while the two graphs of Figure 7 indicate the variation of upper surface velocity V2 ( h).Note that in the case of a low viscous effect (high values of Weissenberg number), the surfaces contact.The axial normal stress generated in the film cannot counteract the constant gravitational force applied.For smaller W e (larger viscous effect), the evolution is almost linear for the first instants and then obtains an asymptotic value at a finite separation.When the Weissenberg number decreases, viscous effects become more important and normal stress is created which can counterbalance the applied load.This asymptotic phenomenon is more pronounced at smaller h0 because viscous effects are stronger for a thinner film.
The evolution of squeeze velocity in Figure 7 is consistent with the gap behavior of Figure 6.The speed is almost constant just after inception, and the generated normal stress may or may not be sufficient to resist a contact, depending of the value of Weissenberg number W e. As above, low W e implies relatively high viscous effects.
The behavior of the dimensionless force (load) generated by fluid forces F ( t) is portrayed in Figure 8.At high Weissenberg number the reaction of the fluid is insignificant.Gradually, the viscous effects appear to create a load that reaches a maximum and then decreases as the squeezing velocity tends to zero -a sort of overshoot phenomenon.For an initial aspect ratio h0 = 10 −3 , this decrease is clearly visible and occurs earlier for a nonzero Weissenberg number, which is in agreement with the velocity curves of Figure 7.
Figure 9 represents the dimensionless shear stress on the top plate, at the outside edge, the radial position r = 1.The shear stresses substantially follow the same trends as those obtained for the load, Figure 8.
The coefficient of friction relative to the Newtonian case f f N is shown in Figure 10.This quantity is practically constant for low values of W e. The presence of elastic effects appears to increase the shear strength.
The Newtonian fluid and the PTT fluid have essentially the same parabolic velocity profile shape as shown in Figure 11.In both cases shown ( h0 = 10 the dimensionless time t = 1, but these represent different instants of physical time and thus different values of maximum velocity.

Conclusion
Study of the behavior of lubrication mechanisms in human joints in transient loading requires sophisticated modeling to properly assess the forces involved, and their effects.The Phan-Thien Tanner continuum model was used which at the microscopic level represents the interaction of a viscous fluid solvent with a polymer macromolecule solute.Using an analytical approach based on reasonable assumptions concerning the fluid flow, we have highlighted some of the complexities inherent in the nature of such viscoelastic fluids.In some cases overshoot loads are obtained whose origin can be attributed to the effect of disturbance of the solvent flow by such macromolecules.Results are organized according to the Weissenberg number which is proportional to the product of shear rate and relaxation time.Among the numerous results, it was observed that at low values of Weissenberg number, viscous forces prevail and load capacity is relatively high.In a simple dynamic system consisting of a disk descending under its gravitational load supported by a film of fluid, the disk reaches an asymptotic value of gap height at small values of Weissenberg number, but contacts at higher values.At higher values of Weissenberg (low viscous effects) we observe that the time to contact is shorter.These results can be interpreted in several ways but we believe can be linked to disturbance of the entanglement of macromolecules, by entrainment in the direction of flow.This analysis can potentially be applied to the behavior of synovial fluid if the behavior of hyaluronic acid can be clearly identified experimentally in healthy and pathological conditions.

h 2 ,
obtained by solving the classical Reynolds equation for the same conditions.