Role of protein concentration on transient film thickness in synovial fluid lubricated joints

A computational model of protein aggregation lubrication has been developed for predicting transient behaviour in lubricated prosthetics. The model uses an advection-diffusion equation to simulate protein transport in order to map concentration changes throughout the contact and inlet zones of an elasto-hydrodynamic contact. Concentration increases lead to exponential increase in fluid viscosity giving rise to lubricating film thicknesses an order of magnitude larger than would be expected using conventional elasto-hydrodynamic theory. The model parameters have been calibrated such that good agreement in transient film thickness is achieved with observed experimental results.


Introduction
Joint replacements have been performed in the UK since the 1960s with the first hip implants.Knee and ankle replacement surgery began in the 1970s and the number of all procedures has increased rapidly ever since.Over 250,000 joint replacements were carried out in the tax year 2018/2019 in the UK, with knee and hip replacements accounting for the vast majority of these numbers [1].As the number of procedures continues to increase, more and more people will be living with prosthetic joints.
A major issue affecting implant success is wear of the implant surfaces [2,3].Wear debris can build up within the joint and cause tissue inflammation and in turn loosening of the implant.Beyond this, there are many other adverse biological responses to wear debris.Metallic ions have been found in the lymphatic system, blood and urine of joint replacement patients and have been linked to lymphadenopathy and necrosis [4].For these reasons, an analysis of the lubrication mechanisms that affect wear are central to research into extending the life of joint prostheses.

Synovial Fluid
The lubricating fluid for healthy joints is synovial fluid, which is a complex water-based mixture of proteins, glycoproteins, hyaluronan, phospholipids and cholesterol and patient specific studies have shown that the concentration of macromolecules in synovial fluid can vary significantly [5,6].Once replacement surgery has been performed, the fluid that occupies the jointperiprosthetic synovial fluidhas differing concentrations of these components and is more akin to diseased synovial fluid as found in patients with arthritis.Periprosthetic synovial fluid has an increased concentration of proteins, higher pH level and decreased bulk viscosity when compared to that of healthy patients [7].In experiments, diluted bovine calf serum is frequently used as a substitute for synovial fluid and has its own protein concentration profile.A summary of each fluid composition can be found in Table 1.
Synovial fluid shows clear shear thinning behaviour in experiments [8].While healthy, rheumatoid and osteoarthritic synovial fluid shows decreasing viscosity at shear rates of 10 3 s − 1 , periprosthetic synovial fluid shows the beginnings of a high shear plateau at these magnitudes.Cases in which high shear rate is encountered are often treated as Newtonian [2,7,9].Further studies by Castellanos et al. [10] show that for bovine calf serum, most often used in experiments, the high shear rate plateau occurs at even lower shear rates of around 10 2 s − 1 .
The concentration of protein content plays an important role in the rheological properties of synovial fluid.Studies by Coleman et al. [11] show that at constant shear rate, the increased content of hyaluronan in a saline solution induces exponentially increasing viscosity.Castellanos et al. [10] show rheometry results for varying concentrations of bovine serum albumin over a range of shear rates from 10 − 2 -10 4 s − 1 showing the effect of protein concentration on the full shear thinning behaviour of the fluid.Results show the onset of a high-shear plateau at shear rates of 10 2 s − 1 and also that the high-shear viscosity of bovine serum albumin increases exponentially with increasing protein concentration.

Protein aggregation lubrication
When pin-on-flat tribometers are employed using diluted bovine calf serum as a synovial fluid analogue between a glass disk and a femoral head implant, interferometer images show a different lubrication mechanism is at work when compared to classical isoviscous elastohydrodynamic contacts.
Studies by Myant, Cann et al. [7,[12][13][14][15][16], show that when using diluted bovine calf serum, which contains macromolecules of albumin and globulin, interferogram images show a distinct inlet region.Fig. 1 shows an example of such an image as obtained by the experiments in 2012 [12].The inlet region appears prominently on the interferogram image because the refraction properties of the fluid are noticeably altered by an increase in protein content in this area.
Another interesting phenomenon noted in the same study is that film thickness is much higher than would be predicted by classical elastohydrodynamics lubrication (EHL) analysis and moreover shows counter-intuitive behaviour.As the velocity of the sliding plate is varied, film thickness decreases as entraining velocity is increased.This is contrary to all conventional, isoviscous EHL predictions that show larger velocities induce larger fluid films.This lubrication method has been named 'Protein Aggregation Lubrication' (PAL) and as the name suggests is based on the hypothesis that proteins present in synovial fluid collect in the inlet to the contact zone.This protein aggregate sits at the inlet and drastically alters the local viscosity allowing much thicker fluid films to form.This aggregate can break down and build back up again giving transient film thickness even under steady conditions and will break down more readily if the entrainment speed increases, causing the phenomenon of decreasing fluid film with increasing speed.
In subsequent experiments also by Myant and Cann in 2014 [14], the transient nature of the protein aggregate region was examined using an inlet disruption test.For this test, the sliding glass was held at a constant velocity and the inlet aggregate allowed to form.The protein matter was then stripped from the contact zone by a passing scratch on the glass surface.
Fig. 2 shows the central film thickness and inlet length (as defined in Fig. 1) results of this inlet disruption test.The data show a steady film thickness that is disturbed by the scratch placed on the glass sliding surface.This removes the protein aggregate from the contact zone and causes the central film thickness to drop to a level within the range predicted by Hooke's eq.(7-28 nm) [17].After the scratch has passed, the protein aggregate starts to rebuild and much higher film thickness is restored.
Fig. 3 shows interferogram images from this experiment that demonstrate the inlet aggregate growth over time up to a steady inlet length.The time indicated is zero at the point that the protein begins to aggregate as indicated by the grey region of Fig. 2.

Aims and Objectives
The aim of this study is to create a computational model that captures the behaviour of protein aggregation lubrication, namely larger than predicted film thickness and the presence of protein aggregation at the inlet of the contact zone.Stokesian dynamics and force coupling methods [18,19] have been used to successfully model the multiphase flows of particulates suspended in fluids, however, these methods are computationally intensive and require a great deal of knowledge of the initial state of the suspension.It has been shown that biological fluids including blood and synovial fluid [20] can be treated as a homogeneous fluid whose properties change based on a concentration of particulates that obey an advection-diffusion-reaction or advection-diffusion eq.
A model for simulating concentration-dependent shear thinning fluid has been proposed by Hron et al. [21].The concentration of protein content is simulated using an advection-diffusion model over a simple 2D wall-bounded domain.The fundamental phenomenon of protein aggregation lubrication however is that the thin film nature of lubricated contacts causes physical blockages to the free flowing protein macromoleculesan effect that has not yet been considered.No such

Methods
The inlet disruption experiment by Myant and Cann [14] discussed in §1.2 -was used as the benchmark for the computational model developed in this study.The benefit of this experiment is that the sliding speed is constant and there are no transient start-up effects.

Geometry and Materials
Fig. 4 shows a schematic of the plane-on-sphere geometry of the experiment, including key parameters.Table 2 gives values for these key parameters, however the exact bulk viscosity of the lubricating fluid is unknown.The fluid used in the experiments was bovine calf serum diluted by water in a ratio of one part serum to three parts water by weight.This gives the resulting solution a protein content of 18 g/l.
The surface roughness of CoCrMo femoral heads is cited frequently as R a, fh = 0.01 μm [6,7,9] and the roughness of the silica-coated glass disc of this tribometer is given as R a, gd = 0.015 μm [16], resulting in a combined equivalent roughness of For a brief period where the film height is at a minimum directly after the protein has been stripped, this likely means that the contact is operating in the mixed lubrication regime, however, as the inlet aggregate builds, the film height rapidly increases to bring the contact into the full fluid-film regime.For this reason, surface roughness effects have been ignored for this study.The values of Table 2 result in a mean Hertzian dry contact pressure of approximately 90 MPa and as such is considered low enough that piezoviscous effects are also ignored [9].The model will therefore use isoviscous elasto-hydrodynamic theory to solve for the geometry and pressure of the contact zone.
For experiments involving film heights in the range of 100 nm and sliding velocities of the order of 10 mm.s − 1 , shear rates are expected in the order of 10 5 s − 1 , in excess of any expected high shear rate plateau of bovine calf serum of around 10 2 s − 1 [10] -as discussed in § 1.1.Shear thinning effects are therefore assumed negligible in this study.
Fig. 5 shows a diagrammatic overview of how the proposed protein aggregation lubrication is constructed.
For each time-step during PAL simulation, the contact geometry and pressure field are solved using the 'Elasto-Hydrodynamic Solver' and the 'Reynolds Equation Super-Solver', the flow rate of the protein matter is then calculated from the flow rate of the lubricating fluid via a geometry-dependent attenuation function and the transport of the protein is simulated using an advection-diffusion equation.The resulting protein concentration field determines the local viscosity using a concentration-dependent constitutive equation and this viscosity field is used for the next time-step.

Theory
To analyse fluid mechanics between deformable surfaces, an EHL model was used.The EHL model takes into account hydrodynamic lubrication using the Reynolds and load balance equations but can also take into account the effects of high pressures in the fluid film that elastically deform the contact surfaces.
In the sphere-on-plane geometry considered, the fluid gap, h, is much smaller than the scale of the contact zone in the x-and y-directions and thus the Reynolds equation can be applied to describe the fluid dynamics in a lubricated contact [22].
The computational model for transient protein aggregation uses an elasto-hydrodynamic solver coupled with an advection-diffusion solver to model protein transport.

Elasto-Hydrodynamic Lubrication
For this study, the surface velocities are limited to a single sliding direction and if the domain is chosen such that the x-axis aligns with this movement the Reynolds equation becomes: ∂h ∂t where p denoted the fluid pressure and η is the fluid viscosity.For fluid lubrication, the cavitation condition must also be considered which limits the pressure from becoming lower than the vapour point of the liquid.Compared to the pressures that can form inside the contact, the vapour pressure of synovial fluid is close to zero and thus the vapour point can be treated as zero.This condition necessitates that p ≥ 0 within the domain and p = ∂p/∂x = 0 at the cavitation boundaryand downstream of this boundarythe location of which is found during the solution [23].
In conjunction with the fluid equation, the EHL equations also solve for the elastic deformation of the surfaces where the deformation at any given point is calculated by integrating the effect of the pressure field and is given for each surface as: where ν denotes the Poisson's ratio and E denotes the Young's modulus of the material, which can be different for both surfaces.Summing the deformation of both surfaces and applying them to parabolic estimate of the undeformed spherical surface results in: where x and y denote the position of the deformation being calculated.
Here, x ′ and y ′ denote the dummy variables over which the integral is taken, signifying the position of the pressure contributing to the deformation.E ′ is an effective modulus that maps the Poisson's ratio and Young's modulus of each surface onto a problem whereby only the upper surface deforms and the lower surface is treated as rigid.It is defined as: The final equation in the system is a load balance to ensure that the pressure generated by the fluid is enough to balance the applied load, w.A simple integral over the pressure field is used: (2.5) Eqs. 2.1, 2.3 and 2.5 can be non-dimensionalised using parameters from the Hertzian solution of a dry contact of a sphere on a plane [24].The solution of such a problem forms a circular contact of radius, a h , where The maximum pressure, p h , and maximum deformation, h h , of the dry (Hertzian) contact solution are, Using these values as references for non-dimensionalisation, the following scaling was applied to the elasto-hydrodynamic equations -Eq's.(2.1), (2.3), (2.5).
Here, ρ 0 and η 0 are the initial density and viscosity respectively.
Using this scaling the non-dimensional Reynolds equation can be written ∂ ∂X where the dimensionless parameter, λ, is given by

.10)
The film thickness equation becomes And finally, the force balance equation becomes

Protein Transport
Local concentration can be calculated using a scalar transport Eq. [25] of the form, ∂c ∂t where c, v and D are the concentration, velocity and diffusion constant respectively.Because of the thin-film nature of the lubrication problem, it is reasonable to assume that the concentration is constant across the thickness of the fluid film and thus, Eq. 2.13 can be integrated in this direction, resulting in the following equation.The assumption of constant concentration along the z-direction allows for many simplifications, resulting in The remaining integrals are of either velocities or velocity gradients, which are are simply volumetric flow rates or gradients thereof respectively.The flow rates, q x and q y , are defined as such reducing Eq. 2,15 further to h ∂c ∂t The remaining integration can be dealt with using the Leibniz integral rule to interchange to integral and the differential operators.The extra terms associated with this rule are all zero due to the no-slip condition at the boundaries.The integral of the ∂w/∂z term is simply the vertical velocities at the surfaces.These velocities are ∂h/∂t and 0 at the upper and lower surfaces respectively, thus the advection-diffusion equation can be written as: In order to solve the scalar transport equation, the pressure solution from the EHL equations is used to calculate volumetric flow rate in each direction as follows: ∂p ∂x ∂p ∂y Non-dimensionalisation of these equations is carried out using the same scaling as used in Eq. 2.8 with the added parameters and by scaling the concentration by an initial value c 0 , such that The resulting non-dimensional equations can be written as where Pe is the Péclet number, defined as Pe = u s a h /D.The Péclet number is assumed to be in the order of 10 4 -10 5 based on available study of the diffusion constant of synovial proteins [26].The study shows that the diffusion constant and thus the Péclet number are concentration dependant, however with advection dominating the equations the effect of concentration on the diffusion is negligible.

Flow Rate
The pressure solution of the Reynolds equation can be used to determine the flow rate of the fluid in the lubricated contact, such that where Q x, ∞ and Q y, ∞ are the fluid flow rate in the x-and y-direction respectively.In the cavitation region, flow rate is conserved such that Q x, ∞ downstream of the cavitation line is constant and equal to the flow rate at the point of cavitation [23].This model does not calculate the reformation point and thus calculated flow rates in this region are not entirely accurate, however flow downstream of the cavitation line is outside the area of interest.
If the protein matter is modelled as moving freely with the fluid, the concentration of protein matter throughout the contact zone would remain uniform over time.Instead, at points in the domain where the fluid film is thin enough, protein flow rate will differ from that of the fluid as the gap becomes too small for the macro particles to flow freely and unobstructed.This effect is achieved by introducing a scaling factor, ϕ that is a function of film height such that (2.25) The conditions of such a scaling function are that ϕ → 0 as H → 0 such that flow rate of protein mass decreases to zero as the dimensionless film height approaches zero, and also that ϕ → 1 as H → ∞ such that the flow of protein returns to the unobstructed Q ∞ as the film height becomes large.
The function proposed in this study is the hyperbolic tangent the behaviour of which can be tuned via the free parameter, a, as depicted in Fig. 6.

Constitutive Equation
As has been discussed in §1.1, fluid viscosity increases exponentially with increasing concentration of protein, as shown in experiments by Coleman et al. [11] and Castellanos et al. [10].It is therefore proposed in this study that the constitutive equation be of the form where η pf is the protein free viscosity i.e. at zero concentration and c 0 is the initial concentration at zero time.
Knowing the initial condition, η = η 0 , c = c 0 , leads to the following relation Substituting back into Eq.2.28 gives such that when non-dimensionalised, the constitutive equation becomes (2.31) In the experiments by Myant and Cann [14], the post-scratch central

Numerical Method
Calculating a solution to Eqs. 2.9-2.11involves an iterative procedure illustrated in Fig. 7.The Reynolds equation is used to solve for pressure which is then passed to the deformation and load balance equations.These equations are, respectively, used to calculate the deformation of the surface, v, and separation value, H 0 .This new height profile is then returned to the Reynolds equation and repeated until the solution converges.
For this, finite difference methods have been utilised such as those developed by Venner [27] to resolve the governing equations.This involves a second order central discretisation for the first two terms of the Reynolds equation and a first order upwind scheme for the wedge term: ∂(ρH)/∂X.The squeeze term, ∂(ρH)/∂T, is also first order, using only the previous film height solution.This term is neglected when calculating the initial condition.The resulting discretisation method for the Reynolds equation is: ) Here, ε = ρH 3 /ηλ and values taken at the intermediate locations are defined in the manner: To discretise the film thickness equation, use is made of a deformation coefficient matrix, K.This matrix stores the Hertzian deformation of every point in the domain subject to a unit point load at every other location in turn.These point load solutions can then be combined linearly to calculate the deformation for any pressure profile.For the study of a sphere on a plane, it is possible to calculate K analytically [27] as follows: where,

.35)
Using this approach, the film thickness equation can be written as: (2.36) The force balance equation can be calculated using any number of discretised numerical integration techniques.For this study, the trapezium rule is applied.
It is common to use a multigrid method for the solution of highly loaded EHL contact problems, however, due to the negligible effects of piezo-viscosity, in this instance a single grid was sufficient.Other than the removal of piezo-viscous modelling and the single grid simplification, the solution method remains the same as that put forward by Venner [27] using a combination of Gauss-Seidel line relaxation and Jacobi distributed line relaxation methods.The solution domain spans − 5 < X < 5, − 5 < Y < 5 with an N × N grid where N = 257.This was deemed sufficient for mesh independence by comparing the steady solution over a range of mesh densities, N = 2 5 + 1 to N = 2 9 + 1.
It is important that the flow rates calculated from the EHL pressure solution accurately conserve mass when used to model the advectiondiffusion of protein content.As the Reynolds equation is itself a conservation of mass equation, it is therefore required that the pressure solution is such that the residual of the Reynolds equation is sufficiently low.It is typically the case that EHL solvers are considered converged when the following conditions on the pressure and force balance are achieved [28].
Here, n denotes the iteration step and the solution scheme is terminated when the changes in pressure solution each iteration become smaller than a prescribed tolerance, as does the force imbalance between the summed pressure and load.The tolerance, ε tol is commonly in the region of 10 − 4 -10 − 6 .
This approach however, does not impose any explicit convergence tolerance on the residual of the Reynolds equation.It is found that although changes in pressure solution maybe be considered converged, the residuals may remain large as mass continuity is based on the second derivative of pressure.If an explicit condition on the residual of the Reynolds equation is introduced with a tolerance of 10 − 6 the EHL solver must run for several orders of magnitude longer to achieve convergence.Instead, the coupled equations are solved until a convergence condition on the film height is achieved and once the film height and force balance are sufficiently converged, it is possible to fix the height and solve only the Reynolds equation to achieve the desired residual convergence.This 'super-solver' for pressure only can use a simple Gauss-Seidel point-wise relaxation scheme with unity relaxation, but must still be solved iteratively to account for cavitation.Terminating once the residual has reached less than 10 − 4 was found to be adequate.The advection-diffusion problem was solved using a implicit Euler scheme.The domain is larger at − 10 < X < 10, − 10 < Y < 10, however because of the inherent computational expense of implicit methods, it is desirable not to increase the number of grid-points when expanding the domain.Instead a new grid N = 257 is used over the larger domain and is thus coarser.The pressure and film heights outside the EHL solution domain (− 5 < X, Y < 5) are zero and the undeformed shape respectively.This is the initial condition of the Reynolds equation 'super-solver' which smooths out any discontinuities in pressure gradient at the boundary of the old and new domain and also globally reduces the Reynolds residual to an acceptable level.Flow rates are then calculated from this solution, scaled with the flow factor, ϕ, and the implicit scheme calculates one time-step, tracking the protein concentration.
The constitutive relation is calculated after each advection-diffusion stage and used as the viscosity field for the next EHL timestep.The viscosity is linearly interpolated to the finer grid used for solving the EHL problem.

Parameter Search
A series of simulations were run spanning a range of values for the two free parametersa, the hyperbolic tangent gradient of Eq. 2.27 and η, the ratio of protein free viscosity to initial viscosity that defines the concentration relation, Eq. 2.31.The ranges are summarised below in Table 3 and give a total of 400 initial simulations.These initial simulations were allowed to run for 72 h before being terminated.This period proved adequate for determining the behaviour of fluid film development with respect to the parameter space.
The initial spread of film thickness curves in Fig. 8 is wide and in order to evaluate which curves best approximate the experimental results two factors were used.The first factor being the stable film height that the simulations approach.This film height plateau can be seen in some but not all of the simulations in Fig. 8.The results of experiment show that the film thickness reaches a stable value of 140 nmbased on the average of the last 20 data points to the nearest nanometre.
To quantify the similarity of the initial delay and subsequent growth rate of the lubricating film, the sum of squared residuals is calculated.These two evaluation methods were adequate to determine an optimal pair of parameters, η and a, whereby the simulation film thickness grows at a similar rate and plateaus at a similar point as experimental results.

Plateau Height
Multiple simulations have similar growth rates but plateau to different stable film thicknesses.An equation for how the parameters affect this stable height is found simply by plotting the end points of the simulations that appear to have reached a stable height.Simulations are deemed to have plateaued by testing the gradient has reduced to an acceptable level of 0.1 nm per timestep (12 nm/s).To ensure that only simulations that had stabilised after a significant growth phase were included, only simulations with a final film thickness of greater that 50 nm were plotted.This cut off line can be seen in Fig. 8 where below this point, film heights are either still increasing or have remained constant.Fig. 9 shows the final film thicknesses of simulations that have reached a stable height and a surface of the form has been fitted using the non-linear least squares regression.The resulting coefficients and fit parameters can be seen in Table 4. Fig. 9 also shows the fitted function contours of constant height, which form straight lines.For the desired stable film thickness of 140 nm this line has the form

Curve Similarity
To quantify the similarity of the initial delay and subsequent growth rate of the lubricating film, the sum of squared residuals is calculated.The sum of squared residuals (SSR) is a measure of curve similarity most often used in least squares regression whereby the sum is minimised to obtain a good model fit to data.
To calculate the SSR between the simulations and experiment, both must have the same sample rate.The simulations are run with a dimensionless time-step of 0.1 (which corresponds to 8.1064 × 10 − 4 s) for 2000 time-steps, however results are stored only every 10 steps such that the sample rate of the simulation results is once every 8.1064 ms with 200 data points excluding the initial solution.The experimental results are interpolated to this same sample rate, so that the SSR can be calculated.
To compare the individual simulations it is desirable to sum over the same number of samples, n.By comparing the simulation length of all 400 simulations in Fig. 10 it can be seen that the majority of simulations terminated between 1/3 and 2/3 of the total 200 recorded samples.There is a distinct group of simulations that do not reach this threshold and it is likely that these simulations became unstable and failed.The SSR is therefore calculated over the first 67 samples equating to the threshold at 1/3 of total simulation time (0.53 s of experiment time).This is an adequate length to analyse the initial growth behaviour of the central film thickness.Simulations that did not reach this threshold were rejected.
Fig. 11 shows the resulting contour plots of the sum of squared residuals of each simulation in the parameter space as compared to experiment.For clarity of the contour plot, the square root of the SSR is plotted.The plot does not have an apparent minimal point but is minimal along a line.This would not be the case if simulations reached the full 2000 time-steps.For full length simulations, the curve similarity test would include the full film thickness development and final film height plateau and would therefore be the only test necessary.Fig. 11 highlights the minimal values along the rows and columns of the parameter grid and for simulations with parameters below η ≈ 0.13 and a ≈ 0.2 fit a line with the equation a = 1.5η.
The two equations that give optimal SSR and stable height plateau are thus a = 1.5η a = − 5.22η + 0.698 (3.3) and can be used to define the optimal parameters for the model.Solving this system of equations gives η = 0.104 and a = 0.156.

Results and Discussion
Using the results of the parameter search, a full length simulation was run and the resulting transient central film thickness, shown in Fig. 12, can be seen to follow the growth rate found in experiments and also reach the stable height of 140 nm.A brief dip in the film thickness within the first 0.2 s is replicated however it is not of the same magnitude.It can be seen in Fig. 8 that none of the simulations within the parameter range considered express the extent of this initial thinning.It is most likely the case that this simplified model does not capture this phenomenon, which could potentially be because the of shear thinning behaviour or because of the assumption of uniform concentration across the contact after the collapse of the stable film height.
The optimal parameters were calibrated for Péclet number equal to 10 4 , however full simulations were also run at an order of magnitude above and below.Fig. 12 shows that at lower Péclet numbers, these parameters no longer give a solution that agrees well with experiments however as the Péclet number increases, changes to the solution are minimal, indicating that the parameters give solutions that agree with experiment for the entire range 10 4 -10 5 .

Table 3
Parameter range and step size for simulations.The unusual step size for η stems from the choice of range for dimensional variable, η pf , being 0.001-0.02Pas⋅s with a step size of 0.001 Pa⋅s.

Concentration and Viscosity
Snapshots of the instantaneous concentration field are shown in Fig. 13.Initially, protein matter aggregates at the inlet of the contact zone and the surface starts to deform noticeably in this region alone due to the increased viscosity and pressure.As the surfaces deform, the larger gap allows the protein aggregate to move slowly into the contact zone and deform the surfaces further downstream until the entire contact zone is flooded with high protein-content fluid.The fluid in the contact zone has a concentration over double that of the bovine calf serum that surrounds the contact region.The relative concentration reaches a maximum of C = 2.58, giving a relative viscosity peak of η = 35.7.This results in a dimensional peak protein concentration of around c = 46.4g/l and peak viscosity of η = 2.68 Pa⋅s.The surface geometry itself shifts from an almost flat contact zone, typical of elastohydrodynamic lubrication, to a relatively steeper converging wedge.
This effect is due to the high viscosity moving the contact zone from the isoviscous-elastic regime towards the isoviscous-rigid regime.To quantify the change in regimes, work by Dowson and Hamrock [29], later revised by Esfahanian and Hamrock [30] mapped the approximate domains of lubrication regimes against a dimensionless elasticity parameter, g E , and dimensionless viscosity parameter, g V .The parameters are defined as where W, U, and G, are the dimensionless load, speed and material parameters respectively as introduced by Dowson and Higginson [31].Fig. 14 shows the diagrammatic regime map of fluid lubrication regimes and an illustration of how increasing viscosity affects the behaviour of the lubricated contact.For the purposes of this figure, the pressureviscosity coefficient of bovine calf serum is taken to be of the same order as water, being roughly α = 1 GPa⋅s − 1 [32].Initially, the bulk viscosity throughout the contact zone is η = 0.075 Pa⋅s, which results in parameters g E and g V that are firmly in the isoviscous-elastic regime.The maximum viscosity reached in the PAL simulation was η = 2.68 Pa⋅s, and using this value the parameters g E and g V have moved into the

Table 4
Coefficients and goodness of fit parameters for final film height response surface.The same contours also showing the points which are minimal in both the x-and y-directions, the line that these points lie upon (dotted) and also the line upon which simulations lead to a plateau of 140 nm (dashed).The intersection of these two lines show the optimal parameters for matching experimental results.isoviscous-rigid regime, whereby the contact behaves like a conventional hydrodynamic converging wedge.

Film Height and Inlet Length
By using the interferometer images of the experiment by Myant and Cann [14], it is possible to compare the transient shape of the film height in the contact zone and also the size and shape of the inlet aggregate.
of height contours between simulation and experiment can be seen in Fig. 15.The deformation of the contact zone shape follows similar behaviour whereby the deformation occurs initially only at the inlet side of the contact zone.Differences between simulation and experiment are visible from 0.8 s onwards as the inlet aggregate forms a channel through the contact zone in the experiment but not in the simulation.As seen in Fig. 13 the protein concentration increases across the width of the contact in the modelled simulation.
The innermost ring, indicating 225 nm film thickness, continues to decrease throughout the simulation due to the continued separation of the two surfaces whereas the channel that can be observed in the experiment leads to an increased central film thickness without this continued global separation.As a result, the innermost ring of experimental images decreases in size less so than that of simulations, but there is a consistent agreement throughout between the location of the innermost contours at the inlet.
If the transient height is viewed along the line y = 0, the simulations agree well with experiment.Fig. 16 shows the full simulation height profile plotted with the contour location from the experimental interferogram images.It can be seen that the behaviour whereby the contact zone deforms first at the inlet, while the central film thickness remains small, is well captured by the PAL model.
The formation of the high concentration channel can also be seen in Fig. 17 whereby high concentration regions have been plotted overlaying interferometer images.The apparent regions of protein aggregate for experiments have been highlighted, and also shown are the regions for which protein concentration has doubled.The factor of two for concentration increase is arbitrary but shows comparable inlet lengths for visualisation.Unlike experiments, the simulation concentration shows a much broader accumulation across the y-direction, leading to a contact zone that is eventually flooded with high protein fluid.
As well as the shape, it is possible to compare the transient development of the inlet length between the interferometer images and that of simulation.Figure 18a shows how the inlet length, s, is defined for interferometer images.For the simulation it is defined in a similar way as the distance between when the dimensionless protein concentration becomes C = 2 and the innermost height contour of h = 225 nm as seen in Fig. 11.The selection of C = 2 is arbitrary, but offers good comparison when analysing the behaviour of the inlet length over time.It can be seen that the inlet develops linearly in experiments, reaching a steady value coinciding with the film height plateau.In simulation, however, the inlet length develops rapidly to begin with and grows only slowly thereafter.The initial jump is due to the time before a value of C = 2 is reached.This disparity in inlet length can also be seen in Fig. 17 whereby the contours of simulation extend further upstream than that of experiments at time-step t = 0.4 s.
The behaviour that the simulations do not manage to capture is the localised nature of the aggregate build up.The flow rate attenuation function, ϕ, does not change over time, but physically the nature of the size distribution of protein matter particulates will indeed change as more macromolecules aggregate.This will lead to an increase of large particles at the inlet as more protein matter is accumulated here.This effect could be modelled with a concentration-dependent function, ϕ(H, C), although the exact nature of aggregate size distribution is more complex than this and will likely also depend on the shearing behaviour of the fluid.

Pressure
Because the concentration initially aggregates exclusively in the inlet region, noticeable changes in viscosity and pressure occur in this region almost immediately.Fig. 19 shows the transient nature of the pressure in the inlet and contact zone over the course of a second.Initially the pressure distribution is typical for a highly loaded elastohydrodynamically lubricated contact but soon develops a second peak

Varying Velocity Conditions
Using the parameters, a and η, it is possible to run simulations of varying sliding velocity.Although there are no transient measurements for the range of sliding velocities, it is noted that one of the phenomena of protein aggregation lubrication is that steady film thicknesses were observed to decrease as sliding velocity increased [12].
Fig. 20 shows that as sliding velocity increases, so too does the final steady central film thickness of the lubricating contact.This is consistent with EHL theory but not of the phenomenon observed in experiments.It is believed that the decrease in film thickness with increased sliding velocity occurs because the higher momentum fluid prevents protein aggregation to the same extent.Work by Pritchard, Wilson et al. [33,34] uses a shear-dependent 'build-up' and 'break-down' term in a scalar transport equation to simulate aggregation behaviour in colloidal systems.This approach combined with a concentration-dependent attenuated flow rate would make for an interesting avenue for further research.

Conclusions
Protein aggregation lubrication can be modelled by using classical elasto-hydrodynamic lubrication theory coupled with a protein concentration-dependent constitutive equation for fluid viscosity.The transport of protein matter in the thin film lubricated contact is subject to an advection-diffusion model using geometry-based attenuation of flow rates.To maintain mass continuity when using the pressure profile within the contact zone to derive fluid flow rates, a strict convergence criterion for the residual of the Reynolds equation when calculating the elasto-hydrodynamic pressure and film thickness is required.This can be achieved with the use of a decoupled Reynolds equation solver.
Simulation results were compared to experiments of a cobalt-chromium-molybdenum alloy femoral head implant tested in a pin-on-flat tribometer using diluted bovine calf serum as a synovial fluid analogue.It was shown that simulated central film thickness results agree well with previous experimentally-observed behaviour by using  an attenuation function to limit protein transport and an exponential, concentration-dependent, local viscosity.The results show high concentration gradients in the lubricated contact zone and care must be taken to ensure that an appropriate advection scheme is used that can model this phenomenon.
By comparing film height profiles with interferogram results of experiment, it can be seen that the geometry of the surfaces develop from an almost flat contact zone, typical of EHL solutions, to a steeper converging wedge.This effect is due to the rapid increase in fluid viscosity at the inlet as protein matter aggregates and represents a move from the isoviscous-elastic regime towards the isoviscous-rigid regime.
Although the transient development of large central film thickness is reproduced well by this computational PAL model, more work must be carried out to better capture the shape of the inlet aggregate and effects of sliding speed.
Further refinement of the model could also be studied, including closer investigation of the flow attenuation function and protein aggregate size distribution.It may also be necessary to include concentration and time dependence in the attenuation function to capture the transient nature of the size of aggregated macromolecules.Other well-studied features such as the effect of surface roughness, surface adsorption and the shear-dependence of aggregation could also be introduced to better capture the behaviour of protein aggregation lubrication and the phenomenon of decreasing film thickness with increasing sliding velocity, which was not captured by the current model.
The implication of protein aggregation lubrication on the design of prosthetic joints is most notably that transient phenomena occur even at steady operating conditions.This study has shown that additions to current elasto-hydrodynamic modelling techniques are needed to capture these phenomena.These results show that for steady operating conditions, contact zones are larger and contact pressures are lower for protein-rich lubricating fluids than EHL analysis would suggest.
The inclusion of protein transport modelling in real geometries and gait cycles is likely to be of increased importance.The repeated build up and break down of inlet aggregation, and therefore protein concentration have been shown to have a significant impact on the film thickness and, by implication, wear of the joint.It is likely that the design of the joint and materials not only influence the fluid film thickness and pressure, but also influence the transport of protein concentration, which in turn has a further influence on the film thickness.The health of the joint and therefore the protein composition of the synovial fluid is also likely to have an impact on the lubrication of the joint and is something which could influence optimum joint design.

Fig. 1 .
Fig. 1.Interferogram of contact zone during constant load and sliding speed experiment.Flow direction is from top to bottom and a clear distinction between bulk fluid and fluid in the inlet region can be seen.The length of the inlet region, s, is also shown.Figure adapted from Myant and Cann [12].

Fig. 2 .
Fig. 2. Results of the inlet disruption test by Myant and Cann[14] showing transient behaviour of the central film thickness, h cen , (black) and inlet aggregate length, s, (red) as a scratch removes the protein matter at the inlet to the contact zone.The subsequent regime whereby the inlet aggregate is reformed is highlighted.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Interferometer images of experiment by Myant and Cann [14] at select time-steps after the beginning of the protein build up regime.

Fig. 4 .
Fig.4.Geometry of a sphere-on-plane lubrication problem in which the lower surface has a velocity, u s , the sphere has an undeformed radius, R, and applied load, w.The film has height h = h(x, y, t) in the z-direction.

Fig. 5 .
Fig. 5. Overview of transient PAL model and individual blocks that create the full transient solver.

Fig. 6 .
Fig. 6.Proposed scaling factor function for limiting flow rate of protein transportation.The gradient parameter, a, can be used for tuning the functionvarious values of a shown in grey.

Fig. 7 .
Fig. 7. Flow diagram of elasto-hydrodynamic lubrication solution method for each transient time step.

Fig. 8 .
Fig. 8. Initial parameter search showing range of simulation solutions (solid grey) as compared to experimental results (black points).Also shown are dotted lines indicating 50 nm, used for determining solutions that have reached a stable plateau and of 0.535 s, before which is the regime in which the curve similarity is calculated.

Fig. 9 .
Fig. 9. (Left) Final film thickness of simulations that have reached a steady plateau and the best fit surface resulting from a non-linear least squares regression.(Right) Contours of constant film height for this surface.

Fig. 10 .
Fig. 10.Comparison of simulation length across all 400 simulations.The majority of simulations terminated between 1/3 and 2/3 of the total 200 recorded samples.

Fig. 11 .
Fig. 11.(Left) Contour plot of the sum of squared residuals (SSR) over the parameter space.(Right)The same contours also showing the points which are minimal in both the x-and y-directions, the line that these points lie upon (dotted) and also the line upon which simulations lead to a plateau of 140 nm (dashed).The intersection of these two lines show the optimal parameters for matching experimental results.

Fig. 12 .
Fig. 12. Transient central film thickness plots for experiment and simulations showing increasing film thickness development over time.The optimal parameters, η = 0.104, a = 0.156, have been run for the Péclet number, Pe = 10 4 , and also for Péclet numbers an order of magnitude higher and lower.

Fig. 13 .
Fig. 13.Contour plots of protein concentration in the inlet and contact regions for the transient build up regime.The central film thickness plots (top) show the current point in time.The slices at y = 0 (middle) show film height development along with moving protein aggregate and the top down view (bottom) shows the full concentration profile over the 2D computational domain.

Fig. 14 .
Fig. 14.Fluid film lubrication regime map indicating the effect of increasing viscosity on the behaviour of the contact zone of the protein aggregation lubrication simulation.

Fig. 15 .
Fig. 15.Height contours of contact zone for simulation (solid red rings) shown overlaid on interferometer images created from data of experiment by Myant and Cann [14].The red rings show height contours of 225 nm matching the quarter wavelength of the interferometer image.Flow direction is from left to right.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 16 .
Fig. 16.Comparison of transient height profile along the line y = 0 between simulation and data from experimental interferogram.The experimental results consist of the height contour locations and the central film thickness measurements.

Fig. 17 .
Fig. 17.Contours of C = 2 for simulation (dotted red) shown overlaid on interferometer images created from data of experiment by Myant and Cann [14].Also highlighted are the apparent protein aggregate regions of experimental images (black dotted).Flow direction is from left to right.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 18 .
Fig. 18.(a) Example of inlet length, s, as defined as the length between the apparent start of the aggregate matter to the inlet of the contact zone.(b) Comparison of measured inlet length in experiment by Myant and Cann [14] to inlet length, s, of simulation using C = 2 as the starting point of inlet aggregate measurement.

Fig. 19 .
Fig. 19.Pressure results (right axis) along y = 0 over the first second of protein build up.Also shown in grey (left axis) is the transient film height profile at y = 0 and velocity profiles of the lubricating fluid.

Fig. 20 .
Fig. 20.Effect of varying sliding velocity on (a) transient central film thickness development and (b) on final steady central film height.Also shown in plot (b) is the resulting steady film thickness of classical EHL solutions.

Table 2
[14]rimental parameters for the inlet disruption study by Myant and Cann[14]on pure sliding contact of a stationary femoral head on 20 mm/s sliding glass.