- Email: [email protected]
Simulation of high tensile Poisson's ratios of articular cartilage with a finite element fibril-reinforced hyperelastic model
Available online at www.sciencedirect.com
Medical Engineering & Physics 30 (2008) 590–598
Simulation of high tensile Poisson’s ratios of articular cartilage with a finite element fibril-reinforced hyperelastic model Jos´e Jaime Garc´ıa ∗ Escuela de Ingenier´ıa Civil y Geom´atica, Universidad del Valle, Calle 13, Carrera 100, Edificio 350, Cali, Colombia Received 7 February 2007; received in revised form 26 June 2007; accepted 27 June 2007
Abstract Analyses with a finite element fibril-reinforced hyperelastic model were undertaken in this study to simulate high tensile Poisson’s ratios that have been consistently documented in experimental studies of articular cartilage. The solid phase was represented by an isotropic matrix reinforced with four sets of fibrils, two of them aligned in orthogonal directions and two oblique fibrils in a symmetric configuration respect to the orthogonal axes. Two distinct hyperelastic functions were used to represent the matrix and the fibrils. Results of the analyses showed that only by considering non-orthogonal fibrils was it possible to represent Poisson’s ratios higher than one. Constrains in the grips and finite deformations played a minor role in the calculated Poisson’s ratio. This study also showed that the model with oblique fibrils at 45◦ was able to represent significant differences in Poisson’s ratios near 1 documented in experimental studies. However, even considering constrains in the grips, this model was not capable to simulate Poisson’s ratios near 2 that have been reported in other studies. The study also confirmed that only with a high relation between the stiffness of the fibers and that of the matrix was it possible to obtain high Poisson’s ratios for the tissue. Results suggest that analytical models with a finite number of fibrils are appropriate to represent main mechanical effects of articular cartilage. © 2007 IPEM. Published by Elsevier Ltd. All rights reserved. Keywords: Anisotropy of cartilage; Tensile Poisson’s ratio; Hyperelastic fibril-reinforced models
1. Introduction Articular cartilage is a highly heterogeneous and anisotropic material covering the end of the bones inside diarthrodial joints. Towards an understanding of the relations between mechanical parameters and osteoarthrosis, many analytical and experimental studies have been developed to determine the mechanical response of articular cartilage. This response is mainly represented by non-linear stress–strain curves under transient loads and equilibrium [1,2], nonlinear tension–compression response [1], intrinsic viscous effects of the solid phase [3,4], and non-linear permeability [5]. Anisotropy of articular cartilage has been further demonstrated in various experimental studies exhibiting tensile Poisson’s ratios higher than one [1,2,6]. A finite element model able of describing most of the independent experiments is a valuable tool to accurately assess the mechanical ∗
Tel.: +57 2 339 2420; fax: +57 2 339 2086. E-mail address: [email protected].
performance of articular cartilage, needed in many practical applications. Within the framework of the biphasic theory [5] and among various approaches that have been proposed to represent the mechanical behavior of the solid phase, fibrilreinforced models have been successful of representing main experimental observations as well as the microstructure of the tissue [7,8,9,10]. Only few of these models have documented Poisson’s ratios under tension. With the assumption of infinitesimal strains Wu and Herzog [11] use homogenization theory to determine the effective elastic constants of a mixture of matrix, orthogonal fibrils and chondrocytes, documenting Poisson’s ratios from the stable models lower than 0.4, which are equally valid for tension and compression since their model does not consider the tension–compression non-linearity. A recent study by Klisch et al. [12] considering finite deformations and a bimodular-orthotropic function documents Poisson’s ratios under tension one order of magnitude lower than those found experimentally. A couple of
1350-4533/$ – see front matter © 2007 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2007.06.004
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
models [4,13] that consider the tension–compression nonlinearity by allowing only tension in the fibrils have been incapable to describe high tensile Poisson’s ratios. One possible explanation for this inconsistency, as suggested by Li et al [4], is that, due to the constrains in the grips, the stress field under a tension test is really two-dimensional, and that the Poisson’s ratio depends on the aspect ratio of the specimen used in the experimental set up. To better represent the architecture of the tissue, Wilson et al. [8] incorporate sets of secondary fibrils at 45◦ but no Poisson’s ratios under tension have been documented with this model. A recent analysis by Ley and Szeri [14] with a fibril-reinforced model under infinitesimal deformations and true uniaxial loading shows that high tensile Poisson’s ratios may be simulated considering non-orthogonal fibril distributions. It has not yet been investigated the capability of models with a finite number of fibrils to reproduce high tensile Poisson’s ratios. Moreover, no study has been conducted to determine the influence of constrains in the grips and finite deformations in the measurement of the tensile Poisson’s ratio of articular cartilage. A finite element fibril-reinforced hyperelastic model was used in this study to simulate high tensile Poisson’s ratios documented in experimental studies by Huang et al. [1], Charlebois et al. [2] and Elliot et al. [6]. The study was specially aimed to quantify the effects in the Poisson’s ratio of constrains in the grips, alignment of the fibrils and finite deformations.
2. Methods For the equilibrium simulation treated here no permeability or intrinsic viscous parameters of our general finite element model [13] were required. In addition, the swelling pressure was not included considering that its main effect is to increase the stiffness of the non-fibrillar matrix, which is
591
still low compared to the stiffness of the fibers. Thus, the solid phase of articular cartilage was represented by the proteoglycan matrix, described by the isotropic hyperelastic function proposed by Holmes and Mow [15], reinforced with four sets of collagen fibers, described by the hyperelastic function proposed by Limbert and Middlelton [16]. The hyperelastic function of the matrix is fully defined by three elastic parameters α0 , α1 , α2 that can be adjusted to fit the Young’s modulus and Poisson’s ratio of the matrix at zero deformation and the degree of non-linearity of the stress–strain curve under compression [17]. The hyperelastic function of the fibers is defined by the Young’s modulus at zero deformation Ef0 and a parameter γ that can be used to adjust the degree of non-linearity of the stress–strain curve under tension. One objective of this study was to determine the capability of a model with a finite number of fibrils to explain experimental results. A model like this may be efficiently implemented in a 3-D geometry. Hence, similar to the secondary set of fibrils used by Wilson et al. [8], four sets of fibers were considered, one in the axial (y) and one in the transverse (x) direction of loading (orthogonal fibers), and two oblique fibrils in a symmetric configuration respect to the orthogonal axes (Fig. 1). The angle of the oblique fibers respect to the transverse axis was changed to see the influence in the Poisson’s ratio of the tissue. Simulations were conducted of the tensile experiments undertaken by Huang et al. [1] with humeral and glenoid human cartilage and Charlebois et al. [2] and Elliot et al. [6] with bovine and human patellar cartilage, respectively. The finite element models consisted of four-node plane-stress elements distributed in a domain equal to a fourth of the experimental samples, according to the divisions shown in Table 1, with symmetry conditions prescribed in the orthogonal axes (Fig. 1). Vertical displacements were prescribed in the upper nodes of the model to complete 12, 7.7 and 16% axial deformation, for the studies of Huang et al. [1],
Fig. 1. Dimensions of the specimens, meshed region, orientation of the oblique fibrils and mesh for the data by Charlebois et al. [2].
592
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Table 1 Geometry and mesh for each experimental study
Charlebois et al. [2] Huang et al. [1] Elliot et al. [6]
Dimensions between grips (mm) length (L) × width (w)
Aspect ratio (L/w)
Number of divisions along each direction
Axial deformation (%)
4×2 5 × 1.5 8.33 × 1.91
2 3.33 4.36
8×4 10 × 3 16 × 4
7.7 12 16
Charlebois et al. [2] and Elliot et al. [6], respectively. In the updated Lagrangian finite element code the fibers were taken into account by calculating at each integration point the residual vector according to the hyperelastic function by Limbert and Middleton [16]. To improve convergence of the iterative scheme, the tangent stiffness matrix for the fibers was calculated and added to the tangent stiffness of the matrix. In order to study the influence of the grips in the Poisson’s ratio, analyses were accomplished considering the transverse constrains in the grips (constrained model) and without constrains in the grips (unconstrained model). The unconstrained model was intended to reproduce a true uniaxial loading condition. Additionally, the Poisson’s ratio for the constrained model (apparent Poisson’s ratio) was calculated with the lateral displacement at the center and with the average of the lateral displacement at the center and at the quarter length of the model, similar to Huang et al. [1]. The Poisson’s ratio was calculated as the negative ratio of transverse strain over axial strain. For the analysis under finite deformations and in order to take into account important variations of geometry respect to the initial configuration, another Poisson’s ratio was calculated, called here the logarithmic Poisson’s ratio, defined as the negative ratio of the transverse logarithmic strain over the axial logarithmic strain. A first analysis was conducted using the Hooke’s law to develop a simple formula (presented in Appendix A and called here the theoretical formula) for calculating the Poisson’s ratio of the whole tissue in terms of the Poisson’s ratio of the matrix and the stress component of the fibers in the transverse direction. The theoretical formula was used with the transverse stress component obtained from finite element analysis to predict a value, called here the theoretical Poisson’s ratio. The transverse stress component of the fibers was normalized (normalized transverse stress) by dividing it over the product of the axial strain and the Young’s modulus of the matrix. The theoretical formula and preliminary finite element results (Fig. 2) showed that the tensile Poisson’s ratio of the tissue was independent of the Poisson’s ratio of the
Fig. 2. Apparent Poisson’s ratio of the humeral surface for two values of the Poisson’s ratio of the matrix and different orientations of the oblique fibers.
matrix for a wide range of fibril orientations. Thus, a value of 1/3 was adopted in this study for the Poisson’s ratio of the matrix, which appears to be more consistent with Poisson’s ratios documented in compression experiments [18]. Hence, following a procedure described by Garc´ıa and Cort´es [13] various sets of material properties were calculated for the matrix and the fibrils using the data reported by Huang et al. [1] for the glenoid and humeral head. From the tensile experiments documented by Charlebois et al. [2] and Elliot et al. [6] it is not possible to determine the elastic constants of the matrix. Thus, the same properties for the matrix used for the data by Huang et al. [1] for humeral cartilage were initially adopted in the analyses of the data by Charlebois et al. [2] and Elliot et al. [6]. Next, other analyses were conducted considering a Young’s modulus for the matrix equal to half and a third that obtained from the study by Huang et al. [1]. For each orientation of the oblique fibers, fibril properties were adjusted to yield the reported experimental Young’s
Table 2 Elastic properties for the matrix used in the simulations of the tensile experiments Assumed νm0 H-s H-m G-s G-m
1/3
Reported HA0 (MPa)
Calculated Em0 (MPa) Em0 = 1.5HA0
Reported n
0.116 0.141 0.138 0.178
0.174 0.212 0.207 0.267
1.15 0.81 1.16 0.88
Calculated α1 = 0 0
Calculated α2 = n/2
Calculated α0 (MPa) α0 = 3Em0 /16α2
0.58 0.41 0.58 0.44
0.0563 0.097 0.0669 0.114
These properties were calculated from the confined compression data reported by Huang et al. [5]. H-s: humeral surface, H-m: humeral middle, G-s: glenoid surface, G-m: glenoid middle.
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
593
modulus E0+ in the axial direction using the well-known theory of composite materials [19] under uniaxial loading and infinitesimal deformations. The Poisson’s ratio obtained from composite materials theory was also used to verify results of the finite element code under infinitesimal deformations and true uniaxial loading. For the calculation of the non-linear stiffness of the material, the parameters A and B of the formula σ = A exp(Bε) used by Huang et al. [1] and Elliot et al. [6] at a relaxation time of 30 min were adjusted with the hyperelastic function for the fibers as explained by Garc´ıa and Cort´es [13]. Since Charlebois et al. [2] report the lateral displacement profile at a relaxation time (30 min) different from the times used to report stiffness (full relaxation and 15 min), the same fibril properties by Elliot et al. [6] were used in the simulation of the data by Charlebois et al. [2]. A summary of material properties is presented in Tables 2 and 3.
3. Results Comparison of lateral displacement profiles (Fig. 3) obtained with various angles for the oblique fibers with those obtained experimentally [2] indicates that oblique fibers at about 60◦ respect to the transverse axis, i.e. 30◦ respect to the axis of loading, give a better description of the lateral profile (Fig. 4). For oblique fibers at 45◦ and for orthogonal fibrils the variation of the lateral profile with the Young’s moduli of the constituents was negligible (Fig. 3). For oblique fibrils at 60◦ , increments in the maximum lateral displacement (and thus of the Poisson’s ratio) of 18 and 30% were obtained for 2- and 3-fold increments, respectively, of the ratio of Young’s moduli of the constituents. The finite element model yielded a good approximation of the infinitesimal-strain closed-form solution under 1% deformation (Fig. 5), except for angles of the oblique fibrils less than 30◦ , where the closed-form solution from composite materials, which does not consider the tension–compression non-linearity, predicts lower values of the Poisson’s ratio. The theoretical formula (Eq. (A.4)) predicts high variations of the Poisson’s ratio with the normalized transverse stress. Mainly, for low values of the normalized transverse stress, the Poisson’s ratio is approximately equal to the Poisson’s ratio of the matrix. On the other hand, for higher values of the normalized transverse stress, the Poisson’s ratio is approximately equal to the normalized transverse stress. Foregoing results are consistent with those obtained with the finite element program (Fig. 2) showing that the Poisson’s ratio of the tissue is approximately equal to that of the matrix for low (<25◦ ) and high (>85◦ ) angles of the oblique fibers and tends to be higher than one for a range of middle angles of the oblique fibers (∼45–75◦ ). Moreover, the logarithmic Poisson’s ratio showed a good correlation with the Poisson’s ratio obtained with the theoretical formula (Fig. 6). Regardless of the effects of the grips and aspect ratio of the
Fig. 3. Profiles of lateral displacement for three angles of the oblique fibrils and different relations between the Young’s modulus of the fibers and that of the matrix. Experimental data by Charlebois et al. [2]. (a) Orthogonal fibers; (b) oblique fibers at 45◦ ; (c) oblique fibers at 60◦ .
594
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Table 3 Tensile stiffness used in the simulations Reported A (MPa)
Reported B
E0+ (MPa)
Non-linear parameter γ
Huang et al. [1], human
H-s H-m G-s G-m
0.652 0.422 0.562 0.296
10.27 8.80 8.49 7.98
6.57 3.55 4.55 2.24
4.6 3.4 3.2 2.8
Elliot et al. [6], human
P-s P-m
0.411a 0.381a
14.15a 5.05a
5.81a 1.92a
8.1 0.5
Charlebois et al. [2], bovine
P-s
0.411a
14.15a
5.81a
8.1
P-s: patellar surface, P-m: patellar middle. a Calculated from the reported data.
specimens (Fig. 7), Poisson’s ratios higher than one could only be obtained considering fiber orientations in the range of ∼45–80◦ , which is consistent with the range of angles displaying higher values of the transverse stress component and with the values predicted with the theoretical formula (Fig. 6). Even specimens with the lowest aspect ratio (2.0) yielded an apparent Poisson’s ratio with maximum differences of about 10% (Fig. 7a) respect to the value predicted by the unconstrained model. For the lowest aspect ratio the appar-
Fig. 5. Poisson’s ratios obtained with the closed-form solution from composite material theory and the finite element program under 1% deformation. This simulation was accomplished for a Poisson’s ratio of the matrix equal to 0.33.
ent Poisson’s ratio obtained with the lateral displacement at the center was better than that obtained with two points, while the opposite trend was observed for the higher aspect ratios (Fig. 7). For a range of angles for the oblique fibers
Fig. 4. Deformed mesh and contours of transverse displacement for the samples of Charlebois et al. [2] with the oblique fibers at 60◦ with respect to the transverse axis. The whole sample was considered in this simulation.
Fig. 6. Comparison of the Poisson’s ratio, the logarithmic Poisson’s ratio, the theoretical Poisson’s ratio and the normalized transverse stress from analyses undertaken with the data by Huang et al. [1] for the humeral surface.
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
595
Fig. 7. Poisson’s ratio calculated with and without constrains in the grips for different values of the aspect ratio of the specimen. With constrain in the grips the Poisson’s ratio was calculated with the average of the lateral displacement at two points and with the lateral displacement at the center (one point). (a) Charlebois et al. [2] L/w = 2, E0+ /Em0 = 33.4. (b) Huang et al [1] L/w = 3.33, E0+ /Em0 = 16.7. (c) Elliot et al. [6] L/w = 4.36, E0+ /Em0 = 11.0.
(∼50–80◦ ) a positive correlation was observed between the apparent Poisson’s ratio and the ratio of the Young’s moduli of the constituents (Fig. 7). A reduction of Poisson’s ratio with strain was observed for fiber orientations of 60◦ and no significant changes were obtained for fiber orientations of 30◦ and 45◦ (Fig. 8a). However, when using the logarithmic formula, slight increases can be observed for fiber orientations of 60◦ and 45◦ , and a good correlation with the Poisson’s ratio calculated with the theoretical formula (Fig. 8b). For the data by Charlebois et al. [2] and the patella surface by Elliot et al. [6] the experimental Poisson’s ratios were about 2-fold those obtained with the finite element model with fibers at 45◦ (Fig. 9). For the other experimental data the
finite element model with the fibers at 45◦ provided a fairly good approximation of the experimental Poisson’s ratio, with percent differences in the range of 8–20% (Fig. 9).
4. Discussion This study showed the ability of fibril-reinforced models under finite deformations to describe high Poisson’s ratios reported in tensile tests of articular cartilage. The result was feasible given the anisotropy of the tissue, mainly characterized by a high relation between the stiffness of the fibrils and that of the matrix. Only considering oblique fibrils with orientations of the order 45–80◦ respect to the transverse axis
596
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Fig. 8. Variations with strain of Poisson’s ratio (P) and the logarithmic (L) Poisson’s ratio for the data by Huang et al. [1] for the humeral surface, and comparison with the theoretical formula (T) documented in Appendix A, for three orientations of the oblique fibers (60◦ , 45◦ and 30◦ ). (a) Poisson’s ratio (P) calculated as the negative ratio of transverse to axial strain. (b) Poisson’s ratio (L) calculated with the logarithmic formula.
Fig. 9. Comparison of the experimental and numerical Poisson’s ratios with the oblique fibers at 45◦ .
of loading was it possible to simulate Poisson’s ratios higher than one, which is in agreement with recent analyses of a reinforced model employing fibril distribution functions under infinitesimal deformations and uniaxial loading [14]. Unlike the use of fibril distributions, the consideration of a limited number of fibrils allows a relatively easy implementation into a non-linear finite element code. The need to consider oblique fibrils for a proper description of the high tensile Poisson’s ratio suggests that articular cartilage contains, even for directions parallel to the split lines, fibril distributions aligned in non-orthogonal directions, which appears to be consistent with the findings by Clark [20] who documents fibrils in the radial zone of cartilage not perfectly straight or parallel. The normalized transverse stress showed to be a good predictor of the Poisson’s ratio of the tissue. In particular, fibers at angles lower than 30◦ are under a negative strain and zero stress and the Poisson’s ratio results to be equal to that of the matrix. For angles higher than 30◦ , tension in the fibers begins to grow, thus increasing the normalized transverse stress (Fig. 6). Since the Young’s modulus of the matrix is significantly lower than that of the fibers, the matrix needs to assume a high strain in order to equilibrate the transverse component of fibril stress, thus increasing the lateral displacement and the Poisson’s ratio of the tissue. On the other hand, fibers aligned near the axial direction (∼80–90◦ ) take an important stress, but have a low component in the transverse direction, which is not capable to significantly alter the Poisson’s ratio of the tissue which, under these conditions, is similar to that of the matrix. This finding was previously obtained with models using orthogonal fibers [4,12,13]. Reduction of the Poisson’s ratio with strain for fibers at 60◦ may be initially attributed to the change of orientation of the fibers which may lead to a reduction of the normalized transverse stress. However, results showed that the transverse stress in the fibers, and therefore, the theoretical Poisson’s ratio, increased with strain (Fig. 8a). The reason for this disagreement appears to be the significant change of dimensions with strain, more marked for the model with the fibers at 60◦ , which displays an important reduction of the transverse dimension. When the change of geometry was taken into account using the logarithmic formula for strain, a slight increase of the Poisson’s ratio was appreciated and a better correlation with the theoretical Poisson’s ratio (Fig. 8b). Moreover, the logarithmic Poisson’s ratio also showed a good correlation with the theoretical Poisson’s ratio for the whole range of fiber orientations (Fig. 6), which supports the conclusion that high Poisson’s ratios are only generated when there is an important transverse stress component and the stiffness of the matrix is much lower than that of the fibers. From the foregoing results it appears that non-orthogonal fibrils should be considered into fibril-reinforced models for a proper description of the mechanical response of articular cartilage. One good candidate is the model proposed by
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Wilson et al. [8] that includes sets of secondary fibrils at 45◦ with respect to the fibrils in orthogonal directions. The model used here, with fibrils at 45◦ , was capable to give a fairly good approximation of most of the experimental Poisson’s ratios. In particular, considering site differences of the experimental stiffness of the matrix and the fibrils, the model was able to reproduce (Fig. 9) significant differences obtained in the study by Huang et al. [1], who document that the Poisson’s ratio for the humeral head is significantly greater than that of the glenoid (p = 0.01) and that the Poisson’s ratio from the superficial zone tended to be greater (p = 0.051) than that from the middle zone. However, this model was unable, even considering large values of the ratio between the stiffness of the fiber and that of the matrix (Fig. 9), to describe Poisson’s ratios near 2, that have been documented in tension tests of articular cartilage. When necessary, refinements should be considered in this model, like the addition of sets of fibrils in directions higher than 45◦ , for a more accurate description of experimental data. The variation of fiber properties with position can also be included in a complete model of articular cartilage in order to properly describe differences with site, like those expected for the deep zone, where the fibrils are perpendicular to the subchondral bone, suggesting that a low in-plane Poisson’s ratio may be found for this location. For all aspect ratios, constrains in the grips made the apparent Poisson’s ratio calculated with the lateral displacement at the center of the specimen rather larger than the true Poisson’s ratio obtained under uniaxial loading. For the lowest aspect ratio, a better estimation of the Poisson’s ratio was accomplished using the lateral displacement at the center only, compared with the Poisson’s ratio calculated with the lateral displacement at two points. This is due to the large variations with position of the lateral displacement of the specimen, significantly affected in the whole length by the grips. Even though it is clear that the models studied here, reinforced with four sets of fibrils, are only a rough description of the real architecture of articular cartilage, which may contain fibrils in different directions [20], it can still be concluded from the analyses of curves like those of Fig. 7 that even with small aspect ratios of the order of two the Poisson’s ratio can be measured with errors less than 10% using the lateral displacement at the center of the specimen. The models used in this study represent a crude approximation of the microstructure of articular cartilage, also characterized by the presence of chondrocytes and fibril capsules around cells [20]. However, these models with a finite set of fibrils are a good compromise between computational efficiency and accuracy to describe main experimental findings. Other mechanical implications arising from the inclusion of oblique fibers must be analyzed in future studies, like the stiffening effect when the tissue is subjected to bidirectional loading, which may be of interest to understand the transient mechanical response of the tissue under unconfined compression [2].
597
Fig. A.1. Free body diagram of one element of the matrix representing the biaxial stress field when the tissue is under uniaxial loading. No shear stresses are generated given the symmetric configurations for the oblique fibers.
Acknowledgments Thanks are due to Universidad del Valle, Cali, Colombia for the time to undertake this study. The author gratefully acknowledges Dr. Carlos Coronado for revising the manuscript.
Appendix A When the whole tissue is under uniaxial stress, the matrix is under a biaxial stress field composed of the axial stress σ my and the transverse stress components of the fibrils, σ fx (Fig. A.1). Thus, using the Hooke’s law, the strain components are σfx νm0 σmy εx = − − (A.1) Em0 Em0 εy =
σmy νm0 σfx + Em0 Em0
(A.2)
where Em0 and νm0 are the Young’s modulus and the Poisson’s ratio of the matrix at zero deformation. By substituting the axial stress from Eq. (A.2) into Eq. (A.1), the transverse strain can be found to be εx = −νm0 εy −
2 ) σfx (1 − νm0 . Em0
(A.3)
Hence, the Poisson’s ratio of the tissue at zero deformation, ν0 , can be calculated as νm0 σfx −εx σfx (A.4) = + νm0 1 − ν0 = εy Em0 εy Em0 εy where σfx /(Em0 εy ) is the normalized transverse stress. From Eq. (A.4) two limiting cases can be considered. If the normalized transverse stress is low compared to 1, the Poisson’s ratio is σfx σfx 1 ⇒ ν0 ≈ + νm0 . (A.5) Em0 εy Em0 εy On the other hand, if the normalized transverse stress is high, the axial strain from Eq. (A.2) is approximately equal
598
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
to νm0 σfx εy ≈ Em0
(A.6)
and the Poisson’s ratio of the tissue becomes approximately equal to the normalized transverse stress, i.e. σfx ν0 ≈ . (A.7) Em0 εy Conflict of interest The author of this study does not have any financial or personal relationships with other people or organizations, which could result in an inappropriate influence of this study.
References [1] Huang C, Stankiewicz A, Ateshian GA, Mow VC. Anisotropy, inhomogeneity, and tension–compression nonlinearity of human glenohumeral cartilage in finite deformation. J Biomech 2005;38:799– 809. [2] Charlebois M, McKee MD, Buschmann MD. Nonlinear tensile properties of bovine articular cartilage and their variation with age and depth. J Biomech Eng 2004;126:129–37. [3] Mak AF. The apparent viscoelastic behavior of articular cartilage—the contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J Biomech Eng 1986;108:123–9. [4] Li LP, Herzog W, Korhonen RK, Jurvelin JS. The role of viscoelasticity of collagen fibers in articular cartilage: axial tension versus compression. Med Eng Phys 2005;27:51–7. [5] Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments. J Biomech Eng 1980;102:73–84. [6] Elliot DM, Narmoneva DA, Setton LA. Direct measurement of the Poisson’s ratio of human patella cartilage in tension. J Biomech Eng 2002;124:223–8.
[7] Li LP, Soulhat J, Buschmann MD, Shirazi-Adl A. Nonlinear analysis of cartilage in unconfined ramp compression using a fibril reinforced poroelastic model. Clin Biomech 1999;14:673–82. [8] Wilson W, van Donkelaar CC, van Rietbergen B, Ito K, Huiskes R. Stresses in the local collagen network of articular cartilage: a poroviscoelastic fibril-reinforced finite element study. J Biomech 2004;37:357–66. [9] Federico S, Grillo A, La Rosa G, Giaquinta G, Herzog W. A transversely isotropic, transversely homogeneous microstructural–statistical model of articular cartilage. J Biomech 2005;38:2008–18. [10] Klisch SM, Lotz JC. Application of a fiber-reinforced continuum theory to multiple deformations of the annulus fibrous. J Biomech 1999;32:1027–36. [11] Wu JZ, Herzog W. Elastic anisotropy of articular cartilage is associated with the microstructures of collagen fibers and chondrocytes. J Biomech 2002;35:931–42. [12] Klisch SM, Sah RL, Davol A. Bimodular-orthotropic-policonvex strain energy functions for the collagen proteoglycan solid matrix of articular cartilage. In: Proceedings of BIO2006, 2006: 157425. [13] Garc´ıa JJ, Cort´es DH. A biphasic viscohyperelastic fibril-reinforced model for articular cartilage: formulation and comparison with experimental data. J Biomech 2007;40:1737–44. [14] Lei F, Szeri AZ. The influence of fibril organization on the mechanical behavior of articular cartilage. Proc R Soc 2006;462:3301– 22. [15] Holmes MH, Mow VC. The non-linear characteristics of soft gels and hydrated connective tissues in ultrafiltration. J Biomech 1990;23:1145–56. [16] Limbert G, Middleton J. A transversely isotropic viscohyperelastic material: application to the modeling of biological soft connective tissues. Int J Solids Struct 2004;41:4237–60. [17] Garc´ıa JJ, Cort´es DH. A nonlinear biphasic viscohyperelastic model for articular cartilage. J Biomech 2006;29:2991–8. [18] Korhonen RK, Laasanen MS, T¨oyr¨as J, Rieppo J, Hirvonen J, Helminen HJ, et al. Comparison of the equilibrium response of articular cartilage in unconfined compression, confined compression and indentation. J Biomech 2002;35:903–9. [19] Christensen RM. Mechanics of composite materials. Malabar: Krieger Publishing Company; 1991. [20] Clark JM. The organization of collagen in cryofractured rabbit articular cartilage: a scanning electron microscopic study. J Orthop Res 1985;3:17–29.
Medical Engineering & Physics 30 (2008) 590–598
Simulation of high tensile Poisson’s ratios of articular cartilage with a finite element fibril-reinforced hyperelastic model Jos´e Jaime Garc´ıa ∗ Escuela de Ingenier´ıa Civil y Geom´atica, Universidad del Valle, Calle 13, Carrera 100, Edificio 350, Cali, Colombia Received 7 February 2007; received in revised form 26 June 2007; accepted 27 June 2007
Abstract Analyses with a finite element fibril-reinforced hyperelastic model were undertaken in this study to simulate high tensile Poisson’s ratios that have been consistently documented in experimental studies of articular cartilage. The solid phase was represented by an isotropic matrix reinforced with four sets of fibrils, two of them aligned in orthogonal directions and two oblique fibrils in a symmetric configuration respect to the orthogonal axes. Two distinct hyperelastic functions were used to represent the matrix and the fibrils. Results of the analyses showed that only by considering non-orthogonal fibrils was it possible to represent Poisson’s ratios higher than one. Constrains in the grips and finite deformations played a minor role in the calculated Poisson’s ratio. This study also showed that the model with oblique fibrils at 45◦ was able to represent significant differences in Poisson’s ratios near 1 documented in experimental studies. However, even considering constrains in the grips, this model was not capable to simulate Poisson’s ratios near 2 that have been reported in other studies. The study also confirmed that only with a high relation between the stiffness of the fibers and that of the matrix was it possible to obtain high Poisson’s ratios for the tissue. Results suggest that analytical models with a finite number of fibrils are appropriate to represent main mechanical effects of articular cartilage. © 2007 IPEM. Published by Elsevier Ltd. All rights reserved. Keywords: Anisotropy of cartilage; Tensile Poisson’s ratio; Hyperelastic fibril-reinforced models
1. Introduction Articular cartilage is a highly heterogeneous and anisotropic material covering the end of the bones inside diarthrodial joints. Towards an understanding of the relations between mechanical parameters and osteoarthrosis, many analytical and experimental studies have been developed to determine the mechanical response of articular cartilage. This response is mainly represented by non-linear stress–strain curves under transient loads and equilibrium [1,2], nonlinear tension–compression response [1], intrinsic viscous effects of the solid phase [3,4], and non-linear permeability [5]. Anisotropy of articular cartilage has been further demonstrated in various experimental studies exhibiting tensile Poisson’s ratios higher than one [1,2,6]. A finite element model able of describing most of the independent experiments is a valuable tool to accurately assess the mechanical ∗
Tel.: +57 2 339 2420; fax: +57 2 339 2086. E-mail address: [email protected].
performance of articular cartilage, needed in many practical applications. Within the framework of the biphasic theory [5] and among various approaches that have been proposed to represent the mechanical behavior of the solid phase, fibrilreinforced models have been successful of representing main experimental observations as well as the microstructure of the tissue [7,8,9,10]. Only few of these models have documented Poisson’s ratios under tension. With the assumption of infinitesimal strains Wu and Herzog [11] use homogenization theory to determine the effective elastic constants of a mixture of matrix, orthogonal fibrils and chondrocytes, documenting Poisson’s ratios from the stable models lower than 0.4, which are equally valid for tension and compression since their model does not consider the tension–compression non-linearity. A recent study by Klisch et al. [12] considering finite deformations and a bimodular-orthotropic function documents Poisson’s ratios under tension one order of magnitude lower than those found experimentally. A couple of
1350-4533/$ – see front matter © 2007 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2007.06.004
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
models [4,13] that consider the tension–compression nonlinearity by allowing only tension in the fibrils have been incapable to describe high tensile Poisson’s ratios. One possible explanation for this inconsistency, as suggested by Li et al [4], is that, due to the constrains in the grips, the stress field under a tension test is really two-dimensional, and that the Poisson’s ratio depends on the aspect ratio of the specimen used in the experimental set up. To better represent the architecture of the tissue, Wilson et al. [8] incorporate sets of secondary fibrils at 45◦ but no Poisson’s ratios under tension have been documented with this model. A recent analysis by Ley and Szeri [14] with a fibril-reinforced model under infinitesimal deformations and true uniaxial loading shows that high tensile Poisson’s ratios may be simulated considering non-orthogonal fibril distributions. It has not yet been investigated the capability of models with a finite number of fibrils to reproduce high tensile Poisson’s ratios. Moreover, no study has been conducted to determine the influence of constrains in the grips and finite deformations in the measurement of the tensile Poisson’s ratio of articular cartilage. A finite element fibril-reinforced hyperelastic model was used in this study to simulate high tensile Poisson’s ratios documented in experimental studies by Huang et al. [1], Charlebois et al. [2] and Elliot et al. [6]. The study was specially aimed to quantify the effects in the Poisson’s ratio of constrains in the grips, alignment of the fibrils and finite deformations.
2. Methods For the equilibrium simulation treated here no permeability or intrinsic viscous parameters of our general finite element model [13] were required. In addition, the swelling pressure was not included considering that its main effect is to increase the stiffness of the non-fibrillar matrix, which is
591
still low compared to the stiffness of the fibers. Thus, the solid phase of articular cartilage was represented by the proteoglycan matrix, described by the isotropic hyperelastic function proposed by Holmes and Mow [15], reinforced with four sets of collagen fibers, described by the hyperelastic function proposed by Limbert and Middlelton [16]. The hyperelastic function of the matrix is fully defined by three elastic parameters α0 , α1 , α2 that can be adjusted to fit the Young’s modulus and Poisson’s ratio of the matrix at zero deformation and the degree of non-linearity of the stress–strain curve under compression [17]. The hyperelastic function of the fibers is defined by the Young’s modulus at zero deformation Ef0 and a parameter γ that can be used to adjust the degree of non-linearity of the stress–strain curve under tension. One objective of this study was to determine the capability of a model with a finite number of fibrils to explain experimental results. A model like this may be efficiently implemented in a 3-D geometry. Hence, similar to the secondary set of fibrils used by Wilson et al. [8], four sets of fibers were considered, one in the axial (y) and one in the transverse (x) direction of loading (orthogonal fibers), and two oblique fibrils in a symmetric configuration respect to the orthogonal axes (Fig. 1). The angle of the oblique fibers respect to the transverse axis was changed to see the influence in the Poisson’s ratio of the tissue. Simulations were conducted of the tensile experiments undertaken by Huang et al. [1] with humeral and glenoid human cartilage and Charlebois et al. [2] and Elliot et al. [6] with bovine and human patellar cartilage, respectively. The finite element models consisted of four-node plane-stress elements distributed in a domain equal to a fourth of the experimental samples, according to the divisions shown in Table 1, with symmetry conditions prescribed in the orthogonal axes (Fig. 1). Vertical displacements were prescribed in the upper nodes of the model to complete 12, 7.7 and 16% axial deformation, for the studies of Huang et al. [1],
Fig. 1. Dimensions of the specimens, meshed region, orientation of the oblique fibrils and mesh for the data by Charlebois et al. [2].
592
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Table 1 Geometry and mesh for each experimental study
Charlebois et al. [2] Huang et al. [1] Elliot et al. [6]
Dimensions between grips (mm) length (L) × width (w)
Aspect ratio (L/w)
Number of divisions along each direction
Axial deformation (%)
4×2 5 × 1.5 8.33 × 1.91
2 3.33 4.36
8×4 10 × 3 16 × 4
7.7 12 16
Charlebois et al. [2] and Elliot et al. [6], respectively. In the updated Lagrangian finite element code the fibers were taken into account by calculating at each integration point the residual vector according to the hyperelastic function by Limbert and Middleton [16]. To improve convergence of the iterative scheme, the tangent stiffness matrix for the fibers was calculated and added to the tangent stiffness of the matrix. In order to study the influence of the grips in the Poisson’s ratio, analyses were accomplished considering the transverse constrains in the grips (constrained model) and without constrains in the grips (unconstrained model). The unconstrained model was intended to reproduce a true uniaxial loading condition. Additionally, the Poisson’s ratio for the constrained model (apparent Poisson’s ratio) was calculated with the lateral displacement at the center and with the average of the lateral displacement at the center and at the quarter length of the model, similar to Huang et al. [1]. The Poisson’s ratio was calculated as the negative ratio of transverse strain over axial strain. For the analysis under finite deformations and in order to take into account important variations of geometry respect to the initial configuration, another Poisson’s ratio was calculated, called here the logarithmic Poisson’s ratio, defined as the negative ratio of the transverse logarithmic strain over the axial logarithmic strain. A first analysis was conducted using the Hooke’s law to develop a simple formula (presented in Appendix A and called here the theoretical formula) for calculating the Poisson’s ratio of the whole tissue in terms of the Poisson’s ratio of the matrix and the stress component of the fibers in the transverse direction. The theoretical formula was used with the transverse stress component obtained from finite element analysis to predict a value, called here the theoretical Poisson’s ratio. The transverse stress component of the fibers was normalized (normalized transverse stress) by dividing it over the product of the axial strain and the Young’s modulus of the matrix. The theoretical formula and preliminary finite element results (Fig. 2) showed that the tensile Poisson’s ratio of the tissue was independent of the Poisson’s ratio of the
Fig. 2. Apparent Poisson’s ratio of the humeral surface for two values of the Poisson’s ratio of the matrix and different orientations of the oblique fibers.
matrix for a wide range of fibril orientations. Thus, a value of 1/3 was adopted in this study for the Poisson’s ratio of the matrix, which appears to be more consistent with Poisson’s ratios documented in compression experiments [18]. Hence, following a procedure described by Garc´ıa and Cort´es [13] various sets of material properties were calculated for the matrix and the fibrils using the data reported by Huang et al. [1] for the glenoid and humeral head. From the tensile experiments documented by Charlebois et al. [2] and Elliot et al. [6] it is not possible to determine the elastic constants of the matrix. Thus, the same properties for the matrix used for the data by Huang et al. [1] for humeral cartilage were initially adopted in the analyses of the data by Charlebois et al. [2] and Elliot et al. [6]. Next, other analyses were conducted considering a Young’s modulus for the matrix equal to half and a third that obtained from the study by Huang et al. [1]. For each orientation of the oblique fibers, fibril properties were adjusted to yield the reported experimental Young’s
Table 2 Elastic properties for the matrix used in the simulations of the tensile experiments Assumed νm0 H-s H-m G-s G-m
1/3
Reported HA0 (MPa)
Calculated Em0 (MPa) Em0 = 1.5HA0
Reported n
0.116 0.141 0.138 0.178
0.174 0.212 0.207 0.267
1.15 0.81 1.16 0.88
Calculated α1 = 0 0
Calculated α2 = n/2
Calculated α0 (MPa) α0 = 3Em0 /16α2
0.58 0.41 0.58 0.44
0.0563 0.097 0.0669 0.114
These properties were calculated from the confined compression data reported by Huang et al. [5]. H-s: humeral surface, H-m: humeral middle, G-s: glenoid surface, G-m: glenoid middle.
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
593
modulus E0+ in the axial direction using the well-known theory of composite materials [19] under uniaxial loading and infinitesimal deformations. The Poisson’s ratio obtained from composite materials theory was also used to verify results of the finite element code under infinitesimal deformations and true uniaxial loading. For the calculation of the non-linear stiffness of the material, the parameters A and B of the formula σ = A exp(Bε) used by Huang et al. [1] and Elliot et al. [6] at a relaxation time of 30 min were adjusted with the hyperelastic function for the fibers as explained by Garc´ıa and Cort´es [13]. Since Charlebois et al. [2] report the lateral displacement profile at a relaxation time (30 min) different from the times used to report stiffness (full relaxation and 15 min), the same fibril properties by Elliot et al. [6] were used in the simulation of the data by Charlebois et al. [2]. A summary of material properties is presented in Tables 2 and 3.
3. Results Comparison of lateral displacement profiles (Fig. 3) obtained with various angles for the oblique fibers with those obtained experimentally [2] indicates that oblique fibers at about 60◦ respect to the transverse axis, i.e. 30◦ respect to the axis of loading, give a better description of the lateral profile (Fig. 4). For oblique fibers at 45◦ and for orthogonal fibrils the variation of the lateral profile with the Young’s moduli of the constituents was negligible (Fig. 3). For oblique fibrils at 60◦ , increments in the maximum lateral displacement (and thus of the Poisson’s ratio) of 18 and 30% were obtained for 2- and 3-fold increments, respectively, of the ratio of Young’s moduli of the constituents. The finite element model yielded a good approximation of the infinitesimal-strain closed-form solution under 1% deformation (Fig. 5), except for angles of the oblique fibrils less than 30◦ , where the closed-form solution from composite materials, which does not consider the tension–compression non-linearity, predicts lower values of the Poisson’s ratio. The theoretical formula (Eq. (A.4)) predicts high variations of the Poisson’s ratio with the normalized transverse stress. Mainly, for low values of the normalized transverse stress, the Poisson’s ratio is approximately equal to the Poisson’s ratio of the matrix. On the other hand, for higher values of the normalized transverse stress, the Poisson’s ratio is approximately equal to the normalized transverse stress. Foregoing results are consistent with those obtained with the finite element program (Fig. 2) showing that the Poisson’s ratio of the tissue is approximately equal to that of the matrix for low (<25◦ ) and high (>85◦ ) angles of the oblique fibers and tends to be higher than one for a range of middle angles of the oblique fibers (∼45–75◦ ). Moreover, the logarithmic Poisson’s ratio showed a good correlation with the Poisson’s ratio obtained with the theoretical formula (Fig. 6). Regardless of the effects of the grips and aspect ratio of the
Fig. 3. Profiles of lateral displacement for three angles of the oblique fibrils and different relations between the Young’s modulus of the fibers and that of the matrix. Experimental data by Charlebois et al. [2]. (a) Orthogonal fibers; (b) oblique fibers at 45◦ ; (c) oblique fibers at 60◦ .
594
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Table 3 Tensile stiffness used in the simulations Reported A (MPa)
Reported B
E0+ (MPa)
Non-linear parameter γ
Huang et al. [1], human
H-s H-m G-s G-m
0.652 0.422 0.562 0.296
10.27 8.80 8.49 7.98
6.57 3.55 4.55 2.24
4.6 3.4 3.2 2.8
Elliot et al. [6], human
P-s P-m
0.411a 0.381a
14.15a 5.05a
5.81a 1.92a
8.1 0.5
Charlebois et al. [2], bovine
P-s
0.411a
14.15a
5.81a
8.1
P-s: patellar surface, P-m: patellar middle. a Calculated from the reported data.
specimens (Fig. 7), Poisson’s ratios higher than one could only be obtained considering fiber orientations in the range of ∼45–80◦ , which is consistent with the range of angles displaying higher values of the transverse stress component and with the values predicted with the theoretical formula (Fig. 6). Even specimens with the lowest aspect ratio (2.0) yielded an apparent Poisson’s ratio with maximum differences of about 10% (Fig. 7a) respect to the value predicted by the unconstrained model. For the lowest aspect ratio the appar-
Fig. 5. Poisson’s ratios obtained with the closed-form solution from composite material theory and the finite element program under 1% deformation. This simulation was accomplished for a Poisson’s ratio of the matrix equal to 0.33.
ent Poisson’s ratio obtained with the lateral displacement at the center was better than that obtained with two points, while the opposite trend was observed for the higher aspect ratios (Fig. 7). For a range of angles for the oblique fibers
Fig. 4. Deformed mesh and contours of transverse displacement for the samples of Charlebois et al. [2] with the oblique fibers at 60◦ with respect to the transverse axis. The whole sample was considered in this simulation.
Fig. 6. Comparison of the Poisson’s ratio, the logarithmic Poisson’s ratio, the theoretical Poisson’s ratio and the normalized transverse stress from analyses undertaken with the data by Huang et al. [1] for the humeral surface.
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
595
Fig. 7. Poisson’s ratio calculated with and without constrains in the grips for different values of the aspect ratio of the specimen. With constrain in the grips the Poisson’s ratio was calculated with the average of the lateral displacement at two points and with the lateral displacement at the center (one point). (a) Charlebois et al. [2] L/w = 2, E0+ /Em0 = 33.4. (b) Huang et al [1] L/w = 3.33, E0+ /Em0 = 16.7. (c) Elliot et al. [6] L/w = 4.36, E0+ /Em0 = 11.0.
(∼50–80◦ ) a positive correlation was observed between the apparent Poisson’s ratio and the ratio of the Young’s moduli of the constituents (Fig. 7). A reduction of Poisson’s ratio with strain was observed for fiber orientations of 60◦ and no significant changes were obtained for fiber orientations of 30◦ and 45◦ (Fig. 8a). However, when using the logarithmic formula, slight increases can be observed for fiber orientations of 60◦ and 45◦ , and a good correlation with the Poisson’s ratio calculated with the theoretical formula (Fig. 8b). For the data by Charlebois et al. [2] and the patella surface by Elliot et al. [6] the experimental Poisson’s ratios were about 2-fold those obtained with the finite element model with fibers at 45◦ (Fig. 9). For the other experimental data the
finite element model with the fibers at 45◦ provided a fairly good approximation of the experimental Poisson’s ratio, with percent differences in the range of 8–20% (Fig. 9).
4. Discussion This study showed the ability of fibril-reinforced models under finite deformations to describe high Poisson’s ratios reported in tensile tests of articular cartilage. The result was feasible given the anisotropy of the tissue, mainly characterized by a high relation between the stiffness of the fibrils and that of the matrix. Only considering oblique fibrils with orientations of the order 45–80◦ respect to the transverse axis
596
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Fig. 8. Variations with strain of Poisson’s ratio (P) and the logarithmic (L) Poisson’s ratio for the data by Huang et al. [1] for the humeral surface, and comparison with the theoretical formula (T) documented in Appendix A, for three orientations of the oblique fibers (60◦ , 45◦ and 30◦ ). (a) Poisson’s ratio (P) calculated as the negative ratio of transverse to axial strain. (b) Poisson’s ratio (L) calculated with the logarithmic formula.
Fig. 9. Comparison of the experimental and numerical Poisson’s ratios with the oblique fibers at 45◦ .
of loading was it possible to simulate Poisson’s ratios higher than one, which is in agreement with recent analyses of a reinforced model employing fibril distribution functions under infinitesimal deformations and uniaxial loading [14]. Unlike the use of fibril distributions, the consideration of a limited number of fibrils allows a relatively easy implementation into a non-linear finite element code. The need to consider oblique fibrils for a proper description of the high tensile Poisson’s ratio suggests that articular cartilage contains, even for directions parallel to the split lines, fibril distributions aligned in non-orthogonal directions, which appears to be consistent with the findings by Clark [20] who documents fibrils in the radial zone of cartilage not perfectly straight or parallel. The normalized transverse stress showed to be a good predictor of the Poisson’s ratio of the tissue. In particular, fibers at angles lower than 30◦ are under a negative strain and zero stress and the Poisson’s ratio results to be equal to that of the matrix. For angles higher than 30◦ , tension in the fibers begins to grow, thus increasing the normalized transverse stress (Fig. 6). Since the Young’s modulus of the matrix is significantly lower than that of the fibers, the matrix needs to assume a high strain in order to equilibrate the transverse component of fibril stress, thus increasing the lateral displacement and the Poisson’s ratio of the tissue. On the other hand, fibers aligned near the axial direction (∼80–90◦ ) take an important stress, but have a low component in the transverse direction, which is not capable to significantly alter the Poisson’s ratio of the tissue which, under these conditions, is similar to that of the matrix. This finding was previously obtained with models using orthogonal fibers [4,12,13]. Reduction of the Poisson’s ratio with strain for fibers at 60◦ may be initially attributed to the change of orientation of the fibers which may lead to a reduction of the normalized transverse stress. However, results showed that the transverse stress in the fibers, and therefore, the theoretical Poisson’s ratio, increased with strain (Fig. 8a). The reason for this disagreement appears to be the significant change of dimensions with strain, more marked for the model with the fibers at 60◦ , which displays an important reduction of the transverse dimension. When the change of geometry was taken into account using the logarithmic formula for strain, a slight increase of the Poisson’s ratio was appreciated and a better correlation with the theoretical Poisson’s ratio (Fig. 8b). Moreover, the logarithmic Poisson’s ratio also showed a good correlation with the theoretical Poisson’s ratio for the whole range of fiber orientations (Fig. 6), which supports the conclusion that high Poisson’s ratios are only generated when there is an important transverse stress component and the stiffness of the matrix is much lower than that of the fibers. From the foregoing results it appears that non-orthogonal fibrils should be considered into fibril-reinforced models for a proper description of the mechanical response of articular cartilage. One good candidate is the model proposed by
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
Wilson et al. [8] that includes sets of secondary fibrils at 45◦ with respect to the fibrils in orthogonal directions. The model used here, with fibrils at 45◦ , was capable to give a fairly good approximation of most of the experimental Poisson’s ratios. In particular, considering site differences of the experimental stiffness of the matrix and the fibrils, the model was able to reproduce (Fig. 9) significant differences obtained in the study by Huang et al. [1], who document that the Poisson’s ratio for the humeral head is significantly greater than that of the glenoid (p = 0.01) and that the Poisson’s ratio from the superficial zone tended to be greater (p = 0.051) than that from the middle zone. However, this model was unable, even considering large values of the ratio between the stiffness of the fiber and that of the matrix (Fig. 9), to describe Poisson’s ratios near 2, that have been documented in tension tests of articular cartilage. When necessary, refinements should be considered in this model, like the addition of sets of fibrils in directions higher than 45◦ , for a more accurate description of experimental data. The variation of fiber properties with position can also be included in a complete model of articular cartilage in order to properly describe differences with site, like those expected for the deep zone, where the fibrils are perpendicular to the subchondral bone, suggesting that a low in-plane Poisson’s ratio may be found for this location. For all aspect ratios, constrains in the grips made the apparent Poisson’s ratio calculated with the lateral displacement at the center of the specimen rather larger than the true Poisson’s ratio obtained under uniaxial loading. For the lowest aspect ratio, a better estimation of the Poisson’s ratio was accomplished using the lateral displacement at the center only, compared with the Poisson’s ratio calculated with the lateral displacement at two points. This is due to the large variations with position of the lateral displacement of the specimen, significantly affected in the whole length by the grips. Even though it is clear that the models studied here, reinforced with four sets of fibrils, are only a rough description of the real architecture of articular cartilage, which may contain fibrils in different directions [20], it can still be concluded from the analyses of curves like those of Fig. 7 that even with small aspect ratios of the order of two the Poisson’s ratio can be measured with errors less than 10% using the lateral displacement at the center of the specimen. The models used in this study represent a crude approximation of the microstructure of articular cartilage, also characterized by the presence of chondrocytes and fibril capsules around cells [20]. However, these models with a finite set of fibrils are a good compromise between computational efficiency and accuracy to describe main experimental findings. Other mechanical implications arising from the inclusion of oblique fibers must be analyzed in future studies, like the stiffening effect when the tissue is subjected to bidirectional loading, which may be of interest to understand the transient mechanical response of the tissue under unconfined compression [2].
597
Fig. A.1. Free body diagram of one element of the matrix representing the biaxial stress field when the tissue is under uniaxial loading. No shear stresses are generated given the symmetric configurations for the oblique fibers.
Acknowledgments Thanks are due to Universidad del Valle, Cali, Colombia for the time to undertake this study. The author gratefully acknowledges Dr. Carlos Coronado for revising the manuscript.
Appendix A When the whole tissue is under uniaxial stress, the matrix is under a biaxial stress field composed of the axial stress σ my and the transverse stress components of the fibrils, σ fx (Fig. A.1). Thus, using the Hooke’s law, the strain components are σfx νm0 σmy εx = − − (A.1) Em0 Em0 εy =
σmy νm0 σfx + Em0 Em0
(A.2)
where Em0 and νm0 are the Young’s modulus and the Poisson’s ratio of the matrix at zero deformation. By substituting the axial stress from Eq. (A.2) into Eq. (A.1), the transverse strain can be found to be εx = −νm0 εy −
2 ) σfx (1 − νm0 . Em0
(A.3)
Hence, the Poisson’s ratio of the tissue at zero deformation, ν0 , can be calculated as νm0 σfx −εx σfx (A.4) = + νm0 1 − ν0 = εy Em0 εy Em0 εy where σfx /(Em0 εy ) is the normalized transverse stress. From Eq. (A.4) two limiting cases can be considered. If the normalized transverse stress is low compared to 1, the Poisson’s ratio is σfx σfx 1 ⇒ ν0 ≈ + νm0 . (A.5) Em0 εy Em0 εy On the other hand, if the normalized transverse stress is high, the axial strain from Eq. (A.2) is approximately equal
598
J.J. Garc´ıa / Medical Engineering & Physics 30 (2008) 590–598
to νm0 σfx εy ≈ Em0
(A.6)
and the Poisson’s ratio of the tissue becomes approximately equal to the normalized transverse stress, i.e. σfx ν0 ≈ . (A.7) Em0 εy Conflict of interest The author of this study does not have any financial or personal relationships with other people or organizations, which could result in an inappropriate influence of this study.
References [1] Huang C, Stankiewicz A, Ateshian GA, Mow VC. Anisotropy, inhomogeneity, and tension–compression nonlinearity of human glenohumeral cartilage in finite deformation. J Biomech 2005;38:799– 809. [2] Charlebois M, McKee MD, Buschmann MD. Nonlinear tensile properties of bovine articular cartilage and their variation with age and depth. J Biomech Eng 2004;126:129–37. [3] Mak AF. The apparent viscoelastic behavior of articular cartilage—the contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J Biomech Eng 1986;108:123–9. [4] Li LP, Herzog W, Korhonen RK, Jurvelin JS. The role of viscoelasticity of collagen fibers in articular cartilage: axial tension versus compression. Med Eng Phys 2005;27:51–7. [5] Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments. J Biomech Eng 1980;102:73–84. [6] Elliot DM, Narmoneva DA, Setton LA. Direct measurement of the Poisson’s ratio of human patella cartilage in tension. J Biomech Eng 2002;124:223–8.
[7] Li LP, Soulhat J, Buschmann MD, Shirazi-Adl A. Nonlinear analysis of cartilage in unconfined ramp compression using a fibril reinforced poroelastic model. Clin Biomech 1999;14:673–82. [8] Wilson W, van Donkelaar CC, van Rietbergen B, Ito K, Huiskes R. Stresses in the local collagen network of articular cartilage: a poroviscoelastic fibril-reinforced finite element study. J Biomech 2004;37:357–66. [9] Federico S, Grillo A, La Rosa G, Giaquinta G, Herzog W. A transversely isotropic, transversely homogeneous microstructural–statistical model of articular cartilage. J Biomech 2005;38:2008–18. [10] Klisch SM, Lotz JC. Application of a fiber-reinforced continuum theory to multiple deformations of the annulus fibrous. J Biomech 1999;32:1027–36. [11] Wu JZ, Herzog W. Elastic anisotropy of articular cartilage is associated with the microstructures of collagen fibers and chondrocytes. J Biomech 2002;35:931–42. [12] Klisch SM, Sah RL, Davol A. Bimodular-orthotropic-policonvex strain energy functions for the collagen proteoglycan solid matrix of articular cartilage. In: Proceedings of BIO2006, 2006: 157425. [13] Garc´ıa JJ, Cort´es DH. A biphasic viscohyperelastic fibril-reinforced model for articular cartilage: formulation and comparison with experimental data. J Biomech 2007;40:1737–44. [14] Lei F, Szeri AZ. The influence of fibril organization on the mechanical behavior of articular cartilage. Proc R Soc 2006;462:3301– 22. [15] Holmes MH, Mow VC. The non-linear characteristics of soft gels and hydrated connective tissues in ultrafiltration. J Biomech 1990;23:1145–56. [16] Limbert G, Middleton J. A transversely isotropic viscohyperelastic material: application to the modeling of biological soft connective tissues. Int J Solids Struct 2004;41:4237–60. [17] Garc´ıa JJ, Cort´es DH. A nonlinear biphasic viscohyperelastic model for articular cartilage. J Biomech 2006;29:2991–8. [18] Korhonen RK, Laasanen MS, T¨oyr¨as J, Rieppo J, Hirvonen J, Helminen HJ, et al. Comparison of the equilibrium response of articular cartilage in unconfined compression, confined compression and indentation. J Biomech 2002;35:903–9. [19] Christensen RM. Mechanics of composite materials. Malabar: Krieger Publishing Company; 1991. [20] Clark JM. The organization of collagen in cryofractured rabbit articular cartilage: a scanning electron microscopic study. J Orthop Res 1985;3:17–29.