J Uribe

J Geringer
email: geringer@emse.fr

L Gremillard

B Reynard

Degradation of alumina and zirconia toughened alumina (ZTA) hip prostheses tested under microseparation conditions in a shock device

This paper considers the degradation of alumina and zirconia toughened alumina vs. alumina for hip implants. The materials are as assumed to be load bearing surfaces subjected to shocks in wet conditions. The load is a peak of force; 9 kN was applied over 15 ms at 2 Hz for 800,000 cycles. The volumetric wear and roughness are lower for ZTA than for alumina. The long ZTA ageing did not seem to have a direct influence on the roughness. The ageing increased the wear volumes of ZTA and it was found to have a higher wear resistance compared to alumina.

I.

Introduction

Total hip arthroplasty is continuously improving. New materials for hip joints must be submitted to gait biomechanics. When the hip joint is replaced, the ligament of the femoral head (round ligament) and the extracapsular ligaments are removed. Studies by video fluoroscopy have shown a separation of several millimeters between the cup and head during the swing phase [START_REF] Lombardi | An in vivo determination of total hip arthroplasty pistoning during activity[END_REF][START_REF] Dennis | In vivo determination of hip joint separation and the forces generated due to impact loading conditions[END_REF]. Thus, when the heel touches the ground, a shock of up to 6 times the body weight occurs [START_REF] Bergmann | Hip joint loading during walking and running, measured in two patients[END_REF][START_REF] Jobard | Contribution à l'étude des impacts transmis à l'articulation de la hanche[END_REF]. The microseparation during the gait cycle is shown in Figure 1 [START_REF] Nevelos | Microseparation of the centers of alumina-alumina artificial hip joints during simulator testing produces clinically relevant wear rates and patterns[END_REF][START_REF] Hausselle | Flaws growth modelling in zirconia acetabular cups of hip prostheses under microseparation[END_REF]. Materials should have good wear resistance and good shock resistance.

The bearing surfaces of artificial hip joints are classified into two types: Hard/Soft and Hard/Hard combinations. Metal-on-polyethylene is the more common of hard and soft pairs. This pair is very resistant to shock but less wear resistant, and debris could lead to an inflammatory response and eventually to aseptic loosening. Hard bearing surfaces such as Metal-on-Metal (MoM) have a high wear resistance, but the main disadvantage is the ion release, which has effects that are not well known. Ceramics, used as an alternative for the hard pairs, were investigated for Total Hip Arthroplasty, THA, in the beginning of 1970's. Boutin was the first orthopedic surgeon who implanted this kind of material in artificial hip joints [START_REF] Boutin | Le frottement alumine-alumine en chirurgie de la hanche. 1205 arthroplasties totales: Avril 1970 -juin 1980[END_REF][START_REF] Boutin | Arthroplastie totale de hanche par prothèse en alumine frittée[END_REF][START_REF] Boutin | The use of dense alumina-alumina ceramic combination in total hip replacement[END_REF].

Hip simulator studies, ISO 14242-1, demonstrate that combinations of ceramic on ceramic have a low wear rate compared to other pairs [START_REF] Essner | Hip simulator wear comparison of metal-on-metal, ceramic-onceramic and crosslinked UHMWPE bearings[END_REF][START_REF] Greenwald | Alternative bearing surfaces: The good, the bad, and the ugly[END_REF], which is an advantage. Nevertheless, ceramics are brittle, and some in vivo fractures have been reported. Specific wear zones were highlighted on the head of hip explants [START_REF] Magnissalis | Wear of retrieved ceramic THA components four matched pairs retrieved after 5-13 years in service[END_REF]. The actual process of degradation is related to the wear of alumina. The Low Temperature Degradation (LTD) of zirconia has been related to an increase of the roughness and the production of debris [START_REF] Chevalier | Low-temperature aging of Y-TZP ceramics[END_REF]. The Desmarquet's problem, which occurs with Prozyr® products, has led to huge problems concerning the implantation of zirconia in the 2000's [START_REF] Chevalier | What future for zirconia as a biomaterial?[END_REF]. Afterwards, alumina was the referenced ceramic material that was implanted as hip implants. Thus, much effort has been dedicated to developing new ceramic materials.

Thus, important improvements have been achieved in the biomaterials field, i.e., new ceramics with alumina matrix composites. Zirconia toughened alumina (ZTA) with small grains of zirconia has been developed with a K IC factor higher than those of common ceramics composed of alumina and zirconia ceramics [START_REF] De Aza | Crack growth resistance of alumina, zirconia and zirconia toughened alumina ceramics for joint prostheses[END_REF][START_REF] Chevalier | Reliability assessment in advanced nanocomposite materials for orthopaedic applications[END_REF]. Nevertheless, the phase transformation of zirconia grains remains a concern [START_REF] Deville | Low-temperature ageing of zirconia-toughened alumina ceramics and its implication in biomedical implants[END_REF][START_REF] Chevalier | On the kinetics and impact of tetragonal to monoclinic transformation in an alumina/zirconia composite for arthroplasty applications[END_REF], and studies on avoiding low temperature degradation and ageing of zirconia have been performed [START_REF] Gutknecht | Key role of processing to avoid low temperature ageing in alumina zirconia composites for orthopaedic application[END_REF].

In this paper, we study the resistance of alumina and ZTA ceramics subjected to shocks, i.e., impact loading of a few milliseconds, with a specific device. The influence of ZTA ageing on the wear resistance is also investigated.

II. Materials and methods

Device

A special device was used (Quiri Hydromecanique). The movement of the head is in the vertical direction, the y-axis, and the cup is constrained to move only in the horizontal direction, the x-axis. The applied load was controlled by displacement, which gave more accurate forces than control by force [START_REF] Uribe | Shock machine for the mechanical behaviour of hip prostheses: a description of performance capabilities[END_REF]. The head is positioned on a cone composed of Ti-6Al-4V, and the cup is in a shell cemented to the metallic support. The head and cup were positioned with an angle of 45° to roughly reproduce the anatomical position (Figure 2).

Materials and conditions of tests

The materials tested were 12 ceramic pairs with a diameter 28 mm: 3 alumina balls against alumina cups and 9 zirconia toughened alumina (ZTA, 10 vol. % of zirconia in alumina matrix [START_REF] Gutknecht | Key role of processing to avoid low temperature ageing in alumina zirconia composites for orthopaedic application[END_REF]) against alumina cups. The composition of alumina is in accordance with ISO 6474-1. The pristine microstructures of ZTA could be observed in [START_REF] Gutknecht | Key role of processing to avoid low temperature ageing in alumina zirconia composites for orthopaedic application[END_REF]. It is significant that the zirconia content is under the percolation rate; zirconia grains are dispersed homogeneously in the alumina matrix. The mechanical properties of pure alumina, pure zirconia and ZTA are presented in Table 1. For ZTA, the effect of ageing of the femoral heads was studied. According to [START_REF] Chevalier | On the kinetics and impact of tetragonal to monoclinic transformation in an alumina/zirconia composite for arthroplasty applications[END_REF][START_REF] Chevalier | Low-temperature aging of Y-TZP ceramics[END_REF] and ISO 6474-2, 10 hours of autoclaving at a temperature of 134°C and a pressure of 2 bars correspond to 40 years of in vivo ageing for ZTA of 17 vol. %. In this study, first, three femoral heads did not undergo ageing. Second, six heads composed of ZTA were artificially aged in an autoclave: 3 heads for 3 hours, i.e., approximately 12 years of actual aging. Finally, 3 heads were aged for 20 hours of actual ageing, i.e., approximately 80 years of ultimate ageing. Tests related to shocks were carried out in wet conditions using newborn calf serum as a lubricant; the composition (calf serum heat inactivated, PAA, triple filtering steps at 0.1 µm, protein concentration of 17.5 g.L-1) is the same than the one suggested in ISO 14242-1 [23].

Material

The test conditions are summarized in Table 2. The 3D roughness, i.e., Sa parameter, measured by an optical profilometer (Bruker nanoscope, ex. Veeco, Wyko NT 9100) before the tests ranged between 28 and 34 nm for balls of alumina or ZTA regardless of the ageing time.

To test the ceramic materials in severe conditions, the load was a peak of force of 9 kN, which is equivalent to 6 times the body weight for a patient of 150 kg. The frequency was of 2 Hz, which is the double of the gait (Figure 3). As shown in Figure 3, the force is increased from 0 to 9 kN over roughly 15 ms.

Materials (head/cup)

Protocol for measuring wear volume and weight loss

The wear of the heads was evaluated using two different techniques.

a) Weight loss

Weight loss is a well known and simple method. The balls were weighed before and after the tests in a balance with a resolution of 0.1 mg after carefully drying and removing any debris. Gravimetric measures were taken as a reference for validating the developed profilometer protocol.

b) Optical profilometry

The second technique used optical profilometry. The scheme of the protocol is shown in Fig 4. This method allows non-destructive analysis and is based on scanning white-light interferometry in which a spherical pattern of fringes is formed; finally, the height according to the peaks is measured.

To measure the depth of the wear zone, images over the whole width of the wear stripe are taken. The depth of worn area is calculated by comparing it to the unworn zone. Six images are taken for each wear stripe. The images are processed using the profilometer software (Vision 4.20 2002-2008 Veeco Instruments, Inc). A filter is then applied for the shape and to practice tilt removal. This filter depends on the topography of the worn surface profile. Each 2D profile is exported to be processed in Matlab®. The unworn profile, i.e., the original profile, of the worn area is reconstructed from the worn profile by fitting a quadratic function. It is possible to calculate the worn area by taking the arithmetic difference of both areas. For each image, the wear volume is calculated as the worn area multiplied by the width of the image, L.

Phase transformation

Phase transformation has been reported for retrieved femoral heads, and it has been related to the time of implantation. This phase transformation, from tetragonal to monoclinic, is studied by Raman spectroscopy (Horiba Jobin Yvon HR800, laser λ=514.5 nm) at latitudes of 90°, 45°, 0°, -45°. The method consists of measuring the intensity of the peaks: tetragonal (I t1 and I t2 ) and monoclinic (I m1 and I m2 ). The content of monoclinic phase is then calculated using the Clarke and Adar formula [START_REF] Clarke | Measurement of the crystallographically transformed zone produced by fracture in ceramics containing tetragonal zirconia[END_REF]

+ + + + = (1)
An example of a Raman spectrum for a ZTA ball is shown in Figure 5, and the content of the monoclinic phase is calculated as: Moreover, scanning electron microscopy, JEOL JSM 6400, was practiced on the worn material.

III. Results and discussion

Wear stripe location, width and roughness

At the end of the test, all femoral heads exhibit two wear stripes (Figure 6). For ZTA, the lower wear stripe is no wider than 1 mm. For the alumina femoral heads, the width of the worn zone reached 3 mm. All stripes appear in the contact zone of the head and the rim of the cup near a latitude of 20°. The ratio between the worn and unworn area roughness is approximately 20 times for alumina and 2-3 times for ZTA.

The average roughness measured by profilometry is shown in Fig. 7. The alumina showed a roughness significantly higher than that of ZTA. Ageing of ZTA seems not to have a significant effect on roughness.

Degradation mechanism/ images of the wear stripe

Using the results of optical profilometry, the 3D images of the worn surfaces are shown in Figure 8. The ZTA surface is highly homogeneous. The worn area of alumina is significantly damaged, and it seems to have pits (Figure 8a). The wear stripes of alumina seem to be rougher than that of ZTA.

The SEM observations of the femoral heads show worn surfaces with grains that are pulled out of the alumina (Figure 9). For ZTA with 20 hours of ageing, the grains seem to be 'piled up', and the worn area appears to be covered by compacted debris. At the micron scale, the boundary between the worn and unworn areas was clearly distinguishable in the alumina heads, whereas, for the ZTA heads, the boundary is less well defined.

Phase transformation

The influence of artificial ageing on the wear resistance of femoral heads is investigated. Three couples were tested: not aged, aged for 3 hours and aged for 20 hours. The influence of shocks is described: it is related to the phase transformation. The Raman spectra are taken before and after the test to determine the percent of the monoclinic phase for each step of ageing. The results are shown in Figure 10.

One way or two-way ANOVA tests (p = 0.05) demonstrate that there was no significant difference in the monoclinic content before and after testing. Moreover, they do not show a significant difference according to the duration of ageing. However, Student's t-test (p = 0.05) yielded comparison results between each sample. Each condition of ageing (not aged, 3 hours and 20 hours) is compared with another one, for a total of 14 combinations. First of all, for the sample after 20 hours of ageing, the percentage of monoclinic phase is significantly different in comparison with the sample aged for 3 hours and not aged. Second, femoral heads that are subjected to 3 hours and 20 hours of ageing show different percentages of monoclinic phase after degradation by shocks than before shocks. The surprising point is that the sample subjected to 20 hours of ageing has a lower percentage of monoclinic phase after shocks than before shocks. The discrepancy is greatest for this sample after shocks and after 20 hours of ageing. For the two other samples, i.e., the sample that was not aged and the sample subjected to 3 hours of ageing, shocks increase the percentage of the monoclinic phase. Thus, one might suggest that the crystallographic transformation is promoted by shocks.

Figure 10. Monoclinic content of ZTA femoral heads

Weight loss and wear rate

The results of the wear volume using both the weight loss and profilometry methods are shown in Figure 11. For alumina, the weight loss measurement is higher than the one related to profilometry. Both methods give similar values of the wear volume for ZTA. For ZTA aged 20 hours, the wear volume calculated by profilometry is much higher than weight loss. The choice of the best filter was studied. For alumina balls, a spherical filter is applied within the worn area. For ZTA, the spherical filter showed increased roughness and did not fit the surface shape. A cylindrical filter was then chosen.

The different filters should explain the different ranking between weight loss and profilometry for the alumina head. This point needs further investigation to elucidate the precise role of the filtering process.

The classic weight loss protocol is the simpler and more widely used method for estimating the wear rate of ceramics. However, this technique presents some disadvantages. When the head is removed from the cone, there could be some weight loss due to wear between the cone and the head. The opposite case is also possible: debris accumulation generated by fretting corrosion might increase the mass of the head. Using optical profilometry, it is possible to measure only the wear due to the friction of both surfaces. The wear of alumina/alumina tested in vitro is reported from 0.02 to 1.84 mm 3 / million cycles.

Reference

Wear rate Stewart et al. 2003 [24] 0.2-1.84 mm 3 /million cycles Essner et al. 2005 [START_REF] Essner | Hip simulator wear comparison of metal-on-metal, ceramic-onceramic and crosslinked UHMWPE bearings[END_REF] 0.02-0.08 mm 3 /million cycles Nevelos et al. 2001 [START_REF] Stewart | Severe wear and fracture of zirconia heads against alumina inserts in hip simulator studies with microseparation[END_REF] 1 mm 3 /year Walter et al. [START_REF] Nevelos | Wear of HIPed and non-HIPed aluminaalumina hip joints under standard and severe simulator testing conditions[END_REF] 0.7 mm 3 /year Table 3. Wear rates for alumina from the literature

Comparison of alumina and ZTA

The results for alumina and ZTA are summarized in Table 4. The duration of ageing does not increase the surface roughness according to the ANOVA test (p < 0.05) for the worn and unworn zones of the head. The roughness of the worn zone of the ZTA composite is approximately 5 times lower than that of alumina. As mentioned for the surface roughness, ageing ZTA for 3 hours has no influence on the wear volume. The wear volume of ZTA is 33% lower than that of alumina. ZTA samples, as they are manufactured, can decrease the wear volume under shock degradation according to the experimental conditions of this study. This composite should be relevant for artificial hip joints in the future. No failure of the head occurred during the typical shock tests. The wear volume of ZTA is drastically lower than that of Ultra High Molecular Weight PolyEthylene, UHMWPE, which is used in Metal on Polymer pairs [START_REF] Walter | Edge loading in third generation alumina ceramic-onceramic bearings: stripe wear[END_REF]. These results depend on the manufacturing process (sintering, polishing, etc.), especially for the cone of the head, which is important for the lifetime of the ceramic head. This type of ceramic could have good performance if this step is considered in manufacturing. No cup composed of ZTA was investigated with shock tests or manufactured in this study. For a Ceramic on Ceramic hip joint, cups should be manufactured and tested with a ZTA head for shock degradations. It would be fruitful to compare typical ceramic-ceramic (alumina/alumina) structures with the composites (ZTA/ZTA) available in the orthopedic field. However, a relevant stage is missing in this study. The interactions between the ceramic debris due to shock degradations and bone cells (osteoblasts and/or osteoclasts) should be investigated. The first key point will be the separation of debris, including small particles of ceramic (the order of magnitude should be a few dozens of nanometers). Other points could be related to behavior of cells with debris. Some questions are under study: the impact of these debris on the apoptosis of cells, the inflammation responses, etc. These points will be investigated in the next experimental tests. Finally, simulations by Finite Element Analysis should be performed to give more accurate information about the kinematics and maximum stresses on new design of heads and cups, for example. Moreover, the use of the friction law in finite element analysis has to be developed. A multi-scale analysis is in progress to characterize degradations of the composite under stresses close to the ones applied during shocks between ceramics [START_REF] Kim | Two-dimensional finite element simulation of fracture and fatigue behaviours of alumina microstructures for hip prosthesis[END_REF].

Alumina

IV. Conclusions

The shock resistance of alumina and ZTA heads against alumina cups was investigated. The following conclusions can be drawn:

• The shock device reproduces the in vivo degradation mechanism. The wear stripe location and shape are comparable to cases reported in the literature. For alumina, the wear rate was found to be similar to the values of retrieved femoral heads. Alumina balls have wider wear stripes and a greater increase in roughness within the worn zone (20 times the initial value, i.e., no-worn zone) than ZTA (2-3 times the initial value). The 3D profilometric technique is used for estimating the wear volumes, and it is an original method. The results are promising. This method should be tested with other bearing surfaces for validation, such as Metal on Metal artificial hip joints. • A phase transformation occurred according to the ageing process (Student's t-test) of ZTA heads.

The ageing process has an impact on the percentage of the monoclinic phase; it increases according to the duration of ageing, and shocks promote the transformation. • For wear induced by shocks, ZTA is more resistant than alumina in terms of wear volume. ZTA suffers degradations under no ageing and 3 hours of artificial ageing. For 20 hours of artificial ageing, the ultimate duration, ZTA seems to show weaknesses. In further investigations, experimental results for composites should be compared with retrieved femoral heads. Numerous topics can be developed and investigated with this shock device: heads and cups composed of composites, comparison of wear volume results with explants (according to the schemes of manufacturers), elucidation of the role of ceramic debris on the bone cells and testing new couples, such as Metal on Metal hip joints.

Figure 1 .

 1 Figure 1. Microseparation; distance (d) between the head and cup during the swing phase leading to a shock

Figure 2 .

 2 Figure 2. Shock machine

Figure 3 .

 3 Figure 3. Scheme of the load applied during the gait and the load applied during the shock experiment; BW: Body Weight.

Figure 4 .

 4 Figure 4. Wear volume measured by optical profilometry.

Figure 5 .

 5 Figure 5. Raman spectrum of a femoral head aged for 20 hours.

Figure 6 .Figure 7 .

 67 Figure 6. Location of wear stripes

Figure 8 .Figure 9 .

 89 Figure 8. 3D profiles of a worn area before and after filtering: a) alumina and b) ZTA Figure 9. SEM images of a) the boundary between worn and unworn areas and b) within the worn surface for alumina and ZTA femoral heads

Figure 11 .

 11 Figure 11. Wear volume for alumina and ZTA femoral heads after 0.8 million cycles

K IC (MPa. K I0 (MPa. Ultimate strength (MPa)

					Hardness	Vickers
					(MPa)
	Alumina	4.2	2.4	400-600	1800-2000
	Zirconia	5.4	3.5	1000	1200-1300
	ZTA	6.0	5.0	600	1800
	12Ce-TZP	7.8	5.1	700	1000-1100
	Table 1. Mechanical properties of alumina, zirconia, zirconia toughened alumina, ZTA, and 12 mol. %
	CeO 2 -doped Tetragonal Zirconia Polycrystals, 12Ce-TZP from [21].

Table 2 .

 2 Experimental conditions

		3 alumina/alumina
		3 ZTA/alumina
		3 ZTA, aged 3 hours/alumina
		3 ZTA, aged 20 hours/alumina
	Test conditions	Wet (calf serum)
	Load	Peak of 9 kN at 2 Hz
	Microseparation	1.3 mm
	Test duration	800,000 cycles

Table 4 .

 4 Results for alumina and ZTA; 1: not aged ZTA; 2: aged 3 hours; 3: aged 20 hours.

Acknowledgments

This study was part of the program OPT-HIP. The authors acknowledge the Agence Nationale pour la Recherche, ANR, for financial support and N. Curt for technical collaboration with the shock device. In addition, the authors would like to acknowledge the Region Rhône-Alpes for the grant 'Explora Pro' at Penn State University, PA USA.


Kyungmok Kim

Jean Geringer

Bernard Forest

Fracture simulation for zirconia toughened alumina microstructure

Introduction

Alumina and zirconia have improved the behaviour of artificial articular joints such as hip prostheses. In addition, their chemical inertness and good tribological properties give them an advantage over metals or polymers for a bearing component such as a head of hip prosthesis.

However, these materials maintain drawbacks such as the body reaction around the joints. Wear debris from the materials induces bone resorption, consequent inflammatory reaction and implant loosening. Bone resorption and aseptic loosening are the major problems limiting the lifetime of current implants. Methods improving durability of orthopaedic implants are being explored by proposing new materials. New materials include zirconia toughened alumina (ZTA) microstructured composites. ZTA composites maintain numerous advantages such as high mechanical properties of yttria-stabilized zirconia, high inertness of alumina without any of their respective drawbacks, high hardness for improved wear behaviour, high toughness and crack propagation threshold for higher reliability, and high chemical inertness for longer durability. These properties will allow the design of more reliable and long-lasting implants for less invasive surgery.

In hip prosthesis, material removal occurs resulting from repeated mechanical shocks [START_REF] Hausselle | Flaws growth modelling in zirconia acetabular cups of hip prostheses under microseparation[END_REF]. Dispersion of stabilized ZrO 2 particles in zirconia toughened alumina composites should lead to the improvement of shock resistance in hip prosthesis. Manufacturing processes for ZTA composites include sintering for ceramic densification [START_REF] Magnani | Effect of the composition and sintering process on mechanical properties and residual stresses in zirconia-alumina composites[END_REF].

A typical sintering process for ZTA composites consists of heating from room temperature to 1500 °C, then a thermal treatment of 1500 °C for 2 hours, and cooling up to room temperature.

The mismatch of thermal expansion coefficient between Al 2 O 3 and ZrO 2 grains induces residual stresses in the microstructure of ZTA composite. Particularly, most residual stress is arisen during cooling, directly affecting fracture behaviour under loading conditions.

Several experiments were conducted in order to investigate failure of hip prosthesis [START_REF] Uribe | Finite element modelling of shockinduced damages on ceramic hip prostheses[END_REF][START_REF] Uribe | Shock machine for the mechanical behaviour of hip prostheses: a description of performance capabilities[END_REF][START_REF] Stewart | Severe wear and fracture of zirconia heads against alumina inserts in hip simulator studies with microseparation[END_REF][START_REF] Aza | Crack growth resistance of alumina, zirconia and zirconia toughened alumina ceramics for joint prostheses[END_REF].

Shock experiments were performed for measuring wear rate of alumina heads or cups of hip prosthesis [START_REF] Uribe | Finite element modelling of shockinduced damages on ceramic hip prostheses[END_REF][START_REF] Uribe | Shock machine for the mechanical behaviour of hip prostheses: a description of performance capabilities[END_REF]. The surfaces of alumina heads were progressively locally worn out with respect to number of cycles. Wear and fracture of zirconia heads against alumina cups were investigated with a hip simulator [START_REF] Stewart | Severe wear and fracture of zirconia heads against alumina inserts in hip simulator studies with microseparation[END_REF]. Wear rate between the head and the cup was measured under microseparation conditions. Crack growth resistance of alumina-zirconia ceramics were experimentally determined [START_REF] Aza | Crack growth resistance of alumina, zirconia and zirconia toughened alumina ceramics for joint prostheses[END_REF]. Crack growth rate of alumina toughened zirconia was lower than those of other ceramic biomaterials such as pure alumina or zirconia.

A variety of methods for investigating failure behaviour are used including fracture mechanics, continuum damage mechanics, and so on. Particularly, a method using a cohesive zone law is useful for fracture behaviour of micro-and macro-structures, since this method enables simulating fracture of interfaces between physical parts and characterizing post-yield softening. Elements used in cohesive zone modelling do not represent any physical material but contain cohesive forces arising when they are pulled apart. If these elements satisfy a pre-defined damage criterion, cohesive forces are removed and no stresses remain between physical parts. This traction-separation response in a cohesive zone law can be expressed with bilinear, exponential, power-law, polynomial, or trapezoidal forms. A cohesive zone model was used for assessing for simulating and predicting the onset and the propagation of the delamination in a composite [START_REF] Lampani | Finite element analysis of delamination of a composite component with the cohesive zone model technique[END_REF]. A numerical tool using a cohesive zone law was developed and was applied to finite element method. Some cohesive models were developed for describing fracture of parts of the human body [START_REF] Clech | A model of tension and compression cracks with cohesive zone at a bone-cement interface[END_REF][START_REF] Yang | Fracture length scales in human cortical bone: The necessity of nonlinear fracture models[END_REF][START_REF] Subit | A micromechanical model to predict damage and failure in biological tissues. Application to the ligament-to-bone attachment in the human knee joint[END_REF]. Tension and compression cracks in bone-cement interface were investigated by using cohesive zones [START_REF] Clech | A model of tension and compression cracks with cohesive zone at a bone-cement interface[END_REF]. A cohesive fracture model was formulated, especially for human femoral cortical bones [START_REF] Yang | Fracture length scales in human cortical bone: The necessity of nonlinear fracture models[END_REF]. A micro-mechanical model using cohesive zone theory was developed [START_REF] Subit | A micromechanical model to predict damage and failure in biological tissues. Application to the ligament-to-bone attachment in the human knee joint[END_REF]. Mechanical damages in the ligament-to-bone attachment of a human knee joint were predicted. The cohesive zone model theory was focused on the development of behaviour laws for crack initiation and propagation at an interface within a fibrous material or at the interface between materials. The role of cohesive properties on intergranular crack propagation in brittle polycrystals was studied [START_REF] Shabir | The role of cohesive properties on intergranular crack propagation in brittle polycrystals[END_REF]. The relation between crack path and cohesive law parameters was investigated. Random realizations of a polycrystalline topology were considered instead of using actual microstructures.

Numerical studies for investigating fracture of ceramic microstructures were carried out using finite element methods. A two-dimensional finite element model of compressive fracture in ceramics was developed [START_REF] Warner | Micromechanical finite element modeling of compressive fracture in confined alumina ceramic[END_REF]. Finite element modelling was performed with pure alumina microstructures generated by Voronoi tessellation of randomly positioned seeds. Micro-and macro-scale boundary element modelling was performed with polycrystalline Al 2 O 3 ceramics [START_REF] Sfantos | Multi-scale boundary element modelling of material degradation and fracture[END_REF]. Three-point bending was applied to a model for inducing crack propagation. A twodimensional finite element model was developed for microstructures of pure Al 2 O 3 used for human hip prostheses [START_REF] Kim | Two-dimensional finite element simulation of fracture and fatigue behaviours of alumina microstructures for hip prosthesis[END_REF]. Image-based extended finite element modelling was employed for thermal barrier coatings [START_REF] Michlik | Image-based extended finite element modeling of thermal barrier coatings[END_REF]. Two-dimensional XFEM analysis of actual microstructure containing cracks was developed under assumption of plane stress condition.

However, a finite element model for zirconia toughened alumina composites has not been developed despite its importance. Meanwhile, it may be possible to induce strain singularities close to triple junctions in a microstructure without fine meshing. An extended finite element method for polycrystals with discontinuous grain boundaries was proposed for avoiding singularities [START_REF] Simone | A generalized finite element method for polycrystals with discontinuous grain boundaries[END_REF][START_REF] Menk | Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals[END_REF]. For the purpose of avoiding singularity problem, special enrichment functions describing singular behaviour were implemented in extended finite element model. The developed finite element model using the enrichment function was used for predicting lifetime of solder joints [START_REF] Menk | Lifetime prediction for solder joints with the extended finite element method[END_REF] and for simulating fracture in electronic package applications [START_REF] Cadge | Advancements in fracture and failure simulation for electronic package applications[END_REF].

In dense Al 2 O 3 -10 vol% ZrO 2 microstructures, Al 2 O 3 and ZrO 2 grains vary in shape and size, thereby leading to difficulties in modelling. Image digitalization software and commercial finite element software (ABAQUS ® 6.8) enables modelling such complex grains and grain boundaries.

In this study, actual Al 2 O 3 and ZrO 2 grains were modelled with a method developed by the authors. A two-dimensional finite element model was developed using a cohesive zone law. In this study, fine meshing was employed to avoid strain singularities that may occur close to the triple junctions in an actual microstructure. Fracture and fatigue behaviours of a microstructure were then investigated. Magnitudes of imposed stresses were higher than those found at the microscopic scale. Body fluid exists in actual hip prostheses, constituted by proteins (Ex. globulin, albumin), salts as chloride sodium, chloride calcium, and so on [20]. But, in this study, dry contact condition without the body fluid was considered between a femoral head and an acetabular cup; dry contact could occur after disruption of the lubrication film during human gait.

Finite element simulation

A bilinear, time-independent cohesive zone law was used for describing fracture of grain boundaries in this study. The cohesive zone law with traction-separation response allows specification of mechanical properties such as stiffness, strength, and relative displacement at failure [START_REF]ABAQUS 6.8 Documentation[END_REF]. The bilinear cohesive response is determined with maximum nominal stresses (normal and shear) and displacements (normal and shear) at failure. Once the maximum displacement of a cohesive element exceeds pre-described failure magnitude, stresses of a cohesive element become fully released. The bilinear cohesive zone law was illustrated in Fig 1 (right) [START_REF] Kim | Two-dimensional finite element simulation of fracture and fatigue behaviours of alumina microstructures for hip prosthesis[END_REF].

[Fig. 1] Fig. 1(left) presents a proposed simulation algorithm for determining damage of a microstructure. A two-dimensional finite element model including grains and cohesive zones is generated by digitalizing an actual microstructure image. Material properties of grains and cohesive zones are then defined. Loading and boundary conditions similar to those found at hip prosthesis are applied to models. Plane strain, implicit analysis is performed in the simulation.

During simulation, stresses and strains of grain elements are calculated along with damage variables of cohesive elements. Fracture of grain boundaries is evaluated in terms of the damage variable. When the damage variable of an element is equal to unity, the element is deleted and the deleted element is considered as crack.

Commercial finite element software (ABAQUS ® 6.8 standard) was used for simulating models. Each model was simulated on a cluster with 10 calculation nodes (Intel Xeon Quad Core 3 GHz, 64 bits with 1 GB). Viscous stabilization was implemented to avoid severe convergence arising from softening behaviour and stiffness degradation, as recommended by ABAQUS [START_REF]ABAQUS 6.8 Documentation[END_REF].

Dissipated energy fraction specified for viscous stabilization was 0.002~0.005 (0.2~0.5%) during calculation.

[Fig. 2] Fig. 2a and2d show the microstructure images of Al 2 O 3 -10 vol% ZrO 2 ceramics [START_REF] Gutknecht | Key role of processing to avoid low temperature ageing in alumina[END_REF]. [Table 1]

[Fig. 3] [Table 2]
In this study, conditions, found at the contact surface between a head and a cup, were taken into account. That is, the selected microstructures were assumed to be placed at contact surface.

Because initial surface roughness was about 0.02 µm, the surface roughness in the microstructures was ignored. Contact stresses between a head and a cup in hip prosthesis were reproduced by applying normal and shear stresses to the top surface of a model. The bottom surface was fixed in all directions while loading was being applied. The left and right sides of a model were not fixed, thus the model was permitted to undergo horizontal and vertical deformations.

Little is known about micro-mechanical properties of ZTA composite. Thus, the following assumptions were made in this study.

• The magnitude of normal stress varies according to location within an actual contact surface due to surface roughness, ranging from 5 MPa to 500 MPa [START_REF] Uribe | Finite element modelling of shockinduced damages on ceramic hip prostheses[END_REF]. In this study, four normal stress magnitudes of 50, 100, 250 and 500 MPa were selected for investigating the effect of contact stress.

• Dry contact condition and room temperature were considered. Coulomb friction coefficient of the material is about 0.4 under these conditions [START_REF] Kalin | Wear and friction behavior of alumina ceramics in aqueous solutions with different pH[END_REF][START_REF] Stachowiak | Fretting wear and friction behaviour of engineering ceramics[END_REF]. Shear stress was determined according to the magnitude of normal stress and friction coefficient (i.e. shear stress = 0.4 × normal stress).

• Grains behaved elastically as shown Table 1.

• The magnitude of elastic modulus was differently defined according to the direction. E x was ten percent lower than a typical value of elastic modulus of the material (374 MPa for Al 2 O 3 and 210 MPa for ZrO 2 ), and E y was ten percent greater than the value.

• Poisson ratio values of Al 2 O 3 and ZrO 2 grains were defined as 0.22 and 0.3, respectively.

• Each grain in Al 2 O 3 -10 vol% ZrO 2 microstructure maintained a differently assigned material orientation, defined by an angle subtended from the horizontal axis as shown in Fig. 4. Each angle was chosen according to the shape of a grain. That is, x axis was placed in the longitudinal direction of a grain.

• The mechanical property of cohesive layers was different according to grain boundaries, since selected microstructures have two sorts of grain boundaries: Al 2 O 3 vs. Al 2 O 3 and ZrO 2 vs. Al 2 O 3 . Mechanical properties of cohesive zones are shown in Table 2.

• Fracture energy release rate of the layer between Al 2 O 3 grains was equal to 0.5 J m -2 in the pure tensile mode. This fracture energy release rate was uniformly distributed to all the cohesive layers between Al 2 O 3 grains. The strength magnitude of a cohesive layer defined here is similar to those found in the literature [START_REF] Warner | Micromechanical finite element modeling of compressive fracture in confined alumina ceramic[END_REF][START_REF] Yasuda | Twist angle dependence of interfacial fracture toughness of (0001) twist boundary of alumina[END_REF]. Fracture energy release rate of the layer between ZrO 2 and Al 2 O 3 grains was higher than the value between Al 2 O 3 grains.

• The nominal shear stress was determined on the basis of the grain boundary shear strength's dependence on the compressive yield strength (3000 MPa) similar to friction coefficient (0.13) of a metallic glass [START_REF] Lund | The Mohr-Coulomb criterion from unit shear processes in metallic glass[END_REF].

Results

Stress distribution in a microstructure after cooling

A sintering process of Al 2 O 3 -10 vol% ZrO 2 leads to ceramic densification, resulting to the increase of mechanical strength. After a sintering process, residual stresses remain due to thermal expansion mismatch between Al 2 O 3 and ZrO 2 . In this simulation, linear thermal expansion coefficients of grains were defined as shown in Table 3.

[Table 3] [Fig. 4]
After cooling a microstructure from 1500 °C to 25 °C, stresses of grains in a microstructure were determined. The magnitudes of stresses vary according to locations of the microstructure due to local stress concentration, ranging from -500 MPa to 700 MPa. Thus, the average stress value of grains was computed as Fig. 4. It was identified that ZrO 2 grains were subjected to tension while Al 2 O 3 ones underwent compression in both models. Small shear stresses of both grains remained after cooling. The average magnitude of the compressive stress in Al 2 O 3 grains was lower than the value (150 ± 50 MPa) found in the literature [START_REF] Bartolome | Neutron diffraction residual stress analysis zirconia toughened alumina (ZTA) composites[END_REF]. The alumina-zirconia composite considered in the literature contained 22 percent zirconia in volume, greater than the amount of zirconia in this study.

Crack distribution at various static loadings

A sintering process plays a role in minimising porosity by densificating a ceramic structure. However, voids may occur unpredictably after manufacturing processes. The existence of voids leads to the decrease of mechanical strength, thereby inducing micro-cracks. The effect of a void was investigated by modelling a microstructure containing a single void. Finite element models with a void and without a void were simulated under the same loading and boundary conditions.

All models were sintered prior to being loaded.

[Fig. 5] [Table 4] Fig. 5 shows crack density versus applied normal stress chart. Crack density was defined as the ratio of total crack length to the total length of grain boundaries in a model. Here, the total crack length was defined as the sum of longitudinal distances of cracks generated in a microstructure.

The total crack length was measured after a loading-unloading cycle. It was apparent that crack density values of microstructures with a single void were higher than those without a void. Models A and B were simulated under four different loading conditions. In both models, the crack density gradually increased with increasing the magnitude of applied stress within range of 50-500 MPa. A notable result was that the difference of crack density existed between two models though crack density rates were similar as shown in Table 4. This difference shows that total crack length is affected by grain shape and arrangement.

[Fig. 6] [Fig. 7] Fig. 6 illustrates crack distribution of model A under different loading conditions. At a normal stress of 100 MPa, a few small cracks (< 0.1 µm) occurred at the upper-left side of the model. At 500 MPa, long cracks were observed at the upper side of a model (Fig. 6b). The model with a void maintained more cracks, compared with that without a void (Fig. 6c). In addition, locations where the cracks occurred were different. At 500 MPa, long horizontal cracks were found at the model with a void, thereby making a grain separated from the model (Fig. 6d). Fig. 7 illustrates crack distribution of model B under two loading conditions. At a normal stress of 100 MPa, small cracks (< 0.2 µm) were observed at a variety of locations. At 500 MPa, long cracks were found at the upper side of the model, but grains were still adhered (Fig. 7b). Fig. 7c and Fig. 7d show that longer cracks were generated due to a void, compared with results without a void. These results show that crack length as well as crack path is somewhat different according to a microstructure.

Crack distribution after cyclic loadings

[Fig. 8] [Fig. 9]
Fatigue behaviour of the sintered models was investigated by applying repeated contact stresses.

During simulation, peak values of the normal stress and the shear stress were 100 MPa and 40 MPa, respectively. Fig. 8 illustrates crack distribution of model A with a void. After the initial fatigue cycle, small cracks occurred at the upper-right side of the model. These cracks grew with increasing number of fatigue cycles. Long cracks were then generated by combining small ones.

After the tenth fatigue cycle, grain boundaries were seriously damaged as a result of progressive crack growth. A few grains were found to be separated. Fig. 9 shows crack distribution of model B with a void. It was also observed that small cracks occurred after the initial cycle and grew with respect to number of cycles. After the fifth cycle, most cohesive elements around an upperright grain were removed. It is demonstrated that this proposed method allows simulating progressive crack growth by applying repeated contact stresses to a microstructure.

[Fig. 10] Fig. 10 shows crack density evolution of the models. The slope of the plot corresponds to crack density rate. Each crack density rate can be expressed as a linear function of the number of cycles. The crack density rate of model B was higher than that of model A. In addition, it was found that a void led to the increase of crack density rate. This result shows that crack density rate varies according to a microstructure.

Discussion

Material loss rate in a model with a void can be approximately evaluated with the amount of separated grains with respect to number of cycles. If the amount of cohesive elements remaining between grains is so small that the grains can be considered to be separated in Fig. 8c, the loss area in model A is approximately 1.0×10 -6 mm 2 after ten cycles (the cross section of the material is 25 µm 2 ). Supposing that grain thickness is 0.001 mm, the loss rate is determined as 0.1×10 -9 mm 3 /cycle, i.e. 10 -4 mm 3 /million of cycles for 25 µm 2 . Meanwhile, it was identified from shock experiments that the actual cross section of a head submitted to shocks was approximately 5 µm x 3 mm [START_REF] Uribe | Shock machine for the mechanical behaviour of hip prostheses: a description of performance capabilities[END_REF]. Thus, the loss rate is equal to 1.5 mm 3 /million of cycles. This wear rate is similar to the one measured by Stewart, i.e. 1.84 mm 3 /million of cycles [START_REF] Stewart | Severe wear and fracture of zirconia heads against alumina inserts in hip simulator studies with microseparation[END_REF]. Meanwhile, in the authors' shock degradation experiment with ZTA composite, it was identified that wear volume is equal to 0.12 ± 0.01 mm 3 /million of cycles [START_REF] Geringer | Degradation of alumina and zirconia-toughened alumina (ZTA) hip prostheses tested under microseparation conditions in a shock device[END_REF]. From the microscopic to the macroscopic scale, our calculation is an order of magnitude. The calculated loss rate corresponds to a local wear rate at contact surfaces of hip prosthesis, since microstructures are used. For more accurate calculation of average wear volume, various actual microstructures found at the contact surface needs to be used. It might be necessary to compare calculated wear volume with clinical observations of zirconia toughened alumina. Recently, insertion of zirconia toughened alumina into patient's body is begun. Thus, clinical observation of zirconia toughened alumina microstructures is very limited. Zirconia toughened alumina for hip prosthesis is expensive compared with pure alumina. This proposed method is also useful in the design and manufacture process of materials. It is possible to perform direct comparison among microstructures in terms of fracture. This developed simulation could decrease the cost of zirconia toughened alumina.

In this simulation, dry contact condition was considered, thus full fluid film lubrication contact condition needs to be modelled. In addition, external stresses were applied to the models for reproducing deformation of a microstructure at the contact surface of a head. Displacement could be imposed to the models instead of the shear stress; this modelling can be achieved after identifying actual relative displacement between a head and a cup.

Conclusions

This paper describes finite element modelling for fracture and fatigue behaviours of zirconia toughened alumina (ZTA) ceramic at the micro level. ZTA ceramic is a promising material for hip prosthesis due to its superb mechanical properties. Nevertheless, modelling of fracture behaviour of the ceramic is rather limited. In this paper, a two-dimensional finite element model was developed for investigating crack initiation and propagation of the microstructure of ZTA ceramic. Two different Al 2 O 3 -10 vol% ZrO 2 microstructures were selected and modelled with commercial finite element software and an image digitalization program developed by the authors. Simulation algorithm from model generation to analysis was detailed, enabling cyclic loadings. A bilinear, time-independent cohesive zone law was implemented for describing fracture behaviour of grain boundaries. The following was carried out with this simulation.

First, residual stress modelling was performed with the microstructure. Residual stresses were determined arising from the mismatch of thermal expansion coefficient between alumina and zirconia grains. Results showed that compressive stresses remained in alumina grains and tensile stresses remained in zirconia ones.

Second, the effect of contact stresses in the microstructure was investigated in terms of the crack density. Normal stress and shear stress were applied to the top surface of the models. The total crack length in a model was then measured. It was identified from results that crack density was gradually increased with increasing the magnitude of the stress. In addition, void effect was identified by removing a single grain in a microstructure. It was observed that a void in a microstructure led to the increase of crack density.

Finally, fatigue behaviour of a microstructure was investigated by performing simulations under cyclic loadings. It was apparently observed that crack density linearly increased with respect to number of cycles. Meanwhile, difference of crack density rates existed between selected microstructures. It was also identified that a void led to the increase of crack density.

In conclusion, this proposed finite element simulation offers an effective method for modelling residual stress and for investigating fracture behaviour of ZTA microstructures. Further investigation will focus on three-dimensional finite element modelling and modelling and consideration of full fluid film lubrication condition. Continuous lines denote curve fits.

Figure captions

T n max T s max o δ n o δ s f δ n f δ s Normal Shear f δ ) ( ) ( D o f max o max f max s max n max δ δ δ δ δ δ δ δ δ - - = + = =         +         1 T T T T 2 max s s 2 max n n Bilinear cohesive zone law N i < N t

  Fig.2a and 2dshow the microstructure images of Al 2 O 3 -10 vol% ZrO 2 ceramics[START_REF] Gutknecht | Key role of processing to avoid low temperature ageing in alumina[END_REF]. Fig.2b

Fig. 1

 1 Fig. 1 Simulation algorithm for determining damage. N i denotes the i th cycle.

Fig. 2

 2 Fig. 2 Images of Al 2 O 3 -10%ZrO 2 microstructures [14] and finite element models. The middle

Fig. 3

 3 Fig. 3 Schematic representation of local coordinates of (a) model A and (b) model B. Grey areas

Fig. 4

 4 Fig. 4 Magnitude of the average stress of Al 2 O 3 and ZrO 2 grains after cooling. (a) Model A and (b) model B. x σ and y σ denote normal stresses in x and y directions shown in Fig.4,

Fig. 5

 5 Fig. 5 Crack density versus applied normal stress chart. Shear stresses were applied with the

Fig. 6

 6 Fig.[START_REF] Aza | Crack growth resistance of alumina, zirconia and zirconia toughened alumina ceramics for joint prostheses[END_REF] Illustrations of crack distribution in model A with respect to magnitude of applied contact

Fig. 7 Fig. 8

 78 Fig. 7 Illustrations of crack distribution in model B with respect to applied contact stress: (a) a

Fig. 9

 9 Fig. 9 Illustrations of crack distribution in model B with a void with respect to number of fatigue

Fig.

  Fig. 10 Evolution of crack density with respect to number of fatigue cycles. Markers denote

Fig. 1

 1 Fig.1Simulation algorithm for determining damage. N i denotes the i th cycle.

Fig. 2 Fig. 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8

 245678 Fig. 2 Images of Al 2 O 3 -10%ZrO 2 microstructures [14] and finite element models. The middle images (b, e) are models generated with the actual images (a, d). For convenience, (b) and (e) are denoted as model A and model B, respectively.A single alumina grain in each model was deleted for investigating the effect of a void (c, f).

Fig. 9 Fig. 10

 910 Fig. 10 Evolution of crack density with respect to number of fatigue cycles. Markers denote calculated data.Continuous lines denote curve fits.

  10 Evolution of crack density with respect to number of fatigue cycles. Markers denote calculated data. Continuous lines denote curve fits.

	Figures
	Digitalize a microstructure
	image
	Generate grain elements	Stress (T)
	Add cohesive elements
	between grains
	Define material properties,
	total number of cycles(N t ),
	loading and boundary
	conditions
	Calculate damage variables
	(D), stresses and strains
	Yes
	i=i+1
	No
	Output results


Alexandr Klimchik
email: alexandr.klimchik@mines-nantes.fr

Yier Wu
email: yier.wu@mines-nantes.fr

Stéphane Caro
email: stephane.caro@irccyn.ec-nantes.fr

Benoît Furet

Anatol Pashkevich
email: anatol.pashkevich@mines-nantes.fr

Benoit Furet
email: benoit.furet@univ-nantes.fr

Advanced robot calibration using partial pose measurements

Keywords: Robot calibration, parameter identification, partial pose measurement, accuracy improvement. I

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

To achieve the desired accuracy, each industrial robot must go through the calibration procedure, which deals with proper parameter tuning of the mathematical model embedded in the robot controller. Because of its importance, the problem of robot calibration has been in the focus of research community for many years and has been studied from different aspects [START_REF] Stone | Kinematic Modeling, Identification, and Central of Robot Forward Manipulators[END_REF][START_REF] Hollerbach | A Survey of Kinematic Calibration[END_REF][START_REF] Mooring | Fundamentals of Manipulator Calibration[END_REF][START_REF] Elatta | An Overview of Robot Calibration[END_REF][5]. In spite of this, the issues of the identification accuracy and calibration error reduction have not found enough attention, only limited number of works directly addressed these important problems [START_REF] Mooring | Fundamentals of Manipulator Calibration[END_REF][5][6]. Generally, two main approaches that allows us to improve the identification accuracy without increasing the number of experiments exist. The first one deals with the preliminary optimization of the manipulator measurement configurations (so-called design of calibration experiments). The second approach consists in enhancing the objective function to be minimized inside the identification algorithm (in order to minimize impact of the measurement noise). As follows from the literature analysis, the first approach has been considered in a number of papers [7-10], while the second one received less attention of the researches. For this reason, taking into account particularities of the measurement system used in our experiments, this paper focuses on the improvement of the second method, which looks rather promising here.

In robot calibration, there exists a number of techniques that differ in the measurement equipment, the nature of the experimental data (position, orientation, distance, etc.) and in the optimization algorithm that produces the desired parameters [11][12][START_REF] Driels | Automated partial pose measurement system for manipulator calibration experiments[END_REF][START_REF] Bai | Experiment study of PUMA Robot Calibration Using a Laser Tracking System[END_REF]. At present, the most popular are the socalled open-loop methods that utilise external measurement devices to obtain either full or partial pose of the end-effector (i.e. position and/or orientation) [12][START_REF] Driels | Automated partial pose measurement system for manipulator calibration experiments[END_REF][START_REF] Bai | Experiment study of PUMA Robot Calibration Using a Laser Tracking System[END_REF]. However, it should be noted that the manipulator end-effector orientation cannot be measured directly, so the orientation angles are calculated using positions of several points around the end-effector centre point (TCP). Relevant example can be found in [START_REF] Dumas | A Methodology for Joint Stiffness Identification of Serial Robots[END_REF], where the end-effector orientation is evaluated using three target points located on the special measurement flange (i.e., three orientation angles are computed from nine Cartesian Project ANR COROUSSO, Project ANR ROBOTEX coordinates provided by the laser tracker). However this approach, which is based on the minimization of the squared sum of the position and orientation residuals, does not allow to minimize the measurement errors impact in the best way (in fact, the position and orientation components of the objective to be minimized are not weighted properly from the statistical point of view). To overcome this difficulty, in this paper it is proposed to use only direct measurement information, i.e. to replace the conventional objective function (squared sum of the position and orientation residuals) by a homogeneous squared sum of position residuals for all measurement points. It is clear that this approach is promising for the identification accuracy improvement, but it requires some revisions of the identification algorithm, which is proposed in this paper.

To address the above defined problem, the remainder of this paper is organized as follows. Section 2 defines the research issue and basic assumptions. In Sections 3, the identification algorithm is presented. Section 4 proposes comparison of the developed identification algorithm with the conventional one. Section 5 presents an application example illustrating benefits of the proposed approach. Finally, Section 6 summarizes the main contributions of the paper.

II. PROBLEM STATMENT

Let us consider a serial robot whose end-effector location (position and orientation φ ) is computed using the following vector function

(, )  tp φ p    (, , ) g  tq θΠ
where (.) g defines the manipulator extended geometric model, q is the vector of actuated coordinates, is the vector of robot elastostatic deflections, and the vector of the parameters 0 is presented as the sum of the nominal component 0 Π and geometrical errors θ   ΠΠ Π Π to be identified via calibration.

In addition to the geometric equation ( 1), let us consider the elastostatic model that allows us to compute the deflections caused by the external loading applied to the manipulator end-effector. Relevant equation can be presented in the following form [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF]

θ F      T θθ ••  θ kJ F
where the matrix θ is the manipulator Jacobian with respect to the elastostatic deflections , and is the manipulator compliance matrix to be identified.

(,, ) / g   θΠ Jq θ θ θ k
Assuming that the values  and are relatively small, equation (1) can be linearized and presented in the form

Π θ      0 θθθ •• • T     tg J JkJ Π F •
where the first term 00 i corresponds to the nominal geometric model (i.e. to the case when ), and the matrix  is the manipulator Jacobian with respect to the geometrical parameters .

(,, ) g  gq 0 Π g  ,  θ 0 Π 0 Π (,, ) /  θΠ Jq Π
The above presented equation ( 3) is the basic expression for the robot calibration that allows user to obtain the desired geometric and elastostatic parameters  and θ k . For this purposes a set of experiments are carried out, in which the end-effector locations

Π   i t
are measured by an external device for several manipulator configurations defined by the vectors of the actuated coordinates  . It is also assumed that the corresponding vectors of the external loading are different to ensure non-singularity of the relevant observation matrix used in the identification procedure. However, there are some difficulties here related to the estimation of the orientation components   i φ of the location vectors i (, ) ii  tp φ . There are two main approaches here that are considered below. Usual approach is based on the straightforward utilization of equation [START_REF] Mooring | Fundamentals of Manipulator Calibration[END_REF], where each configuration i q produces six scalar equations corresponding to the components of the sixdimensional location vector

 i q   i F are
 θ 2 0 θi θθ i , 1 •• • m • m T i ii i i        k Π tg J J k J F Π i n    
However, in practice, the orientation components yi ii x .) (. z  cannot be measured directly, so these angles are computed using excessive number of measurements for the same configuration , ich p oduce Cartesian coordinates

i q wh r   ,) | ,1 , ; ( n n pp 3 xij yij zij pj   
for several target points of the measurement tool attached to the manipulator mounting flange (Figure 1). Hence, instead of using scalar equations, that can be theoretically obtained from the measurement data, this conventional approach uses only 6 scalar equations for the identification. These may obviously lead to some loss of the parameter estimation accuracy. Another difficulty is related to the definition of the vector norm in equation ( 4), where the six-dimensional residuals are not homogeneous. It is clear that the position and orientation components must be normalized before computing the squared sum, but it is a non-trivial step affecting the accuracy (in practice, the normalization factors are usually defined intuitively).

3mn m

To overcome this difficulty it is proposed to reformulate the optimization problem (4) using only data directly available from the measurement system, i.e. the Cartesian coordinates of all reference points (Figure 1). This idea allows us to obtain homogeneous identification equations where each residual has the same unit (mm, for instance), and the optimization problem is rewritten as follows ( , ,

T ij xij yij zij p pp  p ) θ 2 θij θθ i () () ( j , ) 0 11 •• • • ij i mn p pp T ij ij ji           Π k pg J J k J F Π m i n m    
Here, the superscripts "(p)" indicate the position components (three Cartesian coordinates) of the corresponding location vectors, the index "i" defines the manipulator configuration number, while the index "j" denotes the reference point number. An obvious advantage of this formulation is simplicity in the vector norm definition (conventional Euclidian norm can be applied here reasonably, the normalization is not required) and elimination of the problem of the weighting coefficient selection. In fact, under assumption that the measurement noise is presented as a set of i.i.d. random values (similar for all directions x, y, z and for all measurement configurations), the optimal linear estimator should have equal weights. Besides, there are some potential benefits in the identification accuracy here, since the total number of the scalar equations increases from 6 to . m 3mn

To prove advantages of the second approach based on the objective function (5), the following sections will be devoted to the development of the dedicated identification algorithm and the comparison study of the conventional and proposed techniques.

III. IDENTIFICATION ALGORITHM

Let us assume that the measurement tool has reference points (

) that are used to estimate relevant vectors of the Cartesian coordinates for manipulator configurations. Using the homogeneous transformation technique, corresponding geometric model (1) can be presented as the matrix product

n 3 n  (, ,) jj j j ix i y i T zi ppp  p  ( •, , ) • base rob j j ot i i to l io  T θΠ TT q T    
where the vectors j i p are incorporated in the forth column of j i , the matrix defines the robot base location, the matrices

T base T , , j T θ 1 j  tl oo
describe locations of the reference points that are observed by the measurement system (see Figure 1), and the matrix function robot i i describes the manipulator geometry and depends on the current values of the actuated coordinates q , the robot elastostatic deflections , and the vector of the parameters to be estimated. the Cartesian coordinates of the reference points ,1 , j i jn  p corresponding to the configuration can be expressed in the following form i q (,,) (,,)

j j i to base bas o e robot i i b l ase robot i i   θΠ pp R p q RR q p θΠ 
This equation should be accompanied by the elastostatic model

θθ ii i
i that allows us to obtain 3 scalar equations for the identification purposes, where the following vectors/matrices are treated as unknowns: , , ,

T (, •  ,)

 θ kJ q Π F θ j tool p θ k Π • mn base p base R
To simplify computations, it is proposed to split the identification procedure into two steps. The first one deals with the estimation of base , base , tool , which are related to the base and tool transformations (assuming that the manipulator parameters are known). The second step focuses on the estimation of θ k and under assumption that the base and tool components are already identified. To achieve desired accuracy, these steps are repeated iteratively several times. that can also be rewritten in a matrix form as

 ji i i i j too base T robot robot robot base l                    p pp I p R r u     
where and are defined as follows

i robot p i robot R       (,,) ;
(,,) Applying to the linear system (12) the least-square technique, the desired vectors defining the base and tool transformations can be expressed as follows

robot robot i i robot robot i i ii  pp q q θθ R ΠΠ R
      

T

Step 2. On the second step, the remaining parameters and θ , which define the manipulator geometry and the elastostatic properties, are estimated. For this purpose, the principal system (6) is linearized and rewritten in the form

Π k  () () • i jj p i robot i i      pJ Π A p • j p χ      
where the subscript "(p)" denotes the positional components (three first rows) of the corresponding matrices, Π is the vector of geometrical errors, the vector χ collect the unknowns that corresponds to all non-zero elements (by definition) of the compliance matrix θ k (structure of the matrix θ defines in accordance with the assumptions of links\joins compliance before the identification), the matrix i k j  J q is the geometric Jacobian computed for the configuration i with respect to the reference point , and j () j p i  A is derived by relevant transformation of the last term of equation (3) into the vector form:

      11 ,..., jT T i i i i ni ni i jj j j    AJ J J F J F  
In the last expression, 1 i ,..., i j j n JJ denote the vector-columns obtained by splitting of the geometric Jacobian j i  J . For the computational convenience expression ( 17) can be presented in the matrix form

      () () ,• jj p j p iii        Π pJA χ
where t is the residual vector corresponding to the jth reference point for the ith manipulator configuration.

  p pp

Appling to this system the least-square technique, the desired vectors , , defining the manipulator geometric and elastostatic properties, can be expressed as

Π χ  1 () () () 11 11 T mn mn T j pj p j p j ii i i ij ij                 Π BB B p χ          where . () () () , jp jp jp ii i    BJ A
It should be noted that, to achieve the desired accuracy, the steps 1 and 2 should be repeated iteratively.

IV. COMPARISON ANALYSIS

To illustrate the efficiency of the proposed technique, let us compare its accuracy with the conventional one that operates with the full pose information. Their distinctions and particularities can be described as follows:

Approach #1 (conventional): The identification is based on the full pose information x lq l q q yl ql q q q q q zl l q q

q q l q                
where 123 are the actuator coordinates and 123 are the link lengths to be identified. The manipulator kinematics and procedure of the input data preparation for the approaches #1 and #2 are presented in Figure 2. . It was assumed that the geometric errors in the link lengths and the actuated joint offsets are respectively

1 3.0 m l m   , 2 2.0 mm l   , 3 and 1 5.0 l   mm 1.0 deg q   , 2 0.5d q eg   , 3 2.0 deg q   .
Besides, it was also assumed that the measurements errors are i.i.d. random variables with the standard deviation 0.01mm   (which perfectly corresponds to the precision of the equipment used for the experimental validations). The desired parameters were estimated using three measurement configurations, which were generated randomly. To obtain reliable statistics, the calibration experiments have been repeated 1000 times.

Simulation results are summarized in Table I, which presents the standard deviations for the estimates of the desired parameters (corresponding mean values of the estimates are equal to the assigned values with the high accuracy). As follows from these results, the proposed approach ensures the accuracy improvement in the estimation of the link length deviations i by the factor of 2.25 ... 3.83, while the accuracy improvement for the joint offsets estimations i is slightly less, up to 3.33. This confirms advantages of the proposed approach, but it should be mentioned that these numbers are obtained for particular set of the measurement configurations and particular normalization factor utilized in the approach #1. However, as follows from our study, approach #2 always provides essentially better results.

l  q  TABLE I. IDENTIFICATION ACCURACY FOR DIFFERENT IDENTIFICATION APPROACHES

Standard deviation Parameter

Approach #1 Approach #2

Improvement factor

V. APPLICATION EXAMPLE

The developed identification technique has been applied to the elastostatic calibration of industrial robot KR-270 (Figure 3). To take into account the influence of the gravity compensator (which creates the closed-loop) and to apply the virtual joint modeling approach developed for strictly serial robot [START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF], an equivalent non-linear virtual spring is used (its stiffness depends on the joint variable 2 ). However, to implement this idea, it is reasonable to consider a set of the to be identified and leads to the linear system of the identification equations similar to those considered in Section III.

To find the measurement configurations that ensure the best identification accuracy, the design of experiments technique has been applied, which is based on the dedicated industry-oriented performance measure proposed in our previous work [START_REF] Klimchik | Design of Calibration Experiments for Identification of Manipulator Elastostatic Parameters[END_REF]. This yielded 15 optimal measurement configurations with five different angles 2 that are distributed between the joint limits almost uniformly. These optimal configurations have been obtained taking into account physical constraints that are related to the joint limits and the possibility to apply the gravity force (work-cell obstacles and safety reasons). The results of the calibration experiment design are summarized in Table 2.

q At the measurement step, the manipulator was sequentially moved from one configuration to another, where the external loading 250 kg was applied to the special end-effector presented in Figure 4 (it allowed us to generate both external forces and torques). Corresponding experiment setup is shown in Figure 3. To measure the reference point Cartesian coordinates, the laser tracker system Leica AT901 was used. To evaluate manipulator elastostatic deflections, the reference point coordinates have been measured twice, before and after application of the external loading.

Joint angles, [deg]

q 1 q 2 q 3 q 4 q 5 q 6 79. Using these measurement data, the two-step identification procedure has been applied (see Section III). On the first step, the tool and base transformations has been computed; corresponding results are presented in Table 3. On the second step, they have been used for the identification of the elastostatic parameters, which are presented in Table 4

VI. CONCLUSIONS

The paper presents the advanced robot calibration technique that is based on the partial pose information, without explicit computation of the end-effector orientation. In contrast to previous works, the proposed technique uses the Cartesian coordinates measurements only, but several reference points of the end-effector are used. The proposed approach allows us to avoid the problem of non-homogeneity of the least-square objective, which arises in the conventional identification technique, where the full-pose information (both position and orientation) is considered. The technique does not require any normalization of the measurement data and can be efficiently applied to both geometric and elastostatic calibration. As follows from the simulation analysis presented in the paper, the developed technique essentially improves the identification accuracy compared to the conventional one. The obtained theoretical results have been successfully applied for the elastostatic parameters identification of an industrial robot used in aerospace industry.

Figure 1 .

 1 Figure 1. Typical measurement tool used for robot calibration

,

  Π k can we written as follows

Π

  Taking into account that any homogeneous transformation matrix a can be split into the rotational and translational components and presented as b

1 . 3 

 13 For the first step, taking into account that the errors in the base orientation are relatively small, the matrix is presented in the following form base identity matrix, vector base r includes the deviations in the base orientation, and the operator "This leads to the following presentation of equation (8)

  u n are treated as unknowns.

  manipulator end-effector (see Figure1).Approach #2 (proposed):The identification is based on the partial pose information, where three measurement points in the objective function to be minimized by the identification algorithm.Let us consider a 3 d.o.f. serial manipulator whose geometry is described by the following equations 

Figure 2 .

 2 Figure 2. Input data for two identification approaches

  compliance coefficients  2 corresponding to a number of different joint angles  2 that cover relevant joint limits. This yields the extended set of the elastostatic parameters

Figure 3 .

 3 Figure 3. Experimental setup for the identification of the elastostatic parameters

Figure 4 .

 4 Figure 4. Measurement tool with three markers on the end-effector used for elastostatic calibration experiments

  known. It is clear that that corresponding system of linear equations can be solved for  

	manipulator configurations   ,1 , i im  q	m	θ is high enough and the vectors , Π k if the number of 

TABLE II .

 II OPTIMAL MEASUREMENT CONFIGURATIONS

  Borm, C. Menq, "Determination of optimal measurement configurations for robot calibration based on observibility measure," International Journal of Robotics Research, vol. 10(1), 1991, pp. 51-63. Daney, "Optimal measurement configurations for Gough platform calibration," IEEE International Conference on Robotics and Automation (ICRA 2002), 2002, pp. 147-152. [10] A. Klimchik, Yier Wu, S. Caro, A. Pashkevich, "Design of Experiments for Calibration of Planar Anthropomorphic Manipulators," In the 2011 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM 2011), 2011, pp. 576-581 [11] M. Ikits, J.M. Hollerbach, "Kinematic calibration using a plane constraint," Proceedings of 1997 IEEE International Conference on Robotics and Automation (ICRA 1997), vol.4, 1997, pp. 3191 -3196 [12] H. Zhuang, O.Masory, Y. Jiahua, "Kinematic calibration of a Stewart platform using pose measurements obtained by a single theodolite," Proceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems 1995, vol. 2, pp. 329 -334.

	. Here  correspond to five different aggregated 25 compliances of the second actuated joint (which take into values 21 ,..., 
	account impact of gravity compensator and are obtained for different angles of 2 q ) and 36 ,...,   are the compliances of actuated joints 3...6 respectively. These parameters have been
	further used for the compliance error compensation for the
	robotic based milling of a aircraft parts, where essential
	improvement of the precision has been achieved.
	base p , mm	base r , mrad	1 tool p, m m	2 tool p, m m	3 tool p, m m
	34.4 31.9 97.8       	52.8 2.2 15.6     	279.2 16.4 91.9      	279.2 25.2 96.1     	28 130.5 1.8 5.6    	   	[5] J. Hollerbach, W. Khalil, M. Gautier, Springer Handbook of robotics, Springer, 2008, Chapter: Model identification, pp. 321-344. [6] Y. Sun, J.M. Hollerbach, "Active robot calibration algorithm," IEEE International Conference on Robotics and Automation (ICRA 2008),
							2008, pp. 1276-1281.
							[7] W. Khalil, M. Gautier, Ch. Enguehard, "Identifiable parameters and
	TABLE IV.	IDENTIFIED VALUES OF MANIPULATOR ELASTOSTATIC PARAMETERS	optimum configurations for robots calibration," Robotica, vol. 9, 1991, pp. 63-70.
		[8] J. χ i Values, [rad/N×mm]
		χ 21 χ 22 χ 6 χ 5 χ 4 χ 3 χ 25 χ 24 χ 23	0.287 ±0.0003 0.302 ±0.0005 0.293 ±0.0010 0.246 ±0.0007 0.416 ±0.0011 2.786 ±0.0071 3.483 ±0.0120 2.074 ±0.0267 0.277 ±0.0004			[9] D.

ACKNOWLEDGMENT

The work presented in this paper was partially funded by the ANR, France (Projects ANR-2010-SEGI-003-02-COROUSSO and ROBOTEX). The authors also thank research engineers Fabien Truchet and Joachim Marais for their great help with the experiments.