Parameter identification of hyperelastic material properties of the heel pad based on an analytical contact mechanics model of a spherical indentation

Author’s Accepted Manuscript Parameter identification of hyperelastic material properties of the heel pad based on an analytical contact mechanics model of a spherical indentation Ryo Suzuki, Kohta Ito, Taeyong Lee, Naomichi Ogihara www.elsevier.com/locate/jmbbm

PII: DOI: Reference:

S1751-6161(16)30333-2 http://dx.doi.org/10.1016/j.jmbbm.2016.09.027 JMBBM2084

To appear in: Journal of the Mechanical Behavior of Biomedical Materials Received date: 3 April 2016 Revised date: 4 September 2016 Accepted date: 21 September 2016 Cite this article as: Ryo Suzuki, Kohta Ito, Taeyong Lee and Naomichi Ogihara, Parameter identification of hyperelastic material properties of the heel pad based on an analytical contact mechanics model of a spherical indentation, Journal of the Mechanical Behavior of Biomedical Materials, http://dx.doi.org/10.1016/j.jmbbm.2016.09.027 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Parameter identification of hyperelastic material properties of the heel pad based on an analytical contact mechanics model of a spherical indentation Ryo Suzuki1, Kohta Ito1, Taeyong Lee2, Naomichi Ogihara1 1

Department of Mechanical Engineering, Keio University, Yokohama, Japan

2

College of Science & Industry Convergence, Ewha Womans University, Seoul, Korea

*

Corresponding Author. Naomichi Ogihara. Department of Mechanical Engineering, Keio University,

3-14-1, Hiyoshi, Kohoku-ku, Yokohama 223-8522, Japan. [email protected]

Abstract Accurate identification of the material properties of the plantar soft tissue is important for computer-aided analysis of foot pathologies and design of therapeutic footwear interventions based on subject-specific models of the foot. However, parameter identification of the hyperelastic material properties of plantar soft tissues usually requires an inverse finite element analysis due to the lack of a practical contact model of the indentation test. In the present study, we derive an analytical contact model of a spherical indentation test in order to directly estimate the material properties of the plantar soft tissue. Force-displacement curves of the heel pads are obtained through an indentation experiment. The experimental data are fit to the analytical stress-strain solution of the spherical indentation in order to obtain the parameters. A spherical indentation approach successfully predicted the non-linear material properties of the heel pad without iterative finite element calculation. The force-displacement curve obtained in the present study was found to be situated lower than those identified in previous studies. The proposed framework for identifying the hyperelastic material parameters may facilitate the development of subject-specific FE modeling of the foot for possible clinical and ergonomic applications. Keywords foot; plantar soft tissue; finite-element analysis; elasticity

1. Introduction Computer modeling of the foot using a finite element (FE) method has recently been applied to predict movements, deformations, and mechanical stresses in the foot bones and the plantar soft tissues in order to obtain a basic understanding of the biomechanics of the foot (Chen et al., 2001; Cheung and Zhang, 2005, 2008; Gu et al., 2010; Chen et al., 2012; Luximon et al., 2012), computer-aided designs of footwear (Cheung et al., 2006; Yu et al., 2008, 2013), analysis of foot pathologies such as the diabetic foot (Jacob and Patil, 1999; Gefen, 2003; Thomas et al., 2004; Fernandez et al., 2012; Guiotto et al., 2014; Chatzistergos et al., 2014), metatarsal fracture (Gu et al., 2011; Brilakis et al., 2012), plantar fasciitis (Gefen, 2002; Cheung et al., 2006; Cheng et al., 2008), and hallux deformity (Budhabhatti et al., 2007; Isvilanonda et al., 2012; Wong et al., 2014, 2015),

and to study therapeutic footwear intervention (Lemmon et al., 1997; Chen et al., 2003; Cheung et al., 2005; Erdemir et al., 2005; Goske et al., 2006; Luo et al., 2011; Spirka et al., 2014; Chatzistergos et al., 2015; Chen et al., 2015). Dynamic FE analyses of the foot during walking have also been attempted recently in attempts to clarify how the foot complex dynamically and functionally interacts with the ground during walking (Halloran et al., 2011; Qian et al., 2013). In many of these FE models, the soft tissue of the foot was modeled by general hyperelastic material models, such as the generalized Rivlin model (polynomial model) or the Ogden model. The values of the coefficients for the material properties of the tissue were generally determined by referring to literature [e.g., Lemmon et al. (1997) for the generalized Rivlin model (polynomial model) and Erdemir et al. (2006) for the Ogden model]. The accuracy of the FE model prediction depends in part on the accuracy of the material parameters defining the non-linear material properties of the foot. Therefore, researchers try to characterize the material properties of plantar soft tissue for more accurate and objective modeling of the foot (Chokhandre et al., 2012; Fontenella et al., 2012; Chatzistergos et al., 2015). In these parameter identification studies, FE models were used to search for the parameters that best fit an experimentally obtained force-displacement curve. However, the iterative use of the FE model in the optimization of the parameters is time consuming and requires FE simulation software. If we could set up a framework to identify the hyperelastic material parameters based on an analytical model of contact mechanics, the parameter identification process would be greatly facilitated, possibly leading to subject-specific FE modeling of the foot and applications in therapeutic footwear intervention. Quantifying the mechanical properties of foot tissue is also important for early diagnosis and possible interventions for tissue abnormalities such as diabetic foot (Fischer, 1987; Klaesner et al., 2001). Diabetic feet often develop severe ulceration, which frequently leads to lower limb amputation (Boulton et al., 2005; Singh et al., 2005). It has been demonstrated that plantar soft tissue tends to become stiffer in diabetic patients (Gefen et al. 2001; Klaesner et al, 2002; Pai and Ledoux, 2010), indicating that parameter identification of the mechanical properties of the foot tissue is clinically significant for possible early diagnosis and interventions of diabetic foot. However, accurate mechanical characterization has been difficult as mechanical response of the soft tissue is influenced by the size and thickness of the soft tissue, the underlying rigid material (bone), as well as the indenter. Therefore, in the present study, we propose a method by which to identify the hyperelastic material properties of the heel pad using an analytical model of contact mechanics. Specifically, we use a spherical indenter to measure the force-displacement relationship of the heel pads because the contact mechanics between a sphere and a half-space is theoretically better characterized than that of a flat-top cylinder, which may have a complex stress concentration at the indenter edge (Mattice et

al., 2006). The material properties were identified by fitting experimentally obtained force-displacement curves to an analytical model of contact mechanics. Finally, in order to evaluate the feasibility of the present study, we tested whether simulated spherical indentation can actually reproduce the actual experimental results.

2. Materials and Methods 2.1 Experiment The hyperelastic material properties of the heel pads were obtained through an indentation experiment. Three university students (24.2± 0.5years, 168.9±11.1 cm, 62.6±11.0 kg) participated in the experiment. Each participant provided written informed consent prior to the experiment, and the experiment was approved by the Ethics Committee of the Faculty of Science and Technology, Keio University. Each subject sat on a chair, and his/her right foot was fixed to a wooden platform using a Velcro strap (Figure 1). There was space beneath the heel so that the heel pad could be compressed using an aluminum spherical indenter. The diameter of the spherical indenter was chosen to be 30 mm, so that the soft tissue did not protrude behind the spherical indenter. Larger indenters are known to produce more reliable measurements than smaller indenters (Spears et al., 2006). The indenter was controlled by a material testing instrument (EZ-X, Shimadzu, Kyoto, Japan), and the displacement of the indenter and the force applied to the indenter were measured simultaneously at 100 Hz. The accuracy of the force and displacement measurement was approximately 1% of indicated value and 0.01 mm, respectively. Using this system, force displacement curves were obtained. For the measurement of the force displacement curve, the heel of each subject was indented at 0.8 mm/s (Fontanella et al., 2012) to 10 mm, which is approximately 50% of initial thickness of the heel pad (Erdemir et al., 2006; Chatzistergos et al., 2015), and the time changes of the indentation force and the displacement were measured. The initial point of deformation was determined as the point at which a threshold force of 0.01 N was detected. We took five measurements for each curve and subject. We allowed at least one minute between two consecutive trails for the heel pad to recover from deformation.

2.2 Parameter identification In the present study, the experimentally obtained force-displacement curves from indentation were fit using the following model of spherical indentation for hyperelastic materials (Mooney, 1940; Rivlin, 1948, 1951). The strain energy function of the conventional Mooney-Rivlin model can be written as follows:

W = c10 ( I1 - 3) + c01 ( I 2 - 3)

(1)

More generally, the strain energy potential can be represented as a polynomial function (generalized Rivlin model):

W=

3

å c (I ij

1

- 3)i ( I 2 - 3) j

(2)

i + j =1

where W is the strain energy function, I1,2 are the first and second deviatoric strain invariants, and cij are coefficients defining the material characteristics. Under the condition of incompressibility and uniaxial loading, the strain invariants can be written as

I1 = l 2 + 2 / l , I 2 = 2l + 1/ l 2 where l is the deviatoric principle stretch ( l = 1 + e ), and

(3)

e is the strain. Therefore, uniaxial

nominal stress-stretch (σ–λ) equation for this material is as follows (Treloar, 1975):

s=

¶W (l ) ¶l

(4)

If a compressive force and indentation are taken to be positive, then we have

sI = -

¶W (lI ) ¶lI

lI = 1 - e I where

(5) (6)

e I is the indentation strain, and s I is the indentation stress. The relationship between the force applied by the spherical indentation (FI) and the

indentation stress can be written as

FI = p a 2s I

(7)

where a is the contact radius. The relationship between the contact radius a and the indentation displacement δ in hyperelastic contact is modeled as follows (Lin et al., 2009):

a = R p -(d / R ) qd r

(8)

where R is the radius of the spherical indenter, and p, q, and r are coefficients. For hyperelastic materials such as the generalized Rivlin or Ogden models, p, q, and r are found to be 0.483, 0.101, and 0.567, respectively, using the FE method. See Appendix A for the calculation details. In spherical indentation, deformation can be generally represented based on Hertz contact mechanics if the indented substrate can be approximated as a linear elastic material of infinite thickness (Long et al., 2011). However, the thickness of the plantar soft tissue is very small in comparison with the amount of compressive deformation dealt with herein. In order to take the effect of the finite thickness of the substrate into consideration, we assume that a large deformation due to spherical indentation can be represented by a combination of Hertzian contact deformation and volumetric deformation (Figure 2) as

e I = e H + eV

(9)

where

eI

is the indentation strain,

eH

is the Hertz strain, and

eV

is the volumetric strain (Tani

et al, 2009). The Hertz strain of a hyperelastic material is empirically determined (Tabor, 1951; Field and Swain, 1995; Taljat et al., 1998; Herbert et al., 2001) as

e H = 0.2

a R

(10)

On the other hand, the volumetric strain can be estimated based on the assumption that the change in volume of compression due to spherical indentation can be approximated by the change in volume of a paraboloid (Figure 2). Therefore, the relationship between the indentation displacement and the volumetric strain can be derived as

ev =

1 d 2 Rd - d 2 dd h ò0 h + 2 Rd - d 2

(11)

where h is the thickness of the substrate material. See the Appendix B for the calculation details. In this model, if h = 20 mm and R = 5, 10, and 15 mm, the Hertz strain accounts for approx. 89.4, 70.7 and 53.7 %, respectively, of the indentation strain when the indentation depth is R/2. The force applied by the spherical indentation can be represented as a function of δ as follows:

FI = p {a(d )}2 s I (1 - e I (d ))

(12)

When the indentation strain is infinitesimally small, the following linear relationship should be established (Lin et al., 2009):

sI =

20 E0 eI 3p (1 -n 2 )

(13)

where E0 is the initial Young’s modulus, and ν is Poisson’s ratio of the material. Under this condition, the relationships between FI and δ are defined for the hyperelastic strain energy functions (Lin et al., 2009). See the Appendix C for the actual representation of Eq. (12). The parameters cij are estimated by solving the following minimization problem:

å{F

ex

(d ) - FI (d )}2 ® min

(14)

where Fex is the experimentally obtained force-displacement curve. In the present study, Microsoft Excel 2013 Solver was used to solve this optimization problem. However, for a robust FE analysis of the generalized Rivlin material, the parameters must be determine so as to satisfy the condition whereby the strain energy function becomes a convex function (Schroder and Neff, 2003). Therefore, the minimization problem was solved under the inequality constraints c10 , c20 , c30 ³ 0 and the equality constraints cij = 0 (except for c10 , c20 , c30 ) (Yeoh, 1983). Since W is here represented as a

function of only I1, the material parameters can theoretically be obtained from only a uniaxial compression test (Hyun et al., 2012). Similarly, the strain energy function of the incompressible Ogden material model, can be represented as (Ogden, 1972)

W (C , a ) =

2C

a2

(l

a

x

+ lya + lza - 3)

(15)

where λ is the deviatoric principal stretch, and C and α are coefficients. The uniaxial nominal stress-stretch (σ–λ) equation for this material is represented as follows:

sI = -

2C

a

(l

a -1

I

- lI-a /2-1 )

(16)

By using the same procedure described above, we solved for the values of C and α that best fit the experimental data. In order to quantify the fitting error, we calculated the root-mean-square (RMS) error Erms for each curve as

Erms =

1 N ( Fex (d i ) - FI (d i )) 2 å N i =1

(17)

where N is the number of values.

2.3 Finite element modeling In order to test whether the FE method with the experimentally identified hyperelastic material parameters can reproduce the measured force-displacement curve, a schematic FE model of the heel with the spherical indenter was constructed (Figure 3). The heel pad is assumed to be a cylinder with a diameter of 60 mm (Erdemir et al., 2006), which is modeled as an axisymmetric FE model. Axisymmetric, four-node quadrilateral ring element was used for the model construction. Size of each element is 0.1 mm × 0.1 mm. The heel pad was situated on a rigid body representing the calcaneus. The indenter was modeled as a rigid sphere and the heel-indenter contact was assumed to be frictionless as the surface of the spherical indenter is very smooth and the use of gel did not alter the measured force-displacement curve. The simulation conditions were identical to those of the experiments described above. We constructed FE models with the obtained hyperelastic strain energy functions, and the calculated force-displacement curves were compared with the measured curves. Implicit static analysis of the contact between the heel and the spherical indenter was conducted from the indentation displacement = 0 to 10 mm at 0.1 mm interval using Marc ver.2013.0.0 (MSC Software, USA) FE software. In the present study, we considered the (1) Mooney-Rivlin, (2) generalized Rivlin (polynomial), and (3) Ogden models in order to determine how the choice of the hyperelastic model

affects the conformation between measured and simulated data.

3. Results Mean measured force displacement curves for five subjects and the corresponding best-fitting curves are shown in Figure 4. The model parameters estimated using the proposed analytical framework for each of the hyperelastic models are presented in Table 1. As shown in Figure 4, the hyperelastic material models were generally well matched to the experimentally obtained force-displacement curves. The RMS errors between the estimated and measured force curves were 0.16±0.041 N, 0.097±0.024 N, and 0.089±0.025 N (mean ± standard deviation) for the Mooney-Rivlin, generalized Rivlin (polynomial), and Ogden models, respectively. The fitting errors were very small for all of the models. Figure 4 also shows the experimental force-displacement curves and the curves obtained from the FE analyses using the parameters presented in Table 1. A substantial discrepancy between the measured and simulated curves was observed for the Mooney-Rivlin model. However, the force-displacement curves were successfully reproduced by the FE simulation using both the generalized Rivlin (polynomial) and Ogden models. The RMS errors between the computed and measured force curves were 1.83±0.31 N, 0.76±0.37 N, and 0.51±0.17 N (mean ± standard deviation) for the Mooney-Rivlin, generalized Rivlin (polynomial), and Ogden models, respectively. Figure 5 compares the force-displacement curve obtained in the present study with curves reported in the literature. The curve measured in the present study was found to be situated below and behind the curves reported in the literature, indicating that the plantar soft tissue of the feet measured in the present study was softer than in these previous studies.

4. Discussion The present study successfully demonstrated that the hyperelastic material properties of the heel pad can be identified in vivo using the proposed analytical model of spherical contact mechanics and spherical indentation. Previously, inverse FE analysis has been used to identify the material properties of plantar soft tissues (e.g., Erdemir et al., 2006). However, the present study established a framework for identifying parameters without conducting repetitive FE analysis, owing to the introduction of spherical indentation and the proposed contact model, which can deal with the large deformation of a relatively thin hyperelastic substrate. The proposed method is certainly applicable to other locations on the plantar surface, such as the ball of the foot, as well as other body surfaces, such as the thigh or buttock, for possible ergonomic applications. The proposed method is therefore effective for use in the quantification of the material properties of plantar soft tissues, as well as other biological materials. We also demonstrated that the experimentally obtained force-displacement curves could be

accurately represented by the generalized Rivlin (polynomial) and Ogden models, but not by the conventional Mooney-Rivlin model. This may be due to the fact that the Mooney-Rivlin model can only account for relatively small deformation because higher-order terms are not taken into consideration. Therefore, although the model could match the experimentally obtained force-displacement curve where the displacement was small, this was not the case for areas of large deformation. The present study demonstrated that the Mooney-Rivlin model was not sufficient for accurately representing the material properties of the heel pad. For large deformations occurring in plantar soft tissue, the generalized Rivlin model or the Ogden model is more appropriate. The force-displacement curve obtained in the present study was found to be situated lower than those reported in previous studies (Lemmon et al., 1997; Erdemier et al., 2006; Fontanella et al., 2012; Chatzistergos et al., 2015), except for a study by Chokhandre et al. (2012), who identified the material properties based on only one cadaver specimen (Figure 5). This indicates that the stiffness of the heel pad measured in the present study was lower than in previous reports. The reason for this difference may be attributable to the fact that we measured the material properties of in vivo plantar soft tissue in only young living subjects, whereas in the previous studies, the subjects were older (Erdemier et al., 2006; Fontanella et al., 2012; Chatzistergos et al., 2015). Some studies demonstrated that the stiffness of the plantar soft tissue increases with age (Hsu et al., 1998; Kwan et al., 2010). Thus, age appears to be a determinant factor of the properties of the heel pad. Additional studies are needed in order to capture in vivo true variations in the heel pad material properties. An in vitro study by Chokhandre et al. (2012) suggested that the heel pad of cadaver feet may be softer than living feet, but this should also be confirmed through further research. Another possible explanation is that those previous studies, except for Chokhandre et al. (2012), used a flat-top cylinder or ultrasound probe as an indenter, but this study used a spherical indenter to measure the force-displacement relationship of the heel pads. The contact mechanics between a sphere and a half-space is theoretically better characterized than that of a flat-top cylinder or ultrasound probe. The use of non-spherical indenter might have caused more complex stress distribution in the soft tissue in the previous studies, possibly leading to the differences in the estimated material parameters. One possible source of measurement error is the fixation of the foot and leg in space. Since we took measurements from living subjects, perfect immobilization was impossible. However, the foot applied approximately 30 N at maximum to the heel pad, which is actually much less than the gravitational force applied to the leg (ca. 100 N; Winter, 2009). Furthermore, we tried our best to fix the foot to the wooden platform as firmly as possible and used a digital video camera to confirm that there was virtually no vertical displacement of the foot during indentation. Nevertheless, our measurement inevitably incorporated the compliance of the talocrural and subtalar joints in addition to that of the heel pad (Aerts et al., 1995). Isolating the material properties of the heel pad from the

measured data remains as a problem to be solved in future studies.

Acknowledgements We wish to express our sincere gratitude to Dr. Motoharu Tateishi of MSC Software for his continuous support throughout the course of the present study. This study was supported in part by a Grant-in-Aid for Scientific Research (#10252610) from the Japan Society for the Promotion of Science.

Appendix A The three parameters in Eq. (8) were determined using an axisymmetric FE model. The indenter was modeled as a rigid sphere with a diameter of 40 mm and the substrate was modeled as a cylinder with a diameter of 400 mm and a thickness of 200 mm, sufficiently larger than the diameter of the sphere to approximate infinite width and thickness of the substrate. Implicit static analysis of the contact between the sphere and the substrate was conducted from the indentation displacement = 0 to 20 mm at 0.5 mm interval and the relationship between δ and a was obtained. The hyperelastic material parameters were determined based on Erdemir et al. (2006) but we confirmed that essentially the identical results were obtained using the parameters presented by Lemmon et al. (1997), and Chokhandre et al., (2012).

Appendix B When a spherical indenter is pushed into a hyperelastic substrate of finite thickness in a FE simulation, the compressed region takes the approximate form of a paraboloid, and the slope of the tangent plane to the paraboloid at the point at which the sphere and paraboloid intersect (Figure 2; point P) becomes approximately 45 degrees. Therefore, we estimated the change in the volume of the substrate due to indentation by calculating the change in the volume of the paraboloid. The equation of a two-dimensional paraboloid is given as

y = a x2 + b

(1)

where α and β are coefficients. If the paraboloid passes through point P (= a, h) and the slope of the tangent plane to the paraboloid at P is -1 (i.e., 45 deg), then the paraboloid can be represented as

y=-

1 2 a x +h+ 2a 2

(2)

where a is the contact radius, and h is the thickness of the substrate material. Therefore, the volume of the paraboloid V can be analytically obtained as

h h é æ a öù V = ò p x( y ) 2 dy = ò p ê 2a ç h + - y ÷ ú dy 0 0 2 øû ë è

(3)

éæ aö h2 ù = 2p a êç h + ÷ h - ú 2ø 2û ëè The change in volume dV due to the indentation can be approximated as

dV = p a 2 dd where dδ is the instantaneous indentation displacement. Consequently, the volumetric strain can be calculated as d

d

0

0

e v = ò de v = ò

d dV =ò 0 V

p a2 éæ aö h2 ù 2p a êç h + ÷ h - ú 2ø 2û ëè

=

1 d a dd h ò0 h + a

=

1 d 2 Rd - d 2 dd h ò0 h + 2 Rd - d 2

dd

(4)

Appendix C The force applied by the spherical indentation FI can be represented as a function of δ for the Mooney-Rivlin, generalized Rivlin (polynomial), and Ogden models as

ì 20c10 æ e I3 - 3e I2 + 3e I ö 20c01 æ e I3 - 3e I2 + 3e I ö ü ï ï FI = p a 2 í + ç ÷ ÷ý , 2 2 2 ç 3 2 ï ï î 3p (1 -n ) è e I - 2e I + 1 ø 3p (1 -n ) è -e I + 3e I - 3e I + 1 ø þ

æ e I3 - 3e I2 + 3e I FI = p a ç 2 è e I - 2e I + 1 2

(5)

öì 40c20 æ e I3 - 3e I2 ö 40c30 æ e I3 - 3e I2 ö ï 40c10 + ÷í ÷+ ÷ 2 2 ç 2 ç øî ï 3p (1 -n ) 3p (1 -n ) è e I - 1 ø 3p (1 -n ) è e I - 1 ø

, (6)

FI = p a 2

a - -1 é ù 40C 2 - (1 - e I )a -1 ú , e (1 ) ê I 2 3pa (1 -n ) ë û

respectively, where a = R p -(d / R ) qd r and

References

e I = e H + eV

(7)

.

2

ü ï ý ï þ

Aerts, P., Ker, R.F., De Clercq, D., Ilsley, D.W., Alexander, R.M., 1995. The mechanical properties of the human heel pad: a paradox resolved. J. Biomech. 28, 1299–1308. Boulton, A.J.M., Vileikyte, L., Ragnarson-Tennvall, G., Apelqvist, J., 2005. The global burden of diabetic foot disease. Lancet. 366, 1719–1724. Brilakis, E., Kaselouris, E., Xypnitos, F., Provatidis, C.G.., Efstathopoulos, N., 2012. Effects of foot posture on fifth metatarsal fracture healing: A finite element study. J. Foot Ankle Surg. 51, 720– 728. Budhabhatti, S.P., Erdemir, A., Petre, M., Sferra, J., Donley, B., Cavanagh, P.R., 2007. Finite element modeling of the first ray of the foot: A tool for the design of interventions. J. Biomech. Eng. 129, 750–756. Chatzistergos, P.E., Naemi, R., Chockalingam, N., 2015. A method for subject-specific modelling and optimization of the cushioning properties of insole materials used in diabetic footwear. Med. Eng. Phys. 37, 531–538. Chatzistergos, P.E., Naemi, R., Sundar, L., Ramachandran, A., Chockalingam, N., 2014. The relationship between the mechanical properties of heel-pad and common clinical measures associated with foot ulcers in patients with diabetes. J. Diabetes Complicat. 28, 488–493. Chen, W.M., Lee, S.J., Lee, P.V., 2015. Plantar pressure relief under the metatarsal heads: Therapeutic insole design using three-dimensional finite element model of the foot. J. Biomech. 48, 659–665. Chen, W.M., Park, J., Park, S.B., Shim, V.P., Lee, T., 2012. Role of gastrocnemius–soleus muscle in forefoot force transmission at heel rise: A 3D finite element analysis. J. Biomech. 45, 1783– 1789. Chen, W.P., Ju, C.W., Tang, F.T., 2003. Effects of total contact insoles on the plantar stress redistribution: a finite element analysis. Clin. Biomech. 18, S17–S24. Chen, W.P., Tang, F.T., Ju, C.W., 2001. Stress distribution of the foot during mid-stance to push-off in barefoot gait: a 3-D finite element analysis. Clin. Biomech. 16, 614–620. Cheng, H.Y., Lin, C.L., Wang, H.W., Chou, S.W., 2008. Finite element analysis of plantar fascia under stretch-the relative contribution of windlass mechanism and Achilles tendon force. J. Biomech. 41, 1937–1944. Cheung, J.T., An, K.N., Zhang, M., 2006. Consequences of partial and total plantar fascia release: A finite element study. Foot Ankle Int. 27, 125–132. Cheung, J.T., Zhang, M., 2005. A 3-dimensional finite element model of the human foot and ankle for insole design. Arch. Phys. Med. Rehabil. 86, 353–361. Cheung, J.T., Zhang, M., 2008. Parametric design of pressure-relieving foot orthosis using statistics-based finite element method. Med. Eng. Phys. 30, 269–277. Cheung, J.T., Zhang, M., An, K.N., 2006. Effect of Achilles tendon loading on plantar fascia tension

in the standing foot. Clin. Biomech. 21, 194–203. Cheung, J.T., Zhang, M., Leung, A.K., Fan, Y.B., 2005. Three-dimensional finite element analysis of the foot during standing: a material sensitivity study. J. Biomech. 38, 1045–1054. Chokhandre, S., Halloran, J.P., van den Bogert, A.J., Erdemir, A., 2012. A three-dimensional inverse finite element analysis of the heel pad. J. Biomech. Eng. 134, 031002. Erdemir, A., Saucerman, J.J., Lemmon, D., Loppnow, B., Turso, B., Ulbrecht, J.S., Cavanagh, P.R., 2005. Local plantar pressure relief in therapeutic footwear: design guidelines from finite element models. J. Biomech. 38, 1798–1806. Erdemir, A., Viveiros, M.L., Ulbrecht, J.S., Cavanagh, P.R., 2006. An inverse finite-element model of heel-pad indentation. J. Biomech. 39, 1279–1286. Fernandez, J.W., Ul Haque, M.Z., Hunter, P.J., Mithraratne, K., 2012. Mechanics of the foot Part 1: A continuum framework for evaluating soft tissue stiffening in the pathologic foot. Int. J. Numer. Meth. Biomed. Eng. 28, 1056–1070. Field, J.S., Swain, M.V., 1995. Determining the mechanical properties of small volumes of material from submicrometer spherical indentations. J. Mater. Res. 10, 101–112. Fischer, A.A., 1987. Tissue compliance meter for objective, quantitative documentation of soft tissue consistency and pathology. Arch. Phys. Med. Rehabil. 68, 122–25. Fontanella, C.G., Matteoli, S., Carniel, E.L., Wilhjelm, J.E., Virga, A., Corvi, A., Natali, A.N., 2012. Investigation on the load-displacement curves of a human healthy heel pad: In vivo compression data compared to numerical results. Med. Eng. Phys. 34, 1253–1259. Gefen, A., 2002. Stress analysis of the standing foot following surgical plantar fascia release. J. Biomech. 35, 629–637. Gefen, A., 2003. Plantar soft tissue loading under the medial metatarsals in the standing diabetic foot. Med. Eng. Phys. 25, 491–499. Gefen, A., Megido-Ravid, M., Azariah, M., Itzchak, Y., Arcan,M., 2001. Integration of plantar soft tissue stiffness measurements in routine MRI of the diabetic foot. Clin. Biomech. (Bristol, Avon) 16, 921–925. Goske, S., Erdemir, A., Petre, M., Budhabhatti, S., Cavanagh, P.R., 2006. Reduction of plantar heel pressures: Insole design using finite element analysis. J. Biomech. 39, 2363–2370. Gu, Y., Li, J., Ren, X., Lake, M.J., Zeng, Y., 2010. Heel skin stiffness effect on the hind foot biomechanics during heel strike. Skin Res. Technol. 16, 291–296. Gu, Y.D., Ren, X.J., Ruan, G.Q., Zeng, Y.J., Li, J.S., 2011. Foot contact surface effect to the metatarsals loading character during inversion landing. International J. Numer. Meth. Biomed. Eng. 27, 476–484. Guiotto, A., Sawacha, Z., Guarneri, G., Avogaro, A., Cobelli, C., 2014. 3D finite element model of the diabetic neuropathic foot: A gait analysis driven approach. J. Biomech. 47, 3064–3071.

Halloran, J.P., Erdemir, A., 2011. Adaptive surrogate modeling for expedited estimation of nonlinear tissue properties through inverse finite element analysis. Ann. Biomed. Eng. 39, 2388–2397. Herbert, E.G., Pharr, G.M., Oliver, W.C., Lucas, B.N., Hay, J.L., 2001. On the measurement of stress–strain curves by spherical indentation. Thin. Solid. Films. 398–399, 331–335. Hsu, T.C., Wang, C.L., Tsai, W.C., Kuo, J.K., Tang, F.T. 1998. Comparison of the mechanical properties of the heel pad between young and elderly adults. Arch Phys Med Rehabil. 79, 1101– 1104. Hyun, H.C., Lee, J.H., Kim, M., Lee, H., 2012. A spherical indentation technique for property evaluation of hyperelastic rubber. J. Mater. Res. 27, 2677–2690. Isvilanonda, V., Dengler, E., Iaquinto, J.M., Sangeorzan, B.J., Ledoux, W.R., 2012. Finite element analysis of the foot: Model validation and comparison between two common treatments of the clawed hallux deformity. Clin. Biomech. 27, 837–844. Jacob, S., Patil, M.K., 1999. Three-dimensional foot modeling and analysis of stresses in normal and early stage Hansen's disease with muscle paralysis. J. Rehabil. Res. Dev. 36, 252–263. Klaesner, J.W., Commean, P.K., Hastings, M.K., Zou, D., Mueller, M.J., 2001. Accuracy and reliability testing of a portable soft tissue indentor. IEEE Trans. Neural. Syst. Rehabil. Eng. 9, 232–40. Klaesner, J.W., Hastings, M.K., Zou, D., Lewis, C., Mueller, M.J., 2002. Plantar tissue stiffness in patients with diabetes mellitus and peripheral neuropathy. Arch. Phys. Med. Rehabil. 83, 1796– 1801. Kwan, R.L., Zheng, Y.P., Cheing, G.L., 2010. The effect of aging on the biomechanical properties of plantar soft tissues. Clin. Biomech. 25, 601–605. Lemmon, D., Shiang, T.Y., Hashmi, A., Ulbrecht, J.S., Cavanagh, P.R., 1997. The effect of insoles in therapeutic footwear: A finite element approach. J. Biomech. 30, 615–620. Lin, D.C., Shreiber, D.I., Dimitriadis, E.K., Horkay, F., 2009. Spherical indentation of soft matter beyond the Hertzian regime numerical and experimental validation of hyperelastic models. Biomech. Model Mechanobiol. 8, 345–358. Long, R., Hall, M.S., Wu, M., Hui, C.Y., 2011. Effects of gel thickness on microscopic indentation measurements of gel modulus. Biophys. J. 101, 643–650. Luo, G., Houston, V.L., Garbarini, M.A., Beattie, A.C., Thongpop, C., 2011. Finite element analysis of heel pad with insoles. J. Biomech. 44, 1559–1565. Luximon, Y., Luximon, A., Yu, J., Zhang, M., 2012. Biomechanical evaluation of heel elevation on load transfer: experimental measurement and finite element analysis. Acta Mech. Sinica. 28, 232–240. Mattice, J.M., Lau, A., Oyen, M.L., Kent, R.W., Murakami, D., Torigaki, T., 2006. Spherical indentation load-relaxation of soft biological tissues. J. Mater. Res. 21, 2003–2010.

Mooney, M., 1940. A Theory of Large Elastic Deformation. J. Appl. Phys. 11, 582–592. Ogden, R.W., 1972. Large Deformation Isotropic Elasticity: On the correlation of theory and experiment for incompressible rubberlike solids. Proc. R. Soc. Lond. A. 326, 565–584. Pai, S., Ledoux, W.R., 2010. The compressive mechanical properties of diabetic and non-diabetic plantars of ttissue. J. Biomech. 43, 1754–1760. Qian, Z., Ren, L., Ding, Y., Hutchinson, J.R., Ren, L., 2013. A dynamic finite element analysis of human foot complex in the sagittal plane during level walking. PLoS One. 8, e79424. Rivlin, R.S., 1948. Large elastic deformations of isotropic materials. IV. Further developments of the general theory. Philos. T. Roy. Soc. A. 326, 565–584. Rivlin, R.S., 1951. Large elastic deformations of isotropic materials. VII. Experiments on the deformation of rubber. Philos. T. Roy. Soc. A. 243, 251–288. Schröder, J., Neff, P., 2003. Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions. Int. J. Solids. Struct. 40, 401–445. Singh, N., Armstrong, D.G., Lipsky, B.A., 2005. Preventing foot ulcers in patients with diabetes. J. Am. Med. Assoc. 293, 217–228. Spears, I.R., Miller-Young, J.E., 2006. The effect of heel-pad thickness and loading protocol on measured heel-pad stiffness and a standardized protocol for inter-subject comparability. Clin. Biomech. 21, 204–212. Spirka, T.A., Erdemir, A., Ewers, Spaulding, S., Yamane, A., Telfer, S., Cavanagh, P.R., 2014. Simple finite element models for use in the design of therapeutic footwear. J. Biomech. 47, 2948–2955. Tabor, D., 1951. The hardness of metals. Oxford University Press, Oxford. Taljat, B., Zacharia, T., Kosel, F., 1998. New analytical procedure to determine stress–strain curve from spherical indentation data. Int. J. Solids. Struct. 35, 4411–4426. Tani, M., Sakuma, A., Shinomiya, M., 2009. Evaluation of thickness and Young’s modulus of soft materials by using spherical indentation testing (in Japanese). Trans JSME 75, 901-908. Thomas, V.J., Patil, K.M., Radhakrishnan, S., 2004. Three-dimensional stress analysis for the mechanics of plantar ulcers in diabetic neuropathy. Med. Biolog. Eng. Comput. 42, 230–235. Treloar, L.R.G., 1975. The Physics of Rubber Elasticity. Clarendon Press, Oxford. Winter, D.A., 2009. Biomechanics and Motor Control of Human Movement, Fourth Edition. John Wiley & Sons, Hoboken. Wong, D.W., Wang, Y., Zhang, M., Leung, A.K., 2015. Functional restoration and risk of non-union of the first metatarsocuneiform arthrodesis for hallux valgus: A finite element approach. J. Biomech. 48, 3142–3148. Wong, D.W., Zhang, M., Yu, J., Leung, A.K., 2014. Biomechanics of first ray hypermobility: An investigation on joint force during walking using finite element analysis. Med. Eng. Phys. 36,

1388–1393. Yeoh, O.H., 1983. Some form of the strain energy function for rubber. Rubber. Chem. Technol. 66, 754–771. Yu, J., Cheung, J.T., Fan, Y., Zhang, Y., Leung, A.K., Zhang, M., 2008. Development of a finite element model of female foot for high-heeled shoe design. Clin. Biomech. 23, S31–S38. Yu, J., Cheung, J.T., Wong, D.W., Cong, Y., Zhang, M., 2013. Biomechanical simulation of high-heeled shoe donning and walking. J. Biomech. 46, 2067–2074.

Figure legends Figure 1.

Setup for the spherical indentation experiment using heel pads.

Figure 2. Contact between the spherical indenter and a hyperelastic material substrate of finite thickness. The compressed region takes the approximate form of a paraboloid, and the slope of the tangent plane to the paraboloid at the point at which the sphere and paraboloid intersect (point P) becomes approximately 45 degrees. Therefore, in the present study, the change in volume of substrate due to spherical indentation was estimated by analytically calculating the change in volume of the paraboloid, and volumetric strain was computed.

Figure 3. Axisymmetric finite-element model of the spherical indentation of the heel pad.

Figure 4. Mean measured force displacement curves for three subjects and corresponding best fitting curves. The curves obtained from the FE analyses using the identified parameters are also shown for the evaluation of the proposed identification method. The curves obtained from the FE analyses are generally in agreement with the measured curves.

Figure 5. Identified mean stress-strain curves of the heel pads and corresponding curves reported in the literature. The shaded area indicates the standard deviation.

Table 1. Hyperelastic material parameters of the heel pad identified in vivo Model Mooney-Rivlin

generalized Rivlin

Ogden

Parameter

Subject A -4

Subject B -4

Subject C -3

Average 䇲1.12×10-3

c10

䇲8.71×10

c01

4.05×10-3

4.29×10-3

5.75×10-3

4.70×10-3

c10

3.95×10-3

5.72×10-3

4.14×10-3

4.60×10-3

c20

7.03×10-4

0

6.09×10-8

2.34×10-4

c30

1.34×10-3

2.40×10-3

3.18×10-3

2.30×10-3

c

9.06×10-3

12.6×10-3

9.06×10-3

10.2×10-3

α

7.90

7.21

9.01

8.04

4.49×10

䇲2.92×10

PC

Velcro strap㻌 Wooden platform㻌 Spherical indenter

Indentation device

Fig. 1

45㼻

x y

V h a

Heel pad

! ! ddV V P

2R Indenter F

Fig. 2

30 mm

20 mm

Frictionless contact Indentation displacement

Heel pad

15 mm

Rigid indenter

Axis of symmetry Fig. 3

Mooney-Rivlin Subject A

Generalized Rivlin

Ogden Measured Optimized FE FEM

Subject B

Subject C

Fig. 4

Lemmon et al. (1997) Erdemir et al. (2006) Fontanella et al. (2012) Chokhandre et al. (2012) Chatzistergos et al. (2015) Generalized Rivlin Ogden

Fig. 5

Graphical Abstract (for review)