Transient mixed lubrication model of the human knee implant

The human knee implant is computationally modelled in the mixed lubrication regime to investigate the tribological performance of the implant. This model includes the complex geometry of the implant components, unlike elliptical contact models that approximate this geometry. Film thickness and pressure results are presented for an ISO gait cycle to determine the lubrication regime present within the implant during its operation. It was found that it was possible for the lubrication regime to span between elastohydrodynamic, mixed and boundary lubrication depending on the operating conditions of the implant. It was observed that the tribological conditions present in one condyle were not necessarily representative of the other. Multiple points of contact were found within the same condyle, which cannot be computed by the elliptical contact solvers. This model can be used to balance forces in all directions, instead of only the normal loads, as often done in elliptical contact models. This work is an initial step towards understanding the role of the complex geometry in the tribological characteristics of the human knee implant when

poroelastic lubrication and even after I left Leeds, I would try to catch up with Duncan whenever I was in Leeds.I remember talking about becoming a father and Duncan sharing his own experiences and occasional advice with me.
Around 20 years ago, I (Leiming) first heard the big name of Prof Duncan Dowson from books, papers and lectures when I started a MEng course in Tribology at Qingdao China.Prof Yang who supervised my project had a chance to visit Leeds and he discussed with Duncan and Prof Jin about my project on a transient TEHL simulation with surface texturing.I remember how happy I was to have Duncan's comments and approval to co-publish my first paper.The Chinese girl at that time could not imagine that one day she would work at the same building, meet Duncan occasionally at the corridor or Jin's office, visit him regularly on a research project in hip joint non-Newtonian lubrication and even have an afternoon tea at his place trying to guess who was the reviewer of our submitted manuscript from their writing style.I just feel I was so lucky as for many overseas academics in Tribology, Duncan is the one who was standing at the lectern followed by a long queue of people waiting for taking a photo or shaking hands with him.There were so many memorable moments during the time working with Duncan.To me, he is a great idol full of passion for science, and being rigorous in research, and that has encouraged me to keep going on.
I (Gregory de Boer) was introduced to Duncan while working as a Research & Teaching Fellow at the School of Mechanical Engineering, University of Leeds in 2016 by Robert Hewson.Having undertaken my PhD in the multiscale modelling of elastohydrodynamic lubrication at Leeds, I was in awe of his academic career and humbled by his humility when discussing his past.During these meetings, we would talk about the problem of modelling poroelasticity and lubrication mechanisms in human joints, from which the direction our research developed, and his experiences were vital for us in achieving meaningful results.Later in 2018, I collaborated further with Duncan on a new approach to modelling hydrodynamic cavitation.These discussions were a fantastic experience from which I learnt a lot from him about the challenges involved in developing his famous EHL solutions in the 1950s, his research career as an academic at Leeds, his family life, and his opinions on my new ideas for modelling cavitation phenomena.

| Background of knee joint replacements
Degenerative action, whether due to natural ageing, disease or injury, is a common occurrence during the lifetime of a human joint.The number of joint replacement surgeries are increasing rapidly.There are more than 0.1 million total knee replacements surgeries performed per year in the United Kingdom [1].A total of 3.48 million patients per year in the United States alone are expected to undergo knee replacement surgeries by 2030, with hip replacement surgeries expected for 0.57 million patients.Between 2005 and 2030, this translates to an increase of 174% for hip and 673% for knee arthroplasties [2][3][4].
Current designs suffer from wear and loosening, which together account for more than half of knee implant failures [5][6][7][8][9][10].The debris from wear can also cause inflammation due to periprosthetic osteolysis, resulting in implant failures due to infections.Instability of the implants is responsible for roughly 20% of implant failures, with misalignment causing a further 10% of implant failures [7].Misalignment of the implant components is a difficult issue for surgeons to address due to surgical and patient variations, causing complications such as varus osteoarthiritis [11].Revision surgery is then required for correction, causing medical risks, discomfort and psychological trauma associated with such surgeries.From 2005 to 2030, the number of revision hip arthroplasties are expected to increase by 137% in the United States.The number of revision knee arthroplasties are projected to increase by 601% [2,3].
The annual global economic impact of total knee arthroplasty surgeries amounts to $40 billion, an increase of 250% between 2005 and 2015.Revision surgeries cost a further $4 billion annually to carry out, representing an increase of 450% in the same period, with these costs projected to rise in the future [12].Note that these costs do not include hospital and other charges; they only represent the cost of the surgery.
Considering the increase in the number of projected arthroplasty surgeries, the high percentage of implant failures caused by wear, coupled with their medical risks and economic impact, attempts have been made to optimise implants to reduce wear and consequently, revision surgeries.There has been some success in reducing wear [13], but investigation into minimising wear is still ongoing [14].Highly cross-linked polyethylene bearings were introduced to address wear to some extent [15], but wear optimisation of knee implants still remains a challenge.
Experimental investigation has been carried out in joint simulators to help understand wear [10,16].Simulators work by repeating the human gait over millions of cycles, with wear measured post-test.This has provided valuable insights into the tribological behaviour of joint implants but is limited by the nature of the measurement process; that is the transient, realtime behaviour is often not recorded.In addition, this is both a time consuming and expensive process, especially when attempts at optimisation of the implant geometry are undertaken.

| Knee lubrication modelling
Computational modelling of the human knee implant has the potential to overcome these limitations by carrying out resource efficient optimisation of the implant design.This is performed in order to minimise wear and improve implant performance, increasing the reliability of the implant, and is the focus of this work.
These studies suggest that the geometry of the opposing surfaces determines the pressure distribution, highlighting the importance of accurately representing the geometry in lubrication modelling.Studies that include the full, complex geometry of the knee implant have not been carried out to investigate the accuracy of elliptical contact approximations and this provides the motivation for this work.
This work investigates the accuracy of elliptical contact models in capturing the tribological phenomenon present in knee implants during its operation.This is accomplished by performing a comparison with a complex geometry model of the knee implant across the elastohydrodynamic, mixed and boundary lubrication regime, in order to highlight the significance of simulating the complex geometry of the implant.The computational complexity of the complex geometry model presented in this work is also compared with an elliptical contact model.This study can help provide a better understanding of wear and facilitate the optimisation of such implants in order to improve their performance.
The model incorporates statistically described surface roughness and asperity contact mechanics.The results are presented for various gait cycles across a range of physiological viscosities to determine the lubrication regime in which the knee implant operates in.Surface parameters are also varied to investigate their effect on the lubrication regime.

| Numerical method
The three-dimensional (3-D) geometry of the right knee implant is illustrated in Figure 1 and can be obtained in online repositories as a bone-implant model [26].For the purposes of this work, the tibial component was considered stationary and its location is fixed, while the femoral component could rotate and translate freely.
The location of the femoral component was described relative to the tibial insert using three parameters as explained below and shown in Figure 2.
Δz: the displacement in the global z direction between the highest point on the tibial insert and the lowest point on the femoral component.Δx: the displacement along the global x axis from the central point of the femoral component.Only the displacements in the x and z directions were considered.θ x : the counter-clockwise rotation about the centre of the femoral component parallel to the y axis, represents the flexion angle of the knee implant.Internal-external or adduction-abduction rotations are not considered.

| Film thickness
The complex geometry was interpolated onto a 2-D finite difference grid by expressing the height of the 3-D geometry in the z direction as a function of x and y for both implant components.It was assumed that the geometry changes linearly between two successive finite difference mesh points.
The two-dimensional film thickness used in this work was then determined by finding the shortest distance from each discrete finite difference mesh point on the surface of the tibial The orientation of the film thickness in the global coordinate system was found by θ and β represent the orientation of the film thickness direction with the global x and y axes, respectively.x 1 ; y 1 ; z 1 are coordinates of P, while x 2 ; y 2 ; z 2 are coordinates of P 0 , as shown in Figure 3.
If no surface roughness was considered, the resulting smooth film thickness (h s ) was described as follows: where Δh is the deformation of the contacting surfaces.

| Surface velocities
In addition to the surface velocities that arise from the translation of the surfaces (u t , v t , w t ), the global surface velocities (u, v, w) also depend on the rotation of the surfaces.However, to account for the curvature of the surfaces, the surface velocity components used in the Reynolds equation must be local to the coordinate system of the film thickness.Using the values of β and θ calculated in Section 2.1, where θ and β represent the orientation of the film thickness direction with respect to the global x and y axes, respectively, the relationship between the global and local surface velocities was described as follows: The local 'squeeze' term was described as follows:

| Hydrodynamic pressure
The applied load was supported by the fluid film pressure, governed by the Reynolds equation described below [27]: In this equation, p f represents the hydrodynamic pressure, h is the film thickness, η represents the viscosity of the lubricating fluid, and u x and u x represent the average local surface velocity in the x and y directions, respectively.Its derivation from the Navier-Stokes equations can be found in literature [27,28].

| Surface roughness
The effect of surface roughness was stochastically modelled based on the flow factors approach of Patir and Cheng [29].The surface of Cobalt-Chrome femoral component was considered to be significantly smoother than that of the tibial insert [30].As a simplification, surface roughness was assumed to exist on one of the contacting surfaces only and was represented by the height of the rough surface α from the reference smooth film thickness.The surface roughness was described using a Gaussian probability distribution function ϕ r and the height of the rough surface α from the reference smooth film thickness:

| Surface deformation
Deformation was only assumed to occur on the surface of the tibial insert.This was because the Young's modulus of the polyethylene tibial insert (500 MPa to 8.1 GPa) was significantly lower than that of the harder Cobalt-Chrome femoral component (200 GPa) [31].In this work, the modulus of 1 GPa was used and there was no particular reason to favour The finite element method (FEM) was used to calculate deformation.The surface was discretised into a mesh and a point load is applied at a single mesh point, with the resulting deformation throughout the surface calculated using the finite element method [32].This deformation represented a single column of a global deformation matrix, which was computed on the finite element mesh.This process was continued for each mesh point and the deformation the point load caused throughout the surface corresponded to its respective column in this global deformation matrix.This resulted in a dense deformation matrix of size ðnx ⋅ nyÞ � ðnx ⋅ nyÞ for a domain of ðnx � nyÞ finite element mesh points.
After the deformation matrix was calculated on the finite element mesh, linear interpolation was performed to transfer the deformation matrix from the finite element mesh to the finite difference mesh, upon which the two-dimensional lubrication equations were solved.The interpolated matrix K had a size of ðnx ⋅ nyÞ � ðnx ⋅ nyÞ for a finite difference pressure domain P of ðnx � nyÞ hydrodynamic mesh points.Deformation (Δh) was calculated by performing a matrix multiplication of the total pressure with this interpolated deformation matrix (K): It was found that performing the matrix multiplication in Equation ( 9) was computationally intensive, especially when it must be performed repeatedly in the iterative process throughout the gait cycle.Additionally, the storage of the elements that constituted the full deformation matrix was also not feasible for the fine finite difference mesh for the lubrication used in this work.
Consequently, simplifications were made in order to speed up the solution process.Upon observation of the finite element solution, it was found that the most influential coefficients of the deformation matrix were within a few elements of the location of the applied load.This was verified by keeping only the maximum value of each column and setting all other elements of the deformation matrix to zero.Subsequently, the deformation matrix was transformed into a diagonal matrix, which simplified the matrix multiplication to an elementwise multiplication in a reduced order, which was significantly faster to compute.This reduced the number of non-zero elements of the deformation matrix from ðnx ⋅ nyÞ � ðnx ⋅ nyÞ to nx � ny.These non-zero diagonal terms were collected in a simplified deformation matrix (K s ) of the same dimensions as the finite difference mesh.
In this work, the 'simplified' deformation matrix refers to the K s matrix, while the 'full' deformation matrix refers to the deformation matrix K after it was interpolated from the finite element mesh onto the hydrodynamic mesh.
The calculation in Equation ( 9) with the full deformation matrix can be represented as follows: If the full deformation matrix is approximated by the simplified deformation matrix, the calculation in Equation ( 10) can be represented as follows: It is important to note that while the numerical verification performed in Section 2.10 included this simplification, it was specific to the model used in this work; a change in the geometry, material properties or the applied load would require these assumptions to be revisited.

| Asperity contact
Asperity contact was modelled based on the work of Greenwood and Williamson [29] and Greenwood and Tripp [33].Asperity contact was assumed to occur when the asperity height α was larger than the film thickness and is based on a statistical model of surface roughness.The probability of a single asperity making contact was described by a Gaussian probability distribution function ϕ a , as used by Greenwood and Tripp [33].The total local area of contact (A asp ) and the load borne by the asperities (F asp ) was defined if the total number of asperities (N asp ) was known within the local area (dA): It was assumed that the asperities deform according to the Hertzian contact theory, with E0 asp describing the effective Young's Modulus of the asperities and r asp representing the radius of the asperities.The total number of asperities can also be defined by their density (β), for convenience.

-
BUTT ET AL.

| Load balancing
The total applied load F was borne both by the lubricating fluid (F f ) and solid asperities (F s ): Investigating the ratio between F f and F s can give insight into the nature of the lubrication regime.A higher percentage of F f signifies the dominance of elastohydrodynamic lubrication, while the inverse signifies solid contact and boundary lubrication.The complex geometry model presented in this work includes mixed lubrication in order to accurately capture this phenomenon and to investigate the transition between the lubrication regimes.In this work, the ratio F s : F is referred to as the ratio of the total load borne by the solid asperities.

| Gait cycle
Transient loading conditions used in this work are those derived from ISO 14243-3 for a total duration of 1 s [34].The following data was obtained from the gait cycle and shown in Figures 4 and 5 [35].
The displacement and rotation of the components can be used to calculate the surface velocity required for the solution of the Reynolds equation.For the sake of clarity, this gait cycle is referred to as the ISO gait cycle throughout this work.Elliptical contact models often use similar data, as their inability to bear loads in the x and y directions is often compensated by prescribing the displacement in these directions instead.

| Discretised equations
Discretisation of the Reynolds equation for the finite difference method was undertaken using a second order central approximation for the pressure terms and second order upstream for the wedge and transient squeeze terms.
The surface roughness, defined in Equations ( 7) and (8), was described using a Gaussian probability density function, using the mean surface roughness (μ r ) and its standard deviation (σ r ), along with the discretised window (Δ) used in the integration.Consequently, the film thickness terms were expressed as follows: The calculation of asperity contact parameters, defined in Equations ( 12) and ( 13), was also discretised using a Gaussian probability density function, given the standard deviation (σ asp ) and mean (μ asp ) for the asperity heights: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The total pressure was due to both the hydrodynamic pressure and asperity contact in a unit area: Deformation, when using the simplified deformation matrix K s , was calculated by Deformation, when using the full deformation matrix (K), was calculated by Δh The 3-D force balance equations were discretised as follows: n zi;j : p f i;j : The values for Δx, and Δz used in the current iteration (n) are updated using a factor ω f : The gait cycle was repeated until the pressure and film thickness results of the last timestep (time between the discretised times) of the current cycle were within 5% of the previous cycle, as prior studies have revealed that a single sweep through the gait cycle is not enough for the fluid film to reach a periodic solution [17,19].At each timestep, the criteria of numerical convergence for the pressure, film thickness and force balance were set as 1%.

| Mesh independent analysis and verification
The finite element and hydrodynamic meshes were varied in density until the pressure and film thickness results were within 10% of the next finer mesh tested.The finite element analysis of the tibial component involved using an adaptive mesh that was finer around the condyles and coarser everywhere else.It was important that the finite element mesh accurately captured the deformation in the contact area, for which the element size needed to be sufficiently smaller than the contact area.The maximum element size was varied between 1.5 and 0.1 mm.This was done by setting the maximum element size and calculating the deformation matrix from the resulting finite element mesh.The hydrodynamic mesh was varied between 512 � 512 and 4096 � 4096 grid points per condyle.The configurations used in the mesh dependence analysis are shown in Table 1 with the results shown in Figure 6.Configuration C is implemented in this work.
Timestep analysis was conducted by varying the number of timesteps between 50 and 200, for the mesh parameters in configuration C for a complete gait cycle (Figure 7).A value of 100 timesteps was used in this work.
To verify the interpolation of the deformation matrix, as well as the decision to retain only the leading terms of the matrix (i.e. using K s instead of K to calculate deformation), the converged pressure solution for a timestep given by the complex geometry model presented in this work was applied on the tibial component in COMSOL to check whether the maximum deformation matched that of the solver used in this work.For the configuration chosen in this work, the discrepancy was less than 10%.

| RESULTS
The results for the calculated film thickness, pressure and loading are presented in Figure 8 for a viscosity of 0.1 Pas, where the standard deviation of the asperity heights (σ asp ) and 212 - surface roughness (σ r ) was 1 μm with a mean μ asp ¼ μ r ¼ 0. The density of asperities (β) was 50 mm −2 , with a radius (r asp ) of 1 μm.
As seen in Figure 8a,b, decrease in film thickness results in an increase in fluid pressure.The condyle experiencing the global maximum pressure and minimum film thickness changes during the gait cycle, due to a shift in the load borne by each condyle, which can be seen in Figure 8c.The majority of the load was initially borne by the medial condyle until about 25% of the gait cycle, after which it shifts to the lateral condyle.This change is mirrored in Figure 8d, where the heavier-loaded condyle has a higher percentage of the load borne by the solid asperities.Consequently, the lubrication regime present in one condyle does not necessarily reflect the lubrication regime at the other condyle.
Figure 8c also shows that the overall load cannot be prescribed individually to each condyle, as is often done in elliptical contact models.During the gait cycle, the load shifts between the condyles and this can only be accurately captured if the condyles are modelled together with their complex geometries.It is evident from Figure 8d that even with a constant fluid viscosity throughout the gait cycle, the lubrication regime shifts between boundary, mixed and elastohydrodynamic lubrication for the simulated parameters, as the ratio of load borne by the asperities varies between 0% and 80% throughout the cycle.Note that a value of 0% signifies pure elastohydrodynamic lubrication, while 100% signifies purely solid contact.
Figure 9 shows the pressure results mapped onto the tibia at various timesteps during the simulation.Contact points can be characterised as localised regions of high pressure.The number of contact points vary throughout the gait cycle, as seen in Figure 9, depending on the gait parameters.A single contact point can be seen in Figure 9a for the medial condyle.Note that in Figure 9b, there are two points of contact per condyle, while in Figure 9c, there are two points of contacts in the lateral condyle but none in the medial condyle.The location of the contact points also varies throughout the gait cycle.This is true even if the number of contact points for a condyle does not change between timesteps.This is visible in the difference between Figure 9b,c, where the location of the contact points changes significantly.
If multiple high-pressure points exist in close proximity, it was possible for the regions of non-zero pressure around these high-pressure points to overlap, causing individual points to lose definition.In Figure 9b, two local pressure peaks are defined in the medial condyle.However, a relatively lower pressure region of around 5 MPa surrounds the two pressure peaks.This observation is significant because it means that any solution method must consider the possibility of overlapping zones of non-zero pressures between local pressure peaks.A similar observation can be made about such points in the lateral condyle.This is in comparison to Figure 9c, where all the high-pressure points are well defined and isolated, but this is not guaranteed at all points of the gait cycle, as seen in Figure 9b.A cross-sectional view of the pressure profile shown in Figure 9a can be seen in Figure 10, where a single well- defined peak can be seen.However, this is not the case in Figure 11, where multiple local peaks can be seen in each condyle in the y-direction, corresponding to the high-pressure regions seen in Figure 9b.
Regions of high pressure corresponded to regions with lower film thicknesses.This can be seen in Figure 12, which corresponds to the pressure results shown in Figure 9a.
The pressure results shown in Figure 9c represent a special case; as seen in Figure 8d, the condyles are subject to purely elastohydrodynamic lubrication beyond 80% of the gait cycle.This results in a lower peak pressure and the peaks occupy a smaller area in comparison to Figure 9a,b, which represent pressure results in the mixed lubrication regime.Pure boundary lubrication was not experienced in this gait cycle because the load borne by the asperities did not exceed 80%, as seen in Figure 8d.

| DISCUSSION
Elliptical contact models solve for a single point of contact in a condyle.This is due to the limitation of the simulated geometry; an elliptical geometry has a single point of contact with a planar surface [36].The results shown in Figure 9b for the complex geometry model presented in this work indicated that it is possible for multiple points of contact to exist within a condyle.If the contact points were to overlap, this could result in an irregularly defined contact zone that can only be found if the complete geometry of the knee implant is simulated.
A finite element modelling of knee joint contact stress based on in vivo kinematics [37] hints at multiple points of contact during the gait cycle.
The sensitivity of the pressure and film thickness results to the complex geometry of the implant can be further studied by varying the geometric profile of the implants, such as the thickness of the tibial plateau, as well as the properties of the materials.Different combinations of Young's Modulus can be studied to examine the effects of such changes on the results.The deformation of such a complex surface is much different from the simplified half-infinite plane model in the Hertzian contact theory, which highlighted the importance of using the FEM in the calculation of the deformation.
During the solution process of elliptical contact models, the loads on each individual condyle are often prescribed as a predefined ratio of the total applied load.It was found that the load borne by the condyles varied to the extent where these ratios cannot be pre-defined.
Elliptical contact models are unable to balance loads in the non-normal directions because the geometric simplification results in an elliptical geometry contacting with an infinitely flat surface.These loads exist even for the ISO gait cycle, as can be seen in Figure 13, though they do not need to be calculated as part of the solution process.The gait cycle used in this work only considers rotations about the y-axis, that is, the flexion angle of the knee implant.A more realistic model of the knee implant could consider all 3-D rotations and displacements, as well as internal-external rotations.Such a model could also solve the moments arising from the free movement of the implant.In the presented model, these moments are not considered.However, the physiology of the knee is complicated, in that there is a complex layer of tissue surrounding the knee implants, which can affect the extent to which the knee implant is free to move.These restrictions could be incorporated into the model, though the physiology of the knee is highly individualised, which can result in patient-specific restrictions to the movement of the knee implant.The effect of different gait cycles can be studied to investigate the difference in the pressure and film thickness profiles obtained under different loading conditions.
The finite difference and finite element mesh density used in this work can be increased in order to obtain more accurate results, giving additional insight into the sensitivity of the solution to the choice of mesh density, especially if different implant models are used.The mesh sensitivity analysis performed in this work is specific to the implant model and the loading conditions presented in this work.
While the viscosity was varied to investigate its effect on the film thickness and pressure results, the complex geometry model assumes an isoviscous lubricating fluid between the implant component.A more complex pressure- In this work, discretisation of the finite difference method was undertaken using a second order central approximation for the pressure terms and second order upstream for the wedge and transient squeeze terms.The effect of other discretisation techniques on the stability and speed of the solution can also be investigated.Other numerical methods that have been used in elliptical contact models, such as wavelet [38] and homotopy [39] methods, can be employed in the complex geometry model to investigate the accuracy of the results and the rate of solution convergence.
The surface roughness and asperity contact models used in this work are only based on Gaussian probability density functions.The effect of incorporating other probability distributions may be studied in the future.
The simplified deformation matrix (K s ) is an approximation of the full deformation matrix (K), made to facilitate a faster and less resource intensive solution at the expense of accuracy.The effect of this approximation can be studied to improve the accuracy of the calculated deformation by varying the number of terms kept in the deformation matrix.In the presented model, only the diagonal terms of the full deformation matrix were used to calculate deformation, but the effect of improving this interpolation can be considered in future work.
There was an approximation in the calculation of velocity components in local coordinate system of the film thickness which is always normal to the surface.The angles β and θ (representing the orientation of the film thickness direction) do not have any order to them, but the rotation angles used to calculate the projection of coordinate systems in Equation (4) do.An accurate coordinate transformation will be addressed in future work.
The sensitivity of the pressure and film thickness results to the complex geometry of the implant can be further studied by varying the geometric profile of the implants, as well as the properties of the materials.Different combinations of Young's Modulus can be studied in order to examine the effects of such changes on the results.The deformation of such a complex surface is much different from the simplified half-infinite plane The formation of the multiple contact areas of the knee implants is highly dependent on the respective geometries between different implant designs.It includes the material, thickness of tibia components, and conformity of the mating surfaces which is time-dependent in the motion of gait cycles.This work highlights the importance of using real geometry to study the lubrication mechanism of knee implants.
The model presented in this work can be extended as part of a method to calculate wear; the pressure solution obtained from this model can be used to predict wear.This can be incorporated into an iterative optimisation cycle in which the geometry or choice of material can be optimised for minimum wear of the implant for a given gait cycle.

| CONCLUSION
A mixed lubrication model of the human knee implant is presented in this work, which simulates the complex geometry of the implant.Investigation of the presented results revealed that the lubrication regime varies between elastohydrodynamic, mixed and boundary lubrication during the simulation of the gait cycle.It is necessary to model these lubrication regimes together in order to capture this transition accurately.Furthermore, the lubrication regime present in one condyle does not necessarily mirror the regime present in the other condyle.
In addition, this model addresses the following limitations of traditional elliptical contact solvers: a) Multiple points of contact per condyle are found at various points in the gait cycle, which cannot be found with elliptical contact models.An overlap in non-zero pressure regions is also seen if the points of contact are sufficiently close to each other.b) The ratio of load borne by each condyle is found to vary due to the conditions in the gait cycle as well as the geometry, and therefore a pre-defined ratio of load distribution between the condyles cannot be used, as is often used in elliptical contact models.c) An additional method of solving the lubrication equations is presented, where the forces in all directions can be prescribed and the displacement of the components, along with the pressure and film thickness, are solved for in the model.Due to the limitations of point contact, non-normal loads cannot be balanced, so displacements must be prescribed in these directions instead to achieve the solution for pressure and film thickness.
Consequently, simulating the complex geometry of the implants, along with mixed lubrication, results in a more accurate representation of the tribological performance of knee implants at physiological conditions.This can lead to a better optimisation of the geometry in order to reduce wear and increase product life.However, the computational intensity of simulating the complex geometry of the knee implant is significantly higher compared to the conventional elliptical contact models unless approximations to the deformation matrix are made.

F I G U R E 1
The initial three-dimensional placement of the implant components.Dimensions are in millimetres F I G U R E 2 Side-on views of the knee implant components showing (a) Δx and θ x (b) Δz: This is shown in a zoomed-in view of (a) component to the opposing surface of the femoral component (h ab ).

F I G U R E 3
Definition of h ab as used in this work.A, B, C and D represent the vertices of a discretised planar region on a finite difference grid the 1 GPa value over the other lower values such as 660 or 700 MPa, as long as the condition of the tibia component's modulus being significantly lower than that of the femoral component was satisfied.

F I G U R E 4
Axial force (the total applied load in the z direction) and flexion angle (the rotation of the components about the y axis) derived from the ISO gait cycle F I G U R E 5 Displacement in the x direction throughout the gait cycle BUTT ET AL.

T A B L E 1 F I G U R E 6 F I G U R E 7
Configurations used for mesh dependence analysis Error in results normalised against configuration D Difference in results normalised against the results for 200 timesteps BUTT ET AL.

8
Results of simulating the ISO gait cycle for each condyle for η ¼ 0:1 Pas; σ asp ¼ σ r ¼ 1 μm, β ¼ 50 mm −2 , μ asp ¼ μ r ¼ 0: (a) maximum pressure, (b) minimum film thickness, (c) load in z direction, (d) ratio of the total load borne by the solid asperities Cross-sectional slice of the pressure profile across the point of maximum pressure at 2% of the ISO gait cycle for the medial condyle in (a) x-direction and (b) y-direction F I G U R E 9 Pressure distribution (in MPa) mapped onto the tibial insert at (a) 2% of the ISO gait cycle, (b) 48% of the ISO gait cycle, (c) 84% of the ISO gait cycle viscosity model can be implemented to investigate the effect of a non-Newtonian fluid when simulating complex geometries.

F I G U R E 1 2 1
Film thickness profile at 2% of the gait cycle.(a) 3-D contour plot of film thickness.Due to the scaling of the film thickness, only the locations having less than 1.5 times the global minimum film thickness are shown.(b) Cross-sectional film thickness across y at the point of global minimum film thickness Cross-sectional slice of the pressure profile across the point of global maximum pressure at 48% of the ISO gait cycle for each condyle: (a) medial; (b) lateral model in the Hertzian contact theory, which highlighted the importance of using the FEM in the calculation of the deformation.

3
Non-normal loads for the ISO gait cycle at η ¼ 0:1 Pas: (a) load in x-direction, (b) load in y-direction BUTT ET AL.