Mae 493 Indentation Test

ARTICLE IN PRESS

Journal of Biomechanics 39 (2006) 1603–1610 www.elsevier.com/locate/jbiomech www.JBiomech.com

Characterization of the mechanical properties of skin by inverse analysis combined with the indentation test Alexandre Delalleaua,b,, Gwendal Jossea, Jean-Michel Lagardea, Hassan Zahouanib, Jean-Michel Bergheaub a

Laboratoire de Tribologie et Dynamique des Syste`mes, UMR 5513, CNRS, ECL, ENISE, 58 rue Jean Parot, 42023 Saint Etienne Cedex 02, France b Institut de Recherche Pierre Fabre, CERPER Hoˆtel Dieu Saint Jacques, 2 rue Viguerie, BP 3071, 31025 Toulouse Cedex 3, France Accepted 12 May 2005

Abstract This study proposes a new method to determine the mechanical properties of human skin by the use of the indentation test [Pailler-Mattei, 2004. Caracte´risation me´canique et tribologique de la peau humaine in vivo, Ph.D. Thesis, ECL-no. 2004-31; PaillerMattei, Zahouani, 2004. Journal of Adhesion Science and Technology 18, 1739–1758]. The principle of the measurements consists in applying an in vivo compressive stress [Zhang et al., 1994. Proceedings of the Institution of Mechanical Engineers 208, 217–222; Bosboom et al., 2001. Journal of Biomechanics 34, 1365–1368; Oomens et al., 1984. Selected Proceedings of Meetings of European Society of Biomechanics, pp. 227–232; Oomens et al., 1987. Journal of Biomechanics 20(9), 877–885] on the skin tissue of an individual’s forearm. These measurements show an increase in the normal contact force as a function of the indentation depth. The interpretation of such results usually requires a long and tedious phenomenological study. We propose a new method to determine the mechanical parameters which control the response of skin tissue. This method is threefold: experimental, numerical, and comparative. It consists combining experimental results with a numerical finite elements model in order to find out the required parameters. This process uses a scheme of extended Kalman filters (EKF) [Gu et al., 2003. Materials Science and Engineering A345, 223–233; Nakamura et al., 2000. Acta Mater 48, 4293–4306; Leustean and Rosu, 2003. Certifying Kalman filters. RIACS Technical Report 03.02, 27pp. http://gureni.cs.uiuc.edu/
grosu/download/luta+leo.pdf; Welch and Bishop, An introduction to Kalman filter, University of North Carolina at Chapel Hill, 16p. http://www.cs.unc.edu/
welch/kalman/]. The first results presented in this study correspond to a simplified numerical modeling of the global system. The skin is assumed to be a semi-infinite layer with an isotropic linear elastic mechanical behavior [Zhang et al., 1994. Proceedings of the Institution of Mechanical Engineers 208, 217–222] This analysis will be extended to more realistic models in further works. r 2005 Elsevier Ltd. All rights reserved. Keywords: Indentation; Inverse methods; Finite elements; Mechanical properties; Kalman filtering; Human skin

1. Introduction Numerous in vivo measurement tests have been designed to determine the behavior of biological tissues. Corresponding author. Laboratoire de Tribologie et Dynamique des Syste`mes, UMR 5513, CNRS, ECL, ENISE, 58 rue Jean Parot, 42023 Saint Etienne Cedex 02, France. Tel.: +33 4 77 42 75 35; fax: +33 4 77 42 75 39. E-mail address: [email protected] (A. Delalleau).

0021-9290/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2005.05.001

These include the suction test (Agache, 2000), the torsion test (Agache, 2000), the compressive test (Pailler-Mattei, 2004; Pailler-Mattei and Zahouani, 2004; Zhang et al., 1994; Bosboom et al., 2001; Oomens et al., 1984; Oomens et al., 1987), and the longitudinal tensile test (Agache, 2000; Khatyr et al., 2004; Retel et al., 2001). Each of these tests makes use of a specific approach related to the type of stress applied. In this paper, we suggest a new method which can be applied to any one of these tests. We present a general

ARTICLE IN PRESS 1604

A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

process to define the mechanical parameters of biological tissues which do not depend on the type of stress applied or on the hypotheses made. Also standard inverse algorithms require phenomenological studies. What we propose here is a self-regulatory method which does not necessitate such time-consuming studies. This study presents the coupling, based on the inverse method, of the indentation test with its numerical simulation, modeled by finite elements using the SYSTUSs software (SYSTUSs, 2004). We do not take into account the non linear anisotropic viscoelastic mechanical characteristics of real skin. We consider skin simply as a linear elastic isotropic medium (Zhang et al., 1994). This paper demonstrates the validity of the inverse method we propose, applied to the hypotheses of isotropy and elasticity, and further works will deal with more realistic skin behaviors. We present different methods, first for the in vivo experiments, then for numerical simulations and finally those for the inverse problem. We then discuss in more detail the extended Kalman filtering (EKF) method (Leustean and Rosu, 2003; Welch and Bishop). After validating the algorithm relative to the problem, we study an actual case of indentation.

Fig. 1. Indentation system and positioning of the subject’s arm (Pailler-Mattei, 2004).

The in vivo indentation test constitutes the first element in the processing sequence. It is therefore particularly important to perform the experiment under adequate and clearly defined conditions (Agache, 2000). This test is conducted using a AISI52100 steel spherical indenter with a radius of 6.22 mm. The experimental curve shows the compression of the forearm of a twentyyear old female subject (Fig. 1) (Pailler-Mattei, 2004). The indenter is controlled by applied strength. The experiment is conducted at low speeds. The viscous components of the skin—related to its particular structure (Oomens et al., 1984)—are not considered in this study. During the indentation, the variation of the contact strength as a function of the indentation depth is measured. Thus we obtain the strength/indentation depth curve, which is then used to deal with the inverse method. The precision of the strength and that of the displacement sensors are respectively 0.1 mN and 0.1 mm for reliable recordings.

precision to obtain reliable results. First, the geometry of the system faithfully follows the size and shape of the indenter used experimentally. However, in order to simplify the problem, we have run these numerical simulations on a semi-infinite medium (Pailler-Mattei, 2004; Zhang et al., 1994). Then, the indenter is assumed to be rigid. The contact between the indenter and the forearm may be treated as a simple sphere–plane contact which makes calculations possible with an axisymmetrical option. Although some studies have measured the friction coefficient between steel and skin (Asserin et al., 2000; Zahouani et al., 2002), we assume this contact to be frictionless. It is also submitted to a penalties method (SYSTUSs, 2004). The mass inertial effects are neglected and four-node quadrangle elements are used (Fig. 2a). In a first approach, the behavior of the skin is assumed to be isotropic, linear, and elastic (Zhang et al., 1994). It follows the Hooke’s law: sij ¼ 2m
ij þ l
kk dij , where sij is the stress tensor,
ij is the strain tensor, l ¼ Eu=ðð1 þ uÞð1  2uÞÞ, and m ¼ E=ð2ð1 þ uÞÞ are the Lame´’s coefficients. The calculations are performed in a large strains/large displacements option, using an updated Lagrangian formulation. The displacements of the indenter are applied step by step. At each step, nonlinear finite elements equations at equilibrium are solved using a Newton–Raphson method. Each modeling shows approximately 2200 nodes and elements, and the CPU time is about 130 s (PentiumIV 3.4 GHz, RAM 2Go, HD 10000 rpm).

2.2. Numerical simulation

2.3. Kalman filtering

Here numerical simulation (SYSTUSs, 2004) is a tool to interpret the experimental results obtained. The entire modeling process should be conducted with the greatest

Kalman filters are generally used in automated systems or robotics (Leustean and Rosu, 2003; Welch and Bishop) for signal processing or for predicting

2. Methods 2.1. The in vivo indentation test

ARTICLE IN PRESS A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

1605

presents the variations of the simulated indentation strength F sim as a function of three parameters: indentation depth h, Young’s modulus E, and Poisson’s ratio u, see Eq. (1). In most cases, the simulated space shows 100  8  6 points ðh  E  uÞ. F sim ¼ f ðh; E; uÞ.

(1)

The experimental curve (Fig. 3) shows the growth of the measured strength F mes as a function of h. At a given indentation depth, hn , there is a corresponding simulated surface F sim which is a function of E and u (Fig. 4). For each point n of the experimental curve ðF mes n ; hn Þ, Kalman’s process can be used to determine the mechanical parameters. This is done on the simulated surface corresponding to hn , by minimizing the difference existing between F mes and F sim n n . The parameters E n and un are calculated using a cubic Lagrange interpolation of the simulated values. In order to follow the experimental curve, different values of E n and un are determined for different values of hn . To assess the parameters corresponding to the proposed curve, the algorithm makes a prediction (a posteriori parameters are noted þ ) as a function of the prediction made in the previous step (a priori parameters are noted  ). This relationship is described using a correction made by the difference between F mes n and F sim n , weighted by a matrix ½K n  called the Kalman gain matrix, see Eq. (2). 

E

þ

 ¼

u

n

E



u

n

þ ½K n ðF mes  F sim n n Þ.

(2)

The matrix ½K n , whose value depends on the gradients of F sim as a function of E and u, takes into account the predictive nature of the model (Fig. 5). Knowing these gradients enables us to anticipate the variations of F for the next indentation step. Fig. 2. Finite element modeling realized. (a) Mesh and geometry; (b) von Mises equivalent stress (MPa), E ¼ 10 kPa, u ¼ 0:4, indentation depth 800 mm.

Indentation strength (mN)

system movements. We use these predictive capabilities to assess the mechanical parameters of our system (Gu et al., 2003; Nakamura et al., 2000). Given the nonlinear growth of the strength as a function of the indentation depth, we will use a linearized form of Kalman filtering: the EKF (Welch and Bishop). In a case of isotropic linear elasticity, the mechanical parameters to be determined are Young’s modulus E and Poisson’s ratio u. In order to deal with the problem, the numerical simulations are conducted first. The simulated space

12 10 8 6 4 2 0 0

200 400 Indentation depth (µm) Fig. 3. Experimental indentation curve.

600

ARTICLE IN PRESS 1606

A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

Fig. 4. Correspondence between the experimental curve (a) and the simulated space (b).

Fig. 5. Processing algorithm.

2.4. Preliminary tests In order to check the relevance of this algorithm and to assess the influence of the assumptions we made, we use a blind test. It consists of values taken from a simulated curve as if they were an experimental one. Simulated experimental curves will be referred to as experimental*. The programming of this algorithm has been done with MATLABs 6.5 release 13. We have also used internal functions of this software to compute the gradients and to do interpolations. 2.4.1. Filter convergence We have tested different configurations to check that this method works. These tests concern different simulated spaces, different values to be determined and different initial conditions. In order to establish a correspondence between the experimental* and the simulated spaces, we have to consider identical indentation

increments. Thus, the strengths measured during the experiments and the simulations can be interpolated in two ways to correspond to the experimental* displacements:

 

By using a cubic Lagrange interpolation. This has the advantage of following the variations of the curve more closely. By using a quartic polynomial interpolation. This process has the advantage of smoothing out the imperfections of the curve, which facilitates processing but harms precision.

Let us consider the following example. The parameters to be determined are E by ¼ 31 kPa and uby ¼ 0:4. The algorithm is initialized by considering E init ¼ 10 kPa and uinit ¼ 0:3. The measurements include 100 recorded points corresponding to 100 different indentation depths. Fig. 6 shows the fluctuations of the mechanical

ARTICLE IN PRESS A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

parameters E n of the experimental* curve calculated for each point using the two interpolation modes cited above. Fig. 7 does the same with un . As expected, the curves approximated with the Lagrange interpolation are much less regular than those corresponding to a polynomial interpolation. The results we have obtained are E det ¼ 32:9 kPa and det u ¼ 0:3001. Fig. 6 shows a slightly crenellated curve due to Lagrange’s interpolation. Fig. 7 shows a difference between the two values. This is due to the polynomial interpolation in which accuracy is not as good as in Lagrange’s interpolation. These differences are not considered to be significant. Young’s modulus reaches plausible values very rapidly. Poisson’s ratio undergoes a slow change which does not allow us to draw conclusions regarding the study. The low variation in Poisson’s ratio can be explained while observing the strength gradients relative to the two parameters to be determined. The indentation strength presents gradients relative to Young’s modulus which induce variations 100 times greater than those relative to Poisson’s ratios. As the results are far from the required parameters and in order to improve this method, a new procedure has been created, the global stabilization process.

Lagrange's interpolation

Quartic polynomial (4) interpolation

Young's modulus (kPa)

33.6 33.5 33.4 33.3 33.2 33.1 33 32.9 0

20

40

60

80

100

Step #

Fig. 6. Evolution of Young’s modulus for the blind-eyes test. The approximations used are ‘‘cubic Lagrange’’ and ‘‘quartic polynomial’’.

Lagrange's interpolation

2.4.2. The global stabilization process This method consists in an initial calculation with the proposed experimental curve (Fig. 8, step 1). The same curve is then continued in the opposite direction back to the point of departure (Fig. 8, step 2). The parameters determined at the end of this second phase will then be considered in a new calculation which follows the initial experimental curve (Fig. 8, step 3). Each global stabilization process therefore includes an odd number of iterations corresponding to the steps described above. The algorithm adapts in order to stabilize to better values (Fig. 9). We must now determine a value which will stop calculations. This value corresponds to a variation of less than 10n of both parameters at the end of each iteration. It is referred to as C ¼ 10n . For lower values the algorithm is assumed to converge. To test the robustness of the method presented we have carried out two different sets of experiments. The first one (Table 1) shows that if we vary the initial parameters, only the number of iterations necessary for the convergence is modified. In each case-study, the parameters are close to those required (E by ¼ 31 kPa, uby ¼ 0:4). The second test concerns the determination of different values from similar initial parameters (E init ¼ 30 kPa, uinit ¼ 0:35) (Table 2). The values thus defined correspond to the reference values, although a margin error still remains which can be reduced by the convergence parameter. The number of iterations can be linked to the parameters to be determined. If these parameters are located in a low-influence area—i.e. where gradients show low values—convergence takes much longer than if they are located in a high influence area. Table 2 illustrates this phenomenon. To determine a high Young’s modulus, fewer iterations are needed because in our case-study a high Young’s modulus implies high gradients. Characterizing both these parameters may appear a little unusual. The analytical solution developed by Hertz proposes an indentation strength which depends

Quartic polynomial (4) interpolation

20 18 Contact Strength (mN)

300.15 Poisson's ratio (x1000)

1607

300.145 300.14 300.135 300.13 300.125 300.12 300.115

16 14 12 1

10

2

3

8 6 4 2

300.11 0

20

40

60

80

100

Step #

Fig. 7. Evolution of Poisson’s ratio for the blind-eyes test. The approximations used are ‘‘cubic Lagrange’’ and ‘‘quartic polynomial’’.

0 0

50

100

150 200 250 300 Number of considered points

Fig. 8. Global stabilization process.

350

400

ARTICLE IN PRESS A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

1608

0.4

31.6

0.395

31.5

Poisson's ratio

Young's modulus (kPa)

31.7

31.4 31.3 31.2

0.39 0.385 0.38

31.1 31

0.375 100

120

140

160

180

100

200

120

140 160 Loop #

Loop #

180

200

Fig. 9. Global stabilization process for convergence after 201 iterations (E by ¼ 31 kPa, uby ¼ 0:4, E init ¼ 10 kPa, uinit ¼ 0:3).

Table 1 Convergences for blind-eyes test with different initial parameters (E by ¼ 31 kPa, uby ¼ 0:4) Initial parameters

E init ¼ 10 kPa, uinit ¼ 0:3

E init ¼ 70 kPa, uinit ¼ 0:3

E init ¼ 10 kPa, uinit ¼ 0:45

E init ¼ 70 kPa, uinit ¼ 0:45

E init ¼ 50 kPa, uinit ¼ 0:2

E init ¼ 10 kPa, uinit ¼ 0:2

E init ¼ 70 kPa, uinit ¼ 0:2

Number of iterations C ¼ 106 Value of E det (kPa) Value of udet

235

237

187

185

309

301

309

31.02

31.02

30.98

30.98

31.02

31.02

0.399

0.399

0.4

0.4

0.399

0.399

30.98 0.399

Table 2 Convergences for blind-eyes test with different parameters to be determined (E init ¼ 30 kPa, uinit ¼ 0:35) Parameters to be determined

E by ¼ 21 kPa, uby ¼ 0:3

551 Number of iterations C ¼ 106 Value of E det (kPa) 20.99 0.3 Value of udet

E by ¼ 21 kPa, uby ¼ 0:4

E by ¼ 41 kPa, uby ¼ 0:3

E by ¼ 41 kPa, uby ¼ 0:4

E by ¼ 61 kPa, uby ¼ 0:3

E by ¼ 61 kPa, uby ¼ 0:4

385

229

167

187

131

21 0.4

40.97 0.3

41.01 0.4

60.98 0.3

61.01 0.4

only on one and only one parameter, the equivalent Young’s modulus E  ¼ E=ð1  u2 Þ. But this solution can only be used in cases of small strains. It gives the system an infinite number of solutions ðE i ; ui Þ. However, to model the problem, we need a large strains formulation which leads to a single solution ðE; uÞ. The analytical solution for complex systems is not known, but thanks to this method, the mechanical parameters can effectively be determined. Let us now study the predictive qualities of the extended Kalman algorithm for the determination of mechanical parameters regarding the behavior of skin tissue in vivo.

(Pailler-Mattei, 2004). The finite elements models correspond to the calculations presented in Section 2.2. The simulated space shows 8 points for Young’s modulus, 6 points for Poisson’s ratio, and 100 points for the strengths. The final simulated space therefore contains 4800 points, which is assumed to be precise enough. The initial values are E init ¼ 10 kPa and uinit ¼ 0:3. In order to avoid problems related to the initialization of the experimental system—beginning of contact, start of data acquisition—the study domain corresponds to an indentation depth of 50–600 mm. The global stabilization process gives the following values: E det ¼ 5:67 kPa;

3. Experimental results

udet ¼ 0:48.

The experimental curve (Fig. 3) plots 100 recorded measurements for an indentation depth of 600 mm

These are obtained after stabilization of the algorithm (Fig. 10) with a convergence parameter of C ¼ 105 .

ARTICLE IN PRESS A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

1609

In order to determine skin tissue parameters with more precision and draw closer to their real behavior (Khatyr et al., 2004; Retel et al., 2001; Holzapfel, 2000; Silver et al., 2001), this study must incorporate more complex behaviors such as viscoelasticity (Oomens et al., 1984). The anisotropy of the tissue is also a major field of investigation (Khatyr et al., 2004). The most significant problem with compressive tests concerns the stress of the ground substance and the underlying tissues such as tendons, muscles, or bones (Oomens et al., 1987). That is why to take into account each phenomenon, and to characterize only the behavior of the skin, the last point to be studied concerns the thickness of the tissue. The method to determine mechanical parameters will also be adapted to multi-layer systems. Then finite

Running time is about 480 s (PentiumIV 3.4 GHz, RAM 2Go, HD 10 000 rpm). We can also observe the robustness of the algorithm while studying variations on determined parameters with different initial values (Table 3). This table shows that, in spite of different initial values, the computation always reaches the same parameters. The number of iterations seems to be related to Poisson’s ratio. Because of its low gradients, the Poisson’s ratio slows down the whole calculation process. These results should be studied for more experimental tests in order to be confirmed.

4. Discussion Simulated curve

Indentation strength (mN)

This study had two primary goals. It was intended to define a method to obtain mechanical parameters for in vivo experiments. It also tested the quality of the algorithm using a simple numerical model—a semiinfinite isotropic elastic linear medium (Zhang et al., 1994). The final simulated curve, whose parameters are described in Section 3, shows that the results do not correspond exactly to the behavior of the tissue (Fig. 11). It was demonstrated that the method converges very efficiently. The number of iterations can be reduced if the algorithm is improved as it will be presented in a future publication. The differences observed arise from simplifications made in the model.

Experimental curve

12 10 8 6 4 2 0 0

100

200 300 400 Indentation depth (µm)

500

600

Fig. 11. Comparison of simulated and experimental curves for the results obtained in Section 3.

6.4 Poisson's ratio

Young's modulus (kPa)

0.5

6.2 6 5.8

0.45

0.4

0.35

0.3

5.6 0

50

100 150 Iteration #

200

250

0

50

100 150 Iteration #

200

250

Fig. 10. Evolution of Young’s modulus and Poisson’s ratio for an experimental indentation. The approximation made is a bicubic Lagrange approximation and the global stabilization process was used.

Table 3 Convergences for the experimental case for different initial parameters E init ðkPaÞ uinit Number of iterations C ¼ 106 e Value of E det (kPa) Value of udet

10 0.3 232 0.375 5.67 0.48

10 0.4 198 0.375 5.67 0.48

20 0.3 232 0.375 5.67 0.48

20 0.4 198 0.375 5.67 0.48

1 0.3 232 0.375 5.67 0.48

1 0.4 198 0.375 5.67 0.48

ARTICLE IN PRESS 1610

A. Delalleau et al. / Journal of Biomechanics 39 (2006) 1603–1610

element modeling could take into account the real geometry of the skin. The thickness of the dermis may be measured in 20 MHz high-frequency ultrasounds and that of the hypodermis in low-frequency ultrasounds. The Optical Coherence Tomography, OCT, could give us information about the epidermis and the stratum corneum. We are currently studying more realistic systems using multi-layer models. Due to the high number of parameters to be determined, try to reduce the solutions to only one is quite a daunting task. Complex mechanical behaviors are also more difficult to identify. Obviously, such models will reflect skin behavior much more realistically. And it will allow us to understand the phenomena resulting from stresses applied to skin tissue as well as the role played by each layer.

Acknowledgments The authors wish to thank C. Pailler-Mattei for experimental data. References Agache, P., 2000. Physiologie de la peau et explorations fonctionnelles cutane´es, collection explorations fonctionnelles humaines, Editions Me´dicales Internationales, Technique & Documentation, ISSN: 1160-0861, ISBN: 2-7430-0360-X. Asserin, J., Zahouani, H., Humbert, P., Couturaud, V., Mougin, D., 2000. Measurement of the friction coefficient of the human skin in vivo. Quantification of the cutaneous smoothness. Colloids and Surfaces B: Biointerfaces 19, 1–12. Bosboom, E.M.H., Hesselink, M.K.C., Oomens, C.W.J., Bouten, C.V.C., Drost, M.R., Baaijens, F.P.T., 2001. Passive transverse mechanical properties of skeletal muscle under in vivo compression. Journal of Biomechanics 34, 1365–1368. Gu, Y., Nakamura, T., Prchlik, L., Sampath, S., Wallace, J., 2003. Micro indentation and inverse analysis to characterize

elastic-plastic graded materials. Materials Science and Engineering A345, 223–233. Holzapfel, G.A., 2000. Biomechanics of soft tissues, Handbook of Material behavior, 12pp., http://www.biomech.tugraz.at/papers/ report7.pdf Khatyr, F., Imberdis, C., Vescovo, P., Varchon, D., Lagarde, J.M., 2004. Model of the viscoelastic behaviour of skin in vivo and study of anisotropy. Skin Research and Technology 10 (2), 96. Leustean, L., Rosu, G., 2003. Certifying Kalman filters, RIACS Technical Report 03.02, 27pp. http://gureni.cs.uiuc.edu/
grosu/ download/luta+leo.pdf Nakamura, T., Wang, T., Sampath, S., 2000. Determination of properties of graded materials by inverse analysis and instrumented indentation. Acta Mater 48, 4293–4306. Oomens, C.W.J., Van Campen, D.H., Grootenboer, H.J., De Boer, L.J., 1984. Experimental and theoretical compression on porcine skin, Biomechanics: Current Interdisciplinary Research, Selected Proceedings of Meetings of European Society of Biomechanics, pp. 227–232. Oomens, C.W.J., Van Campen, D.H., Grootenboer, H.J., 1987. A mixture approach to the mechanics of skin. Journal of Biomechanics 20 (9), 877–885. Pailler-Mattei, C., 2004. Caracte´risation me´canique et tribologique de la peau humaine in vivo, PhD. Thesis. ECL-no. 2004-31. Pailler-Mattei, C., Zahouani, H., 2004. Study of adhesion strengths and mechanical properties of human skin in vivo. Journal of Adhesion Science and Technology 18, 1739–1758. Retel, V., Vescovo, P., Jacquet, E., Varchon, D., Burtheret, B., 2001. Non-linear model of skin mechanical behavior analysis with finite element method. Skin Research and Technology 7, 152–158. Silver, F.H., Freeman, J.W., DeVore, D., 2001. Viscoelastic properties of human skin and processed dermis. Skin Research and Technology 7, 18–23. SYSTUSs version, 2004. User manual, ESI group. Welch, G., Bishop, G., 2004. An introduction to the Kalman Filter, University of North Carolina at Chapel Hill. 16pp., ohttp:// www.cs.unc.edu/
welch/media/pdf/kalman_intro.pdf4. Zahouani, H., Asserin, J., Humbert, P., 2002. Mechanical properties of the skin during friction assessment. CRC Press LLC, Boca Raton. Zhang, M., Turner-Smith, A.R., Roberts, V.C., 1994. The reaction of skin and soft tissue to shear forces applied externally to the skin surface. Proceedings of the Institution of Mechanical Engineers 208, 217–222.