On the Effect of Friction on Tibiofemoral Joint Kinematics

: The effect of friction on nonlinear dynamics and vibration of total knee arthroplasties is yet to be investigated and understood. This research work aims at studying the influence of friction on nonlinear dynamics, friction ‐ induced vibration, and damage of tibiofemoral joints. For this pur ‐ pose, a spatial dynamic knee model is developed using an asymmetric nonlinear elastic model ac ‐ counting for knee joint ligaments and a penalty contact model to compute normal contact stresses in the joint while contact detection is treated such that the associated computational time is reduced. Several friction models are considered and embedded in the dynamic model to estimate tangential friction forces in the knee joint. External loads and moments, due to the presence of all soft tissues, e.g., muscles and hip ‐ joint reaction forces, applied to the femoral bone are determined using a mus ‐ culoskeletal approach. In the post ‐ processing stage, damage, i.e., wear and creep, are estimated us ‐ ing three wear models and an empirical creep formulation, respectively. In addition, a FFT analysis is performed to evaluate likely friction ‐ induced vibration of tibiofemoral joints. Mesh density anal ‐ ysis is performed and the methodology is assessed against outcomes available in the literature. It can be concluded that friction influences not only the tribology, but also dynamics of the knee joint, and friction ‐ induced vibration is likely to take place when the friction coefficient increases.


Introduction
Friction acts as a resistance to relative motion, which can, for example, help human beings walk and create desired sounds from instruments such as the violin [1]. However, friction can lead to energy dissipation and be harmful to machine elements and biomedical implants due to, for example, material loss (wear) and degradation, affecting their performance and lifetime [2]. Implant-bearing wear due to friction also is believed to play a notable role in the failure of artificial human joints [3]. Moreover, friction can contribute to aseptic loosening of implants due to wear, shear forces, and high frictional torque, leading to bone fracture, instability and falls [4,5]. Friction also influences the relative motion of knee components significantly, which has several important implications. Excessive anterior-posterior displacements or internal-external rotations lead to the contact occurrence at the edges of the plastic body [6], which can be especially damaging to the knee prosthesis, as is reported in retrieval studies [7,8]. Inadequate sliding motion can also impair human body function by increasing the likelihood of flexion contracture, reducing the range of flexion-extension and internal-external rotations, and the over-tensing of the soft tissues [9,10]. Furthermore, artificial human joints can undergo three-dimensional vibration due to friction, leading to an oscillating movement on top of the gross motion of joint components [11]. Such an oscillating phenomenon can influence both the sliding distance and polyethylene wear, which can eventually cause prostheses to wear excessively, consistent with clinical reports. Friction-induced vibration can also lead to squeaking in biomedical joints [11]. In addition, when ultra-high molecular weight polyethylene (UHMWPE) slides against a metallic counterface, changes in the orientation of molecular chains occur. They preferentially become oriented along a specific path, called principal molecular orientation (PMO), owing to friction while resulting in an orientation softening in the perpendicular direction of the PMO. They both contribute to a directional dependency of wear resistance of the polyethylene while depending on the joint trajectory [12][13][14][15]. The relative motion of the knee components can also influence the creep deformation as it is dependent on the applied load and under-force duration, both of which can change due to friction [16][17][18].
Friction is a force developing and acting against the relative motion of solid surfaces and fluid layers that slides against one another. The origin of the word "friction" comes from the Latin verb "fricare" meaning "to rub". It is worth noting that friction is not a fundamental force because it is reducible to underlying phenomena and forces such as atomic force and electromagnetic attraction between particles of in-contact surfaces. In addition, the friction coefficient is not a material property but a means of describing a friction system. In general, the friction coefficient depends on contact stress, relative tangential speed, surface roughness, the material property of contacting pair, surface contamination, temperature, among others [19]. The scientific study of friction was pioneered by Leonardo da Vinci (1452-1519) who developed the basic friction laws in 1493 but did not document them [20]. The first two classical friction laws were rediscovered and recorded in written form by Amontons (1663Amontons ( -1705 in 1699 [21,22]. He established that the friction force is directly proportional to the applied load but independent of the apparent contact area. Later, Leonhard Euler (1707-1783), through several experimental tests, clarified the distinctions between static and kinematic frictions [23]. It is worth noting that the first who used the Greek symbol mu (μ) for the friction coefficient was Euler. The first mathematical relationship to computing friction force was formulated by a French engineer, Charles Augustin Coulomb (1736-1806), stating that friction is independent of the magnitude of relative velocity between sliding bodies [21]. In his formulation, the friction force is set equal to the applied load multiplied by a coefficient that is the so-called Coulomb coefficient, still being used by both academia and industrial engineers. However, such a model cannot address some underlying physics associated with the friction phenomenon.
In the last decades, a great deal of research has been devoted to advancing phenomenological friction models. Such models have aimed at addressing the distinction between the static and kinematic friction, friction phases (stick, stick-slip, and pure slip), velocitydependent property of friction, break-away force (pre-sliding displacement), and frictional lag [24,25]. In addition to the variety of friction models that address the above phenomena, the friction of UHMWPE, of which the tibial insert of total knee arthroplasty is made, demonstrates a dependency on contact pressure. In 2001, Wang et al. performed a series of experimental tests using a ball-and-socket test rig and found that friction decreases with increasing contact stress [26]. Following their work, Saikko in 2006 reported the same finding while performing pin-on-disc friction tests [27]. Both the above research studies found a significant correlation of the friction coefficient to the maximum contact pressure and reported associated empirical friction equations. In 2017, Barcenas-Sanchez et al. used a ball-on-disc test setup to obtain the friction coefficient of total knee arthroplasty during a walking gait cycle. Subsequently, they proposed a regression model to represent experimental data mathematically [28]. In 2018, Garcia-Garcia et al. continued their research work while conducting a series of experiments to develop prediction models of the friction coefficient between the metallic ball and UHMWPE disk lubricated with diluted bovine serum [29]. Regression parameters for a suggested model that fits well in experimental data were determined for the stance and swing phases separately. They later successfully proposed a formulation that represents the experimental data well to address the behavior of the friction coefficient over the transition from the mixed to the full-film lubrication regime during a walking gait cycle [30]. The fitted model relates the friction coefficient to both contact stress and relative velocity. As the frictional events occur at the atomic level and are of multi-physics nature owing to interactions of chemical, electromagnetic, and mechanical processes, some frictional phenomena are yet to be understood and mathematically formulated despite great achievements so far [19,31].
Friction models are to be integrated into algorithms responsible for modeling dynamics of the knee joint. Sathasivam and Walker (1997) investigated the dynamics of total knee replacement under multiple friction coefficients while they used the Coulomb law to compute friction forces [32]. Much of other research work studying the tribology of TKA commonly assumed a fixed friction coefficient, i.e., 0.04, for the whole numerical simulation. Fisher and Dowson reported that the friction coefficient in artificial human joints is in the range of 0.03-0.10 [33]. However, other studies obtained higher ranges of the friction coefficient for artificial human joints, for example, 0.05-0.33 [27], 0.06-0.25 [26], and 0.06-0.17 [28,29]. Focusing on the dynamic modeling of total knee arthroplasty (TKA), some investigations have aimed at developing anatomy-based dynamic models to consider the kinetic and kinematic behavior of TKAs like quasi-static models [34][35][36][37][38][39][40][41], 2D dynamic approaches [42][43][44][45][46][47][48], and sophisticated 3D mathematical dynamic solutions [49][50][51]. One of the very first successful attempts to determine the 3D dynamic solution of the knee joint was carried out by Abdel-Rahman and Hefzy [49]. Their model did not account for the deformation of the articular surfaces and the real geometry of the tibial insert. Later, Caruntu and Hefzy improved that model to include deformable contacts at the articular surfaces [50]. Although much of the available knee joint dynamic models focused on the joint itself rather than whole-body simulation, Piazza and Delp utilized a full-body dynamic model [52]. The drawback of their approach was to employ the rigid-contact theory to simulate the knee joint, incapable of calculating contact stresses. Recently, Askari and Andersen developed a 3D anatomy-based dynamic approach that links an in-detail forward dynamic methodology of the knee articulation with a deformable contact model to a musculoskeletal system. The latter is responsible for determining physiological forces and moments from surrounding tissues, like muscles and hip-joint reaction forces [40].
Given the state of the art in this field, the effects of friction on the TKA's kinematics are yet unclear and remain to be investigated. Therefore, this research work aims at studying friction effects on the joint motion, friction-induced vibration, and the damage of TKAs. Moreover, this research work carries out a comparative study on several phenomenological and empirical friction models integrated into the algorithms accounting for the dynamics of TKA. The computational framework developed in this study also is observable in Figure 1.

An Overview on Friction Models
Friction is generally described as the resistance to motion when two surfaces slide against each other. Friction between solid bodies is an extremely complicated physical phenomenon encompassing elastic and plastic deformations of the articulating bodies, chemical reactions, interactions with wear particles, micro-fractures and cracks, excitation of electrons, and the transfer of particles from one body to another. However, it is possible to represent this complex phenomenon using either simple friction laws with a couple of friction parameters or more sophisticated friction models that commonly require determining a larger number of parameters experimentally. One of the first and simplest friction laws is the Coulomb friction law. Coulomb (1736-1806) determined that the frictional force between two bodies pressed together with normal force ‖ ‖ can be calculated by multiplying the normal force and friction coefficient that can be computed empirically [21]. The relation between the applied load and relative velocity has the following mathematical form: where F is the friction force, the tangential sliding speed, is the Coulomb friction coefficient, and finally is external, tangential force. Using this model to deal with friction in dynamic systems leads to a computational burden due to the discontinuity at zero velocity. The friction force can reach any value between ‖ ‖ and ‖ ‖ due to a sudden alteration in the velocity direction at zero speed. From a physical viewpoint, knowing the external force, imposed on the sliding bodies, can resolve this issue as the friction load is set equal to the externally applied force according to Newton's first law. Although obtaining the external load is experimentally possible at zero velocity, this is not straightforward to determine in numerical simulations. It is worth noting that the Coulomb friction model only needs one input parameter, i.e., , which is why both academic and industrial communities commonly use this model for systems that do not require very precise quantifications. Another issue with the Coulomb friction law is that it does not capture the distinction between static and kinematic friction. It is well-known that the external force required to make bodies in contact move from their rest status is higher than that when they are in relative motion. Therefore, the Coulomb friction law is modified by replacing the Coulomb coefficient in Equation (1) with a static friction coefficient, i.e., ( ), which is sometimes called Coulomb friction with stiction.
The discontinuity issue at zero velocity is still present in this model, making it hard to detect stick and stick-slip phases of friction between bodies in contact. Moreover, the transition from the stick friction phase to the slip one is not smooth and continuous. To resolve these issues, smooth and continuous Coulomb friction models were developed to prevent the computational burden caused by the force discontinuity by introducing a smooth curve around zero velocity. To this end, several continuous functions have been adopted in the literature, like linear, exponential, and trigonometric function types. One such model is considered here that can be obtained by multiplying the hyperbolic tangent function of velocity to the Coulomb friction model, Equation (1) [53]. This modification, which lets the new model be valid for all from a computational viewpoint, is given by in which k is a coefficient that determines how fast the hyperbolic tangent function changes from near 0 to near +1, which can also be represented by ‖ ‖ where depicts a chosen velocity at which one expects the tanh function to take a value close to its maximum. Another model uses a linear frictionvelocity relationship at the proximity of zero velocity to smooth out the discontinuity, which can be expressed as follows [53]: in which is the slope of linear function between zero velocity and 1 ‖ ‖ ⁄ to eliminate the difficulty in determining the friction force at zero sliding speed. Obviously, the friction force at zero relative velocity between the interacting bodies is not zero and is equal to the tangential applied force. Hence, the main drawback of the modified Coulomb friction laws explained above is that they assume zero friction at zero relative velocity, which does not make sense from a physics viewpoint.
With the study of lubrication by Rayleigh and Stokes, and the classical paper on theoretical modeling of fluid-film lubrication, which Osborne Reynolds published [54], the following expression for the viscous friction can be written (5) in which k  is the viscous coefficient. In full-film contacts, the viscous model can be assumed a good representation of behavior of the system. Moreover, the model can act well enough in other contact conditions provided that the viscous coefficient is tuned. Sometimes a better fit to experimental data can be achieved by assuming a nonlinear term of velocity in Equation (5). As such, Equation (5) can be recast as follows [24]: where the power depends on the geometry of bodies in relative motion. The viscous friction term can be superimposed on the other friction models such as the Coulomb friction [53], which results in, for example, the following formulation.
The Coulomb friction law, and its modified versions, do not explain the negative damping characteristic of friction [11], which was observed by Stribeck as a continuous nonlinear velocity-dependent phenomenon called the Stribeck effect [55].This dependency of friction on the relative sliding velocity was also observed by Panovko and Gubanova [56] and Ibrahim [57] experimentally. The general form of the Stribeck friction model can be written as follows: where ‖ ‖ is a friction coefficient function for which multiple formulations are available in the literature [58] and one of the common forms of which is given below.
Here, is the Stribeck velocity, and is an exponent that depends on the geometry of the application. It is worth noting that depicts both the Coulomb friction coefficient and kinetic friction coefficient in this paper while being recognizable according to the context. Friction parameters in Equation (9) can be obtained by performing the friction test for motions with constant velocity. Although the Stribeck friction model is velocitydependent, this model can be considered a steady friction model. Moreover, the reason for decreasing friction with increasing speed in dry sliding metallic bodies was experimentally investigated and was found to be due to the material softening as a result of high temperatures generated at the contact surfaces [59,60]. Moreover, the friction in a lubricated contact decreases with increasing the relative velocity due to a decrease in the number of surface asperities being in contact until full-film lubrication is obtained. Afterward, the friction can gain a value, and either be constant, increase, or decrease with increasing the sliding speed due to viscous and thermal effects. In the case of an increase in friction when full-film lubrication is built, one can add the viscous friction presented in Equation (5) to the relationship Equation (9) to mimic the intuitive physics. The Stribeck model can provide a good representation of the friction between sliding surfaces. The discontinuity problem of the Coulomb friction model is a problem also observed in the Stribeck model. To cope with this issue, the same procedures taken for the Coulomb friction, like using a hyperbolic tangent function, are applicable here as well. Bengisu and Akay suggested a formulation, Equation (10), which is quite different from that in Equation (9), while resolving the aforementioned discontinuity issue [61].
The first part of the above friction coefficient function uses a near-zero continuous curve to avoid divergence of the numerical model. The second term is the Stribeck friction relation. In addition, 0   is the negative slope of the sliding state [62]. After the friction coefficient starts from zero, it increases to peak friction, which Bengisu and Akay [61] referred to as the static friction coefficient, . The friction coefficient then reduces with increasing tangential velocity until the friction finally reaches a steady state, i.e., .
As was discussed previously, the aforementioned models are not efficient to determine friction magnitude when the velocity is zero nor to manage the stiction phenomenon. Karnopp introduced a concept of the dead-velocity band of the zero proximity, i.e., ‖ ‖ , to cope with the difficulties encountered to detect zero velocity and prevent switching between multi-state relationships to transit between stick, stick-slip and pure slip phases of the friction [63]. Within the abovementioned velocity interval, the speed is assumed to be null. The friction phase also is stiction, where the friction force can be equal to either tangential external force or static friction. For the velocities outside of that dead band, the friction force can get the form of Coulomb friction, Stibeck friction, and so forth. The following equation can express this idea.
However, the concept Karnopp used to introduce this friction model does not comply with friction physics. Moreover, this model requires the external load as an input, which is often not available explicitly and is an internal state of the dynamical system.
In addition to the phenomenological friction models, there are empirical friction models proposed for UHMWPE in the literature. Wang et al. [26] reported the following formula derived from experimental data, showing a significant correlation to the contact stress. However, the Coulomb friction model and its modified versions cannot address changes in friction coefficient due to contact stress: in which stands for the friction coefficient and is the contact stress with a unit of MPa. Furthermore, another empirical formulation was suggested by Saikko to relate the friction coefficient to the contact stress as follows [27]: where is a reference pressure (1.1 MPa). The above relationship is correct for 1.1 MPa and below this pressure, the coefficient of friction increases with increasing the contact pressure. Recently, several empirical friction models for tibiofemoral joints were developed by a research group in Mexico [28,29]. Montes-Seguedo et al. [30] proposed a power-law model from which the friction coefficient can be estimated in both mixed lubrication regime and full-film lubrication regime as a function of the maximum pressure ( ), entrainment velocity ( ) in mm/s, and sliding-to-rolling ratio (SRR), which is written as follows.
Models discussed above are not able to capture dynamic friction characteristics such as pre-sliding and frictional lag. Dynamic friction models such as the Dahl approach and LuGre model can address such phenomena. The discussion associated with the dynamic friction models stands out of the scope of the present study. Interested readers are referred to the following references [24,[64][65][66][67] for further information

Contact Solver: Tangential Friction and Normal Contact Forces
As was discussed before, when two surfaces enter into a contact phase and tend to slide against each other, friction develops and acts as a resistance to the relative motion. In this study, the friction forces are taken into account employing several friction models, listed in Table 1 where , is the tangential relative velocity of the node , in contact as the tibial insert is assumed stationary, and , , , is the friction coefficient at the contacting node, which can be considered as a function of either relative velocity or contact pressure and can be substituted from each of the modified Coulomb friction and Stribeck friction models, outlined in Table 1. Moreover, static and dynamic friction coefficients can be dependent upon the contact stress based on empirical models reported by Wang et al. [26] and Saikko [27]. Eventually, seven case scenarios are defined in Table 1, each of which is considered in this study to study kinematics and tribology of TKAs. In Table 1, is a constant value that takes a value greater than one. Moreover, the pressure-dependent friction coefficient is used in the models A2 and B2, illustrated in Figure 2, to show the difference between such models. As can be seen, the friction coefficient is a function of both relative velocity and pressure when model B2 is employed, while it varies merely with contact pressure in model A2. It is worth adding that both models use the pressure-dependent friction relationship Wang et al. proposed in [26]. However, such a relationship is integrated into the Coulomb friction formula for model A2 and the Stribeck one for model B2. It is worth noting that the empirical model Wang et al. [26] suggested does not address the friction coefficient once the contact stress approaches zero. Therefore, the strategy considered here is to prevent the friction coefficient to exceed the value 0.25 that corresponds to the contact stress of 0.37 MPa. In addition, the models A3 and B3 are written based on the relationship recommended by Saikko. His empirical mathematical formulation is valid for the contact stress equal or greater than 1.1 MPa. According to the plot that shows the relationship between the contact stress and friction coefficient in his paper, it can be concluded that the friction coefficient decreases once the contact stress diminishes below 1.1 MPa until it gains an average value of 0.25 once the contact stress reaches zero value. Therefore, a quadratic function is fitted to the data while outliers are removed to address the value of friction once the contact stress gains a value less than 1.1 MPa. The model C1 is the one developed by Montes-Seguedo et al. [30], in which is the velocity of the femoral part at the contact point, which is calculated according to the radius of the functional flexion-extension axis, multiplied by the respective angular velocity and , the speed of the tibial insert at the tibial plateau surface. The latter is computed as the summation of the anterior-posterior linear velocity and the multiplication of the radial distance from the tibial axis to the bottom of the tibial plateau and the internalexternal angular velocity. For further information on how to determine these velocity terms, the reader is referred to the reference [28]. At this stage, a short description of the joint kinematic is presented to determine the velocity of master nodes, including their tangential speed. The relative velocity of a master node like , with respect to its corresponding slave node , locating on body 2, i.e., the tibial insert, can be determined as follows: , , , where , is the velocity vector of the point , with respect to , , and depict the translational velocity of the local coordinate system and the angular velocity of body 1, respectively, , depicts the coordinate vector of the point , , and is the location of the origin of the local coordinate of body 1, i.e., the femoral body, at time tn.
, is the vector normal to the femoral body at the point , at time tn. and , and , also are the normal and tangential relative velocities of the point , . As body 2 is assumed stationary, the components of the angular velocity vector resolved in the fixed basis could be evaluated based on Euler angles as follows [68]: (19) where 0 cos sin cos 0 sin cos cos 1 0 sin (20) in which , , are Euler angles and , , depict their time derivatives. Eventually, the friction force at each element on the femoral components is evaluated knowing the tangential sliding velocities, Equation (18) at all element nodes, the normal load applied to the element, and contact pressure, which are addressed later in this paper. Finally, the resultant friction force tangentially applied to the femoral component can be obtained from the following integration that is performed over the contact area: , , where Ω depicts the contact area, , is friction force at each element and , depicts the size of element area. The present study implements a concept proposed in [69] to evaluate contact stresses, but it is impossible to derive a closed-form mathematical formulation for the contact problem at hand. Askari obtained nonlinear and promising contact models while resolving the discontinuity issue with the KelvinVoigt model at both the beginning and end of the contact process [69]. To this end, he simulated the contact by using a set of independent springs and dashpots, which represent the stiffness and viscous damping of the contact between articulating bodies, respectively. Readers interested are referred to the reference [69][70][71][72] and references therein for further information about the model employed in this study. Such a contact model does not put any kinematic constraints on the dynamical system. It is well-known that the KelvinVoigt model suffers from a discontinuity at both the beginning and end of the contact process due to the form of the damping term in its formulation [73,74]. The concept proposed in the reference [69] showed that simulating a soft contact using a set of independent KelvinVoigt elements spread all over the contact surfaces instead of just one spring-damper element as is used in the KelvinVoigt contact model can resolve the above discontinuity. The penetration depth at each node on the polyethylene surface, , , is known from the minimum distance technique discussed later in the present paper (see Section 5.3). Its time derivative , is also obtained from the normal velocity term given in Equation (17). Therefore, the normal contact pressure, i.e., , , at any node of , can be determined based on the following KelvinVoigt contact relationship: in which is the contact stiffness and the relaxation time, is obtained from Λ/Ξ where Λ 1 / 1 1 2 and Ξ is the tibia thickness [40], E and are Young's modulus and Poisson's ratio of the polyethylene that are, respectively, 463 MPa and 0.46 as are reported in [75][76][77]. Such a contact model belongs to the category of the penalty non-adhesive contact approaches and penalizes the indentation of the master body into the slave one. As independent KelvinVoigt elements are spread all over the contact area, the normal contact force applied to the femoral components should be computed by performing the following integration over the contact area: (23) where is the normal unit vector of an element with the area on the femur-bearing surface. It is worth noting that the following criterion is used to assess whether a node is in contact or not.
The resultant contact area can also be obtained according to the summation of the area associated with all elements engaged in contact area having nodes with non-zero contact pressure.

Patient-Specific Musculoskeletal Modeling
To estimate patient-specific TKA knee joint loads, an already validated, patient-specific musculoskeletal model is applied [35]. This model is based on data for one male subject (age: 86, height: 1.80 m and mass: 75 kg) with an instrumented, posterior cruciateretaining TKA prosthesis from the 5th Grand Challenge Competition to Predict In Vivo Knee Loads [78][79][80]. For more details of the implant shape and design, the reader is referred to the Grand Challenge knee dataset [78]. Among others, the dataset contains preand post-operative Computed Tomography (CT) scans, trajectories of skin markers and ground reactions during standing reference trials and dynamic movement trials, including gait at a self-selected speed as well as measurements from the instrumented knee prosthesis. Additionally, Stereolithography (STL) 3D geometries of the femoral component, tibial tray and insert, patellar button, and segmentation of the post-operative CTs of the partial pelvis, femur, patella, tibia, fibula, partial talus and partial calcaneus are also included in the dataset [81]. For the dynamic analysis, the one gait trial of level walking at self-selected speed (PS_ngait_og_ss1) is utilized [78]. The musculoskeletal model is developed using the AnyBody Modeling System (AMS) v. 7.1 (AnyBody Technology A/S, Aalborg, Denmark) and based on the human model from the AnyBody Managed Model Repository (AMMR) v. 1.6. From this model, we extract all the knee flexion-extension angles as well as all loads associated with muscle, hip joint, gravitational, inertial, etc., which are imposed on the femur. The tibiofemoral joint reaction forces and knee ligament loads are included in the joint model described below and, therefore, are not part of that model's boundary conditions.

Dynamics Solver
In this section, a forward dynamic model of the tibiofemoral joint is presented. We focus on an assembly consisting of the tibia and femur as the tibiofemoral joint's components along with its ligaments to develop the forward dynamic model to which damage prediction approaches are later added. In the following, the development process of such a model illustrated in Figure 1 is described in detail. A forward dynamics model is first derived, which requires the evaluation of contact stresses and ligament forces and moments before being solved. To determine contact stresses, triangulated surfaces of bearing components are smoothed and contact detection is performed.

Dynamics Modeling
The governing equations of the system motion are derived based on the free body diagram of the tibiofemoral joint represented in the second right-hand panel of Figure 1. The resulting forces and moments applied to the femoral bone due to the existence of all surrounding soft tissues such as muscles and hip-joint reactions are computed by means of musculoskeletal modeling. Such loadings and moments are moved to the origin of the femur coordinate system and written in the global coordinate system, illustrated in the second right-hand panel of Figure 1, as Fext and Mext. The loads and moments owing to the contact, friction, and ligaments are also determined and summed up to those external ones. The governing equations of the translational motion of the femoral component, therefore, can be formulated using Newton's second law of motion as follows [69,82]: , (25) in which the mass of the femoral bone is designated by M. The translational acceleration vector of the femur is also depicted by . Furthermore, stands for the vectorial summation of forces imposed on the bone, which is mathematically formulated in Equation (25). In this equation, the external forces are depoicted by , resultant contact forces , friction force , and finally ligaments loads are designated by . In this solution setup, three coordinate systems are defined, one of which is the global coordinate system (XYZ) attached to the underneath of the tibial insert and on the top of the tibial tray. It is worth mentioning that the tibia is considered stationary, while the femur has six degrees of freedom to translate and rotate. One degree of freedom of this dynamic system is, later, constrained while knowing the flexion-extension angle of the knee joint over a level walking gait. The other coordinate system is a local system (xyz) attached to the STL center of the femoral component, rotating and translating with the femur. However, the third coordinate system, , is considered with the same origin as the femoral coordinate system, although it does not rotate with the femur and it has bases aligned with those of the global coordinate system. The angular-momentum equation can now be written in the body coordinate system of the femur, which governs the rotational motion of the femur as follows [69]: * * * * * * where * is the summation of not only external moments but also those resulted from the contact, friction and ligaments, which are computed with respect to the body coordinate system of the femoral part (xyz). Compared to the computation of moments in the coordinate system xyz, it is much easier to obtain moments in the third coordinate system discussed above. The summation of moments in the third coordinate system is designated by and is related to * by the following relationship: * in which is the transpose of the rotation tensor of the femur, * and * appeared in Equation (26) and stand for the tensor of inertia and the angular velocity vector of the femoral component in the body coordinate system, respectively. In the dynamic model developed in this study, the rotation tensor is defined based on the Euler angles of the zx-y sequence, given by the array of Euler angles , , to characterize the rotation of the femur as well as the location of nodes on the femoral-bearing surface. The respective rotation tensor can be cast by: ≡ cos cos sin sin sin sin cos cos sin sin sin cos sin cos cos sin sin cos cos sin sin cos sin cos cos sin sin cos cos (28) where the flexion-extension angle, which is the amount of the planar rotation of the femur around the axis z, is named by . From the time derivatives of Euler angles, the angular velocity vector in the body coordinate system can be given by: * * where * cos sin cos 0 sin 0 1 cos cos sin 0 The tangent operator * is tangent to the rotation manifold. There is a possibility to simplify Equation (26) while setting the body coordinate system such that it coincides with the principal axes of the mass moment of inertia tensor, i.e., * . Hence, the moment of inertia tensor is reduced to a diagonal form, and the governing equations turns to the so-called Euler's equations for the angular motion of the dynamic system as follows. * * * * * * * * * * * * * * * * * * * * * Characterizing of the femoral bone along with surrounding tissues of the male subject, the femoral mass, i.e., M, is 7.5 kg and its moment of inertia tensor is obtained as (0.4516 0 0; 0 0.0213 0; 0 0 0.4516) [35,78]. The adaptive RungeKuttaFehlberg approach is employed to integrate Equations (25) and (31) over time of interest [83]. A maximum error is set and the error magnitude at each time step of the numerical simulation is assessed while comparing outcomes from the above method of different orders 4 and 5. Once the error acquired is greater than the set maximum error, the time step is halved and the computational process is done again until the accuracy and stability of outcomes become guaranteed. The present study also considers a minimum value of the integration time-step size 5 × 10−5s along with an integration tolerance 5 × 10−6 [84].

Ligament Modeling: Ligament Forces and Moments
In this study, ligament forces are modeled using an asymmetric nonlinear elastic model. The force, , each ligament, i.e., PCL, LCL, and MCL, imposed on the femur can be evaluated from the following formula, which also includes a slack region [85].
Here, is the ligament stiffness and depicts the strain while is a constant with value of 0.03, which is associated with the transition phase between linear and nonlinear regions of the force-strain curve. In order to determine the ligament strain, the slack length, l0, is obtained as below: where is the reference strain and lr stands for the bundle length estimated when the leg is fully extended. Finally, the strain, , at any time is computed from the following equation: (34) in which l is the instantaneous bundle length of the ligament, which is determined according to the motion of the femur over the walking activity. Wrapping surfaces are also employed to prevent ligaments from penetrating the bones and the implant. Stiffness and reference strain assigned to each ligament bundle are extracted from [35,86], which are adapted from the literature [85]. To compute ligament forces, the insertion points of every single ligament are assigned before the numerical simulation starts. Since attachment sites could not be determined from the dataset, they are estimated according to descriptions found in the literature [35,[87][88]. The motion of the femoral component at each time step is acquired using motion equations, given by Equations (25) and (31), and the rotation tensor, Equation (28). Therefore, we can now determine the current locations of insertion points in the global coordinate system. Moreover, each ligament is discretized into several elements to calculate their time-dependent length according to the femoral motion. The magnitude of each force is, in turn, obtained while solving Equations (32)- (34) and the corresponding vector direction is considered to be along the ligament tangent at its insertion point. The resulting ligament force and moment imposed on the femur are computed by ,

Contact Detection
STL 3D geometry files of the knee components are available from the Grand Challenge knee dataset [78]. After removing facets and respective vertices that are not likely to get into contact, the triangulated surfaces of the contacting pair are smoothed using a Laplacian smoothing technique. As the tibial insert is assumed stationary, its bearing surface is constructed using the Non-Uniform Rational B-Spline (NURBS) surface methodology [40,[89][90] and, in turn, is meshed using a uniform mesh size in the first place. A bounding box scheme is applied to the meshed surface of the tibia to reduce the computational time of searching spatial contact. All boxes, susceptible to the contact incident while deemed based on a one-by-one surface function-based scheme, are retained. In the next step, each box becomes open and each node inside is assessed for whether it is in contact with the femoral component or not using the minimum distance technique described as follows. The minimum distance of a node such as , on the slave-bearing surface with respect to the femoral body can be obtained as is described as follows. A facet with number i, belonging to the triangulated surface of the body 1, is considered and a vector is drawn from the first vertex of the facet i, i.e., , to the node of interest (see Figure 3). This vector is, in turn, projected on the vector normal to the facet under consideration, which is designated by , ⃗ and its magnitude is equal to its ℓ norm , ⃗ . Then, the projection of the node , on the plane of the facet i is obtained and called . . By this time, the question is whether that projected point, i.e., . , places inside the domain of the facet i. To check it out, the summation of angles between each two of the three following vectors, , ⃗ , , ⃗ , and , ⃗ , as illustrated in Figure 3, is acquired. The ℓ norm of each of those three vectors is also computed. If either of the following cases takes place, it is confirmed that the point , corresponds to the likely contact point on the master body and the minimum distance of the node , with respect to the femoral surface is then acquired, i.e., , , ⃗ , and stored as the penetration depth , in the penetration matrix D: (i) the summation of the above three angles is equal to 2 ; (ii) the ℓ norm of one of those vectors is zero, so , and the respective vertex are coincident.

Damage Formulation
Now that the dynamic model is ready to be used, the tribology of TKA is formulated to predict wear. Archard wear model is commonly used by the tribology community to describe adhesive and abrasive wear mechanisms, although it is often adopted for a wide range of applications as a result of its efficiency and simplicity [36]. Employing Archard wear law, the linear wear rate can be computed using the following expression [91,92]: ℎ (36) in which the wear depth and sliding distance are presented by h and s, respectively. The parameter kW stands for the wear factor that has the following unit mm 3 N −1 m −1 . Although the wear factor has been considered constant by several research studies, Equation (37), Kang et al. observed from experimental data that the wear factor is dependent on the contact pressure and cross-shear ratio. Therefore, they suggested a new formulation for the wear factor that varies with not only the cross-shear ratio, ℑ, but also the contact pressure. They proposed a relationship to compute the wear factor of the UHMWPE tibial insert [93] accounting for such dependencies given in Equation (38) in which the average contact stress, i.e., , for a given element is computed by averaging contact pressure over one gait cycle, and the definition of the cross-shear ratio and the procedure to compute it is welldescribed in [93].The described wear techniques listed in Table 2 are used in this study to estimate wear in total knee arthroplasties Wear factor given in Equation (38) is dependent on the contact pressure and is not easily implemented into computational wear modeling. Abdelgaied et al. [77] developed a wear model based on the idea that volumetric wear (W) is proportional to the contact area (A) and sliding distance (S) and a non-dimensional wear coefficient (C) determined experimentally to complete the wear formulation, Equation (39). The cross-shear ratio is CS, and parameters a, b and c are constant and determined from the experimental measurements of a multi-directional pin-on-disk wear test [77] as a = 8.5173e-65, b = 9.3652e-60, and c = −6.7454. .

Cross-Shear Ratio
Under the sliding motion of a metallic part on ultra-high molecular weight polyethylene, the polymeric chains align in the direction of the maximum frictional work according to the experimental findings [15,94], which is the so-called principal molecular orientation (PMO) [12,15]. The cross-shear ratio can be quantified as the quotient of the frictional work done perpendicular to the PMO direction ( ) to the total frictional work ( ) as is formulated below [95][96][97].
The frictional work can be calculated based on is the incremental sliding distance that can be computed by ∆ , ∆ at each element. The direction ( ) of incremental motion along the contact point trajectory is calculated to determine the vector of the incremental sliding distance, ∆ . Moreover, considering a test PMO direction (axis) with an angle , frictional force and incremental sliding distance components parallel to and perpendicular to that axis can be calculated for each increment [12,[95][96]. Assuming that the friction coefficient does not vary during a cycle, the cross-shear can be formulated as follows.
The dynamic model developed in the previous sections obtains the sliding distance and contact force with time, which allows calculating Equation (41) for each angle , from 0 to π, until the principal molecular orientation, , is reached where the cross-shear ratio is minimized.

Creep Estimation
In addition to the wear phenomenon, creep deformation contributes to the surface damage of the UHMWPE tibial insert due to the viscoelasticity of UHMWPE, where deformation under a constant load varies with time. Performing experimental creep tests, Lee and Pienkowski [16] presented a nice formulation to compute the creep deformation, which is a function of both pressure and a logarithmic timescale as follows [92]: where , is the linear damage at an element (k, j) due to creep, N is the total number of cycles in service, Γ , is the thickness of element (k, j) of the acetabular cup and ∆ , associated with an element (k, j) of the cup surface is non-zero just when the corresponding i is in the set of time values where the contact stress is non-zero. The unit of time and pressure are minute and MPa, respectively, according to Lee and Pienkowski experiment.

Convergence and Mesh Analysis
To guarantee the accuracy and convergence of dynamics and contact simulations, a mesh density analysis is performed [98]. Three different element sizes to discretize the domain of the bearing components projected on the XZ plane are considered, i.e., 0.2, 0.4, and 0.6 mm. Multiple system parameters that can be influenced by the mesh size are used to assess the mesh density, such as contact forces and moments on both medial and lateral condyles, and maximum contact pressures on both condyles. According to the outcomes we obtained at different percentages of the gait cycle, it has been concluded that mesh sizes 0.2 and 0.4 produce accurate results. The latter is used in this study to perform the computational analysis. The respective outcomes corresponding to two pressure peaks in a gait walking cycle, for example, are listed in Table 3. Moreover, a block size of 4 × 3 is considered and a one-by-one surface function scheme is incorporated for contact detection in order to accelerate the computational modeling. Finally, the total computational time consumed to solve the motion equations and to acquire wear values for the present method is less than 4 hours while the dynamic model developed in this study has been implemented in MATLAB (R2017a) and run on a 2.7 GHz personal computer with Intel ® Core(TM) i7-6820HQ CPU.

Dynamics and Tribology of the Joint
The developed model is employed to compute linear and volumetric wear rates for one million cycles while using three available wear techniques that are listed in Table 2.
The outcomes that are obtained for all seven friction models listed in Table 1 are reported  in Table 5. The highest values of both linear and volumetric wear rates are reported employing Method i with a constant wear factor of 2.2 × 10 −16 mm 3 N −1 m [36]. Using the wear factor suggested by Kang et al. [93], that is called Method ii, gives less wear depth on the medial condyle than that on the lateral one, which opposes the outcomes associated with the other two methodologies. The wear values resulted from Method iii, associated with the concept presented by Abdelgaied et al. [77] in which wear coefficient is considered instead of wear factor, stand between those from Method i and Method ii. As can be seen from Table 5, the highest wear values are obtained using the friction models A1 and C1, while the lowest amounts are related to those acquired from A3. Considerable variations in obtained wear rates are observed while using these three models. Hence, it can be concluded that care should be taken to estimate wear in terms of the wear model that one employs. The wear map and distribution obtained from the following four friction models, i.e., A1, A3, B1, and C1, considered in this article are illustrated in Figure 4. It can be observed that the maximum wear on the lateral condyle is shifted from the tibial plateau to the lateral periphery of the joint when using model B1. Moreover, the area undergoing high material losses on the medial condyle splits into two sub-areas in Figure 4b while lower linear wear rates are observed when using the model B1 compared to the others. The maximum medial condyle wear rates are obtained once model A1 is utilized as can be observed from the intense red color of the wear map. The loci of the motion of the femur center with respect to the tibia on the transverse plane for the following five friction models, i.e., A1, B1, B2, B3, and C1, are illustrated in Figure 5. Static and dynamic friction coefficients are 0.1 and 0.085 while and are equal to 1 mm/s and 10, respectively. To estimate such friction coefficients for the sake of computational analysis, the average magnitude of friction coefficient from empirical models, i.e., Saikko model and Wang et al.'s, are obtained over the range of contact pressure the knee joint undergoes over a walking gait cycle. The trajectories show the same trends to some extent, although some variations are observed; for example, the one corresponding to model A1 is more influenced by friction than that acquired from model C1 as can be observed in the left-hand side plot of Figure 5. The external-internal rotation angle is also plotted against the abduction-adduction angle over a gait cycle in Figure 5. The models A1 and B1 lead to the very similar shapes of the rotational path, while the same can be observed in the case of models B2 and C1. The ranges of variation of both angles when using model B3 decrease. Phase portrait diagrams of the dynamic motion of the tibiofemoral joint are also presented in Figure 6 for the friction models A1, B1, and C1. It can be seen that the range of velocity increases once model C1 is utilized, while model A1 produces the highest acceleration magnitudes and the lowest maximum velocities among others. The norm of position vector varies in the range of 44.26-46.30 mm. Moreover, the highest speed occurs during the swing phase at around 70% of the walking gait cycle in Figure 6a,c,e. Phase portrait diagrams can help one study the type of dynamic response observed in the knee joint while employing different friction models. Figure 6b,d,f show chaotic behaviors of the tibiofemoral motion, which can be considered sensitive to the used friction model, and densely filled in the phase portrait diagram [103]. The chaotic response suggests that the vibration can occur due to friction, the so-called friction-induced vibration, and impact/contact between the femoral part and the tibial insert. This finding is in agreement with the acceleration level found for the frequency at which the joint is likely to vibrate, as is shown in Figure 7, where the highest acceleration is obtained once the model A1 is used. A FFT frequency analysis is performed to characterize the friction-induced vibration of the tibiofemoral joint. It is found that the tibiofemoral joint does not vibrate when the friction is below a critical value as the joint does not oscillate using Coulomb or Stribeck friction models with friction coefficient 0.05. It has also tested that the frictionless joint does not undergo oscillation. However, increasing the friction coefficient from 0.05 to 0.1, the joint vibrates as can be observed from Figure 7c

Coulomb and Stribeck Friction Models
In this section, an in-detail comparison between the two well-known friction models, Coulomb and Stribeck friction laws, is conducted. The effect of friction coefficient on linear and volumetric wear rates is investigated while using Coulomb friction law and acquired outcomes are listed in Table 6 using Method iii [77]. It can be observed that the maximum linear wear rates are obtained in the case study with μc = 0.1. The lateral linear wear with friction coefficient 0.2 is lower than the others, which means it is not a direct relationship between the wear values and friction coefficient. The interesting outcome is that the maximum volumetric wear rate occurs for the case study with μc = 0.05 while the minimum volumetric and linear wear rates are associated with the frictionless joint and the one with a friction coefficient of 0.2. When the friction is higher, greater linear wear can be expected as the worn area decreases due to the changes in the trajectory of the femoral part. However, there is another important contributing factor to wear occurrence, which can prevent the wear to increase. When UHMWPE slides against a metallic counterface, molecular chains preferentially become oriented and hardened along the direction of the PMO, resulting in higher wear resistance of polyethylene [14,104]. The effect of parameters in the Stribeck friction model on linear and volumetric wear rates are also studied and presented in Table 7. It can be observed that the higher the magnitude of is, the greater the linear wear rates are. Moreover, increasing the amount of leads to an increase in linear wear rates. The higher friction coefficients become, the higher the linear wear rate on the medial condyle is, while the less the value of friction coefficient, the greater the linear wear rate on the lateral condyle would be. The maximum volumetric wear rate also takes place once the friction coefficients 0.2/0.13 are considered. The effect of friction coefficients in the Coulomb and Stribeck friction models on the loci of the tibiofemoral joint on the transverse plane is also investigated and illustrated in Figure 8. It can be seen that increasing the friction coefficient leads to a significant change in the trajectory of the femur center in the lateral-medial direction in particular. However, the friction influence is smaller in the posterior-anterior direction, although some variations have been observed for the case study with friction coefficient 0.05 and 0.2 for Coulomb friction law and 0.2/0.125 for Stribeck model, in particular, compared to others. This is also observed in that the trajectory loops of all case studies are closed. Although friction and viscous contact lead to energy loss in the system, energy is continuously introduced in the system due to the rheonomically constrained flexion-extension rotation and boundary conditions, e.g., loads and moments imposed on the femoral bone from muscles, hip joints, among others, which is why the femur center motion converges to a steady-state periodic motion. It can be inferred that the trajectory obtained while using the Stribeck friction model is less influenced by friction force due to the difference existing between the dynamic and static friction coefficient and the velocity-dependent characteristic of friction once the Stribeck model is utilized. The damage map and distribution using the Coulomb friction model with four different friction coefficients after five million cycles in service are considered and illustrated in Figure 9. It can be observed that the distribution of damage in both condyles expands in the anterior-posterior direction once the friction coefficient increases. Moreover, increasing friction gives rise to damage values on both condyles as can be observed in Figure  9 and the more intense color of the damage map becoming red. It is worth mentioning that the damaged area on the lateral condyle of the knee joint is seen to increase significantly with increasing the friction coefficient. The other finding is the maximum damage on the lateral side of the joint occurs close to the middle plateau of the tibia when μc = 0.05, while it is shifted to the anterior part of the lateral condyle with increasing friction as is observable in Figure 9b,d. The damaged area on the medial condyle with high damage increases in the case study with μc = 0.2, as seen in Figure 9d. The same study is performed to assess the effects of both static and dynamic friction coefficients used in the Stribeck friction model on the damage of tibiofemoral joint five-year post-operation, illustrated in Figure 10. It can be observed that increasing friction coefficients leads to a decrease in damage magnitude on both condyles, while a greater area of lateral condyle gets engaged in wear occurrence along the anterior-posterior direction, in particular. The highest damage magnitudes also happen in the case of the frictionless joint as can be seen in Figure  10a.

Conclusions and Future Research Directions
This study developed a spatial dynamic methodology to investigate the effect of friction on nonlinear dynamics, vibration, wear, and creep of tibiofemoral joints. A mesh density analysis was performed to end up with a fine mesh size that guarantees the accuracy of obtained outcomes. Comparing acquired results to those available in the literature also made it possible to evaluate the developed model. In conclusion, friction can considerably influence the motion and damage of TKAs. From the phase portrait diagrams, chaotic behaviors were observed in the tibiofemoral movement. It was also inferred that frictioninduced vibration takes place when the friction coefficient increases. Moreover, using an empirical friction model resulted in vibration occurrence. Making any changes in the magnitudes of parameters in the Stribeck friction model led to the alteration of wear rates on both condyles. Furthermore, it cannot be expected that increasing the friction coefficient can always cause wear to rise owing to the strong nonlinearities of the wear mechanism. The three wear models employed in this study to predict the wear of TKAs also produced outcomes with significant differences, which should be taken into account by designers and engineers in their designing process. Eventually, it was illustrated that friction models influenced both translational and rotational loci of tibiofemoral motion, and consequently, the damage distribution.
As future research directions, the proposed model can be extended to account for the dynamics of the patella-femoral joint in addition to the tibiofemoral joint as well as to take fluid lubrication into account while handling fluidstructure interaction [105][106][107]. Moreover, the suggested approach can be utilized to optimize the design and material properties of TKAs to improve their lifetime and tribology performance, using optimization approaches such as genetic algorithm (GA) and particle swarm optimization (PSO) algorithms [108]. The friction models considered in this article are not able to capture dynamic characteristics of the friction phenomenon, e.g., pre-sliding and frictional lag. The dynamic-based friction approaches like Dahl and LuGre are to be employed accounting for these characteristics, which can lead to much more precise outcomes than those obtained in our study. As was illustrated, increasing friction leads to friction-induced vibration that should be considered in more detail experimentally as well. In addition, the viscoelasticity of polyethylene tibial insert can be considered using more sophisticated models, namely the generalized Maxwell model, Wagner model, and Prony series [109].
The present study assumed that the geometry evolution of the polyethylene tibia owing to material loss and creep deformation does not affect the contact stresses and knee trajectory, whereby the wear magnitude varies linearly during the wear analysis. Surface evolution leads to a reduction in contact stresses over time, resulting in decreasing wear rates and creep deformation [36]. Therefore, it can be deduced that the present study overestimated linear wear rates and creep deformation, both of which are dependent on contact pressure. The presented methodology can be extended to include the surface variations of the plastic tibial insert due to the wear and creep phenomena. The improved model will address the effect of the geometry update on both contact stresses and contact point trajectory, and eventually, wear and creep. Furthermore, damage distribution was found to vary with changes in the friction coefficient of knee components. Therefore, a greater area of the polyethylene tibia can get damaged when friction increases. A future direction of this research is to consider how any alteration in the damage map can influence the performance and lifetime of TKAs. The knowledge gained through such a study can contribute to the design and material development of TKAs.