- Email: [email protected]
Constitutive formulation and analysis of heel pad tissues mechanics
Medical Engineering & Physics 32 (2010) 516–522
Contents lists available at ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Constitutive formulation and analysis of heel pad tissues mechanics A.N. Natali ∗ , C.G. Fontanella, E.L. Carniel University of Padova, Centre of Mechanics of Biological Materials, Via F. Marzolo 9, I-35131 Padova, Italy
a r t i c l e
i n f o
Article history: Received 28 May 2009 Received in revised form 11 February 2010 Accepted 18 February 2010 Keywords: Heel pad Tissue mechanics Constitutive model Visco-elasticity Mechanical tests
a b s t r a c t This paper presents a visco-hyperelastic constitutive model developed to describe the biomechanical response of heel pad tissues. The model takes into account the typical features of the mechanical response such as large displacement, strain phenomena, and non-linear elasticity together with time-dependent effects. The constitutive model was formulated, starting from the analysis of the complex structural and micro-structural configuration of the tissues, to evaluate the relationship between tissue histology and mechanical properties. To define the constitutive model, experimental data from mechanical tests were analyzed. To obtain information about the mechanical response of the tissue so that the constitutive parameters could be established, data from both in vitro and in vivo tests were investigated. Specifically, the first evaluation of the constitutive parameters was performed by a coupled deterministic and stochastic optimization method, accounting for data from in vitro tests. The comparison of constitutive model results and experimental data confirmed the model’s capability to describe the compression behaviour of the heel pad tissues, regarding both constant strain rate and stress relaxation tests. Based on the data from additional experimental tests, some of the constitutive parameters were modified in order to interpret the in vivo mechanical response of the heel pad tissues. This approach made it possible to interpret the actual mechanical function of the tissues. © 2010 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction The configuration of the foot has been investigated with regard to weight-bearing capacity during static standing and ambulation [1]. Specific loose connective tissues, here defined as plantar tissues, surround bony structures in the foot plantar region. Their conformation and mechanical properties strongly influence body posture during static standing and the gait cycle. Specifically, the mechanical role of plantar tissues is to perform load transmission phenomena, also to dampen shocks generated during the gait or running cycle [2], making it possible to achieve a smooth distribution of plantar pressure [3]. A specific knowledge of plantar tissues’ mechanical behaviour contributes to the investigation of the mechanical functionality of the foot, considering the interaction phenomena between the foot and surrounding structures, such as footwear products [4–6]. The mechanics of plantar tissues must be investigated by both an experimental and a computational approach. Experimental activities are necessary to acquire data about the tissues’ mechanical response. These data will represent the necessary input for the development and validation of computational models and will give an appropriate constitutive formulation for plantar tissues.
∗ Corresponding author. Tel.: +39 049 827 5598; fax: +39 049 827 5604. E-mail address: [email protected] (A.N. Natali).
Detailed information about the structural and micro-structural configuration of the tissues, together with data from mechanical tests, is necessary to define the proper constitutive formulation [7–9]. The plantar region is a composite material of fatty and connective tissues located between skin dermis, fascia and bony segments. Within the plantar region, particular attention is paid to the heel pad, which is the portion of the plantar region interposed between calcaneum and skin that plays a fundamental role in foot mechanics (Fig. 1a). The heel pad is mainly organized according to a honeycomb configuration [10,11]. Fat tissue chambers are embedded and separated from each other by connective septa (Fig. 1b and c). In different locations within the heel pad structure, chambers and septa are characterized by specific dimensions and orientations [12]. These aspects are heavily influenced by the specific local loading conditions. The mechanical properties of the heel pad were investigated using in vivo or ex vivo techniques. Different loading conditions were applied over the heel region and the corresponding displacement or strain was measured [6,13–16], leading to a comprehensive characterization of the mechanical response of the plantar region. The specific stress–strain behaviour of the tissue was analyzed using in vitro tests based on more simple geometry specimens, boundary conditions and loading protocols. Experimental data from in vitro tests can be found in the literature. Miller-Young et al. [17] developed compression tests on cylindrical specimens from
1350-4533/$ – see front matter © 2010 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2010.02.018
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
517
a formulation capable of interpreting the actual response of heel pad living tissues, constitutive parameters were subsequently rearranged accounting for in vivo behaviour. This task was approached by analyzing data from in vivo tests, paying special attention to the relationship between tissue stiffness and strain rate. This relationship was analyzed by the proposed constitutive formulation in order to properly modify the specific constitutive parameters. 2. Materials and methods 2.1. Morphometry and histology The heel pad is the portion of plantar tissues interposed between the calcaneum and the skin (Fig. 1a). The thickness ranges between 14.4 and 24.5 mm with an average value of 18 mm [26]. The distribution of overall heel pad thickness varies according to local loading conditions. For example, the highest thickness values are found where severe interactions develop, such as loading by impacts [27]. The structure of the heel pad is complex and hierarchically organized [28–30]. As reported in Fig. 1b, the heel pad mainly consists of adipose tissue chambers, separated by connective tissue septa [27]. The adipose chambers contain closely packed fat cells with diameters ranging between 100 and 200 m [15] and arranged in lobules within a framework of connective bundles. Larger groups of lobules are enclosed by further bundles of collagen and elastin fibres, which thus form a supporting network (Fig. 1c). This honeycomb configuration of septa entails an overall closed-cells structure [10]. Adipose chambers dimension and orientation depend on the specific location within the heel pad (Fig. 1a and b). Vertically oriented chambers are located in the central portion of the heel pad. The chambers are small in the region close to dermis, progressively increasing towards the calcaneal perichondrium. In the lateral and posterior regions, chambers are small and more transversally oriented [27]. The connective septa configuration varies according to the geometries of the adjacent adipose chambers. Collagen fibres within septa are spirally arranged and superiorly fixed to bone or other septa, inferiorly to other septa or dermis [30]. Compressive stiffness is caused by the almost incompressible behaviour of fat tissue bounded by collagen network [31]. The spiral distribution of collagen fibres within the fatty chambers and connective septa is suitable to limit the bulging of chambers themselves [27]. Fig. 1. Morphometric and histological configuration of the heel pad: magnetic resonance image according to sagittal section (a), longitudinal section by a schematic representation (b) and optical microscope picture of a transversal section (c), from Weissengruber et al. [30].
the human heel pad under different loading conditions and evaluated stress relaxation phenomena that were also studied by Ledoux and Blevins [18] on similar geometry specimens. Different constitutive formulations are reported in literature for heel pad tissues, with particular attention to the hyperelastic and visco-elastic response [18–23]. The main goal of the proposed work was to define the constitutive model and evaluate the associated parameters, aiming at a reliable and efficient characterization of the mechanical behaviour of heel pad tissues. A specific visco-hyperelastic constitutive model was formulated. The proposed model was capable of accounting for the typical response of heel pad tissues, such as large displacements and strains, an almost incompressible behaviour and time-dependent effects. The constitutive parameters were determined by a specific procedure, consisting of optimization methods capable of minimizing the discrepancy between experimental and model results [24,25]. Constitutive parameters were preliminarily evaluated based on experimental data from in vitro tests. In order to come up with
2.2. Mechanical tests The formulation of a proper constitutive model for heel pad tissues requires the analysis of data from experimental tests [16–18,32,33], to correctly interpret the mechanical response of the tissues. More specifically, the constitutive parameters are evaluated by taking into account the experimental results of MillerYoung et al. [17], because this data is both accurate and complete. Mechanical tests were performed on tissue specimens from human cadavers according to unconfined compression loading conditions. Twenty cadaveric feet were obtained from 6 males and 4 females. The fat pad was removed from the surrounding tissues. Cylindrically shaped samples were cut from the fat pad with a constant diameter and height of about 8 and 10 mm. Five samples were removed from 10 feet in the anterior, central, posterior, medial and lateral regions of the heel pad, in directions perpendicular to the skin surface. A sample was cut from the central region of the heel pad of the other 10 feet, with the long axis parallel to the skin surface in the medial-lateral direction. Tests were developed up to 50% strain, according to different strain rates of 0.1, 0.01, 2100 and 4200%/s. Stress relaxation tests were performed by
518
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
loading the specimen up to 40% strain and holding constant for 60 s.
Table 1 Hyperelastic parameters for heel pad tissues. Kv (MPa)
Experimental data from the mechanical tests showed the strong non-linearity of the mechanical behaviour of heel pad tissues. The stress–strain behaviour of the tissue was influenced by strain rate [17]. Stress relaxation phenomena develop under constant strain conditions [18]. The constitutive analysis of heel pad tissues was consequently performed by a visco-hyperelastic formulation, leading to the following relationship between stress and strain history [25,34]: n
qi (t, ∞ , i , i )
(1)
i=1
where S was the second Piola-Kirchhoff stress tensor, S∞ was stress measure defining the equilibrium response of the material, as the stress condition when viscous phenomena was completely relaxed, C was the right Cauchy-Green strain tensor, qi were internal variables defining the non-equilibrated stresses because of viscous phenomena [35,36], n was the number of viscous branches adopted [37] and t was time. The equilibrium response of the material was specified according to a potential W∞ , which depended on the strain state C and hyperelastic parameters that evaluated tissue stiffness properties at equilibrium. Specifically, constitutive parameters Kv , C1 were related respectively to the initial bulk modulus and shear stiffness, while r and ˛1 were parameters that specified volumetric and deviatoric stiffness evolution with strain, respectively, because of non-linearity of the material response [25]. The evolution of viscous variables was defined according to convolution integral formulations [36,38], depending on stress–strain history and viscous parameters, as relative stiffness ∞ , i and relaxation time i [37]. More details about the constitutive model formulation are reported in Appendix A.
2.4. Constitutive parameters evaluation Constitutive parameters were evaluated by comparing experimental data and results from the constitutive formulation [24]. A specific formulation of the model was developed, accounting for the configuration of mechanical tests. Because nominal stress and tissue stretch were experimentally measured, the mechanical model had to provide for a relationship between the first Piola-Kirchhoff stress tensor, as P = FS, F being the deformation gradient, and principal stretch components, as 1 , 2 , 3 . Assuming direction 1 as the loading direction, the model had to provide the stress compomod for an experimentally imposed stretch exp . Within the nent P11 1 developed constitutive formulation, the stress component P11 was a function of all the principal stretches. Only 1 was experimentally determined and 2 , 3 had be evaluated by analytical methods [39]. Null value could be assumed for lateral stress components P22 , P33 , according to the uniaxial loading condition. Equations P22 (1 , 2 , 3 ) = 0 and P33 (1 , 2 , 3 ) = 0 then made it possible to evaluate 2 , 3 for prescribed values of 1 . The comparison between model results and experimental data, exp as P11 , was performed by a cost function that depended on consti-
C1 (MPa)
˛1 −3
1.06 × 10
4.73 × 10
2.3. Constitutive model
S(C, qi ) = S ∞ (C, Kv , r, C1 , ˛1 ) +
r −2
4.65 × 10
1
1.19 × 100
tutive parameters [24], as:
˝() =
P mod (, 1 1 1 − 11 exp ned P ned
exp
)
2 1/2 (2)
11
i=1
where was the set of constitutive parameters and ned the number of experimental data. The function ˝ was a measurement of the overall difference between experimental and model results when constitutive parameters were adopted. The optimization problem entailed evaluation of the set of constitutive parameters opt that minimizes ˝. The optimization procedure had to account for specific limitations on the parameters domain in order to ensure material stability requirements, as the tendency of hyperelastic strain energy function was to increase strictly with strain and the tangent Poisson ratio needed to be defined within an admissible range [40]. These conditions might be difficult to define by specific boundaries on parameters domain. Material stability requirements were consequently enforced by penalty contributions to the cost function [41]: ˝p () = ˝() +
k ()
(3)
k
Penalty terms k assumed high values when constitutive parameters and/or the associated model results did not satisfy prescribed conditions, leading to the minimum of the cost function to be located within an admissible domain of parameters. If the adopted constitutive model was strongly non-linear, the cost function was characterised by multimodal behaviour and the function presented a global minimum and further local minima. Solving the optimization problem by deterministic methods may result in the definition of only one of the local minima, without generating the optimal solution [42]. The transition out of local minima was possible, from an operational point of view, by adopting stochastic algorithms, as the simulated annealing, making it possible to move towards the global minimum [43]. Nevertheless, stochastic techniques did not guarantee an exact minimum, but only moved closer to the minimum itself. In order to be able to find the minimum, an additional deterministic step had to be introduced. The problem, in agreement with the notes reported in [24,25,39], was consequently approached, using a specific coupled stochastic–deterministic optimization procedure which made it possible to identify and analyze the different minima and provide an efficient and reliable identification of the global minimum [44]. 3. Results The optimization procedure was used in vitro experimental results [17] which lead to the constitutive parameters reported in Tables 1 and 2, with regard to hyperelastic and viscous response, respectively. The hyperelastic parameters entailed the equilibrium initial shear stiffness G∞ to be 9.30 × 10−3 MPa. The shear stiffness parameter was assumed as a reference for evaluating tissue stiffness
Table 2 Viscous parameters for heel pad tissues. 1
1 (s)
2
2 (s)
3
3 (s)
4
4 (s)
8.95 × 10−1
6.23 × 10−4
4.29 × 10−5
1.55 × 10−2
2.76 × 10−3
9.88 × 104
3.12 × 10−3
9.82 × 105
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
Fig. 2. Comparison of experimental data (empty dots) and model results (continuous lines) for unconfined compression tests on tissue samples at different strain rates: 0.01%/s (a) and 0.1%/s (b). The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
properties in the light of the relationship with hyperelastic parameters, as G∞ = 2C1 . The relative stiffness at equilibrium, defined as = 1 − 1 − 2 − 3 , was about 9.96 × 10−2 . It was associated with the mechanical response of the material when the strain rate approached zero [37]. The instantaneous initial shear stiffness, as the shear stiffness when strain rate approaches infinitum, could be evaluated as G0 = G∞ / ∞ . Considering a time-dependent response, the initial shear stiffness ranged between 9.30 × 10−3 and 9.33 × 10−2 MPa for strain rates included between zero and infinitum, respectively. The achieved values were in agreement with the data reported in literature [22]. Both experimental (empty circles) and model (continuous lines) results were reported in Figs. 2 and 3, corresponding to constant strain rate tests developed for different strain rates, as 0.1 (Fig. 2a), 0.01 (Fig. 2b), 2100 (Fig. 3a) and 4200%/s (Fig. 3b). The quality of the fit was reported as the value assumed by the cost function. The results from stress relaxation test are reported in Fig. 4a, while in Fig. 4b the same results were proposed according to normalized stress values. The agreement between experimental data and model results was good for both constant strain rate and stress relaxation tests. The proposed constitutive formulation together with the constitutive parameters was adopted to analyze experimental data from relaxation tests performed by Ledoux and Blevins [18]. Specifically, tests were performed on specimens from human cadaveric feet. All feet were taken from donors without any influencing pathology. Specimens of plantar soft tissues were acquired from different locations, as subcalcaneal, subhallucal, submetatarsal and lateral submidfoot. In particular, data on specimens from the subcalcaneal region, as the heel pad, were compared with results from the proposed constitutive formulation, as reported in Fig. 5. The experimental magnitude of relaxation phenomena was lower than model
519
Fig. 3. Comparison of experimental data (empty dots) and model results (continuous lines) for unconfined compression tests on tissue samples at different strain rates: 2100%/s (a) and 4200%/s (b). The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
Fig. 4. Comparison of experimental data (empty dots) and model results (continuous lines) for stress relaxation tests (40% compression strain): nominal stress vs. time (a) and normalized stress vs. time (b). The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
520
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522 Table 3 Constitutive parameters adopted for the analysis of data from in vivo tests. Kv (MPa)
C1 (MPa) −1
1.25 × 10
Fig. 5. Experimental data (empty dots) from stress relaxation test developed by Ledoux and Blevins [18] together with model prediction (continuous line). Terms a and b evaluated the discrepancy between experimental and model results with regard to instantaneous and equilibrium response, respectively. The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
1 −2
1.23 × 10
2 −1
7.17 × 10
3 −1
1.55 × 10
4 −2
6.52 × 10
6.26 × 10−2
matical formulation was developed to analyze such relationship, as reported in Appendix B. In Fig. 6, model results were reported by a continuous line together with experimental values. It was necessary to rearrange the constitutive parameters to be able to interpret the in vivo responses. Specifically, the hyperelastic parameters Kv , C1 and viscous relative stiffness were properly modified, as reported in Table 3. The modification of hyperelastic and viscous parameters made it possible to account for the influence of in vivo tissue configuration on equilibrium and time-dependent response of the material, respectively. The results showed the capability of the model to interpret the mechanical response of the tissue, with regard to the different experimental situations considered. 4. Discussion
predictions (indicated with term b) while the experimental instantaneous stress was greater than model values (indicated with term a). The different responses were interpreted as a consequence of the age of the donors. The data from Ledoux and Blevins [18] is significantly lower than the data adopted for the constitutive parameters evaluation [17]. In fact, different authors reported the influence of ageing phenomena on tissue stiffness properties [13,45,46], showing results in agreement with the interpretation of the analyzed experimental data. Relevant differences in the mechanical properties of heel pad tissues with regard to in vitro and in vivo conditions have been found by many authors [12,17,18]. Zheng et al. [4] analyzed the mechanical properties of the heel pad by performing indentation tests on healthy young people (21–24 years). A similar technique was adopted by Erdemir et al. [47] and Erdemir et al. [48] on healthy people with average ages 47 and 52 years, respectively. A different approach was adopted by Gefen et al. [49]. In this approach the mechanical properties of the heel pad were investigated by measuring the heel pad tissue deformation and the heel-ground contact pressure during the stance phase of the gait. The analysis of available data made it possible to obtain information about the in vivo relationship between tissue stiffness and strain rate. Specifically, in Fig. 6 the experimental relationship between initial shear stiffness and strain rate was represented by empty dots. Starting from the proposed constitutive model, a specific mathe-
Fig. 6. Influence of strain rate on initial shear modulus for heel pad tissue: model results (continuous line) and data from in vivo experimental tests (empty dots). Experimental data were taken from different authors: (1) Zheng et al. [4], (2) Erdemir et al. [48], (3) Erdemir et al. [47], (4) Gefen et al. [49]. The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
A visco-hyperelastic constitutive model has been developed to interpret the mechanical response of heel pad tissues, accounting for data from morphometric and histological analyses and mechanical tests. In agreement with experimental evidence, the constitutive formulation was capable of accounting for typical features of heel pad tissues mechanical behaviour, such as geometric non-linearity, an almost incompressible response, and time-dependent effects. The constitutive parameters were evaluated by analyzing data from in vitro experimental tests, accounting for different compression strain rates and stress relaxation conditions. Tests in vitro were mainly considered because of the reliable and extensive field of data available. Information about the living system composed of skin, heel pad tissues and bone, could be achieved by in vivo tests. Constitutive parameters were evaluated analyzing data from a specific set of experimental tests, showing good agreement between data and model results. The results confirmed the mechanical coherence between the proposed constitutive formulation and the mechanics of heel pad tissues, by adopting a constitutive model according to a phenomenological approach. The reliability of the achieved parameters was confirmed by the correspondence between experimental and model results. It was also confirmed by the fact that the evaluation accounted for data from tests developed by different authors, according to different loading conditions, and analyzing information from both in vitro and in vivo situations. The ability of the model to interpret the different experimental data represented a relevant aspect of the proposed approach to heel pad mechanics with regard to previous important contributions [18–23]. The coupled computational–experimental approach proposed here was part of a more general investigation addressed to the mechanical response of different tissues [8,24,25,34,39,50].Nevertheless, further efforts are necessary for a comprehensive characterization of plantar tissue mechanics. For example, more experimental data are required from in vitro, ex vivo and in vivo tests, accounting for mechanical properties of both adipose tissues and skin, together with their coupling. The preliminary characterization of the heel pad and skin mechanical response can be developed by analyzing data from in vitro tests. Validation should be performed by analyzing data from ex vivo and in vivo tests [51]. This action should improve the results of the procedure for identification of the constitutive parameters, aiming at a better configuration of numerical models. Such models should be adopted to investigate the biomechanical functionality of the foot, accounting also for interaction phenomena between the foot and surrounding elements, as already reported in literature
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
[52–56]. In addition, when interpreting foot mechanics addressed to direct support of clinical and therapeutic activity, the conditions of ageing and disease must be evaluated [13,16,45] for their effects on the mechanical properties of the plantar tissues. Appendix A. The visco-hyperelastic constitutive model is defined by specifying the Helmholtz free energy function as it follows [25,34]: (C, qi ) = W ∞ (C) +
n
t 0
i=1
1 i ˙ q (s) : Cds 2
(A.1)
where the former term W ∞ is the hyperelastic potential associated to the equilibrium mechanical response of the tissue, while the latter term specified the viscous contribution. The second PiolaKirchhoff stress tensor S is computed according to principles of thermodynamics [35,38], as: S(C, qi ) = 2
(C, qi )
∂
Leading to Eq. (1). S∞ is the stress measure defining the equilibrium response of the material, as the stress condition when viscous phenomena have completely relaxed. The equilibrium stress is evaluated by the equilibrium hyperelastic potential W∞ : S ∞ (C) = 2
∂W ∞ (C) ∂C
(A.3)
Almost incompressible behaviour of the tissue is assumed and the strain energy function can be split into volumetric U∞ and iso˜ ∞ contributions [8]: volumetric W ˜ ∞ (˜I1 ) W ∞ (˜I1 , J) = U ∞ (J) + W
(A.4)
where J was the deformation Jacobian, as J = det(C), while ˜I1 is the first iso-volumetric invariant of the right Cauchy-Green strain tensor, as ˜I1 = tr(J −2/3 C). Because of the strong non-linearity of tissue response, specific polynomial and exponential formulations are provided for, as: Kv [(J − 1)2 + J −r + rJ − (r + 1)] U (J) = 2 + r(r + 1)
(A.5)
C ˜ ∞ (˜I1 ) = 1 {exp[˛1 (˜I1 − 3)] − 1} W ˛1
(A.6)
∞
According to the theory of visco-elasticity, evolution laws of viscous variables are defined by integral formulations [38], as: qi (t) =
i ∞i
t
t − s
exp − 0
i
S ∞ ds
(A.7)
Appendix B. With regard to a pure shear stress configuration, the viscohyperelastic relationship between stress and strain history (1) can be re-written according to the following formulation [38]:
t
(, t) =
g(t − s) 0
d 0 ds ds
(B.1)
where and are shear stress and strain, respectively, 0 specified the instantaneous shear response of the material, as the material response when strain rate approaches infinitum, while g is a standard relaxation function, as: g(t) = ∞ +
n i=1
According to the proposed hyperelastic formulation, the instantaneous stress can be computed as: 0 () = 2C10 exp[˛1 ]
i exp
−t i
(B.2)
(B.3)
where C10 is an instantaneous hyperelastic parameter, evaluated as C10 = C1 / ∞ . The aim pertains to the evaluation of the initial shear modulus, as the tissue shear stiffness corresponding to small strain values. The previous relationship can be linearized by a Taylor series for small shear strain values: 0 = 2C10 →0
(B.4)
During a constant strain rate test, shear strain linearly increases with time according to strain rate k: ␥(t) = kt
(B.5)
Accounting for Eqs. (B.2), (B.4) and (B.5), convolution integral (B.1) leads to the following formulation [57]:
(A.2)
∂C
521
␥→0 = ∞ 2C10 ␥ +
n
i 2C10 k i
1 − exp
−␥
i=1
k i
(B.6)
According to standard solid mechanics definitions, the initial shear modulus can be evaluated as: G=
−␥ d = ∞ 2C10 + i 2C10 exp d␥ ␥→0 k i n
(B.7)
i=1
The average value of shear stiffness for strain rate k can be defined applying the integral mean value theorem over the shear strain range [0, L ], where L is the strain limit of almost-linear shear response, leading to the following relationship: 1 i 0 i 2C1 k ␥L n
G(k) = ∞ 2C10 +
i=1
1 − exp
−␥ L
k i
(B.8)
Conflict of interest statement We, the authors, declare that there are not any financial and personal relationships with other people or organizations that could inappropriately influence (bias) our work. References [1] Jahss MH, Michelson JD, Kummer F. Investigations into the fat pads of the sole of the foot: heel pressure studies. Foot Ankle Int 1992;13:233–42. [2] Whittle MW. Generation and attenuation of transient impulsive forces beneath the foot: a review. Gait Posture 1999;10:264–75. [3] Ledoux WR, Hillstrom HJ. The distributed plantar vertical force of neutrally aligned and pes planus feet. Gait Posture 2002;15:1–9. [4] Zheng YP, Choi YKC, Wong K, Chan S, Mak AFT. Biomechanical assessment of plantar foot tissue in diabetic patients using an ultrasound indentation system. Ultrasound Med Biol 2000;26:451–6. [5] Gefen A, Megido-Ravid M, Azariah M, Itzchak Y, Arcan M. Integration of plantar soft tissue stiffness measurements in routine MRI of diabetic foot. Clin Biomech 2001;16:921–1025. [6] Rome K, Webb P, Unsworth A, Haslock I. Heel pad stiffness in runners with plantar heel pain. Clin Biomech 2001;16:901–5. [7] Maurel W, Wu Y, Magnenat Thalmann N, Thalmann D. Biomechanical models for soft tissue simulation. Berlin: Springer; 1998. [8] Natali AN, Pavan PG, Carniel EL, Lucisano ME, Taglialavoro G. Anisotropic elastodamage constitutive model for the biomechanical analysis of tendons. Med Eng Phys 2005;27:209–14. [9] Kroon M, Holzapfel GA. A new constitutive model for multi-layered collagenous tissues. J Biomech 2008;41:2766–71. [10] Jahss MH, Michelson JD, Desai P, Kaye R, Kummer F, Buschman W, et al. Investigations into the fat pads of the sole of the foot: anatomy and histology. Foot Ankle Int 1992;13:227–32. [11] Snow SW, Bohne WHO. Observations on the fibrous retinacula of the heel pad. Foot Ankle Int 2006;27:632–5.
522
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
[12] Hsu C-C, Tsai W-C, Wang C-L, Pao S-H, Shau Y-W, Chuan Y-S. Microchambers and macrochambers in heel pads: are they functionally different? J Appl Physiol 2007;102:2227–31. [13] Kinoshita H, Francio PR, Murase T, Kawai S, Ogawa T. The mechanical properties of the heel pad in elderly adults. Eur J Appl Physiol 1996;73:404–9. [14] Cavanagh PR. Plantar soft tissue thickness during ground contact in walking. J Biomech 1999;32:623–8. [15] Tong J, Lim CS, Goh OL. Technique to study the biomechanical properties of the human calcaneal heel pad. Foot 2003;13:83–91. [16] Hsu C-C, Tsai W-C, Chen CP-C, Shau Y-W, Wang C-L, Chen MJ-L, et al. Effects of aging on the plantar soft tissue properties under the metatarsal heads at different impact velocities. Ultrasound Med Biol 2005;31:1423–9. [17] Miller-Young JE, Duncan NA, Baroud G. Material properties of the human calcaneal fat pad in compression: experiment and theory. J Biomech 2002;35:1523–31. [18] Ledoux WR, Blevins JJ. The compressive material properties of the plantar soft tissue. J Biomech 2007;40:2975–81. [19] Ledoux WR, Meaney DF, Hillstrom HJ. A quasi-linear, viscoelastic, structural model of the plantar soft tissue with frequency-sensitive damping properties. J Biomech Eng 2004;126:831–7. [20] Wu JZ, Dong RG, Schopper AW. Analysis of effects of friction on the deformation behaviour of soft tissue in unconfined compression tests. J Biomech 2004;37:147–55. [21] Spears IR, Miller-Young JE, Watrìers M, Rome K. The effects of loading condition on stress in the barefooted heel pad. Med Sci Sports Exerc 2005;37: 1030–6. [22] Spears IR, Miller-Young JE. The effect of heel-pad thickness and loading protocol on measured heel-pad stiffness and a standardized protocol for inter-subject comparability. Clin Biomech 2006;21:204–12. [23] Ling HY, Choi PC, Zheng YPS, Lau KT. Extraction of mechanical properties of foot plantar tissue using ultrasound indentation associated with genetic algorithm. J Mater Sci Mater Med 2007;18:1579–86. [24] Natali AN, Carniel EL, Pavan PG, Dario P, Izzo I, Menciassi A. Hyperelastic models for the analysis of soft tissue mechanics: definition of constitutive parameters. In: IEEE/RAS-EMBS international conference on biomedical robotics and biomechatronics. 2006. p. 1–8. [25] Natali AN, Carniel EL, Pavan PG, Sander FG, Dorow C, Geiger M. A constitutive model for the analysis of the biomechanical response of the periodontal ligament. J Biomech Eng 2008;130:031004. [26] Uzel M, Cetinus E, Ekerbicer HC, Karaoguz A. Heel pad thickness and athletic activity in healthy young adults: a sonographic study. J Clin Ultrasound 2006;34:231–6. [27] Rome K. Mechanical properties of the heel pad: current theory and review of the literature. Foot 1998;8:179–85. [28] Fritsch H. Sectional anatomy of connective tissue structures in the hindfoot of the newborn child and the adult. Anat Rec 1996;246:147–54. [29] Muehleman C, Li J, Zhong Z, Brankov JG, Wernick MN. Multiple-image radiography for human soft tissue. J Anat 2006;208:115–24. [30] Weissengruber GE, Egger GF, Hutchinson JR, Groenewald HB, Elsässer L, Famini D, et al. The structure of the cushions in the feet of African elephants (Loxodonta africana). J Anat 2006;209:781–92. [31] Torp-Pedersen ST, Matteoli S, Wilhjelm JE, Amris K, Bech JI, Christensen R, et al. Diagnostic accuracy of heel pad palpation – a phantom study. J Forensic Leg Med 2008;15:437–42. [32] Aerts P, Ker RF, De Clercq D, Ilsley DW. The effects of isolation on the mechanics of the human heel pad. J Anat 1996;188:417–23. [33] Ker RF. The time-dependent mechanical properties of the human heel pad in the context of locomotion. J Exp Biol 1996;199:1501–8.
[34] Natali AN, Pavan PG, Carniel EL, Dorow C. Visco-elastic response of the periodontal ligament: an experimental–numerical approach. J Connect Tissue Res 2004;45:222–30. [35] Holzapfel GA. Non linear solid mechanics – a continuum approach for engineering. New York: John Wiley & Sons Ltd.; 2000. ˜ E, Pena ˜ JA, Doblaré M. On modelling nonlinear viscoelastic effects in lig[36] Pena aments. J Biomech 2008;41:2659–66. [37] McCrum NG, Buckley CP, Bucknall CB. Principles of polymer engineering. New York: Oxford Science Publication; 1997. [38] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [39] Natali AN, Carniel EL, Pavan PG, Bourauel C, Ziegler A, Keilig L. Experimental–numerical analysis of minipig’s multi-rooted teeth. J Biomech 2007;40:1701–8. [40] Schröder J, Neff P, Balzani D. A variational approach for materially stable anisotropic hyperelasticity. Int J Solids Struct 2005;42:4352–71. [41] Boukari D, Fiacco AV. Survey of penalty, exact-penalty and multiplier methods from 1968 to 1993. Optimization 1995;32:301–34. [42] Siarry P, Berthiau G, Durbin F, Haussy J. Enhanced simulated annealing for globally minimizing functions of many-continuous variables. ACM Trans Math Software 1997;23:209–28. [43] Corana A, Marchesi M, Martini C, Ridella S. Minimizing multimodal functions of continuous variables with the simulated annealing algorithm. ACM Trans Math Software 1987;13:262–80. [44] Sylos Labini M, Covitti A, Delvecchio G, Medoro N. A stochastic–deterministic optimization method for eliminating the environmental impact produced by a grounding system. In: Prooceedings of ICREPQ 2003. 2003. [45] Hsu TC, Wang CL, Tsai WC, Kuo JK, Tang FT. Comparison of the mechanical properties of the heel pad between young and elderly adults. Arch Phys Med Rehabil 1998;79:1101–4. [46] Natali AN, Pavan PG, Carniel EL, Dario P, Izzo I. Characterization of soft tissue mechanics with aging. IEEE Eng Med Biol 2008;27:15–22. [47] Erdemir A, Viveiros ML, Cavanagh PR. A numerical-experimental approach for characterizing subject specific hyperelastic properties of the heel pad. In: Summer bioengineering conference, sonesta beach resort in key Biscayne. 2003. [48] Erdemir A, Viveiros ML, Ulbrecht JS, Cavanagh PR. An inverse finite-element model of heel-pad indentation. J Biomech 2006;39:1279–86. [49] Gefen A, Megido-Ravid M, Itzchak Y. In vivo biomechanical behaviour of the human heel pad during the stance phase of gait. J Biomech 2001;34:1661–5. [50] Natali AN, Carniel EL, Gregersen H. Biomechanical behaviour of oesophageal tissues: material and structural configuration, experimental data and constitutive analysis. Med Eng Phys 2009, doi:10.1016/j.medengphy.2009.07.003. [51] Miller-Young JE. Factors affecting human heel pad mechanics: a finite element study. Doctoral Dissertation, Department of Mechanical Engineering, University of Calgary, Canada; 2003. [52] Lemmon D, Shiang TY, Hashmi A, Ulbrecht JS, Cavanagh PR. The effect of insole in therapeutic footwear – a finite element approach. J Biomech 1997;30:615–20. [53] Chen W-P, Ju C-W, Tang F-T. Effects of total contact insoles on the plantar stress redistribution: a finite element analysis. Clin Biomech 2003;18:S17– 24. [54] Cheung JT-M, Zhang M. A 3-dimensional finite element model of the human foot and ankle for insole design. Arch Phys Med Rehabil 2005;86:353–8. [55] Cheung JT-M, Nigg BM. Clinical applications of computational simulation of foot and ankle. Sport Orthop Traumatol 2007;23:264–71. [56] Cheung JT-M, Zhang M. Parametric design of pressure-relieving foot orthosis using statistics-based finite element method. Med Eng Phys 2008;30:269–77. [57] Natali AN, Carniel EL, Pavan PG. Constitutive modelling of inelastic behaviour of cortical bone. Med Eng Phys 2008;30:905–12.
Contents lists available at ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Constitutive formulation and analysis of heel pad tissues mechanics A.N. Natali ∗ , C.G. Fontanella, E.L. Carniel University of Padova, Centre of Mechanics of Biological Materials, Via F. Marzolo 9, I-35131 Padova, Italy
a r t i c l e
i n f o
Article history: Received 28 May 2009 Received in revised form 11 February 2010 Accepted 18 February 2010 Keywords: Heel pad Tissue mechanics Constitutive model Visco-elasticity Mechanical tests
a b s t r a c t This paper presents a visco-hyperelastic constitutive model developed to describe the biomechanical response of heel pad tissues. The model takes into account the typical features of the mechanical response such as large displacement, strain phenomena, and non-linear elasticity together with time-dependent effects. The constitutive model was formulated, starting from the analysis of the complex structural and micro-structural configuration of the tissues, to evaluate the relationship between tissue histology and mechanical properties. To define the constitutive model, experimental data from mechanical tests were analyzed. To obtain information about the mechanical response of the tissue so that the constitutive parameters could be established, data from both in vitro and in vivo tests were investigated. Specifically, the first evaluation of the constitutive parameters was performed by a coupled deterministic and stochastic optimization method, accounting for data from in vitro tests. The comparison of constitutive model results and experimental data confirmed the model’s capability to describe the compression behaviour of the heel pad tissues, regarding both constant strain rate and stress relaxation tests. Based on the data from additional experimental tests, some of the constitutive parameters were modified in order to interpret the in vivo mechanical response of the heel pad tissues. This approach made it possible to interpret the actual mechanical function of the tissues. © 2010 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction The configuration of the foot has been investigated with regard to weight-bearing capacity during static standing and ambulation [1]. Specific loose connective tissues, here defined as plantar tissues, surround bony structures in the foot plantar region. Their conformation and mechanical properties strongly influence body posture during static standing and the gait cycle. Specifically, the mechanical role of plantar tissues is to perform load transmission phenomena, also to dampen shocks generated during the gait or running cycle [2], making it possible to achieve a smooth distribution of plantar pressure [3]. A specific knowledge of plantar tissues’ mechanical behaviour contributes to the investigation of the mechanical functionality of the foot, considering the interaction phenomena between the foot and surrounding structures, such as footwear products [4–6]. The mechanics of plantar tissues must be investigated by both an experimental and a computational approach. Experimental activities are necessary to acquire data about the tissues’ mechanical response. These data will represent the necessary input for the development and validation of computational models and will give an appropriate constitutive formulation for plantar tissues.
∗ Corresponding author. Tel.: +39 049 827 5598; fax: +39 049 827 5604. E-mail address: [email protected] (A.N. Natali).
Detailed information about the structural and micro-structural configuration of the tissues, together with data from mechanical tests, is necessary to define the proper constitutive formulation [7–9]. The plantar region is a composite material of fatty and connective tissues located between skin dermis, fascia and bony segments. Within the plantar region, particular attention is paid to the heel pad, which is the portion of the plantar region interposed between calcaneum and skin that plays a fundamental role in foot mechanics (Fig. 1a). The heel pad is mainly organized according to a honeycomb configuration [10,11]. Fat tissue chambers are embedded and separated from each other by connective septa (Fig. 1b and c). In different locations within the heel pad structure, chambers and septa are characterized by specific dimensions and orientations [12]. These aspects are heavily influenced by the specific local loading conditions. The mechanical properties of the heel pad were investigated using in vivo or ex vivo techniques. Different loading conditions were applied over the heel region and the corresponding displacement or strain was measured [6,13–16], leading to a comprehensive characterization of the mechanical response of the plantar region. The specific stress–strain behaviour of the tissue was analyzed using in vitro tests based on more simple geometry specimens, boundary conditions and loading protocols. Experimental data from in vitro tests can be found in the literature. Miller-Young et al. [17] developed compression tests on cylindrical specimens from
1350-4533/$ – see front matter © 2010 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2010.02.018
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
517
a formulation capable of interpreting the actual response of heel pad living tissues, constitutive parameters were subsequently rearranged accounting for in vivo behaviour. This task was approached by analyzing data from in vivo tests, paying special attention to the relationship between tissue stiffness and strain rate. This relationship was analyzed by the proposed constitutive formulation in order to properly modify the specific constitutive parameters. 2. Materials and methods 2.1. Morphometry and histology The heel pad is the portion of plantar tissues interposed between the calcaneum and the skin (Fig. 1a). The thickness ranges between 14.4 and 24.5 mm with an average value of 18 mm [26]. The distribution of overall heel pad thickness varies according to local loading conditions. For example, the highest thickness values are found where severe interactions develop, such as loading by impacts [27]. The structure of the heel pad is complex and hierarchically organized [28–30]. As reported in Fig. 1b, the heel pad mainly consists of adipose tissue chambers, separated by connective tissue septa [27]. The adipose chambers contain closely packed fat cells with diameters ranging between 100 and 200 m [15] and arranged in lobules within a framework of connective bundles. Larger groups of lobules are enclosed by further bundles of collagen and elastin fibres, which thus form a supporting network (Fig. 1c). This honeycomb configuration of septa entails an overall closed-cells structure [10]. Adipose chambers dimension and orientation depend on the specific location within the heel pad (Fig. 1a and b). Vertically oriented chambers are located in the central portion of the heel pad. The chambers are small in the region close to dermis, progressively increasing towards the calcaneal perichondrium. In the lateral and posterior regions, chambers are small and more transversally oriented [27]. The connective septa configuration varies according to the geometries of the adjacent adipose chambers. Collagen fibres within septa are spirally arranged and superiorly fixed to bone or other septa, inferiorly to other septa or dermis [30]. Compressive stiffness is caused by the almost incompressible behaviour of fat tissue bounded by collagen network [31]. The spiral distribution of collagen fibres within the fatty chambers and connective septa is suitable to limit the bulging of chambers themselves [27]. Fig. 1. Morphometric and histological configuration of the heel pad: magnetic resonance image according to sagittal section (a), longitudinal section by a schematic representation (b) and optical microscope picture of a transversal section (c), from Weissengruber et al. [30].
the human heel pad under different loading conditions and evaluated stress relaxation phenomena that were also studied by Ledoux and Blevins [18] on similar geometry specimens. Different constitutive formulations are reported in literature for heel pad tissues, with particular attention to the hyperelastic and visco-elastic response [18–23]. The main goal of the proposed work was to define the constitutive model and evaluate the associated parameters, aiming at a reliable and efficient characterization of the mechanical behaviour of heel pad tissues. A specific visco-hyperelastic constitutive model was formulated. The proposed model was capable of accounting for the typical response of heel pad tissues, such as large displacements and strains, an almost incompressible behaviour and time-dependent effects. The constitutive parameters were determined by a specific procedure, consisting of optimization methods capable of minimizing the discrepancy between experimental and model results [24,25]. Constitutive parameters were preliminarily evaluated based on experimental data from in vitro tests. In order to come up with
2.2. Mechanical tests The formulation of a proper constitutive model for heel pad tissues requires the analysis of data from experimental tests [16–18,32,33], to correctly interpret the mechanical response of the tissues. More specifically, the constitutive parameters are evaluated by taking into account the experimental results of MillerYoung et al. [17], because this data is both accurate and complete. Mechanical tests were performed on tissue specimens from human cadavers according to unconfined compression loading conditions. Twenty cadaveric feet were obtained from 6 males and 4 females. The fat pad was removed from the surrounding tissues. Cylindrically shaped samples were cut from the fat pad with a constant diameter and height of about 8 and 10 mm. Five samples were removed from 10 feet in the anterior, central, posterior, medial and lateral regions of the heel pad, in directions perpendicular to the skin surface. A sample was cut from the central region of the heel pad of the other 10 feet, with the long axis parallel to the skin surface in the medial-lateral direction. Tests were developed up to 50% strain, according to different strain rates of 0.1, 0.01, 2100 and 4200%/s. Stress relaxation tests were performed by
518
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
loading the specimen up to 40% strain and holding constant for 60 s.
Table 1 Hyperelastic parameters for heel pad tissues. Kv (MPa)
Experimental data from the mechanical tests showed the strong non-linearity of the mechanical behaviour of heel pad tissues. The stress–strain behaviour of the tissue was influenced by strain rate [17]. Stress relaxation phenomena develop under constant strain conditions [18]. The constitutive analysis of heel pad tissues was consequently performed by a visco-hyperelastic formulation, leading to the following relationship between stress and strain history [25,34]: n
qi (t, ∞ , i , i )
(1)
i=1
where S was the second Piola-Kirchhoff stress tensor, S∞ was stress measure defining the equilibrium response of the material, as the stress condition when viscous phenomena was completely relaxed, C was the right Cauchy-Green strain tensor, qi were internal variables defining the non-equilibrated stresses because of viscous phenomena [35,36], n was the number of viscous branches adopted [37] and t was time. The equilibrium response of the material was specified according to a potential W∞ , which depended on the strain state C and hyperelastic parameters that evaluated tissue stiffness properties at equilibrium. Specifically, constitutive parameters Kv , C1 were related respectively to the initial bulk modulus and shear stiffness, while r and ˛1 were parameters that specified volumetric and deviatoric stiffness evolution with strain, respectively, because of non-linearity of the material response [25]. The evolution of viscous variables was defined according to convolution integral formulations [36,38], depending on stress–strain history and viscous parameters, as relative stiffness ∞ , i and relaxation time i [37]. More details about the constitutive model formulation are reported in Appendix A.
2.4. Constitutive parameters evaluation Constitutive parameters were evaluated by comparing experimental data and results from the constitutive formulation [24]. A specific formulation of the model was developed, accounting for the configuration of mechanical tests. Because nominal stress and tissue stretch were experimentally measured, the mechanical model had to provide for a relationship between the first Piola-Kirchhoff stress tensor, as P = FS, F being the deformation gradient, and principal stretch components, as 1 , 2 , 3 . Assuming direction 1 as the loading direction, the model had to provide the stress compomod for an experimentally imposed stretch exp . Within the nent P11 1 developed constitutive formulation, the stress component P11 was a function of all the principal stretches. Only 1 was experimentally determined and 2 , 3 had be evaluated by analytical methods [39]. Null value could be assumed for lateral stress components P22 , P33 , according to the uniaxial loading condition. Equations P22 (1 , 2 , 3 ) = 0 and P33 (1 , 2 , 3 ) = 0 then made it possible to evaluate 2 , 3 for prescribed values of 1 . The comparison between model results and experimental data, exp as P11 , was performed by a cost function that depended on consti-
C1 (MPa)
˛1 −3
1.06 × 10
4.73 × 10
2.3. Constitutive model
S(C, qi ) = S ∞ (C, Kv , r, C1 , ˛1 ) +
r −2
4.65 × 10
1
1.19 × 100
tutive parameters [24], as:
˝() =
P mod (, 1 1 1 − 11 exp ned P ned
exp
)
2 1/2 (2)
11
i=1
where was the set of constitutive parameters and ned the number of experimental data. The function ˝ was a measurement of the overall difference between experimental and model results when constitutive parameters were adopted. The optimization problem entailed evaluation of the set of constitutive parameters opt that minimizes ˝. The optimization procedure had to account for specific limitations on the parameters domain in order to ensure material stability requirements, as the tendency of hyperelastic strain energy function was to increase strictly with strain and the tangent Poisson ratio needed to be defined within an admissible range [40]. These conditions might be difficult to define by specific boundaries on parameters domain. Material stability requirements were consequently enforced by penalty contributions to the cost function [41]: ˝p () = ˝() +
k ()
(3)
k
Penalty terms k assumed high values when constitutive parameters and/or the associated model results did not satisfy prescribed conditions, leading to the minimum of the cost function to be located within an admissible domain of parameters. If the adopted constitutive model was strongly non-linear, the cost function was characterised by multimodal behaviour and the function presented a global minimum and further local minima. Solving the optimization problem by deterministic methods may result in the definition of only one of the local minima, without generating the optimal solution [42]. The transition out of local minima was possible, from an operational point of view, by adopting stochastic algorithms, as the simulated annealing, making it possible to move towards the global minimum [43]. Nevertheless, stochastic techniques did not guarantee an exact minimum, but only moved closer to the minimum itself. In order to be able to find the minimum, an additional deterministic step had to be introduced. The problem, in agreement with the notes reported in [24,25,39], was consequently approached, using a specific coupled stochastic–deterministic optimization procedure which made it possible to identify and analyze the different minima and provide an efficient and reliable identification of the global minimum [44]. 3. Results The optimization procedure was used in vitro experimental results [17] which lead to the constitutive parameters reported in Tables 1 and 2, with regard to hyperelastic and viscous response, respectively. The hyperelastic parameters entailed the equilibrium initial shear stiffness G∞ to be 9.30 × 10−3 MPa. The shear stiffness parameter was assumed as a reference for evaluating tissue stiffness
Table 2 Viscous parameters for heel pad tissues. 1
1 (s)
2
2 (s)
3
3 (s)
4
4 (s)
8.95 × 10−1
6.23 × 10−4
4.29 × 10−5
1.55 × 10−2
2.76 × 10−3
9.88 × 104
3.12 × 10−3
9.82 × 105
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
Fig. 2. Comparison of experimental data (empty dots) and model results (continuous lines) for unconfined compression tests on tissue samples at different strain rates: 0.01%/s (a) and 0.1%/s (b). The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
properties in the light of the relationship with hyperelastic parameters, as G∞ = 2C1 . The relative stiffness at equilibrium, defined as = 1 − 1 − 2 − 3 , was about 9.96 × 10−2 . It was associated with the mechanical response of the material when the strain rate approached zero [37]. The instantaneous initial shear stiffness, as the shear stiffness when strain rate approaches infinitum, could be evaluated as G0 = G∞ / ∞ . Considering a time-dependent response, the initial shear stiffness ranged between 9.30 × 10−3 and 9.33 × 10−2 MPa for strain rates included between zero and infinitum, respectively. The achieved values were in agreement with the data reported in literature [22]. Both experimental (empty circles) and model (continuous lines) results were reported in Figs. 2 and 3, corresponding to constant strain rate tests developed for different strain rates, as 0.1 (Fig. 2a), 0.01 (Fig. 2b), 2100 (Fig. 3a) and 4200%/s (Fig. 3b). The quality of the fit was reported as the value assumed by the cost function. The results from stress relaxation test are reported in Fig. 4a, while in Fig. 4b the same results were proposed according to normalized stress values. The agreement between experimental data and model results was good for both constant strain rate and stress relaxation tests. The proposed constitutive formulation together with the constitutive parameters was adopted to analyze experimental data from relaxation tests performed by Ledoux and Blevins [18]. Specifically, tests were performed on specimens from human cadaveric feet. All feet were taken from donors without any influencing pathology. Specimens of plantar soft tissues were acquired from different locations, as subcalcaneal, subhallucal, submetatarsal and lateral submidfoot. In particular, data on specimens from the subcalcaneal region, as the heel pad, were compared with results from the proposed constitutive formulation, as reported in Fig. 5. The experimental magnitude of relaxation phenomena was lower than model
519
Fig. 3. Comparison of experimental data (empty dots) and model results (continuous lines) for unconfined compression tests on tissue samples at different strain rates: 2100%/s (a) and 4200%/s (b). The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
Fig. 4. Comparison of experimental data (empty dots) and model results (continuous lines) for stress relaxation tests (40% compression strain): nominal stress vs. time (a) and normalized stress vs. time (b). The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
520
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522 Table 3 Constitutive parameters adopted for the analysis of data from in vivo tests. Kv (MPa)
C1 (MPa) −1
1.25 × 10
Fig. 5. Experimental data (empty dots) from stress relaxation test developed by Ledoux and Blevins [18] together with model prediction (continuous line). Terms a and b evaluated the discrepancy between experimental and model results with regard to instantaneous and equilibrium response, respectively. The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
1 −2
1.23 × 10
2 −1
7.17 × 10
3 −1
1.55 × 10
4 −2
6.52 × 10
6.26 × 10−2
matical formulation was developed to analyze such relationship, as reported in Appendix B. In Fig. 6, model results were reported by a continuous line together with experimental values. It was necessary to rearrange the constitutive parameters to be able to interpret the in vivo responses. Specifically, the hyperelastic parameters Kv , C1 and viscous relative stiffness were properly modified, as reported in Table 3. The modification of hyperelastic and viscous parameters made it possible to account for the influence of in vivo tissue configuration on equilibrium and time-dependent response of the material, respectively. The results showed the capability of the model to interpret the mechanical response of the tissue, with regard to the different experimental situations considered. 4. Discussion
predictions (indicated with term b) while the experimental instantaneous stress was greater than model values (indicated with term a). The different responses were interpreted as a consequence of the age of the donors. The data from Ledoux and Blevins [18] is significantly lower than the data adopted for the constitutive parameters evaluation [17]. In fact, different authors reported the influence of ageing phenomena on tissue stiffness properties [13,45,46], showing results in agreement with the interpretation of the analyzed experimental data. Relevant differences in the mechanical properties of heel pad tissues with regard to in vitro and in vivo conditions have been found by many authors [12,17,18]. Zheng et al. [4] analyzed the mechanical properties of the heel pad by performing indentation tests on healthy young people (21–24 years). A similar technique was adopted by Erdemir et al. [47] and Erdemir et al. [48] on healthy people with average ages 47 and 52 years, respectively. A different approach was adopted by Gefen et al. [49]. In this approach the mechanical properties of the heel pad were investigated by measuring the heel pad tissue deformation and the heel-ground contact pressure during the stance phase of the gait. The analysis of available data made it possible to obtain information about the in vivo relationship between tissue stiffness and strain rate. Specifically, in Fig. 6 the experimental relationship between initial shear stiffness and strain rate was represented by empty dots. Starting from the proposed constitutive model, a specific mathe-
Fig. 6. Influence of strain rate on initial shear modulus for heel pad tissue: model results (continuous line) and data from in vivo experimental tests (empty dots). Experimental data were taken from different authors: (1) Zheng et al. [4], (2) Erdemir et al. [48], (3) Erdemir et al. [47], (4) Gefen et al. [49]. The value reported for ˝ represented the discrepancy between experimental data and model results as the value assumed by the cost function.
A visco-hyperelastic constitutive model has been developed to interpret the mechanical response of heel pad tissues, accounting for data from morphometric and histological analyses and mechanical tests. In agreement with experimental evidence, the constitutive formulation was capable of accounting for typical features of heel pad tissues mechanical behaviour, such as geometric non-linearity, an almost incompressible response, and time-dependent effects. The constitutive parameters were evaluated by analyzing data from in vitro experimental tests, accounting for different compression strain rates and stress relaxation conditions. Tests in vitro were mainly considered because of the reliable and extensive field of data available. Information about the living system composed of skin, heel pad tissues and bone, could be achieved by in vivo tests. Constitutive parameters were evaluated analyzing data from a specific set of experimental tests, showing good agreement between data and model results. The results confirmed the mechanical coherence between the proposed constitutive formulation and the mechanics of heel pad tissues, by adopting a constitutive model according to a phenomenological approach. The reliability of the achieved parameters was confirmed by the correspondence between experimental and model results. It was also confirmed by the fact that the evaluation accounted for data from tests developed by different authors, according to different loading conditions, and analyzing information from both in vitro and in vivo situations. The ability of the model to interpret the different experimental data represented a relevant aspect of the proposed approach to heel pad mechanics with regard to previous important contributions [18–23]. The coupled computational–experimental approach proposed here was part of a more general investigation addressed to the mechanical response of different tissues [8,24,25,34,39,50].Nevertheless, further efforts are necessary for a comprehensive characterization of plantar tissue mechanics. For example, more experimental data are required from in vitro, ex vivo and in vivo tests, accounting for mechanical properties of both adipose tissues and skin, together with their coupling. The preliminary characterization of the heel pad and skin mechanical response can be developed by analyzing data from in vitro tests. Validation should be performed by analyzing data from ex vivo and in vivo tests [51]. This action should improve the results of the procedure for identification of the constitutive parameters, aiming at a better configuration of numerical models. Such models should be adopted to investigate the biomechanical functionality of the foot, accounting also for interaction phenomena between the foot and surrounding elements, as already reported in literature
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
[52–56]. In addition, when interpreting foot mechanics addressed to direct support of clinical and therapeutic activity, the conditions of ageing and disease must be evaluated [13,16,45] for their effects on the mechanical properties of the plantar tissues. Appendix A. The visco-hyperelastic constitutive model is defined by specifying the Helmholtz free energy function as it follows [25,34]: (C, qi ) = W ∞ (C) +
n
t 0
i=1
1 i ˙ q (s) : Cds 2
(A.1)
where the former term W ∞ is the hyperelastic potential associated to the equilibrium mechanical response of the tissue, while the latter term specified the viscous contribution. The second PiolaKirchhoff stress tensor S is computed according to principles of thermodynamics [35,38], as: S(C, qi ) = 2
(C, qi )
∂
Leading to Eq. (1). S∞ is the stress measure defining the equilibrium response of the material, as the stress condition when viscous phenomena have completely relaxed. The equilibrium stress is evaluated by the equilibrium hyperelastic potential W∞ : S ∞ (C) = 2
∂W ∞ (C) ∂C
(A.3)
Almost incompressible behaviour of the tissue is assumed and the strain energy function can be split into volumetric U∞ and iso˜ ∞ contributions [8]: volumetric W ˜ ∞ (˜I1 ) W ∞ (˜I1 , J) = U ∞ (J) + W
(A.4)
where J was the deformation Jacobian, as J = det(C), while ˜I1 is the first iso-volumetric invariant of the right Cauchy-Green strain tensor, as ˜I1 = tr(J −2/3 C). Because of the strong non-linearity of tissue response, specific polynomial and exponential formulations are provided for, as: Kv [(J − 1)2 + J −r + rJ − (r + 1)] U (J) = 2 + r(r + 1)
(A.5)
C ˜ ∞ (˜I1 ) = 1 {exp[˛1 (˜I1 − 3)] − 1} W ˛1
(A.6)
∞
According to the theory of visco-elasticity, evolution laws of viscous variables are defined by integral formulations [38], as: qi (t) =
i ∞i
t
t − s
exp − 0
i
S ∞ ds
(A.7)
Appendix B. With regard to a pure shear stress configuration, the viscohyperelastic relationship between stress and strain history (1) can be re-written according to the following formulation [38]:
t
(, t) =
g(t − s) 0
d 0 ds ds
(B.1)
where and are shear stress and strain, respectively, 0 specified the instantaneous shear response of the material, as the material response when strain rate approaches infinitum, while g is a standard relaxation function, as: g(t) = ∞ +
n i=1
According to the proposed hyperelastic formulation, the instantaneous stress can be computed as: 0 () = 2C10 exp[˛1 ]
i exp
−t i
(B.2)
(B.3)
where C10 is an instantaneous hyperelastic parameter, evaluated as C10 = C1 / ∞ . The aim pertains to the evaluation of the initial shear modulus, as the tissue shear stiffness corresponding to small strain values. The previous relationship can be linearized by a Taylor series for small shear strain values: 0 = 2C10 →0
(B.4)
During a constant strain rate test, shear strain linearly increases with time according to strain rate k: ␥(t) = kt
(B.5)
Accounting for Eqs. (B.2), (B.4) and (B.5), convolution integral (B.1) leads to the following formulation [57]:
(A.2)
∂C
521
␥→0 = ∞ 2C10 ␥ +
n
i 2C10 k i
1 − exp
−␥
i=1
k i
(B.6)
According to standard solid mechanics definitions, the initial shear modulus can be evaluated as: G=
−␥ d = ∞ 2C10 + i 2C10 exp d␥ ␥→0 k i n
(B.7)
i=1
The average value of shear stiffness for strain rate k can be defined applying the integral mean value theorem over the shear strain range [0, L ], where L is the strain limit of almost-linear shear response, leading to the following relationship: 1 i 0 i 2C1 k ␥L n
G(k) = ∞ 2C10 +
i=1
1 − exp
−␥ L
k i
(B.8)
Conflict of interest statement We, the authors, declare that there are not any financial and personal relationships with other people or organizations that could inappropriately influence (bias) our work. References [1] Jahss MH, Michelson JD, Kummer F. Investigations into the fat pads of the sole of the foot: heel pressure studies. Foot Ankle Int 1992;13:233–42. [2] Whittle MW. Generation and attenuation of transient impulsive forces beneath the foot: a review. Gait Posture 1999;10:264–75. [3] Ledoux WR, Hillstrom HJ. The distributed plantar vertical force of neutrally aligned and pes planus feet. Gait Posture 2002;15:1–9. [4] Zheng YP, Choi YKC, Wong K, Chan S, Mak AFT. Biomechanical assessment of plantar foot tissue in diabetic patients using an ultrasound indentation system. Ultrasound Med Biol 2000;26:451–6. [5] Gefen A, Megido-Ravid M, Azariah M, Itzchak Y, Arcan M. Integration of plantar soft tissue stiffness measurements in routine MRI of diabetic foot. Clin Biomech 2001;16:921–1025. [6] Rome K, Webb P, Unsworth A, Haslock I. Heel pad stiffness in runners with plantar heel pain. Clin Biomech 2001;16:901–5. [7] Maurel W, Wu Y, Magnenat Thalmann N, Thalmann D. Biomechanical models for soft tissue simulation. Berlin: Springer; 1998. [8] Natali AN, Pavan PG, Carniel EL, Lucisano ME, Taglialavoro G. Anisotropic elastodamage constitutive model for the biomechanical analysis of tendons. Med Eng Phys 2005;27:209–14. [9] Kroon M, Holzapfel GA. A new constitutive model for multi-layered collagenous tissues. J Biomech 2008;41:2766–71. [10] Jahss MH, Michelson JD, Desai P, Kaye R, Kummer F, Buschman W, et al. Investigations into the fat pads of the sole of the foot: anatomy and histology. Foot Ankle Int 1992;13:227–32. [11] Snow SW, Bohne WHO. Observations on the fibrous retinacula of the heel pad. Foot Ankle Int 2006;27:632–5.
522
A.N. Natali et al. / Medical Engineering & Physics 32 (2010) 516–522
[12] Hsu C-C, Tsai W-C, Wang C-L, Pao S-H, Shau Y-W, Chuan Y-S. Microchambers and macrochambers in heel pads: are they functionally different? J Appl Physiol 2007;102:2227–31. [13] Kinoshita H, Francio PR, Murase T, Kawai S, Ogawa T. The mechanical properties of the heel pad in elderly adults. Eur J Appl Physiol 1996;73:404–9. [14] Cavanagh PR. Plantar soft tissue thickness during ground contact in walking. J Biomech 1999;32:623–8. [15] Tong J, Lim CS, Goh OL. Technique to study the biomechanical properties of the human calcaneal heel pad. Foot 2003;13:83–91. [16] Hsu C-C, Tsai W-C, Chen CP-C, Shau Y-W, Wang C-L, Chen MJ-L, et al. Effects of aging on the plantar soft tissue properties under the metatarsal heads at different impact velocities. Ultrasound Med Biol 2005;31:1423–9. [17] Miller-Young JE, Duncan NA, Baroud G. Material properties of the human calcaneal fat pad in compression: experiment and theory. J Biomech 2002;35:1523–31. [18] Ledoux WR, Blevins JJ. The compressive material properties of the plantar soft tissue. J Biomech 2007;40:2975–81. [19] Ledoux WR, Meaney DF, Hillstrom HJ. A quasi-linear, viscoelastic, structural model of the plantar soft tissue with frequency-sensitive damping properties. J Biomech Eng 2004;126:831–7. [20] Wu JZ, Dong RG, Schopper AW. Analysis of effects of friction on the deformation behaviour of soft tissue in unconfined compression tests. J Biomech 2004;37:147–55. [21] Spears IR, Miller-Young JE, Watrìers M, Rome K. The effects of loading condition on stress in the barefooted heel pad. Med Sci Sports Exerc 2005;37: 1030–6. [22] Spears IR, Miller-Young JE. The effect of heel-pad thickness and loading protocol on measured heel-pad stiffness and a standardized protocol for inter-subject comparability. Clin Biomech 2006;21:204–12. [23] Ling HY, Choi PC, Zheng YPS, Lau KT. Extraction of mechanical properties of foot plantar tissue using ultrasound indentation associated with genetic algorithm. J Mater Sci Mater Med 2007;18:1579–86. [24] Natali AN, Carniel EL, Pavan PG, Dario P, Izzo I, Menciassi A. Hyperelastic models for the analysis of soft tissue mechanics: definition of constitutive parameters. In: IEEE/RAS-EMBS international conference on biomedical robotics and biomechatronics. 2006. p. 1–8. [25] Natali AN, Carniel EL, Pavan PG, Sander FG, Dorow C, Geiger M. A constitutive model for the analysis of the biomechanical response of the periodontal ligament. J Biomech Eng 2008;130:031004. [26] Uzel M, Cetinus E, Ekerbicer HC, Karaoguz A. Heel pad thickness and athletic activity in healthy young adults: a sonographic study. J Clin Ultrasound 2006;34:231–6. [27] Rome K. Mechanical properties of the heel pad: current theory and review of the literature. Foot 1998;8:179–85. [28] Fritsch H. Sectional anatomy of connective tissue structures in the hindfoot of the newborn child and the adult. Anat Rec 1996;246:147–54. [29] Muehleman C, Li J, Zhong Z, Brankov JG, Wernick MN. Multiple-image radiography for human soft tissue. J Anat 2006;208:115–24. [30] Weissengruber GE, Egger GF, Hutchinson JR, Groenewald HB, Elsässer L, Famini D, et al. The structure of the cushions in the feet of African elephants (Loxodonta africana). J Anat 2006;209:781–92. [31] Torp-Pedersen ST, Matteoli S, Wilhjelm JE, Amris K, Bech JI, Christensen R, et al. Diagnostic accuracy of heel pad palpation – a phantom study. J Forensic Leg Med 2008;15:437–42. [32] Aerts P, Ker RF, De Clercq D, Ilsley DW. The effects of isolation on the mechanics of the human heel pad. J Anat 1996;188:417–23. [33] Ker RF. The time-dependent mechanical properties of the human heel pad in the context of locomotion. J Exp Biol 1996;199:1501–8.
[34] Natali AN, Pavan PG, Carniel EL, Dorow C. Visco-elastic response of the periodontal ligament: an experimental–numerical approach. J Connect Tissue Res 2004;45:222–30. [35] Holzapfel GA. Non linear solid mechanics – a continuum approach for engineering. New York: John Wiley & Sons Ltd.; 2000. ˜ E, Pena ˜ JA, Doblaré M. On modelling nonlinear viscoelastic effects in lig[36] Pena aments. J Biomech 2008;41:2659–66. [37] McCrum NG, Buckley CP, Bucknall CB. Principles of polymer engineering. New York: Oxford Science Publication; 1997. [38] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [39] Natali AN, Carniel EL, Pavan PG, Bourauel C, Ziegler A, Keilig L. Experimental–numerical analysis of minipig’s multi-rooted teeth. J Biomech 2007;40:1701–8. [40] Schröder J, Neff P, Balzani D. A variational approach for materially stable anisotropic hyperelasticity. Int J Solids Struct 2005;42:4352–71. [41] Boukari D, Fiacco AV. Survey of penalty, exact-penalty and multiplier methods from 1968 to 1993. Optimization 1995;32:301–34. [42] Siarry P, Berthiau G, Durbin F, Haussy J. Enhanced simulated annealing for globally minimizing functions of many-continuous variables. ACM Trans Math Software 1997;23:209–28. [43] Corana A, Marchesi M, Martini C, Ridella S. Minimizing multimodal functions of continuous variables with the simulated annealing algorithm. ACM Trans Math Software 1987;13:262–80. [44] Sylos Labini M, Covitti A, Delvecchio G, Medoro N. A stochastic–deterministic optimization method for eliminating the environmental impact produced by a grounding system. In: Prooceedings of ICREPQ 2003. 2003. [45] Hsu TC, Wang CL, Tsai WC, Kuo JK, Tang FT. Comparison of the mechanical properties of the heel pad between young and elderly adults. Arch Phys Med Rehabil 1998;79:1101–4. [46] Natali AN, Pavan PG, Carniel EL, Dario P, Izzo I. Characterization of soft tissue mechanics with aging. IEEE Eng Med Biol 2008;27:15–22. [47] Erdemir A, Viveiros ML, Cavanagh PR. A numerical-experimental approach for characterizing subject specific hyperelastic properties of the heel pad. In: Summer bioengineering conference, sonesta beach resort in key Biscayne. 2003. [48] Erdemir A, Viveiros ML, Ulbrecht JS, Cavanagh PR. An inverse finite-element model of heel-pad indentation. J Biomech 2006;39:1279–86. [49] Gefen A, Megido-Ravid M, Itzchak Y. In vivo biomechanical behaviour of the human heel pad during the stance phase of gait. J Biomech 2001;34:1661–5. [50] Natali AN, Carniel EL, Gregersen H. Biomechanical behaviour of oesophageal tissues: material and structural configuration, experimental data and constitutive analysis. Med Eng Phys 2009, doi:10.1016/j.medengphy.2009.07.003. [51] Miller-Young JE. Factors affecting human heel pad mechanics: a finite element study. Doctoral Dissertation, Department of Mechanical Engineering, University of Calgary, Canada; 2003. [52] Lemmon D, Shiang TY, Hashmi A, Ulbrecht JS, Cavanagh PR. The effect of insole in therapeutic footwear – a finite element approach. J Biomech 1997;30:615–20. [53] Chen W-P, Ju C-W, Tang F-T. Effects of total contact insoles on the plantar stress redistribution: a finite element analysis. Clin Biomech 2003;18:S17– 24. [54] Cheung JT-M, Zhang M. A 3-dimensional finite element model of the human foot and ankle for insole design. Arch Phys Med Rehabil 2005;86:353–8. [55] Cheung JT-M, Nigg BM. Clinical applications of computational simulation of foot and ankle. Sport Orthop Traumatol 2007;23:264–71. [56] Cheung JT-M, Zhang M. Parametric design of pressure-relieving foot orthosis using statistics-based finite element method. Med Eng Phys 2008;30:269–77. [57] Natali AN, Carniel EL, Pavan PG. Constitutive modelling of inelastic behaviour of cortical bone. Med Eng Phys 2008;30:905–12.