1. Introduction

Femoral acetabular impingement (FAI) is a medical condition in which abnormal bony protrusions on the femoral head or the rim of the acetabulum produce stress concentrations in the hip joint during motion. Damage incurred in the articular cartilage due to the increased stress in impinged areas is thought to be a precursor to cartilage degeneration and osteoarthritis [1, 40]. A preoperative simulation of real patient hip morphology that could be manipulated during the planning of procedures would be of great benefit to surgeons and could also conceivably reduce the necessity for total hip arthroplasties if the FAI is indentified early.

Articular cartilage is an essential component to any such simulation despite a thickness not greater than 4 mm [1, 29]. It is a biphasic soft-solid tissue comprised of an interstitial fluid phase and a solid matrix that interact to distribute and dampen the load carried by a joint. Significant progress has been made towards the accurate simulation of articular cartilage with computational methods, beginning with the biphasic theory of Mow et al. [29]. Mak et al. [25] extended this to the biphasic poroviscoelastic (BPVE) model with fluid flow-dependent and independent parameters to simulate viscous effects. Both the transverse isotropic model [1, 9] and the conewise elasticity model [38] provide means of including directionality in loading response. Soulhat et al. [39] proposed the more rigorous approach to the simulation of cartilage fibrils, which was subsequently extended by Li et al. [21, 23] to the biphasic fibril-reinforced poroelastic model. It is beyond the scope of this paper to provide a thorough review of the field of cartilage simulation, but other significant works include those by DiSilvestro [9], Korhonen [18], Wilson [42], Garcia [13], and Seifzadeh [33], to name only a few.

The finite element method (FEM) is the traditional means of modeling articular cartilage. While FEM is considered to be the gold standard for accuracy in computational methods, it is an unrealistic expectation that with the current technology included in a reasonable desktop PC that FEM could be made to perform interactively in a preoperative simulation considering the many nonlinear properties of articular cartilage. Other methods have been proposed for the simulation of soft tissue, perhaps the predominant of which is known as mass-spring [36, 1, 28]. However, mass-spring is impractical for the intended FAI simulations, since it is a non-volumetric method where spring elements connecting point masses must be tuned for each desired scenario [1, 44].

Smoothed particle hydrodynamics (SPH) is a Lagrangian particle technique that employs a nearest neighbour search to reduce the computational complexity of a simulation. While SPH was originally developed for the simulation of exploding celestial objects [14], considerable work has since been performed with fluid dynamics simulations of the Navier-Stokes equations [30, 6]. SPH has more sparingly been applied to solids, beginning with the work of Desbrun and Cani [11], and extended to bending, fracturing, and high impact tests [15, 7, 5]. Solenthaler et al. [37], provided an approach to derive the deformation field of an elastic solid, which was improved upon by Becker et al. [4] by a corotation method to prevent unrealistic forces in rotations. Porosity in an elastic solid material using SPH was proposed by Lenaerts et al. [20]. The recent work by Akinci et al. provided means to combine the simulation of solid and fluid SPH [1]. The extension of SPH to soft tissue is relatively sparse, with a few examples of fluid SPH confined by mass-spring meshes [31, 32], and a simulation of a virtual liver by Hieber et al. [17] that was limited to small deformations in two dimensions.

The research presented here uses SPH to simulate articular cartilage as a biphasic fibril-reinforced poroviscoelastic material employing a combination of components drawn from previously published research and newly introduced methods. The objectives were to implement such a model and validate it using indentation and unconfined compression tests while considering the necessity of computational efficiency as applicable to a future real-time hip simulation.

2. Materials and Methods

The fundamental SPH integral interpolant is given by Gingold and Monagahan [14] in its generalized form as:

(1)

where is the function of the spatial coordinates , is the smoothing kernel, and is a differential volume element. The integral is advanced in time for each particle using an integration technique, which in the current work is the commonly used leapfrog method. The continuum of a material can be represented by Lagrangian SPH particles with physical equations based upon the above interpolant.

In the current procedure, the behaviour of articular cartilage in response to compression scenarios is considered to be due to the interactions of the stresses in each of the components:

(2)

where is the total stress from all components, while , and are respectively the stress contributions from the non-fibrillar solid matrix, the collagen fibril network, and interstitial fluid support.

2.1. Non-fibrillar Matrix

The model used to represent the non-fibrillar matrix component of cartilage is based upon the works by Solenthaler [37] and Becker [4]. This model was extended to perform all computations in parallel using Nvidia’s CUDA framework and was tested to verify that it exhibited correct rotational behaviour. Since the underlying method remains the same as proposed in the aforementioned papers, for brevity it will not be repeated here except where specifically relevant to new components (in particular please refer to these papers for force evaluation in the solid matrix). One particular such case is that of the hookean constitutive material model, which has been replaced in the current research by the non-linear elastic neo-hookean model suggested by Wilson et al. [43]:

(3)

where is the determinant of the deformation tensor . The deformation tensor follows directly from the calculations to obtain the strain in the elastic solid equations of Solenthaler et al. [37], with the addition of the identity tensor, yielding the relation:

(4)

The bulk and shear modulus are given by:

(5)

(6)

where is the Young’s modulus and the Poisson’s ratio.

2.2. Collagen Fibril Network

Collagen fibre placement within the material is simplified with the assumption that a fibre is created between every particle and each of its neighbours, which in the case of the current tests are arranged in a regular lattice. There is a precedent for establishing fibre location according to element spacing [16]. This allows for the stress response of fibres to be modelled by the deformations and displacements affecting particle neighbourhoods rather than explicitly creating particles for their representation. The number of fibres differs depending on the particle position in the object volume.

The strain in each fibre is determined by:

(7)

where and represent the positions of particle i and its neighbour j, while the superscripts of 0 indicate the initial underformed positions. It is assumed that fibres only strengthen the material in tension so only positive strains are considered. Fibres are modelled as linear springs with stiffness in parallel with a nonlinear spring with stiffness as in the work by Li et al. [21], where and are the modulus and strain of the fibre, respectively. Even without viscoelastic properties this is assumed sufficient to capture the behaviour of the fibril network for the current research. Viscoelastic effects in fibres could not be implemented at this time due to numerical precision constraints.

In real cartilage tissue, fibre orientation is depth dependent, which as is typically assumed in FEM simulations to be 0° to the articular surface in the superficial zone (0% - 10% of depth), 45° in the middle zone (10% - 70% of depth), and 90° in the deep zone (70% - 100%). This is the depth dependency assumed for the current simulations. Individual fibre contribution to the stress profile can be found via a relationship between the particle orientations and the assumed fibre angle. Because the particle neighbourhoods and their associated fibres may have changed orientation, the new assumed fibre angle for the current particle’s zone is found by:

(8)

where is the current assumed fibre angle, is the original undeformed assumed fibre angle, while and are the current and undeformed angles between the particle and its neighbour. This leads to a multiplication factor M defined by:

(9)

The resulting equation for the stress contribution of each fibre is:

(10)

which includes the unit vector defining the direction between particles to assign directionality to the three-dimensional stress-strain relation. It should be noted that the contribution of the above equation is negative because it is defined based on the neighbour-to-current particle direction.

2.3. Interstitial Fluid

The behaviour of interstitial fluid through the porous structure of the non-fibrillar matrix is based upon the work by Lenaerts et al. [20] but has been significantly adapted to conform with the needs of the current research. Each SPH particle begins the simulation with an associated fluid mass that can then be exchanged between neighbouring particles and the external environment. The expression for the saturation of the material is given by Lenaerts et al. as:

(11)

where is the absorbed mass of the fluid in particle , is the fluid density, is the particle porosity, and the particle’s current volume, which is calculated as by Solenthaler et al. [37]. In the work by Lenaerts et al., the porosity was calculated based on the changing fluid density, but it is suggested here that it is beneficial to base the porosity on the changing fluid volume:

(12)

where is the initial porosity of the current particle, is the current porosity, and is the initial particle volume. Now if the ratio of the current volume to the initial volume is reduced to equal the initial porosity of the material (0.5 for example), the porosity will also become equal to zero, simulating all the fluid mass exiting from that porous particle. This is essential for the simulation of cartilage since when all the fluid is exuded from the material the load should be entirely borne by the non-fibrillar matrix and collagen fibrils.

Depth dependency of porosity is incorporated into the current procedure by:

(13)

where is the depth dependent porosity, is a material constant, is the height of the sample and is the current particle depth.

The absorbed fluid mass is dependent on the fluid density by the relation:

(14)

The density of the bulk material must be updated to reflect the additional fluid component to the solid component for use in the equations defining the non-fibrillar matrix:

(15)

Figure 1 provides an illustration of possible fluid forces acting on the fluid mass contained within an SPH particle during the simulation. The capillary potential acts to draw fluid into a porous material via surface tension, while the pore pressure can cause fluid to flow into or out of the material depending on the saturation.

Figure 1. Example fluid mass exchange between SPH particles. Black arrows represent pore (blue) pressure, white arrows represent capillary pressure, grey arrows represent resultant fluid flow between effective solid particle volumes (brown)

Each SPH particle should begin the simulation fully saturated, which leads to the next revision to the work by Lenaerts et al., in that fluid mass should be capable of being exchanged between fully saturated SPH particles. Expanding and contracting particle neighbourhoods to maintain a 100% saturation maximum was found to destabilize the simulation. Instead, the saturation of the material is allowed to exceed 100%, which provides for the necessary fluid behaviour. The capillary potential is then restricted by a boundary condition to prevent undefined behaviour when the saturation exceeds 100%:

(15)

The pore pressure is given by the state equation as in Lenaerts et al.:

(16)

where is again a user-controlled constant, and are the current and original densities, and is a user-controlled constant to enforce incompressibility.

The pore velocity , the rate at which the fluid mass is diffusing, is changed from the work of Lenaerts et al. so that it is more intuitively based on the difference in capilllary and pore pressure between a particle and its neighbours. This allows for the relationship to be reduced to a single equation:

(17)

where is the permeability of the material, is the viscosity of the fluid, and is the SPH kernel employed in the work by Solenthaler [37]. The permeability can be made strain dependent by the relation [43]:

(18)

where is the initial permeability, is a positive constant, and and are the initial and current porosity.

Fluid SPH particles are never explicitly created in the simulation, so fluid mass is exchanged in the material in a Eulerian manner defined by the differential:

(19)

The coefficient is a scalar that is the result of a scalar product between the pore velocity and the unit vector of the neighbour-to-current particle direction, ensuring appropriate contributions to the mass exchange by each neighbouring particle:

(20)

Here is a user-defined constant to control diffusion, noting that in the current procedure this will have an additional effect to that originally intended by Lenaerts et al. in that it will increase diffusion with saturations greater than 1. An explicit Euler integration is performed during the update section of each time step to diffuse the fluid mass:

(21)

Now that the fluid mass can diffuse between particles, the stress contribution during deformation is calculated based on the pore pressure of the current particle:

(22)

This calculation differs from the work by Lenaerts et al. in that it uses the pore pressure of Equation 16 instead of a constant value. This adjustment has been made to more accurately simulate increasing pore stress with increasing deformation of the cartilage material, an effect which can be seen in many FEM cartilage simulations. is a user-defined constant, as in Lenaerts et al.

Since the concern in the current research is only with the cartilage behaviour, it is not necessary to explicitly create SPH particles representing exterior fluid mass. This was accomplished through the creation of so-called void particles, which are essentially containers of infinite size with which the mass from the cartilage particles can be exchanged. They are instantiated in space according to the original particle positions as shown in Figure 2.

Void particles are only considered in two neighbour loops within the simulation. The first is during the calculation of equation 17, which for void particle neighbours becomes:

(23)

Figure 2. Void particle instantiation: The current particle position (blue) and the positions of its neighbouring particles (green) are checked to determine where void particles should be placed (white)

This is the same as equation 17 with the exception that the pore and capillary pressure of a void particle are assumed to be zero, and the volume of the void particle is assumed to be the same as the current particle (while recalling that void particles allow for infinite fluid flow, and so can be considered infinite for that purpose). The fluid mass differential equation is similarly updated with the fluid mass and volume of the current particle used instead:

(24)

Void particles are updated in the integration stage of the simulation to maintain their positions relative to their associated SPH particles.

2.4. Rigid Body Boundary Enforcement

Cartilage compression simulations require the creation of rigid materials to represent indenters and plates that can be used to cause deformation in the material, or for bone in the case of a whole hip joint simulation. This is accomplished here by the virtual particle method [24], which employs boundary particles along the surface of any rigid object with two layers of virtual particles for neighbourhood support. Pressure forces exerted by boundary particles on any cartilage particles that enter their sphere of influence are calculated according to the SPH fluid pressure force calculations [27] using Desbrun’s spiky kernel [10]:

(25)

which is equal to 0 for all other values of r, and where h is the smoothing length of the particle. Pressure is calculated with an incompressible equation of state [27]. The result is insignificant penetration by any cartilage particle into the spherical influence radius of a boundary particle.

2.5. Simulation Parameters

Simulations were performed on an Intel Core i5-3450 CPU with 8GB of RAM and an NVIDIA GeForce 660Ti with 2GB of memory. The C++ program created to perform the simulations was divided into an initialization section with a combination of host and CUDA code, and a main loop that was run each time step with all associated calculations performed in parallel using CUDA. Visualizations were created in the Maya 2013 API through direct correlations to the output of particle positions from the C++ program. Four indentation tests and one compression test were performed to confirm validity. A minimum of two simulations were run for each test, corresponding to a perfectly elastic material model and the proposed biphasic fibril-reinforced poroviscoelastic material model. Indenter tests were prepared to match the complete femoral knee compartment experiments of the female subjects of Vidal-Lesso et al. [41] as closely as possible (for more complete subject details refer to that paper). The compression test was similarly based upon the highly regarded FEM work by Li et al. [22]. These experiments and simulations were chosen due to restrictions on stability within the simulation. Because the stability is based upon the deformation of the non-fibrillar matrix, the fluid mass is currently not able to support in excess of 20% of the load. However, by increasing the permeability of the material (3 x 10^-14 m⁴/Ns versus 1.522 x 10^-15 m⁴/Ns of Wilson et al. [43]), short time relaxation to match the published results can be achieved via faster diffusion (i.e. the solid matrix component is stiffer but fluid diffuses more quickly, thereby decreasing the stiffness of the bulk material at a rate approximating the literature at short time spans). A Poisson’s ratio of 0.45 was selected to simulate near-incompressibility.

Indentation tests 1 and 2 were created to match specimen 1 of the work by Vidal-Lesso et al. [41], with a cartilage depth of 2 mm corresponding to a 5 particle thickness. Indentation tests 3 and 4 were created to match specimen 5 of the work by Vidal-Lesso et al. [41], with a cartilage depth of 3 mm corresponding to a 7 particle thickness. The cartilage thickness is constrained to the 1/2 mm initial particle spacing. Tests 1 and 3 were created as circular plugs 8 mm in diameter, while tests 2 and 4 were squares with 12 mm to each side. Creating both circular plugs and larger squares for similarly constrained simulations - with the larger squares having a greater similarity in surface area to the cartilage experiments of Vidal-Lesso et al. - allows for analysis of the effect of surface area on the simulations. In all indentation tests a circular indenter of 3 mm diameter applied a constant strain rate of 0.21 mm/s (as given in the above work; 0.07%/s strain rate equivalent) to achieve 0.5 mm displacement, followed by a relaxation period of 1.5s, and then an unloading period at the same rate. The compression test was similarly created to match the FEM work of Li et al. [22], where unconfined compression creates a 15% strain over a 15s period on a 3 mm diameter plug followed by 15s of relaxation. In this case, the simulations were not intended to be capable of mimicking the relaxation phenomena of the longer time span of that simulation. Rather, the concern was instead that the reaction force match their results at peak strain, thereby validating the simulation procedure for a second type of scenario in addition to the indenter tests.

Table 1. Material Parameters and Values Used in the Simulations

In a real case, a surgeon is not able to determine cartilage properties from a patient until that cartilage has been extracted during surgery. Therefore, for a preoperative case a set of predefined material parameters should be determined and used for all simulations. With this in mind, the material parameters were tuned by trial and error to indentation test 1 so that the biphasic simulation matched the published experimental results of Vidal-Lesso et al, while still remaining reliably stable for all testing scenarios. These same parameters can be and were used for subsequent simulations. This avoids the necessity of tuning for each one (as would be necessary for a mass-spring simulation).

Again, in the ideal case where the fluid mass is able to support a greater proportion of the load, the material parameters, and Young’s Modulus in particular, which in the biphasic simulations is an order of magnitude greater than the literature, could be made to reflect those values determined from experimentation. However, the material parameters used in these simulations need only be determined once for all subsequent simulations and they produce biomechanical cartilage behaviour that falls within the expected range in all cases (as will be seen in the following section), making them very useful for cases of preoperative evaluation where real tissue testing is not an option. Table 1 summarizes the material parameters used in the simulations. Gravity is ignored.

3. Results and Discussion

Figure 3 provides a visualization of Indenter test 4, which is the largest of the simulations with 6696 particles, and as such provides the most clarity for evaluation. Stress is represented as the average of the magnitudes of the principal stresses with increased brightness representing increased stress. At maximum strain (Figure 3b) the stress is concentrated around the area of the indenter, which over the following 1.5s spreads outward from that region due to the diffusion of fluid mass (Figure 3c). When the indenter is removed from the cartilage (Figure 3d), the material exhibits an area of low stress where the indenter had been because of lower fluid content, yielding a corresponding negative pore pressure.

Figure 3. Visualization of stresses in the particle system representing the cartilage in indenter test 4: a) t = 0s, b) t = 2.38 s, c) t = 3.88 s, d) t = 6.26 s

The snapshots shown in Figure 3 were taken during the simulation procedure as outlined in section 2. The only time when the material is in equilibrium is in Figure 3a. At the investigated rates of deformation there were no visible oscillations or vibrations exhibited by the particles in response to indentation. Upon numerical examination of individual particle fluid mass content, it becomes apparent that the stability of the system is highly dependent on the rate of deformation. At low rates there is very little fluctuation in the fluid mass. This fluctuation increases dramatically with increasing strain rate, eventually becoming unstable as the fluid mass attempts to balance between adjacent particles. This is the primary reason for the stability limit in these simulations, as it was found that a higher proportion of load could be borne by the fluid if rates of deformation were decreased. Lower deformation rates also allowed for a lowering of the solid matrix Young’s modulus closer to values presented in the literature, since more of the stress could be transferred to the fluid mass. It should be noted that the deformation rates in the presented simulations were reproduced from the work of Vidal-Lesso et al., where it was stated that the experimental parameters closely describe the rates experienced during walking.

Figure 4 illustrates the component stresses in the material at peak stain. Figure 4a shows that the fibril stress is largely constrained to the top layer, and to a lesser extent the layer beneath, with additional stress occurring at the bottom corners of the sample. This agrees with what would be expected, as the indenter only causes tensile strain in the horizontal direction so that only fibres aligned to the cartilage surface (superficial layer) and those oriented at 45 degrees (middle layer) would contribute. The stresses in the bottom corners are the result of the indenter causing deformations of the bulk material, yielding stress concentrations. In contrast, the pore stress visualization in Figure 4b is more evenly distributed throughout the layers, since the porosity decreases from a maximum of 0.7 at the cartilage surface to 0.55 in the deep layer. 0.7 is the maximum advisable porosity for the material due to the aforementioned stability constraints.

Figure 4. Visualization of component stresses at peak strain in indenter test 4: a) fibril stress, b) pore stress

The one-sided nature of the collagen fibre stress distribution in Figure 4a is due to the inability of the custom scripts created for Maya 2013 to create a perfectly circular indenter, which is due to the small scale of the simulations. It is thought that the behaviour of fibres is correct considering the uneven indenter geometry and the lack of friction in the simulation. This appears to indicate that the method of establishing fibres between particles in a neighbourhood, while heavily simplified, is sufficient for simulating cartilage to the desired level of accuracy (that accuracy being within the normal range of cartilage behaviour reported in the literature).The reaction force results of the four indenter tests are presented in Figure 5. Reaction forces are calculated by the summation of the unequal forces experienced by indented particles according to Newton’s 2nd law. The elastic simulation of Figure 5a attains strong correlation to the experimental results of Vidal-Lesso et al., but is not as strongly correlated as the biphasic simulation, which also matches the relaxation response. Similar trends can be seen through all four indenter tests. The differences between Figure 5c and Figure 5d emphasize the necessity of matching the surface area of the simulations to the experiments.

Figure 5. Indenter reaction force for biphasic and elastic simulations versus experimental [41]: a) Indentation test 1, b) Indentation test 2, c) Indentation test 3, d) Indentation test 4

Two important points can be noted based on Figure 5. The first is that the larger surface areas of Figure 5b and Figure 5d provide stronger correlation to the experimental. The original experiments by Vidal-Lesso et al. were performed on complete compartments of specimens rather than on cartilage plugs, and so the larger surface area more accurately represents those experiments. The second point is that the elastic simulations of indenter tests 3 and 4, where there is a higher cartilage depth, significantly undershoot the anticipated force responses. The biphasic simulations are more strongly correlated, particularly in the case of indenter test 4 where the biphasic simulation is capable of reproducing the experimental results with a higher degree of accuracy.

The reaction force results of the unconfined compression test are shown in Figure 6a. The longer time frame of this simulation allows the fluid mass more time to diffuse from the area of the indenter before the peak strain is reached so that less relaxation phenomena occurs in the biphasic simulations. The biphasic simulation that does not employ void particles to allow fluid mass to escape from the cartilage plug indicates a higher peak reaction force. This is expected since more fluid remains in the material to support the load. The biphasic simulation with void particles is approaching the elastic simulation force as the simulation continues, which is again to be expected since the solid components should support the entirety of the load in equilibrium.

The stress contributions of individual component stresses are provided in Figure 6b. The total material stress (axial) is calculated by averaging the reaction force over the surface area. The component stresses are numerical results experienced by a sample compressed particle, and their summation is shown to approximately equal the calculated total material stress. The results are what should be attained given the restrictions of stability on the load bearing capacity of the fluid, with the pore stress an order of magnitude lower than the matrix stress. It should be noted that the contribution of the fibril stress is nearly absent (~2 kPa maximum). Unlike in the indenter tests, the strain is evenly distributed across the surface of cartilage, and expansion in the material in the horizontal direction appears to only occur in the deeper layers that are not restrained by horizontal fibres. However, this behaviour may be the result of the small simulation scale.

Figure 6. Unconfined compression simulation results: a) Reaction force of elastic and biphasic (with and without void particles) simulations versus FEM results of Li et al. (“experimental”) [22], b) Component stress contributions and the total material stress

Computational results are summarized in Table 2. Pearson’s correlation coefficient was calculated using the statistical package of Microsoft Excel 2007. Correlation was not calculated for the compression tests since the relaxation of the peak force does not match the results of Li et al. There is a mean drop in frame rate from elastic to biphasic simulations of 29.03% with standard deviation of 0.62. Correlation was calculated to be stronger for biphasic simulations except in the case of indentation test 3.

Table 2. Computational Results of Simulations EXPERIMENT SIMULATION PARTICLE # BOUNDARY # FRAMES/S CORRELATION-r Indentation 1 Elastic 1624 246 98.871 0.976742 Indentation 1 Biphasic 1624 246 69.2433 0.98481 Indentation 2 Elastic 4774 286 46.1606 0.993154 Indentation 2 Biphasic 4774 286 32.8984 0.994124 Indentation 3 Elastic 3206 246 62.7095 0.98623 Indentation 3 Biphasic 3206 246 44.7148 0.992574 Indentation 4 Elastic 6696 286 32.6161 0.98132 Indentation 4 Biphasic 6696 286 23.2458 0.995798 Compression Elastic 98 352 172.403 Compression Biphasic 98 352 138.116

Pearson’s correlation coefficient is provided as a rough conceptualization of the relationship between experimental (and in the case of Li et al., FEM) and simulated results at all levels of strain. However, this statistic is scale-invariant and so is inappropriate for determination of the viability of simulation results without an accompanying visual examination. From an examination of Figure 5c, the biphasic simulation does in fact produce superior results for this test, despite having a lower calculated correlation.

It is worth mentioning again that the unconfined compression simulation (and indeed, all the rest) used parameters that were tuned to Indentation test 1, and so were not tuned specifically to this compression test. The high degree of accuracy of the results of this simulation in attempting to reproduce Li et al’s. FEM results is a testament to the applicability of this method to a great range of cartilage simulations. It should be noted that solid material simulation in SPH requires a greater volume of calculations in comparison to fluid SPH simulations (see for example [6, 35]) and so computational speeds are not directly comparable.

4. Conclusions

A new method has been presented to simulate articular cartilage in indenter and unconfined compression scenarios using SPH. Biphasic fibril-reinforced poroviscoelastic simulations produced superior results to Hookean material models in all cases. The decrease in frame rate of 29.03% (with standard deviation of 0.62%) is considered to be acceptable for the attained increase in accuracy. This may also be the first solid SPH simulation that is capable of producing accurate reaction force results to deformations.

The stability requirements of the simulations placed restrictions on the available research that could be used for comparison, since the vast majority of published results are focused on long-time compression rather than the short-time scenarios presented here. However, short time compression is likely more relevant to FAI since peak stresses occur during motion where articular cartilage does not demonstrate long term relaxation phenomena.

While the tuned parameters from indentation test 1 produced strongly correlated results in all subsequent tests, in the ideal case the stiffness and permeability parameters should be based on results determined through investigation of real cartilage samples. Increasing the stability of the simulations, in particular the fluid mass exchange, is therefore a primary concern with this procedure, but this cannot be accomplished by limiting particular variables since the fluid load support needs to be as high as possible to achieve closest approximations to experimental results. To alleviate this problem, an alternative to the simple leapfrog approach to integrating the particles of the solid portion of the material is suggested to be the next avenue of research, although it would have to be determined how to propagate the fluid flow in the material in this case, since it is currently performed in a Eulerian manner.

Unlike in the case of the indentation tests, it was unclear how to represent localized stresses in the unconfined compression test other than with the method of division over the surface area of the material that was used by Li et al. This is largely due to the current lack of a means to reduce all of the extremely complex interactions of the components to a single stress value that is relevant throughout the bulk material.

The lack of friction in the simulation is a concern that needs to be addressed. Specifically, it was necessary to restrain the bottom layer of the plugs from lateral movement otherwise the entire pieces would slide out from beneath the indenter without causing deformation. However, it should be noted that due to boundary lubrication and other mechanical effects, the real interaction between cartilage surfaces in a joint is nearly frictionless.

The simulations were validated against knee compartment indentation experiments and FEM compression simulations found in the literature. In the future, “short time” experiments (as opposed to the 100 ~ 1000 second experiments that comprise the majority of published work) on articular cartilage extracted from the hip joint to determine material parameters and behaviour would provide a more suitable basis for validation. However, it is thought that the studies used for comparison in this work provide sufficient accuracy.

The four tuned parameters as indicated in Table 1 above were found to be essential to the stability of the simulations. They were determined by a rigorous trial and error methodology that arrived at the values given that allow for the full range of material behaviour presented in the works of Vidal-Lesso et al. experiments and Li et al. simulations to be reproduced. Therefore it is highly recommended to use these same values in the reproduction of this work until the stability issues can be resolved. Nevertheless, although these parameters are not analogous to their real material equivalents, the material behaviour, suggested here to be a far more relevant consideration to a preoperative simulation of impingements, matches that expected of cartilage.

CUDA was found to have a 1.5X to 2.4X increase in computational speed over host code in these simulations. However, although the frame rates presented in Table 2 could be argued to be interactive, they do not represent a real-time simulation. It was found to be necessary due to the high material stiffness to decrease the time step corresponding to each frame to an extremely small size. It is thought that an improved integration method could also allow for a significant increase in time step size.

There is currently a 60,000 particle limit in the simulation, which was sufficient for the current simulations. However, a simulation of a complete hip joint for preoperative purposes may require significantly more particles than this. Since this limit was found to be directly relatable to the amount of onboard memory of the GPU (an investigation not shown in this paper), this limit could be increased by the use of a more powerful video card than the 2GB card used for the current simulations. Alternatively, it may be possible to devise a method of sending the particle data in batches to the video card to be operated upon rather than in bulk.

An important adjustment to the visualization of the simulations that needs to be made in the future is to provide a stress magnitude correlation to the displayed colour. Currently, the colourization is only provided for analysis of general trends, and actual stress values are determined from internal investigation of the resulting numerical values of individual particles. Such an adjustment would make the simulations much more pleasant to analyze. Nevertheless, the visualization in its current state does not detract from the accuracy of the reaction force results.

It is suggested here that the current research has reached a level of maturity to warrant the creation of an MRI-based simulation of a hip joint for further testing.