Gaetan Boyer
email: gaetan.boyer@peritesco.com

Jérôme Molimard
email: molimard@emse.fr

Mohamed Ben Tkaya
email: mohamed.bentkaya@peritesco.com

Hassan Zahouani
email: hassan.zahouani@ec-lyon.fr

Marc Pericoi
email: marc.pericoi@peritesco.com

Stéphane Avril
email: avril@emse.fr

) -J Molimard

- M Ben Tkaya

Keywords: extensometer, full field measurement, skin, mechanical properties

Assessment of the in-plane biomechanical properties of human skin using a finite element model updating approach combined with an optical full-field measurement on a new tensile device

Introduction

The skin is the heaviest and the vastest organ of the human body, and certainly one of the most important part. Its function of interface between the body and the external environment is vital. Skin must have very particular properties to permit exchange, for example perspiration, but also to protect the body against different types of attacks (mechanical, chemical or biological). The skin is composed of three main layers, the epidermis, the dermis and the hypodermis with different properties. At a mechanical point of view, skin is a very complex structure: it is an heterogeneous, visco-elastic, anisotropic, adhesive material showing a non linear stress-strain relationship (Agache et al., 2000). Moreover, the skin has a non zero natural tension. The measurement and the knowledge of these mechanical properties in vivo is essential for various domains, for example for clinical (assessment of the evolution of a pathology) or cosmetic (assessment of the efficiency of a product) research. Many different devices, using various techniques, have been developed to evaluate the mechanical state of the skin, like suction [START_REF] Diridollou | In vivo model of the mechanical properties of the human skin under suction[END_REF], contact [START_REF] Jachowicz | Indentometric analysis of in vivo skin and comparison with artificial skin models[END_REF], (Pailler et al., 2008), [START_REF] Boyer | Dynamic indentation on human skin in vivo: ageing effects[END_REF], or contactless indentation (Boyer et al., 2012), torsion (De Rigal, 2002) or extension [START_REF] Thacker | In vivo extensometer for measurement of the biomechanical properties of human skin[END_REF], [START_REF] Vescovo | In vivo tensile tests on human skin: the extensometers[END_REF], (Lim et al., 2008) tests. Extension test principle consists in attaching two pads on the skin and recording the forces induced by a tangential movement of one or both pads.

These previous studies showed that the skin is a visco-hyper-elastic and anisotropic material, with a complex initial tension. A lot of work have addressed in the past hyperelasticity [START_REF] Manan | Determining hyperelastic parameters of human skin using 2D finite element modelling and simulation[END_REF] or viscosity [START_REF] Boyer | Dynamic indentation on human skin in vivo: ageing effects[END_REF], but few work investigated the anisotropy. If considering the whole in-plane behaviour, with a simplification of the anisotropy to orthotropy, the number of variables should be 14 making very difficult a proper characterization. This paper aims at a specific characterization of anisotropy under the hypothesis of linear elastic orthotropy. This is made possible by using a small strain change, representative of skin strain variations during everyday life (10 %). Consequently results will be apparent elastic moduli. As the sought behaviour is anisotropic, the complete evaluation of its properties with this test considering only forces assessment requires several measurements with different directions. Last, in order to have the best contactless evaluation of strain and to improve the understanding of the skin behaviour, it can be very interesting to evaluate the strain field between the pads in addition to the forces during extension test [START_REF] Marcellier | Optical analysis of displacement and strain fields on human skin[END_REF].This paper presents a new device specially developed to this aim, the whole method used to assess the displacement field of the stressed area, and a complete analysis of the results using inverse approach (Finite Element Model Updating).

Methods

Device

The device is shown in figure 1. It is divided into two parts, the "mechanical" part, which aim is to move the pad and to measure forces and displacement, and the "optical" part, whose aim is to visualize and to record images of the deformed area. The whole acquisition chain is described in part 2.3. {figure 1 should be placed here}

Mechanical part

The extensometer is composed of a static and a movable pad attached to the skin by double-sided adhesive. The choice of a static pad permits to measure the force due only to the stressed area.

Indeed, the movable pad measures the force of the stressed area in addition to the force due to the surrounding area. The size of the pads can be 10×10, 20×20 or 20×50 mm 2 . The gap between the pads ranges from 2 to 27 mm. Extension or compression tests can be performed. Each pad is mounted on a piece fixed onto a tangential force sensor, itself fixed onto a normal force sensor (all force sensors are load cells with full Wheatstone bridge). The normal force sensors allow to control the pressure applied by the pads on the skin. This point has been discussed by [START_REF] Lim | New extensometer to measure in vivo uniaxial mechanical properties of human skin[END_REF] and seems to be very important to improve the accuracy of the measurement. Ideally, tests are performed with a normal force oscillating around zero. Remaining oscillations are related to the patient movements (breathing, hearth beat) . The full scale of normal sensors is 3N and the full scale of tangential sensors is 10N. The normal sensor of the static pad is bonded onto the frame, and the normal sensor of the movable pad is bonded onto the movable part of a micrometric motorised translation stage (Polytech Instrument) controlled by computer. This translation stage allows a movement with a velocity ranging from 0.1 to 1.5 mm/s with an accuracy less than 1 µm. The static part of this translation stage is bonded onto the frame. The frame is mounted on a manual translation stage, which permits to control the normal pressure applied by the pads. This manual translation stage is itself mounted on a complete movable frame allowing X,Y,Z movement and rotation using a ball pivot. This frame and a picture of in vivo measurement on the arm are shown in figure 2. {Figure 2 should be placed here}

Acquisition chain

Signals of the force sensors are converted into voltage using two bridge conditioners PMI-301B (PM Instrumentation) and recorded using a 16 bits acquisition card NI 6221-M (National Instrument). Images are recorded through USB-2 port. The translation stage is controlled by a supplier card PI C-843 (Polytech Instrument). In order to have the best synchronisation possible between forces and displacement signals, the position of the stage is directly recorded using internal counter of the acquisition card. A program developed under LabView environment drives the whole system.

Optical part

The image acquisition of the stressed area is performed by an 8 bits camera (The Imaging Source) placed over the pads with high quality optics (deformation less than 0.1%). A micro metric manual translation stage allows to set up the focal length. The field of view is 1280×960 pixels 2 , corresponding to an area of 20×15 mm 2 (15.6 µm/pixel). During the test, images are recorded at a frequency of 10 frames/s. The lighting of the area is very important. It is performed by two groups of three blue (λ= 470 nm) LEDS in order to minimize the diffusion of the light into the skin. Indeed, the lower the wavelength of the light, the less the penetration. LEDS technology gives an illumination with a negligible heating compared with other technologies. Last, the bulk volume size is very low, leading to a very compact system. A laser projects a line on the stressed area in order to measure a profile. This function is not described in this paper.

Full field measurement method

An example of an image recorded by the camera is shown figure 3. It is used to find the displacement field between the first and the last image of the test. The method used is based on Digital Image Correlation (DIC) [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications[END_REF], and has been presented extensively in [START_REF] Molimard | Frequency-based image analysis of random patterns: an alternative way to classical stereocorrelation[END_REF]. So far, only the basic principle is given here. {figure 3 should be placed here}

2D Digital image correlation

Assuming a reference image im0 described by f(r,s), a deformed image im1 of im0 after a small deformation is described by g(r,s) by the following equation:

(1)

where δx and δy are the components of the displacement of im1 and b(r,s) the noise measurement.

A way to find δx ans δy is to maximize function h defined by:

(2)

where * denotes the cross-correlation product. Found values of r and s found correspond to the maximal probability of displacement (δx, δy). This method can be applied in Fourier space using Fast Fourier Transform function, noted FFT2D. Equation 1 becomes:

(3) where the over-line denotes the complex conjugate. Equation 1 and 3 are used at a local scale, on a typical 32×32 zone of interest (ZOI). The procedure is repeated all over the map, giving a final displacement field. Now, classical digital image correlation performs the cross-product either in Fourier or real space.

More refined approaches exist by the way: for example, the signal could be normalized with respect to the mean local intensity, and/or the mean local contrast (Bornert, 2009). The cross-correlation peak is commonly interpolated in order to reach a sub-pixel displacement accuracy. The interpolation function can be a Gaussian or polynomial function. Basically, this choice does not have a strong theoretical basis, and authors generally have an empirical approach. We propose hereafter an alternative to this weak point in the image correlation approach.

Sub-Pixel algorithm

The sub-pixel algorithm is based on the local analysis of ZOI. The random intensity pattern can be described as a sum of periodic functions using a classical Fourier expansion. Then, the random intensity pattern (the speckle signature) can be decomposed into a sum of period patterns (fringes); finally, this approach is a generalization of the grid method, well adapted to small displacement estimation. Because Fourier Transform requires continuity, the ZOI is windowed using a bitriangular function. So far, the algorithm is an extension of the Windowed Fourier Transform (WFT) algorithm proposed by Surrel (Surrel, 1996). As shown Fig. 4, in the frequency domain, each couple of frequencies is characterized by an amplitude and a phase. This phase is proportional to the displacement normal to the corresponding fringe direction. Note that the phase is defined only if the amplitude is higher than zero. Then, in absence of any phase jump, displacements can be related to any defined phase using the relationship:

(4) where is a vector containing all the phase differences, is the displacement vector, and A is the sensitivity matrix. A is a rectangular matrix with 2 rows and as many lines as defined phases, up to 255 in usual conditions.

Displacements can then be derived from eq. 2 using the pseudo-inverse of A. This operation is possible if . In practice, this means that at least two phases along two different directions exist.

(

) 5
Practical problem of this approach is that the signal to noise ratio is weak for each couple of frequencies. Then, the quality of the measurement is obtained by averaging all the available information throw the pseudo-inverse function. A detection of erroneous phases is implemented in order to increase the system resolution.

One should note also that a phase jump can occur in the Fourier domain. Even if a specific treatment should be developed, it sounds better to use first a pixel correlation algorithm. This ensures that the two ZOIs will be as superimposed as possible and that the fringe order is zero for any point and any couple of frequencies. No deformation of the ZOI is proposed here, considering that target applications will be in the small transformation domain.

{figure 4 should be placed here}

Application to the skin

The information used for DIC on human skin is contained in the surface texture. Therefore, the ZOI has to be adjusted to the skin signature [START_REF] Bornert | Assessment of Digital Image Correlation Measurement Errors: Methodology and Results[END_REF], Badulescu et al., 2013]. Considering the global framework of the Nyquist-Shannon sampling theorem, ZOI size should be strictly higher than 2 times the speckle size. Here, ZOI size has been set to 32, which corresponds to a spatial resolution of 8 pixels i.e. 125 µm (80 % radius of the autocorrelation function). Indeed, it isn't our purpose to detect movements at the wrinkle scale but at a more global level. So far, the phenomenon is low-frequency, so a minimization of the spatial resolution is not critical. A repeatability test showed a resolution of 0.03 px, i.e. 0.5 µm.

During the test, the texture changes, with appearance of new furrows in the direction of the stress, and disappearance of furrows perpendicular to this direction, as shown in figure 5.

Identification of mechanical parameters

Modelling skin mechanical behaviour

Results obtained through the experimental set-up are forces and displacement fields: information potentially lead to mechanical properties of the skin, if an appropriate identification strategy is developed. As a matter of fact, because the stress cannot be deduced from simple assumptions, the mechanical test is not a simple one, and an inverse procedure has to be adopted. In the recent years, identification from optical full field technique has been widely used for inert materials, following different kind of behaviour law, for example anisotropic elasticity (Grediac et al., 1998), damage propagation [START_REF] Claire | A finite element formulation to identify damage fields: the equilibrium gap method[END_REF] or plasticity (Avril et al., 2008b). Methods could be adapted for each case. In the present study, we followed first works using Finite Element Model Updating [START_REF] Molimard | Identification of the four orthotropic plate stiffnesses using a single open hole tensile test[END_REF], later developed and optimized by Silva [START_REF] Silva | Identification of Material Properties using FEMU: Application to the Open Hole Tensile Test Applied Mechanics and Materials[END_REF] because of its versatility.

Here, two main problems have to be addressed: the skin behaviour law and the mechanical model.

Skin behaviour law should be described as a visco-elastic, anisotropic, non linear model and/or with pre-tension (Hendricks, 2001). In this work, if noting that test speed should only exhibit elasticity without viscous effects, and considering that the applied load is small enough not to change the instantaneous stiffness, the behaviour law could be linear elastic. It is proposed here to use an orthotropic linear elastic model for the identification procedure. The resulting value can be seen as a secant definition of the elastic parameters, even if the real elastic behaviour corresponds to a more general law. Note that the skin layers are not represented in this description: the chosen test is purely in-plane, thus is not sensitive to changes across the skin thickness, as found using the Classical Laminate Theory (Tsai, 1980).

Finite Element Model

In the FE model, the skin is idealized as a plate (figure 6). Pads are represented by two areas with a controlled displacement, one being zero, the other corresponding to the displacement in the experiment. This implies that we assume zero sliding under the pads; this assumption could be controlled with the image correlation system.

The plate is clamped on its boundaries in order to take into account the load transmission around the static pad. Consequently, the problem is to find the right dimensions around the area of interest (between the pads). Here, the idealized skin is supposed to be a square in order to minimize the number of unknown variables. Then the size of the square is tuned in such a way that the force ratio between the movable and the static pad were as close as possible to the experiment. Using a 100 mm square plate, 1mm thick, a material behaviour corresponding to commonly accepted values for skin (E 1 = 1 MPa, E 2 = 0.5 MPa, G 12 = 0.4 MPa, ν 12 = 0.4) and the angle between the orthotropic frame of reference and the forearm axis α = 20°, the obtained ratio is close to 2.02 when the mean experimental value is 2.08, minimum and maximum values being respectively 1.7 and 2.7 (see "Results" section).

Last, the plate is meshed using quadratic elements outside the zone of interest, and triangular elements inside the zone of interest. This choice is made to optimize the number of elements (and the calculation time) but also to have the same interpolation procedure and meshing elements both in the experimental and in the numerical procedure. Mean element size within the zone of interest is 0.67 mm. {Figure 6 should be placed here}

Identification procedure

The numerical model should be controlled either on displacement or force. If it were controlled on the displacement the cost function would have involved a distance between the calculated force and the experimental force in the cost function, the two distance terms (on displacement and on force) balanced using a penalty parameter. First tests using such a cost function showed the importance of the penalty term, and the great difficulties in tuning this parameter. Consequently, this approach has been discarded, and the control on force is directly ensured by the finite element code.

The cost function is established for N experimental points for each of the for m load cases used in this study:

(6)

where refers to the experimental displacements, and refers to the numerical displacements.

The optimization algorithm used here is a boundary Levenberg Marquardt method, developed by Le Riche [START_REF] Guyon | Least squares parameter estimation and the Levenberg-Marquardt algorithm: Deterministic analysis, sensitivities and numerical experiments[END_REF]. This gradient method uses the Hessian matrix H, the identification parameters being denoted x 1 , … x n .

(7)

The study of this Hessian matrix gives information on the identifiability of the different parameters.

We chose here to identify the compliance matrix, following previous study by Silva (Silva, 2009), but results will be given using engineering constants (E 1 , E 2 , G 12 , ν 12 ) and the angle α for the sake of simplicity. This study clearly shows that the test developed here is insufficient to identify all the target parameters. Table 1 present an indicator of identifiability (ratio of the highest eigenvalue by the lowest eigenvalue of the Hessian matrix) for 4 different cases (traction along 0° direction, 0° and 90°; 0°, 45° and 90°; 0°, 45°, 90° and 135°).

{Table 1 should be placed here} Finally, the starting point influence has been checked. Again, the number of test-cases is important:

it is necessary to use at least 3 test-cases to achieve a successful identification for any starting point.

This results shows that it is necessary to use a rich enough test to identify multiple parameters.

Because of the poorness of the test used here, it is necessary to multiply the test conditions. In the following, 4 directions will be used.

Results

Validation

Force and displacement validation

Different tests have been performed to validate the device and the method. Calibration on force sensors has been performed on the whole device. We found sensitivity of 2.988 ±0.0012 and 2.830 ±0.0013 N/V for the static and movable normal sensors, and 0.492 ±0.0014 and 0.484 ±0.0016 N/V for the static and movable tangential sensors (evaluated with 5 loads). No hysteresis and no drift of these sensitivities have been observed.

Tests using a laser velocimeter have shown a minimal reproducible displacement of 1 µm, and a difference between the theoretical maximal speed is less than 5 % (normal load being 2 N).

Full field measurement method validation

Uncertainty sources are displacement errors due to noise propagation from the camera, algorithm accuracy, camera calibration i.e. pixel size and out-of-plane displacement. Noise propagation and algorithm accuracy have been evaluated through tests on simulated images of the skin with rigid body motion, simulated images of speckle with sine displacement of amplitude 0.2 and 0.02 pixels [START_REF] Bornert | Assessment of Digital Image Correlation Measurement Errors: Methodology and Results[END_REF], and skin images which simulate an extension with a static part at the bottom of the image. The resolution found was 0.03 pixels. Camera calibration is performed with the procedure proposed by (Bouquet, 2000). Last, out-of plane displacement is measured using the laser line. Expected vertical accuracy is 10 µm. Resulting uncertainty on the in-plane displacement field is 0.05 %. Finally, the overall error estimation without out-of-plane displacement has been evaluated as 0.12 % of the measured displacement.

Last, especially for in vivo tests, movements of the subject can also disturb the measurement. One should note also that surface optical signature changes when tension is applied. This point has been discussed previously in "Methods" section.

Experimental results

The in vivo test have been performed on the forearm of 2 groups of healthy young and old women, respectively 23.2 +/-1.6 (10 subjects) and 60.4 +/6 2.4 (group 2, 10 subjects). Only one typical result is presented here for the seek of clarity (26 years old woman). The size of the pads was 10×10 mm 2 , the initial gap was 10 mm, the displacement of the movable pad was an extension of 1 mm with a velocity of 1 mm/s. The reference axis is the axis of the forearm, and the angle between and the axis of the tension test is noted . 10 images were recorded during the test.

Parameters used for the correlation were a size of subset images of 32×32 pixels 2 and a pitch of 4 pixels. In order to have coherent results, the mesh used for the first displacement reconstruction has been used for the others. Displacement fields and maximum load on static and movable pads for each direction are shown in figure 7 and figure 8 using 5% strain (0.5 mm of displacement on movable pad). Furthermore, we can also observe a difference of 41.6 %, 101.8 % and 86.3 % between the forces on static and movable pads respectively for 0°, 45°, 90° and 135° direction which confirm the influence of the surrounding area for the force measurement.

{Figure 7 should be placed here} {Figure 8 should be placed here}

Identification results

Identification has been performed for each test-case, then results have been summarized as a mean value and a variation (at 95 % probability level). Anisotropy has been clearly outlined in Table 2, even if the variability is high: E 2 is 3.7 times higher than E 1 . On the other hand, even if differences exist between the two groups, these differences are not significant using ANOVA procedure.

{Table 2 should be placed here} Now, it is possible to draw the apparent stiffness according to the direction using Classical Laminate Theory (Tsai, 1980). The result is shown figure 9.

{figure 9 should be placed here}

Discussion

Displacement fields are very different depending of the stress direction (figure 7), and usually asymmetric, showing qualitatively a global anisotropy. The punctual force measurement results confirm the anisotropic behaviour of the skin (figure 8). Differences can be observed also between the two groups (young or old), with a higher anisotropy for young people, but also a difference in the orientation of force frames. Only 8 on 32 test-cases are significantly different using ANOVA procedure. Indeed, the differences are too low compared to the variability among the groups.

Finally, only the identification procedure gave quantitative estimate of the stiffness variations with the direction (Table 2).

The anisotropic behaviour has been taken into account using a simple orthotropic model at low strain. The results obtained expressed in term of Young's modulus can be linked to the whole behaviour observed with ageing. Indeed at low strain level the response of the tissue is mainly due to the elastin fibres network which is responsible to the natural tension state of the skin. This tension decreases with ageing. Then, the decrease of E1 and E2 for elder group should indicates a decrease in the pre-tension in both directions with the age: the higher the pre-tension is, the stiffer is the skin. This explanation is in good agreement with previous studies performed with other methods (Boyer et al., 2012). Note that this conclusion is valid only at low strain for which the response is mainly due to elastin network, higher strain inducing non linear behaviour with collagen fibres response.

Based on this observation, one should imagine as a future work, that the anisotropy can be modelled with a non uniform pre-tension field associated with elastic non linear law instead of an orthotropic behaviour model.

Last, the orientation of the stiffer axis with the arm has to be studied in the future. As a matter of fact, the relative angular position of the hand and the elbow is very important in the definition of a reference mechanical state of the skin. Here, the use of a complete positioning system ensures that the position is the same for all the patients (Figure 10), but one should note that forearm is twisted compared to the position at rest. The stiffer axis would probably be aligned with the forearm axis in the rest position.

Compared to previous studies [START_REF] Khatyr | Model of the viscoelastic behaviour of skin in vivo and study of anisotropy[END_REF], the richness of displacement fields makes possible a more realistic mechanical description, in particular by taking into account the skin reaction around the pad. Finally, the approach permits a full identification of all the parameters of the chosen orthotropic behaviour, including in particular the Poisson ratio value.

{Figure 10 should be placed here}

Conclusion

Measuring anisotropy of skin is a key feature to understand the injury recovery and scare formation.

A complete device and method for the in vivo assessment of the displacement fields in addition to the forces during tension test has been presented. Results have been obtained on the forearm with 4 different traction orientations and with two groups of patients, selected according to their age; they clearly show the anisotropy of the skin by analysing 1/ the displacement fields 2/ the forces and reaction forces.

In a second step, this anisotropy has been evaluated using a complete identification procedure based on the Finite Element Updating Method. Given results are apparent elastic moduli, obtained for a 10% variation above the initial natural skin tension. Major stiffness is found to be 3.7 higher than the minor stiffness. The stiffer direction is mainly across the forearm, but this orientation would have be studied better in the future because of the major influence of the arm position. Even if the data are too scattered to give a clear trend, the stiffness is higher for the group of young patients than for the elder one. It has been proposed to explain this phenomenon with the variation of pretension in the skin with the age, in conjunction with the hyper-elastic behaviour of the skin. This assumption will be checked in the near future.

  {Figure 5 should be placed here} In this case, DIC cannot be performed, because information in the reference and deformed images are too different. The way to solve this problem is to perform correlation between successive images recorded throughout the test. If N is the number of images, we obtain then N-1 displacement maps. The problem is to reconstruct the complete displacement fields with the first image as the reference. We use the finite element method as described by(Avril et al., 2008a). The principle of the method is to mesh the area of interest for each map and to search the reconstructed fields and which minimize an error function F between values at the nodes of the mesh and the experimental data. The reconstructed displacement fields can then be interpolated between each node of the mesh and summed to obtain the displacement fields between the first and the final image.

Figure 2 .

 2 Figure 2. Whole system withmounted on a movable frame (left). Measuring head can be rotated to fit with human body necessities, as shown on the right and (in vivo test on the arm (right).

Figure 3 .

 3 Figure 3. View of the camera: static (a) and movable (b) pads, current frame of reference (x,y). Histogram (full scale: 256 grey levels) and Zone Of Interest (ZOI) have been represented on the upper and lower left angles.

Figure 4 .

 4 Figure 4. Basic principle of sub-pixel algorithm.

Figure 5 .

 5 Figure 5. Reference (left) and deformed (right) image of the same area.

Figure 6 .

 6 Figure 6. Numerical model.

Figure 7 .

 7 Figure 7. Displacement fields along and with various tension direction on a woman forearm (young group).

Figure 8 .

 8 Figure 8. Tangential force F (in N) on static and movable pad (mean values for each group).

Figure 9 .

 9 Figure9. Apparent skin stiffness on the forearm according to the direction.

Figure 10 .

 10 Figure 10. Forearm positioning apparatus.

Figure 1 .

 1 Figure 1. Developed device.

Figure 2 .

 2 Figure 2. Whole system mounted on a movable frame (left). Measuring head can be rotated to fit with human body necessities, as shown on the right (in vivo test on the arm ).

Figure 3 .

 3 Figure 3. View of the camera: static (a) and movable (b) pads, current frame of reference (x,y). Histogram (full scale: 256 grey levels) and Zone Of Interest (ZOI) have been represented on the upper and lower left angles.

Figure 4 .

 4 Figure 4. Basic principle of sub-pixel algorithm.

Figure 5 .

 5 Figure 5. Reference (left) and deformed (right) image of the same area.

Figure 6 .

 6 Figure 6. Numerical model.

Figure 7 .

 7 Figure 7. Displacement fields along

Figure 8 .

 8 Figure 8. Tangential force F (in N) on static and movable pad (mean values for each group).

Figure 9 .

 9 Figure 9. Apparent skin stiffness on the forearm according to the direction.

Figure 10 .

 10 Figure 10. Forearm positioning apparatus.

Table 1 .

 1 Identifiability of the 5 parameters for different test-cases, from 1 loading direction to 4 loading directions. Identifiability coefficient corresponds to the capacity to separate the effects from each material parameter.

Table 2 .

 2 Identified results for young and old groups ~ Figure 1. Two views of the dDeveloped device.

Table 1 . Identifiability of the 5 parameters for different test-cases, from 1 loading direction to 4 loading directions. Identifiability coefficient corresponds to the capacity to separate the effects from each material parameter.

 1

	1 test	2 tests	3 tests	4 tests
	258.8	37.49	29.47	37.06

Table 2 . Identified results for young and old groups

 2


published or not. The documents may come

Strain field measurement in orthogonal machining of a titanium alloy

Madalina Calamaz, Dominique Coupard, Franck

Science Arts & Métiers (SAM)

is an open access repository that collects the work of Arts et Métiers ParisTech researchers and makes it freely available over the web where possible. Any correspondence concerning this service should be sent to the repository Administrator : archiveouverte@ensam.eu


Blaise Bourdin

Jean-Jacques Marigo

Corrado Maurini

Paul Sicsic

Morphogenesis and propagation of complex cracks induced by thermal shocks

Keywords: numbers: 46, 15, Cc 62, 20, mt

We study the genesis and the selective propagation of complex crack networks induced by thermal shock or drying of brittle materials. We use a quasi-static gradient damage model to perform large scale numerical simulations showing that the propagation of fully developed cracks follows Griffith criterion and depends only on the fracture toughness, while crack morphogenesis is driven by the material's internal length. Our numerical simulations feature networks of parallel cracks and selective arrest in two dimensions and hexagonal columnar joints in three dimensions, without any hypotheses on cracks geometry and are in good agreement with available experimental results.

Complex crack patterns are ubiquitous in nature and in technology applications. Yet the theoretical understanding and predictive numerical simulation of how and when complex crack patterns arise (nucleation) and how they evolve (crack propagation) is fraught with challenges. Although approaches based on phase fields [1] or variational regularizations [2] have led to significant advance in the numerical simulation of complex crack patterns, short of introducing initial flaws at the structural scale [3], prescribing ad-hoc stress criteria [4], or accepting global energy minimization arguments whose physical relevance is debated [5][6][7], the predictive understanding of crack nucleation is still an elusive goal.

It is well-accepted that while Griffith-like models are appropriate for crack propagation at the scale of a structure, they are inadequate for the modeling of crack nucleation in brittle materials. Arguably, finer models, where a microscopic (material ) length scale plays a fundamental role, are necessary to determine the critical load and crack geometry at the onset, especially in situations where complex crack patterns arise straight from the nucleation. The consistent combined modeling and numerical simulation of crack nucleation and propagation from the material to the structural length-scale is a challenging and largely open issue.

In this Letter, we study the morphogenesis and the selective growth of complex crack patterns induced by material shrinking under thermal shock. We report unprecedented quantitative agreement between numerical simulations, a theoretical model, and experiments at scales spanning from the material internal length to the structural length-scale. Our numerical simulations predict key features of fracture patterns observed in experiments, such as the formation of periodic patterns and the scaling laws governing their selective propagation in two and three dimensions, and do not require any a priori hypotheses on cracks geometry. The method we use leverages recent progress in the understanding of the links between damage models [10,11] and the variational approach to fracture [2,12]. It is based on a rate-independent gradient damage model with stress-softening based on two material parameters: the fracture toughness which rules the evolution of fully developed cracks, and the material's internal length which controls the initial stages of crack nucleation.

We investigate the thermal shock of a brittle ceramics, a now classical experimental setup [3,13,14] where a sample initially at a uniform temperature T 0 is quenched in a cold bath at temperature T 0 -∆T . We consider a rectangular slab Ω exposed to the thermal shock through its thin faces only. We focus first on very thin slabs, which we represent by a two-dimensional body in plane stress. We assume that within the range of temperatures involved, the material properties remain constant. Denoting by u the displacement field and ε = (∇u+∇ T u)/2 the linear strain tensor, we consider for the sound material a linear elastic behavior of energy density

ψ t (ε) = A 0 (ε -ε th t ) • (ε -ε th t )/2
where A 0 is the isotropic elastic stiffness tensor. The inelastic deformation induced by the time-dependent temperature field T t is ε th t = β(T t -T 0 )I, where I is the identity matrix. The index t is meant to highlight the dependence on time. We neglect the cracks influence on heat transfer so that the temperature field T t solves the heat equation ∂ t T tk c ∇ 2 T t = 0 on Ω. Phase changes, non-uniform convection and other nonlinear aspects of the heat exchange between the fluid and the sample are neglected by assuming that the temperature of the domain boundary exposed to the thermal shock is constant and equal to that of the water bath, i.e. T t = T 0 -∆T . Inertial effects are not considered because the diffusion velocity of the temperature field is much slower than the wave speed in the material at the relevant scales in time and space. This hypothesis is universally accepted in the literature on thermal shock problems [3,4,9,14,15]. We model material failure using a gradient damage model characterized by the energy function

E t (u, α) = Ω ψ t (ε) s(α) + G c 4c w w(α) ℓ + ℓ |∇α| 2 dx, ( 1
)
where α is a scalar damage field varying between 0 (sound material) and 1 (fully damaged material), G c is the material's fracture toughness, ℓ an internal length, and

c w = 1 0
w(s) ds a normalization constant. In a timediscrete setting, the quasi-static evolution is obtained by solving at time t i the following minimization problem min u,α≥αi-1 E ti (u, α), where the unilateral constraint on α enforces the irreversibility condition on the damage. The compliance function s and the energy dissipation function w should be chosen such that (1) converges as ℓ → 0 to a Griffith-like energy Ω\Γ ψ t (ε) dx + G c S(Γ), where S is the surface measure of the crack Γ [2,16,[START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF]. In this model, material interpenetration in the fully damaged area is possible. In all the simulations presented here, it can be checked a posteriori that this issue does not present itself. Here, we use s(α) = 1/(1 -α) 2 and w(α) = α, a choice motivated by the convenience of its numerical implementation and specific analytical studies [11,18]. With this choice the damage model has a stress-softening behavior and remains purely elastic without damage until the stress reaches the critical value:

σ c := G c E w ′ (0) 2 c w ℓ s ′ (0) = 3G c E 8ℓ . (2)
The relation above may be used to determine the numerical value of the internal length for a specific material from the knowledge of its elastic limit σ c , Young modulus E, and fracture toughness G c [18]. The present model is in many aspects similar to the phase-field models of fracture developed independently [19]. Those with singlewell dissipation potentials [1,20] are in the form of (1) with

w(α) = c(1 -g(1 -α)), where g(φ) = 4φ 3 -φ 4 .
One significant difference is that while phase-field models typically involve some form of viscous regularizations, our formulation is rate-independent. In addition, the current literature based on phase-field models is concerned only with the propagation of a pre-existing cracks and does not consider the initiation problem.

The dimensional analysis of the energy (1) highlights three characteristic lengths: the geometric dimension of the domain L, the internal length ℓ and the Griffith length ℓ 0 = G c / Eβ 2 ∆T 2 . Using the material's internal length as the reference unit, the problem can be reformulated in terms of two dimensionless parameters, the dimension of the structure L/ℓ (a geometric parameter) and the intensity of the thermal shock ℓ 0 /ℓ (a loading parameter). This is a significant departure from the classical Griffith setting where the only relevant parameter is L/ℓ 0 [4,5,9]. Figure 1 compares the experiment from [8, FIG. 5(d)] (1 mm × 9.8 mm × 50 mm ceramic slab, ∆T = 380 • C) with the damage field from a numerical solution of the gradient damage model. The material properties, communicated by the authors of [8] are E = 340 GPa, ν = 0.22, G c = 42.47 Jm -2 , σ c = 342.2 MPa, and β = 8 × 10 -6 K -1 , which using (2) gives ℓ = 46 µm and ℓ 0 = 14 µm. As our model is rate independent its solution are independent of k c , up to a change of time scale. The numerical results are obtained through a finite element discretization and the approach of [2,18,21]. The main technical difficulties are the constrained minimization of a non-convex energy, and the need for a spatial discretization adapted to the material length-scale scale ℓ. Cracks correspond to the localized bands where α goes from 0 to 1 and back to 0. The qualitative agreement between experiments and simulation is very good. In particular, our simulations reproduce the key phenomenon: the emergence of a periodic array of parallel short cracks at the initiation and their selective propagation toward the interior of the slab. spacing d as a function of the distance a to the edge exposed to the thermal shock for the final configuration, the agreement is striking. Note that in the experimental results shorter cracks are probably filtered out by the adopted experimental crack detection methods [18]. In a first regime, very short equi-distributed cracks nucleate (the plateaus of the crack spacing for short depth in the numerical experiments), followed by selective arrest and period doubling [22]. In the central region of the plot, we can compare experimental and simulation data with a scaling law obtained in [9] through linear fracture mechanics calculations by imposing a bifurcation condition between crack propagation modes with period doubling or not (solid line). For larger values of a, we observe the final crack arrest caused by the finite size of the sample, again in very good agreement with the experiments. Whereas classical theories can be applied in the second and third regimes consisting of fully developed cracks, they cannot properly account for the nucleation phenomenon observed here without preexisting flaws. Our simulations are initialized with a null damage field, an homogenous material, and an unflawed geometry. The crack nucleation is due to the softening character of the material behavior.

The second series of simulations focuses precisely on the crack nucleation process and hence on short times. In this setting, one can assume that the domain is semiinfinite so that the geometric parameter L/ℓ is infinite and the only parameter is the intensity of the thermal shock ℓ 0 /ℓ. For an undamaged material, the stress is uniaxial and reaches its maximum value σ max = Eβ∆T at the surface of the thermal shock. Since (σ max /σ c ) 2 = 3ℓ/8ℓ 0 , for mild-enough thermal shocks (ℓ 0 /ℓ > 8/3), the critical stress is never reached and the solution remains elastic at all time. If ℓ 0 /ℓ < 8/3, damage takes place at t = 0, is homogeneous in the horizontal direction, and non-null in a band of finite thickness D, which penetrates progressively inside the body until a critical time t * . At t = t * the horizontally-homogeneous solution becomes unstable and the damage field develops oscillations of periodicity λ * (FIG. 2(a)).

An analytical solution for the damage field in the first stage of the evolution and its bifurcation and stability analysis is reported in [23], providing semi-analytical results for the periodicity λ * , the damage penetration D * , and the time t * at the bifurcation. Here we perform several simulations varying ℓ 0 /ℓ and detect the critical parameter at the bifurcation. In Figure 2(b) the numerical simulations (dots) are compared to [23] (solid lines). The good agreement provides an excellent verification of our numerical model. For severe shocks (ℓ 0 ≪ ℓ), the results disclose a welldefinite asymptotic behavior with λ * ∼ √ ℓ 0 ℓ, D * ∼ ℓ, and t * ∼ ℓ 0 ℓ/k c . In this regime we observe numerically that all oscillations at the bifurcation develop in fully formed cracks (max α = 1), which is not the case for milder shocks (ℓ 0 ∼ ℓ). However, the full post-bifurcation analysis remains an open problem at this time.

Experimental studies show that in three-dimensions cracks delineate cells with coarsening polygonal crosssections [START_REF] Goehring | Proc. Natl. Acad. Sci. USA[END_REF]. Because of the complexity of the problem, the few available theoretical and numerical studies are based either on simplified two-dimensional models [START_REF] Jagla | [END_REF]26] or on strong assumptions on the crack geometry [27]. The numerical simulation and analysis of the full three-dimensional problem is a major challenge for classical fracture mechanics tools and remains therefore largely unexplored. the computational cost, the computation is performed on a fragment of width 5 mm and height 1 mm of the domain, and the temperature is assumed constant throughout the thickness of the sample. We observe a transition from transverse cracks to three-dimensional fracture patterns delimiting polygonal cells between 480 • C and 580 • C, which is consistent with [8,Fig 5]. Another series of simulations performed at constant temperature for increasing sample thickness (not shown here) highlight the same behavior: transverse cracks for thin domains, transitioning to three-dimensional cracks for thicknesses between 1 mm and 2 mm. This is also consistent with the observations in [8, Fig 5] and justifies the use of a two-dimensional model, a posteriori. Figure 4 shows a fully tri-dimensional crack pattern obtained for a domain of dimension 150ℓ × 150ℓ × 20ℓ for ℓ 0 = 0.05 ℓ. During the simulation, a disordered pattern of small cells nucleates in the first time steps and propagates quasi-statically inside the domain. A selection mechanism leading to honeycomb patterns with increasingly large and regular cell arises from energy minimization. Tracking the propagation of three-dimensional crack front of Figure 4 using a classical Griffith-based model requiring an explicit description of the crack surface and its propagation criterion would be prohibitively complex. Instead, our three-dimensional computations are performed through a straightforward extension of the discretization and minimization algorithm for the energy (1). Obtaining an accurate scaling law for the cell diameter as a function of the depth as the one of the two-dimensional case of Figure 1 would require simulations on larger domains which would rapidly become computationally prohibitive. Hence, we perform a series of numerical experiments by fixing the loading ℓ 0 and by varying at the same time the internal length ℓ and the domain size L so as to keep their ratio equal to L/ℓ = 40. The number of element is kept constant with a mesh size h = ℓ/5. For each computation, we compute the average cell diameter d as a function of distance from the bottom edge a using a post-processing software. This process does not involve any adjustable parameter, yet our results match the two-dimensional scaling law of [9] over several orders of magnitude, leading us to conclude that the scale selection mechanism in two and three dimensions are identical. In addition, while the initial phase of the evolution depends strongly on ℓ, later time evolution of fully developed cracks at the structural is unaffected by this parameter, matching the general scaling law for a Griffith-based model. This finding is consistent with the properties of the energy functional E t which is known to lead to a Griffith-type propagation criterion [START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF]20,28,29].

Our simulations show that a purely quasi-static model based on energy minimization can fully explain the formation of imperfect polygonal patterns and their selective coarsening as "maturation" mechanism during propagation, a phenomenon sometimes attributed to nonequilibrium processes [START_REF] Goehring | Proc. Natl. Acad. Sci. USA[END_REF]. We show that a carefully chosen gradient damage model can be used to account simultaneously for the nucleation of complex crack patterns and their propagation following Griffith criterion. Further works will be carried on to perform a careful statistical analysis of the geometry of the 3D crack patterns and further comparisons to experimental results. The present modeling framework has a general validity and can be applied to other domains including for example the formation of basalt columns with uniform cross-sectional diameters through the solidification of lava fronts [START_REF] Goehring | Proc. Natl. Acad. Sci. USA[END_REF] or shaping of biological systems as observed of the scales on the heads of crocodiles [30]. We are also considering stronger thermo-mechanical coupling including the effect of cracks on heat transfer as in [15,31].

FIG. 1 :

 1 FIG. 1: (color online). Full scale numerical simulation of a ceramic slab submitted to a thermal shock. (a) Damage field from the numerical simulation (blue α = 0, red α = 1). (b) Experimental results from [8, FIG. 5(d)]. (c) Average crack spacing d as a function of their depth a for (a) and (b).The solid line is an approximate scaling law obtained in[9] by imposing a period doubling condition on a Griffith model. Here ℓ = 46 µm is the material internal length, ℓ0 = Gc/ Eβ 2 ∆T 2 = 14 µm the Griffith length (loading parameter), 2 L = 9.8 mm the total depth of the slab. See FIG.2for the meaning of the distributed damage zone.

FIG. 2 :

 2 Figure1compares the experiment from [8, FIG.5(d)](1 mm × 9.8 mm × 50 mm ceramic slab, ∆T = 380 • C) with the damage field from a numerical solution of the gradient damage model. The material properties, communicated by the authors of[8] are E = 340 GPa, ν = 0.22, G c = 42.47 Jm -2 , σ c = 342.2 MPa, and β = 8 × 10 -6 K -1 , which using (2) gives ℓ = 46 µm and ℓ 0 = 14 µm. As our model is rate independent its solution are independent of k c , up to a change of time scale. The numerical results are obtained through a finite element discretization and the approach of[2, 18,21]. The main technical difficulties are the constrained minimization of a non-convex energy, and the need for a spatial discretization adapted to the material length-scale scale ℓ. Cracks correspond to the localized bands where α goes from 0 to 1 and back to 0. The qualitative agreement between experiments and simulation is very good. In particular, our simulations reproduce the key phenomenon: the emergence of a periodic array of parallel short cracks at the initiation and their selective propagation toward the interior of the slab. Figure1(c) shows a quantitative comparison between the numerical simulation of FIG.1(a) and experimental data from[8] by plotting the average crack

FIG. 3 :

 3 FIG. 3: (color online). Three-dimensional version of the experiment from FIG. 1(b) showing the transition from two to three-dimensional crack patterns. The simulations are performed on a subdomain of dimension 5 mm × 1 mm × 1 mm and temperature contrast (a) 380 • C (ℓ0 = .27 ℓ); (b) 480 • C (ℓ0 = 0.17 ℓ); (c) 580 • C (ℓ0 = 0.12 ℓ); (d) 680 • C (ℓ0 = 0.08 ℓ).

Figure 3

 3 Figure 3 is a three-dimensional version of the simulation from FIG. 1 on a plate of thickness 1 mm, for increasing values of ∆T . The fracture geometry is represented by the level surface α = 0.95. In order to reduce

FIG. 4 :

 4 FIG. 4: (color online). Complex fracture pattern for ℓ0 = Gc/ Eβ 2 ∆T 2 = 0.05 ℓ in a domain of size 150ℓ × 150ℓ × 20ℓ color-coded by distance from the bottom surface where the thermal shock is applied. The problem was discretized in 44M linear finite elements in space (mesh size h = ℓ/5) and 100 time steps. The computation was performed on 1536 cores of the NSF-XSEDE cluster Stampede at Texas Advanced Computing Center in 10 h.

FIG. 5 :

 5 FIG. 5: (color online). Average size d (square root of the average cross-sectional area) of the fracture-delimited cells as a function of the depth a (distance to the exposed face) in cubic domains with edge length L ranging from 2 mm to 2 m compared with the two-dimensional scaling law from [9] (solid line). All the simulations are with ℓ0 = 14 µm and ℓ = L/40. Inset: top view to the crack patterns for (a) L = 2 mm; (b) L = 63.2 mm; (c) L = 2 m.

The authors wish to thank Yingfeng Shao for providing the experimental data used in FIG. 1. B.B. work was supported in part by the National Science Foundation grant DMS-0909267. Some numerical experiments were performed using resources of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575 under the Resource Allocation TG-DMS060014. J.J.M. and C.M. gratefully acknowledge the funding of the ANR program T-Shock OTP J11R087.