An analytical calculation of the fluid load support fraction in a biphasic material: an alternative technique.

BACKGROUND: The fluid load support fraction (W(F)/W(T)) can be used to define the mechanical contribution of the interstitial fluid (W(F)) to the total force (W(T)) in the deformation of cartilage. Traditionally, W(F)/W(T) is calculated using complex experimental setups or time-consuming micromechanical poroelastic Finite Element (FE) simulations. AIM: To define and validate a fast and efficient technique to predict W(F)/W(T) using an analytical approach that can be applied without micromechanical detail or experimental measurement.
METHODOLOGY: Poroelastic FE simulations defined accurate values of W(F)/W(T) for a range of loading configurations and were used to validate subsequent predictions. The analytical prediction of W(F)/W(T) used elastic contact mechanics to calculate W(F), and viscoelastic FE representation to calculate W(T). Subsequently, these independent calculations of W(F) and W(T) provided values of W(F)/W(T) that were compared with the poroelastic FE calculations. RESULTS AND DISCUSSION: The analytical prediction of W(F)/W(T) proved effective and suitably accurate (mean difference S<0.05). This technique demonstrated how W(F) and W(T) can be determined independently, without a biphasic constitutive model. Here we used viscoelasticity to calculate W(T) as an example, however, W(T) could be measured experimentally or predicted computationally.

Introduction

Cartilage is a biphasic material with two constitutive components that support load, a solid phase and a fluid phase. The fluid load support fraction (W/W) describes the fraction of total load (W) that is supported by the interstitial fluid pore pressure (W) [Ateshian et al., (2003); Ateshian (2009)]. Specifically, W is determined from the integration of the interstitial fluid pore pressure p over the contact area A,and W from the integration of the contact pressure, where n is defined as the unit outward normal on dA of the stress σ, over the contact area A,

Traditionally, experiments or micromechanical numerical analyses are used to quantify W/W, i.e. a counter-sunk pressure sensor has been used in experimental studies [Basalo et al., (2004); Krishnan et al., (2004); Park et al., (2003)], whereas micromechanical numerical studies [Soltz and Ateshian (1998)] have used poroelastic finite element (FE) simulations [Park et al., (2003); Pawaskar et al.]. Yet, as experimental studies generally measure isolated cartilage specimens, these techniques cannot measure W/W during in-vitro physiological loads in whole-joints. Furthermore, given that poroelastic FE simulations are typically based on a quasi-static constitutive foundation, using an implicit solution method, these numerical calculations cannot model the dynamic conditions present in physiological loads [Stops et al., (2011)]. Thus, there is no constitutive framework that enables W/W to be calculated in cartilage during dynamic analyses.

In light of the above, and considering the need to calculate W/W in (1) a dynamic load condition and (2) in a whole joint: an explicit formulation of FE enables both a dynamic framework for these conditions. Though the quasi-static framework that underpins poroelasticity is not compatible with explicit FE, due to the problems in implementing inertial effects, explicit FE can model time-dependent effects using a phenomenological approach, e.g. viscoelasticity. The implementation of viscoelasticity does not allow for the micromechanical features of cartilage to be analysed, as poroelasticity allows, but it does provide a reliable representation of the macroscopic mechanics. Furthermore, a viscoelastic representation offers versatility in terms of modelling options, i.e. it can be used in a FE solution, and also a multibody dynamics solution (which is the primary technique for musculoskeletal modelling). Thus, in whole-joint analyses and in physiological loading conditions, viscoelasticity offers an ideal method to represent the time-dependent effects of cartilage.

Hence, a computational methodology that can predict W/W alongside a viscoelastic representation would provide the opportunity to understand cartilage mechanics in a range of settings, from whole-joint analyses to dynamic loading conditions. Consequently, this paper aims to validate an analytical prediction of W/W which complements a viscoelastic representation of cartilage, without the need for poroelastic simulations or complex experiments.

Methodology

Poroelastic W/W calculation

Poroelastic FE simulations were employed to calculate ‘gold standard’ values of W/W; such an approach has been shown to be accurate [Pawaskar et al.]. Axisymmetric models with contact-dependent fluid flow at the contact interface, free flow at the edges, and restricted flow at the bone–cartilage interface were generated. Uniaxial loads were applied through an indenter, perpendicular to the unconfined material, which was bonded to a rigid foundation (Table 1 and Fig. 1). Nine poroelastic FE simulations were performed: three indenter geometries included (1) a rigid planar-surface, (2) a large rigid sphere and (3) a small rigid sphere, all with (a) ramp loads, (b) constant strains (stress relaxation) and (c) sinusoidal strains (Table 1). In all cases, the mesh size was adjusted until convergence (Table 1) was achieved. For these nine conditions, W/W was calculated based on Eqs. (1) and (2).

Table 1

Material Properties and Loading Configurations.

Material formulation	Equilibrium elastic modulus E∞ (MPa)	Poisson's ratio ν∞	Permeability (m4/N s)	Fluid content (%)
Poroelastic	4	0.15	1x10−15	75
Material formulation	Instantaneous elastic modulus E0 (MPa)	Equilibrium elastic modulus E∞ (MPa)	Time decay τ (sec)	Poisson's ratio ν0
Viscoelastic	5.19	4	Planar 677	0.49
Large sphere 469
Small sphere 14

⁎ The viscoelastic properties were implemented in a generalised form, where the generalised relaxation modulus was defined as , where t is time and E is the ratio of the equilibrium and the instantaneous modulus, respectively. Note the three loading configurations resulted in viscoelastic parameters with equivalent elastic properties, but with varying time decay constants. These parameters are a direct result of the contact-dependent fluid flow boundary at the contact interface: the planar surface prohibited flow over the entire upper surface, whereas the small spherical indenter prohibited flow for less than 1% of the upper surface, thus producing very different boundary conditions between each loading configuration.

Fig. 1

An illustration of the load setup for (a) a planar-surface indenter, (b) a large spherical indenter and (c) a small spherical indenter; note the indenter dimensions are for illustrative purposes only, please refer to Table 1 for accurate dimensions. The dimensions of the FE meshes were 12×2 mm, wherein the element size of each mesh (for each loading configuration) was optimised through a sensitivity analysis.

Analytical W/W prediction

Using the W/W values from the poroelastic simulations as a ‘gold standard’, in much the same way as experimental data is used to validate a computational simulation, the subsequent analytical predictions were compared. These analytical predictions were based on the concept that the two constituent phases of cartilage, the interstitial fluid and the deformable solid matrix, can be determined independently and that the mechanical contribution of the fluid can be approximated from the difference between the load supported by the solid matrix (W) and the total load supported by the material (W), [Bonnevie et al., (2011)]. Assuming the deformable solid matrix to be elastic, the value of W was calculated using elastic contact mechanics.

Here, three different equations were used to represent the different indenter geometries, as outlined below. In each case E∞ is the equilibrium elastic modulus, h is the thickness of the cartilage, v∞ is the equilibrium Poisson's ratio. The penetration depth δ and the total load W were determined from the respective viscoelastic simulation.

For the planar-surface configuration, W/W was analytically predicted from the equation developed by Johnson (1985):where with a1 being the contact radius. Note, Eq. (3) is used to calculate the force in the solid phase of the material, and hence, all properties relate to the equilibrium state.

For the large spherical indenter, where the contact radius is large relative to the cartilage thickness, W/W was calculated from:where R is the radius of the sphere, a2 is the contact radius calculated by [Jaffar (1997)].

And finally, for the small spherical indenter, where the contact radius is small relative to the cartilage thickness, W/W was calculated from:where D is calculated as and a3 is the contact radius calculated by ; wherein a Poisson's ratio of 0.15 was assumed [Jaffar (2008)]. It should be noted that alternative asymptotic solutions to Eq. (5) also exist [e.g. Argatov, (2010); Argatov, (2011)].

In Eqs. (3–5), values of W were calculated from viscoelastic FE models whose properties (Table 1) were calibrated to the force–deformation response of the corresponding poroelastic model in the stress-relaxation setup. Both the instantaneous and the equilibrium response of the viscoelastic models proved accurate in representing the equivalent poroelastic models, as can be seen in Table 1 (Viscoelastic Calibration Error). Note, the mesh/configuration used for each loading scenario was identical in the respective poroelastic and viscoelastic models to ensure that there were no compounding influences on the solutions.

Statistical comparison

The nine analytical W/W predictions were compared to the ‘gold standards’ using the mean difference (error) S,where and are the ith increment from the ‘gold standard’ calculations and the analytical, respectively, while graphical illustrations are presented in terms of .

Results

Planar-surface indenter: W/W

For the ramp load, the analytical W/W prediction proved accurate, as can be defined by S=0.01 (Fig. 2a). Likewise, for the constant 2.5% strain, the predicted values were close to the ‘gold standard’, S=0.04 (Fig. 2b), and for the sinusoidal strain, the predicted values again followed the trend of the standard values, with S=0.07, but with discrepancies when the sine curve prescribed a turning point from unloading to loading, i.e. at time points T=2, 4, 6 and 8 s (Fig. 2c).

Fig. 2

Values of the fluid load support fraction (W/W) calculated by a poroelastic FE simulation “gold standard” and predicted by an analytical calculation of the elastic contact mechanics (shown with errors defined as ) during loads with a planar surface: (a) shows a linear ramp load up to 2.5% compressive strain (equivalent to 325 MPa load), (b) shows a constant 2.5% strain and, (c) shows a sinusoidal strain varying from 0.5 to 2.5% strain.

Large spherical indenter: W/W

The predicted values of W/W demonstrated S=0.05 during the ramp loading (Fig. 3a), and S=0.06 for the constant 2.5% strain (Fig. 3b), while the sinusoidal loading demonstrated values in line with the ‘gold standard’, but with discrepancies at the unloading-loading turning points, S=0.05 (Fig. 3c).

Fig. 3

Values of the fluid load support fraction (W/W) calculated by a poroelastic FE simulation “gold standard” and predicted by an analytical calculation of the elastic contact mechanics (shown with errors defined as ) during a series of loads with a large spherical indenter: (a) shows a linear ramp load up to 2.5% compressive strain (equivalent to 16 MPa load), (b) shows a constant 2.5% strain and, (c) shows a sinusoidal strain varying from 0.5 to 2.5% strain.

Small spherical indenter: W/W

The predicted values of W/W produced S=0.05 during the ramp loading (Fig. 4a), and S=0.01 for the 2.5% strain (Fig. 4b). Again, the W/W for the sinusoidal strain proved similar to the ‘gold standard’, yet with discrepancies at the unloading-loading time points, S=0.16 (Fig. 4c).

Fig. 4

Values of the fluid load support fraction (W/W) calculated by a poroelastic FE simulation “gold standard” and predicted by an analytical calculation of the elastic contact mechanics (shown with errors defined as) during a series of loads with a small spherical indenter: (a) shows a linear ramp load up to 2.5% compressive strain (equivalent to 0.02 MPa load), (b) shows a constant 2.5% strain and, (c) shows a sinusoidal strain varying from 0.5 to 2.5% strain. Note, the negative W/W values shown in Fig. 4c, are a result of suction at the contact surface, caused by negative pore pressures during cyclic loading [Krishnan et al., 2005].

Discussion

The analytical prediction of W/W proved, in general, accurate. The planar-surface configuration demonstrated the most desirable results (Fig. 2), while the spherical indenters also proved acceptable (Figs. 3 and 4). However, important discrepancies occurred for the sinusoidal loading configuration for all indenters, and in particular, at the instances when unloading turned to loading (t=2, 4, 6 and 8 s in Figs. 2c, 3c and 4c). It is interesting to note, that this analytical prediction was also capable of capturing the negative W/W values during cyclic loading (Fig. 4c), which as Krishnan et al. (2005) has suggested, is an intriguing and counter-intuitive response by which cartilage supports load.

Limitations

Interestingly, the level of accuracy provided by the analytical approach varied for the three indenter types. This variation most probably arose from the accuracy (or more precisely, the correct use of) the elastic contact equations. For example, eq. (4) assumes a set of boundary conditions that require a rigid sphere loaded on a thin layered material that is bonded to a rigid foundation, where and 0≤v∞≤0.5, [Jaffar (1997)]. Given that the large spherical indenter exhibited a range of a/h ratios during the ramp loading, from ∼0.001 at the instant of load application to a maximum of 6.32 at 2.5% strain, the relatively high absolute errors (S=0.14) in Fig. 3a may arise from a non-conformance of the boundary conditions required by Eq. (4), i.e. at the instant of load application a/h→0 and a conflict between the boundary conditions and the mechanics arose. This predicament is likely to occur in many scenarios when, for example, the bearing surfaces do not fully conform to either the idealised dimensions or the perfect contact, defined in the given contact equation. Thus, it is prudent to note the importance of implementing the correct elastic contact model in order for this analytical W/W calculation to be effective, and consequently, it could be beneficial to derive the elastic contact mechanics from configuration-specific force-deformation data. For scenarios that involve complex geometries, such as a human joint, an elastic FE simulation could provide such information.

A limitation of the viscoelastic simulations was the need to determine the time-constant parameter for each specific geometry (Table 1); as the contact area varied for each indenter geometry, the boundary flow conditions also varied (the contact interface had effectively a no-flow condition due to the impermeable representation of the indenter), and thus, the cartilage exhibited different rates of interstitial fluid pressure reduction for each indenter geometry. However, the instantaneous and the equilibrium moduli were both consistent across all indenters. Thus, when applying this technique to an in-vivo setting, the time constant is the only parameter that may prove difficult to ascertain; the moduli can be readily obtained from experimental studies or published sources. Given the complex boundary conditions present in-vivo (the biphasic jump condition), it may be worthwhile undertaking a computational poroelastic study to acquire this data (an example can be found in [Pawaskar et al., 2010)]. However, we believe this inconvenience is overridden by the benefit of being able to calculate W/W in a dynamic setting.

Applications

The analytical W/W prediction can be employed in many scenarios where the total load W and the penetration depth δ are known. Here, a viscoelastic FE simulation was employed to offer an example, and intended, application. However, an experiment that measures W and δ can apply this technique, and likewise, an explicit FE simulation that uses time-dependent mechanics can also employ such a method. What is important, is that the proposed technique measures/calculates W and δ: future calculations of W/W do not need to arise solely from counter-sunk pressure sensors in expensive experiments, or from time-consuming poroelastic FE simulations.

The proposed analytical method is intended for the use in explicit FE and multibody simulations. It is hoped that this technique will provide essential information on material mechanics within a whole joint during physiological loading conditions. The ability to understand W/W in a dynamic loading environment is imperative to understand how degenerative conditions progress and for the design of replacement therapies (medical implants).

Conclusions

The analytical technique provided a suitable means for first approximation of W/W, and in situations where micromechanical continuum poroelastic models cannot be employed, this technique offers a calculation with a reasonable level of confidence. This verified technique is not intended to replace current experimental or poroelastic FE modelling techniques, but rather to offer an alternative approach for when these traditional calculation techniques are impractical.

Conflict of interest statement

There are no conflicts of interest.