Paul Leplay

Olivier Lafforgue

François Hild

Analysis of asymmetrical creep of a ceramic at 1350°C by Digital Image Correlation

The glass industry requires the use of innovative ceramics that enable for long lifetimes. At very high temperatures, one of the key parameters for ceramics is their creep resistance. The characterization of the creep behavior, usually assessed through flexural tests, can be complex when an asymmetry appears between tension and compression. To detect and quantify such asymmetrical behaviors, the use of Digital Image Correlation (DIC) is proposed. First, several challenges are to be tackled for DIC at high temperature, namely, the random pattern stability, the radiation filtering and the heat haze. They are exacerbated by the limited possibilities to heat ceramics, the nonuniform strain fields and their low levels. Beyond several experimental developments, the strain uncertainties are decreased thanks to the use of two global approaches of DIC based on ad hoc finite-element kinematics. Last, the proposed methodology is applied to the analysis of asymmetric creep at 1350°C of an industrial zircon ceramic designed for its high creep resistance.

Introduction

In the glass industry, the ceramics used in furnaces or for forming are subjected to extreme conditions in terms of temperature, corrosion and creep. The development of innovative ceramics with optimal chemical and mechanical properties is a scientific and industrial challenge for the suppliers. Among the previously listed properties, creep is one of the most prevalent design parameters to ensure the lifetime of ceramic structures at high temperatures. 1, 2 , 3 For a reliable analysis of the final creep resistance, a rigorous approach consists of first performing mechanical creep tests at high temperatures, then identifying the parameters of a suitable constitutive law and finally running numerical simulations to predict the lifetime in various applications.

The first step is devoted to creep tests. For ceramic materials, 3-point or 4-point flexure is one of the most common mechanical tests thanks to its easy implementation. [START_REF] Quinn | Design Data for Engineering Ceramics: A Review of the Flexure Test[END_REF] However, its main drawback is the nonuniform strain/stress field, which implies to be able to distinguish tensile from compressive responses. Most of the time, the tensile and compressive behaviors are assumed to be symmetric so that the vertical deflection measurement is sufficient to calculate the flexural creep strain assuming perfect boundary conditions. But many ceramics, especially at high temperatures, may have an asymmetric creep behavior. [START_REF] Massard | Etude du fluage de réfractaires électrofondus du système alumine-zircone-silice[END_REF][START_REF] Dusserre | Applicability of image correlation techniques to characterise asymmetric refractory creep during bending tests[END_REF][START_REF] Leplay | Identification of asymmetric constitutive laws at high temperature based on Digital Image Correlation[END_REF] For such ceramics, uniaxial tensile tests may be performed, yet they are very delicate. [START_REF] Hild | Stress Heterogeneity Effect on the Strength of Silicon Nitride[END_REF] Compression and flexural tests have been considered, [START_REF] Blond | Modeling of high temperature asymmetric creep behavior of ceramics[END_REF] even though the identification of the tensile behavior is not straightforward. A final route is to only perform flexural tests, but with a more complex strain measurement technique. For instance, several strain gauges can be used, [START_REF] Laws | Derivation of the tensile stress-strain curve from bending data[END_REF] although the real boundary conditions of the test remain partially assessed. An alternative consists of resorting to full-field measurements not only to assess the actual boundary conditions, but also to measure the strain field of the sample. Since high temperatures are reached, non-contact techniques are desirable. Among them, [START_REF]Full-Field Measurements and Identification in Solid Mechanics[END_REF] Digital Image Correlation (DIC) has been chosen for its high potential for direct link with numerical simulations.

DIC has already been applied to analyze high temperature experiments. Lyons et al. [START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF] have shown its feasibility up to 650°C when samples were subjected to translations, thermal expansion and tensile loads. Then creep crack growth was studied in the vicinity of the crack tip for experiments on a nickel-based super-alloy up to 700°C [START_REF] Liu | Experimental Characterization of Crack Tip Deformation Fields in Alloy 718 at High Temperatures[END_REF] . Pan et al. [START_REF] Pan | High-temperature deformation field measurement by combining transient aerodynamic heating simulation system and reliability-guided digital image correlation[END_REF] showed that ordinary equipment can be used to evaluate coefficients of thermal expansion of various metallic materials. Leplay et al. [START_REF] Leplay | Identification of asymmetric constitutive laws at high temperature based on Digital Image Correlation[END_REF] studied the damage behavior of aluminum titanate by performing four-point flexural tests up to 900°C. The feasibility of reliable DIC measurements up to 1200°C has been shown by Grant et al. [START_REF] Grant | High-temperature strain field measurement using digital image correlation[END_REF] , Pan et al. [START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200°C[END_REF] and Hammer et al. [START_REF] Hammer | Strain Measurement at Temperatures Up to 800°C Utilizing Digital Image Correlation[END_REF] . Last, Novak and Zok [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF] achieved temperatures up to 1500°C by controlling the illumination and filtering, and minimizing heat haze effects. Moreover they proposed speckling procedures that allowed for reliable DIC analyses. It is worth noting that stereo DIC is seldom used in these high temperature analyses [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF] . De Strycker et al. [START_REF] Strycker | Measuring the thermal expansion coefficient of tubular steel specimens with digital image correlation techniques[END_REF] evaluated the coefficient of thermal expansion of stainless steels during a cooling from 600°C outside the furnace.

Among all the previous cited works, three main challenges for DIC at high temperature are systematically listed:

 The first issue is related to the surface random texture required for DIC. Up to 300°C, it was shown that regular paint was sufficient. [START_REF] Pommier | Analyse d'une éprouvette de fatigue thermomécanique par corrélation d'images[END_REF] So-called high temperature paints can be used up to 600°C. However, this solution becomes not always feasible at higher temperatures. [START_REF] Hammer | Strain Measurement at Temperatures Up to 800°C Utilizing Digital Image Correlation[END_REF][START_REF] Strycker | Measuring the thermal expansion coefficient of tubular steel specimens with digital image correlation techniques[END_REF] For instance, boron nitride combined with aluminum oxide was used by Lyons et al. [START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF] . Pan et al. [START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200°C[END_REF] propose to mix cobalt oxide with commercial inorganic adhesive to reach temperatures up to 1200°C.  As pointed out early on by Lyons et al. [START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF] , special care should be exercised at high temperature since the so-called black body radiations will lead to violations of brightness conservation. One remedy is to use more powerful lighting devices. However, there is a limit when temperatures of the order of 800°C are reached. Another solution is to use band pass filters to reach higher temperatures. [START_REF] Pan | High-temperature deformation field measurement by combining transient aerodynamic heating simulation system and reliability-guided digital image correlation[END_REF][START_REF] Grant | High-temperature strain field measurement using digital image correlation[END_REF][START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200°C[END_REF] It is also possible to work with ultra-violet light to avoid any radiation. 21  The last issue is related to the existence of a temperature gradient between the sample and the camera lens that induces a variation of the refractive index of the air. Even a constant temperature gradient generates an additional artifact for the optical path. However the main issue is related to varying thermal gradients in the air due to convection at high temperatures.

One solution is to use an air flux to homogenize the heat haze effects. [START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF][START_REF] Hammer | Strain Measurement at Temperatures Up to 800°C Utilizing Digital Image Correlation[END_REF][START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF] Other authors simply neglect these phenomena. [START_REF] Pan | High-temperature deformation field measurement by combining transient aerodynamic heating simulation system and reliability-guided digital image correlation[END_REF][START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200°C[END_REF] The heat haze influence is very dependent on the experimental setup and the way the sample is heated. The worst case is to have a traditional furnace since a large amount of air is then subjected to heat haze between the sample, the window and the lens of the camera. [START_REF] Leplay | Identification of asymmetric constitutive laws at high temperature based on Digital Image Correlation[END_REF][START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF][START_REF] Hammer | Strain Measurement at Temperatures Up to 800°C Utilizing Digital Image Correlation[END_REF] An alternative route is to resort to induction, 15 laser 18 or infrared [START_REF] Pan | High-temperature deformation field measurement by combining transient aerodynamic heating simulation system and reliability-guided digital image correlation[END_REF][START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200°C[END_REF] heating. These solutions limit the convection effects related to the air flow but cannot be applied to every class of materials, in particular to ceramics.

To have a better overview of the existing state of art and thus a good understanding of the objectives of the present paper, the main applications of DIC at high temperature are summarized and compared in Table 1. Several remarks can be made. First, most of the studies deal with metallic materials. It means that the heat haze can be minimized by the several types of possible heating devices and the final strain level is greater than 1-2 %. Second, the main goal is often to measure mean strains over a large region of interest so that the accuracy of the measured final value is improved. Third, the final strain uncertainties are seldom explicitly estimated. 1), it requires a stable paint pattern for ceramics at this temperature, a control of the important heat haze generated by the use of a traditional furnace, the measurement of nonuniform strain fields and very low uncertainties, down to at least 0.01 % given the low strain levels of ceramics. To reach these objectives, the paper is organized as follows. In Section 2, the two DIC methods used herein are presented. Section 3 deals with the experimental set-up, the sample preparation, the correction of the radiation and heat haze effects. A thorough uncertainty analysis is then performed. The results are given in Section 4 for an industrial zircon ceramic in terms of its asymmetrical creep behavior at 1350°C.

Digital Image Correlation methods 2.1 Local DIC

In the sequel, DIC will be used to measure strain fields at elevated temperatures. In all the above-cited works, except Ref. [START_REF] Leplay | Identification of asymmetric constitutive laws at high temperature based on Digital Image Correlation[END_REF] , local approaches to DIC were considered. In the present case, an alternative route will be followed since small strain levels are sought in conjunction with a good spatial resolution, namely, global approaches to DIC. [START_REF] Hild | Digital Image Correlation[END_REF]

Global DIC

Global approaches to DIC consist of performing the registration over the whole Region Of Interest (ROI). One advantage lies in the fact that the continuity of the displacement field is automatically satisfied. Let us consider the gray level picture f in the reference configuration, and g the picture in the deformed configuration. The conservation of gray level reads

)) ( ( ) ( x u x x   g f (1)
where x is any pixel within the ROI. To measure the displacement field u(x), a weak form of the previous equation is considered, namely, the sum of squared differences

       ROI c g f 2 )) ( ( ) ( x u x x (2)
is minimized with respect to the chosen kinematic parameterization of the experiment at hand. In the present study two different kinematic bases are considered. First (Section 2.2.1), a general purpose description corresponds to the choice of a finite-element based DIC approach. [START_REF] Besnard | Finite-element'' displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF] Second (Section 2.2.2), as already proposed when flexural experiments are performed, [START_REF] Hild | Displacement Measurement Technique for Beam Kinematics[END_REF][START_REF] Leplay | Damage law identification of a quasi brittle ceramic from a bending test using Digital Image Correlation[END_REF] a beam kinematics will also be implemented.

To validate the choice of any kinematic basis and thus the measured displacement field, the correlation residual  is defined as the gray level difference for any pixel of the region of interest in the reference configuration f(x) and that of the picture in deformed configuration corrected by the measured displacement g(x+u(x))

)) ( ( ) ( ) ( x u x x x    g f



(3) which should be as close as possible to the acquisition noise.

Finite-element kinematics: Q4-DIC

In the present experiment, it is expected that the displacement field is continuous. One of the simplest kinematic bases is associated with 4-noded quadrilaterals. In the element frame (1 ≤ x ≤ 1, 1 ≤ y ≤ 1), the shape functions

N i (i = 1,4) read ) 1 )( 1 ( 4
1 ) , ( y x y x N i    (4)
so that the two components of the elementary displacement field are expressed as

   4 , 1 ) , ( ) , ( i i i e x y x N y x u  and    4 , 1 ) , ( ) , ( i i i e y y x N y x u  (5)
where  i and  i are the nodal displacements, along two perpendicular directions of the considered element e.

Beam kinematics: Beam-DIC

Since a four-point flexural test is analyzed, a reduced kinematics is also devised. It is based on an Euler-Bernoulli hypothesis. [START_REF] Timoshenko | Strength of Materials[END_REF] Five displacement fields are considered. [START_REF] Leplay | Damage law identification of a quasi brittle ceramic from a bending test using Digital Image Correlation[END_REF] First, in the central part of the beam submitted to four-point flexure, only three kinematic degrees of freedom are needed, namely, the (constant) curvature , a rigid body rotation  and vertical translation . In the horizontal direction, two degrees of freedom  0 and  1 are needed to describe a rigid body translation and tension/compression. The latter may be induced by the friction on the various supports of the flexural setup. The displacement of any point of dimensionless coordinates (x,y) reads

                                        1 0 2 2 0 0 1 2 1 0 ) , ( ) , (      L L x x x x y xy y x u y x u b y b x (6)
where L denotes the beam length. The advantage of this last choice is that the number of unknowns has been drastically reduced so that the measurement uncertainties are lower. [START_REF] Hild | Displacement Measurement Technique for Beam Kinematics[END_REF] Furthermore, the evaluation of the rotation and curvature does not require any post-processing as they are directly related to the measured degrees of freedom and the differentiations of Equation (6).

Experimental setup 3.1 Description of the material and the test

For the experimental work, an isopressed zirconium silicate product (ZS-ISO) made by Saint-Gobain SEFPRO is considered. This material is manufactured by cold isopressing followed by pressureless sintering. The ZS-ISO samples have a density of about 4000 kg/m [START_REF] Wilkinson | Creep Mechanisms in Multiphase Ceramic Materials[END_REF] and an open porosity of about 5 %. In service, this ceramic undergoes extreme conditions, namely, temperature, corrosion and creep.

To evaluate the creep properties of such materials, a parallelepipedic sample (h × w × L = 15 × 10 × 150 mm 3 ) is tested in a 4-point flexural device (outer span, D 1 = 120 mm, inner span D 2 = 40 mm) with a dead load of 129 N. Assuming a linear elastic behavior, this load level generates a tensile/compressive stress of about ± 7 MPa. In front of the furnace equipped with a window, a digital single-lens reflex (DSLR) camera is chosen for its large sensor size, its high definition and low cost (Canon 700D, APS-C sensor, 18 RGB Mpixels). To obtain a small enough angle of view, a 300-mm macro telephoto lens is used with a 2 × focal extender (Figure 1). One image is shot every 15 minutes for 45 hours. The physical size of the gray level pixel is 16.2 µm.

Measuring the response in flexure of ZS-ISO is challenging since this ceramic is designed to have a high creep resistance even at these temperature levels. It is possible neither to speed up creep, nor to wait too long. The final creep strain of ZS-ISO remains very low (i.e., about 0.5-1.0 %), which is an issue for DIC measurements at such high temperatures. The strain uncertainties should be at least one to two orders of magnitude lower than the expected final strain levels. Therefore, specific experimental features must be considered for DIC at 1350°C.

Specific experimental features for DIC at 1350°C

Paint pattern stability and surface contrast

Before the test begins, a random pattern must be deposited onto the sample surface for DIC purposes. The first challenge consists of finding a paint that is stable up to 1350°C for several days on this ceramic. An in-house solution was developed with a non-oxide powder (silicon carbide). The stability of this specific paint was checked for up to 200 hours at 1350°C on ZS-ISO samples. This pattern has a mean correlation radius of about 3 pixels, which is favorable for DIC analyses. [START_REF] Besnard | Finite-element'' displacement fields analysis from digital images: Application to Portevin-Le Châtelier bands[END_REF][START_REF] Triconnet | Parameter choice for optimized digital image correlation[END_REF]

Blackbody radiations and image contrast

The painted sample is heated up to 1350°C. Above 800°C, the blackbody radiations become critical for the sensors of DSLR cameras, despite their inherent band pass filter to remove UV radiations (< 410 nm wavelengths) and IR radiations (> 710 nm wavelengths). At 1350°C, the intensity of the red channel is so high that the final image has no contrast at all, even with a 1000 W Xenon lamp used to lighten the sample. To be able to run DIC analyses at these temperature levels, the benefit of utilizing a blue filter has already been demonstrated. [START_REF] Grant | High-temperature strain field measurement using digital image correlation[END_REF][START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200°C[END_REF]21 Such a blue band pass filter prevents the blackbody radiations from being captured by the sensor. For a constant temperature test, the radiation level may not vanish provided it is constant. Consequently, a filter (i.e., 470 ± 20 nm) is preferred since the green channels also record some information (Figure 1). It is worth remembering that the green sensors have the advantage of being twice as numerous as blue and red sensors in the Bayer matrix in DSLR cameras. The blue and green channels are here both used, so that the final texture and dynamic range of gray level images is slightly enhanced. To still increase the final contrast, the image histogram is stretched directly when the 14-bit *.cr2 files are converted to 8-bit *.tif images.

Heat haze and image stability

The sample is heated inside a furnace with several electrical resistances. This type of heating enables for a homogenous thermal loading of ceramic materials (± 3°C variation in the present case inside the furnace). However, compared to other heating systems suited to metallic samples, 16,18,15 21 such a furnace generates a large amount of heat haze due to the varying refractive index of the hot air in front of the sample. This heat haze is problematic since DIC assumes that the differences from one image to another only come from a physical motion, and not from an optical artifact. A first solution consists of subtracting the heat haze influence to the DIC results, provided an accurate evaluation of all the dynamic motions of the air between the sample and the lens are quantified. [START_REF] Delmas | Shape distortions induced by convective effect on hot object in visible, near infrared and infrared bands[END_REF] A second route is to minimize the heat haze as much as possible. [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF] First, as proposed by Lyons et al. [START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF] , a double glazing made of sapphire was used to close the furnace window to enable for the sample visualization without additional thermal heterogeneity and heat haze. Then, a fan was added to regularize the flow between the window and the lens. Although such a fan was supposed to minimize the heat haze, [START_REF] Lyons | High-temperature deformation measurements using digital-image correlation[END_REF][START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF] no significant benefit was observed in the present case. It would also be possible to add a fan directly inside the furnace to reduce the optical distortions due to temperature gradients. [START_REF] Hammer | Strain Measurement at Temperatures Up to 800°C Utilizing Digital Image Correlation[END_REF] Once the heat haze is minimized as much as possible, it is possible to 'average' it. Instead of taking a lot of images during the test and averaging a posteriori the measured displacement fields, [START_REF] Strycker | Measuring the thermal expansion coefficient of tubular steel specimens with digital image correlation techniques[END_REF] it is proposed to average a priori the images and the heat haze influence as well. A neutral density opaque filter is used to strongly reduce the amount of light captured by the camera (Figure 1). Hence, the integration time of the camera is increased up to several seconds to take one single image with a correct exposure. The heat haze, which is assumed to be random, is continuously averaged over this period of time. Provided the exposure duration is long enough with respect to the heat haze frequency, the images are slightly blurred, but very stable. In that sense, to make an a priori average with 10-s image is equivalent to average a posteriori the results of one thousand 0.01s images. Fast creep rates in comparison with the time exposure may induce artifacts, but this is not the case in this work. The efficiency of such an approach with long exposure images is proven in the next part.

DIC strain uncertainties

To evaluate the optimal exposure duration, several configurations are tested at 1350°C. Seven series of 20 images are shot with an exposure duration ranging from 1/250 s to 30 s thanks to different neutral density filters. No load is applied so that no creep is expected, namely, no strains are expected. Therefore, the variations of the local strains for the 20 acquired images are considered as representative of all the artifacts and uncertainties associated with the experimental setup (e.g., DIC algorithm, sensor noise, heat haze). In the following analyses, the standard uncertainties are evaluated (i.e., they correspond to the standard deviation of the reported quantity). The uncertainty analysis is carried out in two steps:

 First, a Gaussian white noise is added to each of the seven reference images. The added noise has a 3 gray level standard deviation, corresponding to the quality of the camera sensor at ISO100. The reference and noisy images are then registered with the two different DIC codes (Figure 2). The present standard strain uncertainties are about one order of magnitude higher than the usual uncertainties due to the subpixel interpolation. Again, the strain uncertainties only depend on the element size, but not on the exposure duration. Beam-DIC is much less sensitive to noise thanks to its low number of degrees of freedom in conjunction with the large size of its unique element. Last, both results are consistent since the extrapolated results of Q4-DIC tend to those of Beam-DIC. This a priori analysis shows that to reach very low standard strain uncertainties, Beam-DIC is desirable, provided the beam kinematics is consistent with the analyzed experiment.

Figure 2 -A priori standard strain uncertainties due to numerical noise representative of the used camera

 Second, all the images of the seven series are correlated with respect to their own reference image (Figure 3). For the images shot at 1/250 s without neutral density filters, the standard strain uncertainties are very high ( 0.4 % for Q4-DIC with 64 × 64 pixel elements). The reason for such high levels is understood when analyzing the measured displacement field, which reveals the effect of heat haze (Figure 4 left). Very coarse meshes help, but even with one large element and using Beam-DIC, the standard strain uncertainty is still of the order of 0.04 %. The standard strain uncertainties are reduced by one new order of magnitude when the exposure duration is greater than a few seconds. For 30-s images, the standard strain uncertainty is of the order of 0.05 % for Q4-DIC with 64 × 64 pixels, and 0.002 % with Beam-DIC with one element. This result means that the exposure duration is long enough with respect to the frequency of the heat haze. For Q4-DIC, a plateau is reached with the longest exposure durations. Others sources of uncertainties such as the previously studied numerical noise limit the accuracy of a local measurement (Figure 4 The fact that the experimentally determined strain uncertainties (Figure 5) do not reach the a priori estimates (Figure 4) is an indication that the experimental configuration is very challenging for reliable DIC analyses. However, by increasing the time integration, the achieved uncertainties reach very low levels. The latter ones have been obtained by adapting the time integration and the spatial regularization provided by Q4-and Beam-DIC.

Asymmetrical creep behaviour of ZS-ISO at 1350°C 4.1 DIC general results with Q4-DIC and Beam-DIC

During the 45-hour creep test, 180 images with a 30-s exposure duration are shot. The reference image is acquired just after the weight is applied (t = 0 h). On the reference image, only the central zone between the two upper supports where the curvature is assumed to be constant is analyzed (Figure 5). This central zone is either divided into 504 elements of size 64 × 64 pixels for Q4-DIC, or used with one beam element of size 936 × 2345 pixels for Beam-DIC. For both approaches the measured displacement fields are corrected thanks to a small calibration sample put just below the tested sample. The former enables all rigid-body motions due to the furnace, the camera, and the tripod to be excluded from the displacement field measurements. First, the consistency of the measurement results is analyzed in terms of the final deformed shape (Figure 6). For Q4-DIC, which is obtained without any assumption on the flexural kinematics, the global shape of a bent sample is easily recognizable, no noise is visible and the horizontal displacement amplitude is consistent with the results of Beam-DIC (i.e., 9.3 to 11.2 pixels or 150 to 181 µm). The consistency of the results is also checked with the correlation residuals. They are initially quite low (about 2% of the dynamic range) but end at a rather high level (6-9% of the dynamic range) after 45h. With respect to standard room temperature tests, these levels should be considered as problematic. However it is believed that this gradual increase reveals the slow degradation of the non-oxide powder pattern.

The vertical displacements are analyzed via the deflection profiles (Figure 7). The latter is calculated as the mean vertical displacement for each vertical section of the beam (i.e., all vertical nodes for Q4-DIC, all vertical pixels for Beam-DIC). As expected from Figure 8, the final deflection profile corresponds to a 2 nd order polynomial. The final maximum vertical deflection is equal to 1493 µm (or 92.2 pixels) near the center of the sample. This deflection level corresponds to what would have been recorded by a displacement sensor as classically used in flexural tests without DIC. On these deflection profiles, a slight unbalance between the left and right sides is noticed. It means that this 4-point flexural test is not perfectly balanced, namely, the total weight was probably not correctly put in the middle of the two upper supports. Moreover, the two DIC results are fully consistent.

Figure 7 -Final vertical deflection profiles for Q4-DIC and Beam-DIC

From the strain fields shown in next part, it is possible to calculate that the strain levels are on average 10% higher on the left part of the sample. It confirms that this four-point flexural test was not fully balanced between the left and right. For a better analysis of the test imperfection, it is worth noting that it would be possible to use an identification approach based on DIC results taking into account these effects. [START_REF] Bouterf | Damage law identification from full-field displacement measurement: Application to four-point bending test for plasterboard[END_REF]

Quantification of the asymmetrical creep behavior

The longitudinal strain fields better reveal the differences between the two DIC approaches (Figure 8). In Beam-DIC (right), the horizontal field varies linearly with the vertical coordinate since it is directly computed from analytical functions (Equation ( 6)). The minimum and maximum final strains are directly available, namely, ε comp =  0.619 % and ε tens = + 0.877 % after 45h. It corresponds to a ratio of about 1.42 between the tensile and compressive strains. For Q4-DIC, the strain amplitude ranges from  0.76 % at the bottom and + 1.10 % at the top of the beam. A tension/compression asymmetry is also detected, but with a lower accuracy. The strains obtained with Q4-DIC have an uncertainty of the order of 0.09 %, which is consistent with the uncertainty analysis performed previously (Section 3.3). With Q4-DIC, to obtain reliable strain values in tension and compression, a posteriori assumptions have to be added (Figure 9). Thanks to the constant curvature in the central zone of four-point flexure, the local strains are horizontally averaged. These mean values are then linearly interpolated and extrapolated to reach the sample edges. This post-processing yields a compressive strain of ε comp  0.620 % and a tensile strain of ε tens = + 0.883 %, which is very close to the results obtained directly via Beam-DIC. The change of the tensile and compressive strains is plotted in Figure 10. These strains are compared to an analytical value calculated from the vertical deflection measurement (i.e., a 'virtual' displacement sensor). The maximum deflection is either corrected by the motions of the reference sample (Figure 7), or not. Several remarks are drawn:  Since DIC is started after the initial load has been applied, the measured strain field corresponds to the current total strain minus the initial elastic strain. [START_REF] Dusserre | Applicability of image correlation techniques to characterise asymmetric refractory creep during bending tests[END_REF]  The strain values obtained by Q4-DIC and Beam-DIC are equal within ± 0.01 % margin. In the present case, adding the flexural assumption a posteriori or a priori does not make a big difference.

 No significant primary creep is observed for this ZS-ISO ceramic, namely, the creep rate is constant only a few hours after the test begins.

 The mean linear creep rate is 2 × 10 -4 / h in tension, and 1.35 × 10 -4 /h in compression; it corresponds to a tension/compression ratio about 1.5. The first 10 hours are excluded from this analysis in case primary creep occurred.

 When the raw deflection is not corrected with respect to the reference sample, the resulting strain lies between the tensile and compressive DIC strains, but its value is noisy.

 When the raw deflection is corrected by removing the spurious motions, the strain value is much more regular and is now equal to the average between the tensile and compressive strains (i.e., the classical link between deflection and strains in flexure with a symmetric response).

Figure 10 -History of the tensile and compressive creep strains

It is worth noting that the uncertainties for the tensile and compressive strains are extremely low. The distribution of the difference between the strain profiles and linear creep rates is Gaussian and its standard variation is as low as 0.007%. Adding a priori or a posteriori hypotheses on the flexural kinematics enables the strain uncertainty level to be reduced by one order of magnitude with respect to the local values of Q4-DIC (Figure 8). To the authors' knowledge, this uncertainty level has not been reached at such high temperature, in particular above 1000°C.

Last, with these extreme strain values in flexure, ε tens and ε comp , it is possible to calculate the location of the neutral axis. With respect to the sample mid-height, the distance h n is defined as the axis where the creep strain vanishes (7) The neutral axis location is erratic at the beginning of the test when the strain levels are very low (Figure 11). As the strain levels increase, the neutral axis location is more and more accurately determined, and its final location is estimated to 1.3 ± 0.1 mm for the two DIC approaches. Since the approach with the deflection measurement assumes intrinsically a symmetric behavior, the neutral axis does not move from the mid-height.

Figure 11 -Location of the neutral axis, h n , during the creep test

The two approaches developed herein, namely, Q4-DIC and Beam-DIC, have yielded consistent quantifications of the asymmetrical creep behavior of the studied ZS-ISO ceramic at 1350°C in terms of vertical deflection measurements, minimum and maximum flexural strains, creep rates and neutral axis location.

Conclusion

To develop some innovative ceramics for glass industry, the issue of reliable creep characterization must be tackled, especially in case of asymmetry between tension and compression. DIC has been used herein to monitor a creep flexural test at 1350°C on an industrial zircon ceramic.

To be able to measure very low strains by DIC at such high temperatures, several experimental challenges have been addressed, especially the effects of heat haze that have been continuously averaged thanks to neutral density filters. Compared with the initial configuration, one order of magnitude on the standard strain uncertainties has been gained thanks to long exposure images. To decrease the strain uncertainties by a second order of magnitude, two global DIC approaches based on the known specific kinematics in flexure are used to measure the tensile and compressive strains. Standard strain uncertainties of the order of 0.007 % have been achieved during the entire creep test. Under these conditions, it has been shown that the creep strain rate is about 1.5 times higher in tension than in compression. Such an asymmetry would have been impossible to detect with the classical measurement devices of the vertical deflection. Some test imperfections, such as left/right unbalance in four-point flexure, have been observed. They should be recognized and quantified to avoid any wrong interpretation of the test. This is possible thanks to the measurement of displacement fields.

Thanks to this innovative methodology to measure accurately the full creep strain field at such high temperature, the next steps, namely, identification of constitutive parameters and numerical simulations, become possible to predict the lifetime of a structure made of this ceramic in its various industrial applications.

Figure 1 -

 1 Figure 1 -Overview of the experimental set-up for creep tests and DIC at 1350°C

  right).

Figure 3 -Figure 4 -

 34 Figure 3 -A posteriori standard strain uncertainties mainly due to the heat haze and acquisition noise

Figure 5 -

 5 Figure 5 -Reference image at 1350°C, region of interest and mesh for Q4-DIC (left) and Beam-DIC (right).

Figure 6 -

 6 Figure 6 -Final deformed mesh (amplified by 10) for Q4-DIC (left) and Beam-DIC (right). The longitudinal displacement field is also shown and expressed in pixels

Figure 8 -

 8 Figure 8 -Horizontal strain fields (expressed in %) as obtained with Q4-DIC (left) and Beam-DIC (right)

Figure 9 -

 9 Figure 9 -Post-processing Q4-DIC results to estimate the horizontal strain profile

Year / Reference Temperature (°C) Objective Material Heating Max. strain Std. strain resolution Pattern

	1996 12	[25-650]	Mean CTE/tensile strains	Superalloys	Furnace	0.8%	0.02% at 650°C	Boron nitride / Aluminum oxide
	2009 15	[25-1000]	Mean CTE/tensile strain	Superalloy	Ohmic	2%	NA	abraded surface
	2011 16	[25-1200]	Mean CTE strain	Stainless steel	Radiation	2%	NA	Cobalt oxide / inorganic adhesive
	2011 18	[25-1500]	Mean CTE and strain field	Superalloy / C/SiC composite	Laser	2%	0.04% at 800°C	Alumina/Zirconia
	2012 7	[600-900]	Strains in flexure	Ceramic	Furnace	0.5%	0.05% at 900°C	NA
	2013 22	[25-1550]	Mean CTE and strain field	Stainless steel C/SiC composite	Radiation	4%	0.5% at 1300°C	Pre-oxidized / ceramic coating
	2014 21	[100-1150]	Mean strain in tension / CTE	Superalloy	Induction	2%	0.5% at 725°C	Pre-oxidized / VHT paint
	2014 17	600	Mean strain in tension	Ti-6Al-4V	Furnace	25%	NA	Rustoleum paint
	2014 23	2600	Mean strain in tension	Carbon	Ohmic	4%	NA	Sprayed tungsten
	2014 24	720	Necking in tension	Steel	Induction	30%	NA	Electrochemical

Table 1 -

 1 Overview of the existing works on DIC at high temperature

	The goal of this paper is to address asymmetrical creep in a flexural test at 1350°C. Compared with
	existing works (Table

Corresponding email address: paul.leplay@saint-gobain.com

Acknowledgements

The authors would like to thank Saint-Gobain SEFPRO for providing the ZS-ISO samples and supporting this work carried out at CREE (Centre de Recherches et d'Etudes Européen) thanks to the experimental facilities of the Laboratoire d'Application Physique.


José Antonio

Reglero Ruiz

Michel Vincent

Noëlle Billon

Plane-compression properties of microcellular

Keywords: microcellular injection molding, chemical blowing agents, polypropylene, mechanical properties

come

-INTRODUCTION

The commercial and industrial interest in polymeric structural foams has been increasing in the last years, due to the combination of good mechanical properties together with weight reduction. Especially in injection molded parts, the foaming process provides several advantages, such as good surface aspect, good dimensional stability and the reduction of the raw material costs [START_REF] Xu | Microcellular injection molding[END_REF]. The injection foaming process leads to core-shell foamed structures with closed outer skins and a foamed core [START_REF] Mechraoui | Mechanical properties of polypropylene structural foams with fiber-reinforced skins[END_REF][START_REF] Ruiz | Micro and nano cellular amorphous polymers (PMMA,PS) in supercritical CO 2 assisted by nanostructured CO 2 -philic block copolymersone step foaming process[END_REF]. This allows obtaining equivalent mechanical properties to a neat polymer injected part, but with a weight reduction.

One of the main injection processes is based on using chemical blowing agents (CBA) [START_REF] Behravesh | Extrusion of polypropylene foams with hydrocerol and isopentane[END_REF][START_REF] Bociaga | The influence of injection molding parameters and blowing agent addition on selected properties, surface state and structure of HDPE parts[END_REF][START_REF] Bledzki | Effects of the chemical foaming agents, injection parameters and melt-flow index on the microstructure and mechanical properties of microcellular injection-molded wood-fiber/polypropylene composites[END_REF]. The foaming agent is added in the form of pellets to the feeder of the injection machine, in a proportion between 1 and 2 % in weight. The foaming agent usually is a low density polyethylene charged with a percentage of inorganic compounds, such as citric acid or sodium bicarbonate, which decompose at temperatures close to the melting point of the polymer, generating the gas, mainly CO 2 and water vapor.

The core-back expansion process which is investigated in this paper derives from the classical injection molding process, where a command module is added to control the movement of the mobile part of the mold. Once the polymer containing a foaming gas fills up the mold, the mold is accurately opened to quickly increase the cavity volume and rapidly reduce pressure. The sudden pressure drop enhances bubble nucleation and achieves a fine cell structure within the polymer foam. Thus, core-back foam injection molding is effective in achieving a high expansion ratio with uniform fine cell structure. Moreover, compared to classical microcellular injection processes, the core-back foam injection molding process allows increasing the thickness of the final part, thus reducing the foam density.

Although microcellular injection molding process has many advantages, the resulting products usually have severe surface defects such as swirl marks and lack of smoothness [START_REF] Cha | Microcellular injection molding with improved surface quality[END_REF][START_REF] Guo | Ocurrence of surface defects in TPO. Injected foam parts[END_REF]. The introduction of a gas counter pressure (GCP) in microcellular injection molding process has been recently analyzed to improve the surface quality of molded parts, because it can prevent the foaming reaction in the front flow during injection step [START_REF] Lee | A novel method for improving the surface quality of microcellular injection molded parts[END_REF][START_REF] Turng | Development of a hybrid solid-microcellular co-injection molding process[END_REF]. In all these studies, the efficiency of the introduction of the GCP has been demonstrated to reduce the surface defects drastically when microcellular foaming process is employed.

The present work is focused on the analysis of the mechanical properties of the structural foams obtained using this process. Plane compression tests have been carried out, to obtain the Young's modulus of the foamed core. After, a theoretical approach has been used to estimate the Young's modulus of the whole sample, including both the core and the solid outer skin. The influence of several parameters, such as cell size and foam density has been also investigated. The foamed core density has revealed as the most important parameter to determine the mechanical properties, with a small influence of the cell size. The results obtained show an interesting way to estimate the Young's modulus in compression tests of the foamed core, without removing the solid skin. Moreover, simple models have been employed to predict the mechanical properties of these materials.

-MATERIALS AND EXPERIMENTAL PROCEDURES

-Raw Materials

A pure polypropylene (PP) with a melt flow index of 65 g/10 min (ISO R1133), and a density of 0.91 g/cm 3 was employed in all the experiments. Two different endothermic chemical blowing agents referred as CBA-1 and CBA-2 have been used [START_REF] Ruiz | Polymer foaming with chemical blowing agents: experiment and modeling[END_REF]. These foaming agents are PE-based compounds with reactive elements (citric acid, sodium bicarbonate or a mix of both components). CBA-1 contains 70 wt. % of citric acid and CBA-2 contains 70 wt. % of sodium bicarbonate, (wt. % respect to PE matrix). In the following, CBA refers to the compound and not to the reactive elements only.

-Injection tests

Injection tests were carried out using an ENGEL machine, with a maximal injection capacity of 150 cm 3 . It is equipped with a shut-off nozzle. The screw diameter D was 32 mm, with L/D = 24 (L is the screw length). Polymer pellets were introduced in the barrel with 2 wt.% of CBA, and melted at temperatures that can reach 230°C. Mold temperature was fixed to 40 °C. A back pressure is applied at the plasticization unit which values between 6 and 8 MPa. The residence time of the CBA-PP system before injection shot is about 6 minutes. The decomposition reaction of the CBA produces both CO 2 and H 2 0 vapor, which is expected that after foaming and cooling of the samples remains as condensed water vapor. Then the polymer and the dissolved gas (CO 2 + H 2 0 vapor) are injected into the mold. A gas counter pressure (GCP) is applied in the cavity during the filling step to prevent the development of the foaming reaction. In a second step, the mobile part of the mold moves backwards slightly. This movement causes a rapid pressure decrease in the cavity so that it allows the blowing agent to expand and to create a cellular structure in the core of the injected part. No other injection parameters were varied (injection temperature, flow rate or core-back opening course).

Rectangular plates were molded. The final thickness of the sample is 2.5 mm after foaming. The total volume of the sample including sprue and gate was about 36 cm 3 , with a weight of 25 g. The mold temperature was fixed at 40°C. The injection temperature was 230°C. Gas Counter Pressure was fixed to 15 bar, using nitrogen gas, connecting the mold to a highpressure N 2 system using and automatic valve. Two groups of samples were injected using different CBA's. A total number of 15 plates was injected in each group. A minimum of three plates of each test were analyzed to assure a good reproducibility of the morphology results. For the mechanical tests, five plates of each group were analyzed. Figure 1 presents a scheme of the core-back expansion molding process.

-Morphological characterization

Scanning Electron Microscopy images (SEM) were taking using a Philips equipment (Model Philips XL-30). For the morphological observations, samples were fractured at ambient temperature, to assure a good contrast between solid and foamed regions. All the observations were carried out extracting three samples from the center of the plate, in the perpendicular direction to the injection flow.

The morphological parameters were quantified using the ImageJ® software [START_REF] Ruiz | Polymer foaming with chemical blowing agents: experiment and modeling[END_REF]. This software accounts for the number of bubbles in each image and the average radius. Cell density N c in the sample was calculated using equation (1) [START_REF] Gosselin | Cell morphology analysis of high density polymer foams[END_REF], which accounts for the 3D extrapolation starting from a 2D image:

( ) (1)
F  represents the foam core density, and P  the solid polymer density (0.91 g/cm 3 ). V f represents the volume of each extracted sample, in our case 2.5x10x10 mm 3 . The average radius is calculated from equation ( 2):

̅ ∑ ∑ (2)
, where N represents the bubble count. The minimum quantifiable bubble radius is 1 m, and a three different morphological observations were determined using this software to minimize the deviation derived from the extrapolation from 2D to 3D. The software calculates for each bubble i the radius R i , assuming that each bubble is a sphere, which is a first approximation.

Then, it calculates the average bubble radius ̅ of the foamed core in the micrograph by using equation [START_REF] Mechraoui | Mechanical properties of polypropylene structural foams with fiber-reinforced skins[END_REF]. Other geometrical approximations can be taken into account (for example, ellipsoids), but considering spherical-shape bubbles was found to be the most accurate calculation.

-Plane-Compression tests

Low velocity plane compression tests were carried out in five samples extracted of each plate, with dimensions 2.5x10x10 mm 3 , at a strain rate of 0.4 s -1 and at room temperature.

The external edges of the samples were cut to avoid side effects due to injection skin. Samples were tested using an INSTRON type machine, with a 10 KN load cell, up to 75% of stress, reaching the densification area. Young's modulus was calculated by fitting the initial linear part of the strain-stress curve of each sample. All the results of the five samples were averaged.

-RESULTS AND DISCUSSION

-Morphology

Figure 2 illustrates the morphology of the structural foamed samples. First, it is clear that using CBA-1, based on citric acid, leads to lower cell sizes than using CBA-2, which is based on sodium bicarbonate. It is observed also that both samples present an homogeneous structure, with a solid outer skin with total thickness about 700 m. The density values of the foamed core of the samples varied between 0.49 g/cm 3 and 0.71 g/cm 3 . In our work, samples obtained with CBA-2 presented lower density values. It was observed that the solid outer skin total thickness did not vary because the mold temperature was fixed in all the injection tests. Moreover, the foamed core density presents a very small variation because is directly related to the course during core-back opening, which was also fixed in all the experiments. Ii is expected that varying the core-back displacement has a direct influence in the microcellular density. The morphological parameters obtained using the ImageJ® software are presented in table 1.

CBA

Plane-compression tests

Figure 3 presents an example of a strain-stress (,) curve obtained from the plane compression test. This curve corresponds to a typical stress-strain curve of a cellular material [START_REF] Gibson | Cellular solids: structure and properties, Cambridge solid state science series[END_REF]. There are 4 regions clearly defined. Region 1 corresponds to the elastic regime of the material, and it is employed to obtain the elastic Young's modulus. The region number 2 is defined by the yield stress  y , and corresponds to the Plateau regime, and it is related to the energy absorption capability of the material. In an ideal cellular material,  y = p . In the region number 3, an important quantity of the gas inside the bubble escapes due to buckling and rupture of the cell walls, leading to a densification regime in which strain increases greatly. Finally, in region 4 there is no gas remaining in the cells, and the polymer presents a compression behavior similar to a solid. In figure 3,  D and  D represent the stress and strain at which densification process begins. In an ideal behaviour, the slope of the stress-strain curve in the plateau region is null, and no decrease of the slope is produced at the beginning of the densification region. However, in real samples, this small change in the slope can be produced due to inhomogeneities of the cellular structure Values of Young's modulus E f vary between 204 MPa and 222 MPa, about 10% of the value obtained for the solid PP (E s =2.2 GPa). This value corresponds to a high-impact resistance polypropylene, and was obtained experimentally and provided by the raw material supplier.

The elastic limit strain  y employed to calculate the values of Young's modulus was about 5% in all the cases. Some mechanical tests were employed to validate this limit, and after testing only up to the linear regime, samples recovered their initial thickness. On the other hand, densification strain values  y are in the range of 40%. This is a low value compared to low density polyolefin foams [START_REF] Rodriguez-Perez | Mechanical properties of closed cell polyolefin foams[END_REF] which present a higher value of the densification strain (75%), making them suitable for energy absorbing applications. It was observed that during the plane compression tests, the foamed core was compressed while the solid outer skin can be considered as an additional compression plate, and has no influence in the mechanical properties measured. Thus, the Young's modulus measured during compression tests corresponded only to the foamed core.

A semi empirical expression relating the Young's modulus of cellular materials with the Young's modulus of the solid polymer is given by Gibson and Ashby [START_REF] Gibson | Cellular solids: structure and properties, Cambridge solid state science series[END_REF]:

( ) (3)
E f is the Young's modulus of the cellular material, E S is the Young's modulus of the solid material,  f is the foam density and  s is the solid polymer density.

The values of E S (2.2 GPa) and  s (0.91 g/cm 3 ) are known. It is possible to plot the relative Young's modulus versus the relative foam density to determine the parameters C and n. The results are presented in figure 4, for all the experimental samples obtained, in a log-log graph, together with the theoretical correlation to determine C and n.

Figure 4 -Relation between relative young's modulus and relative density (Equation 3)

The obtained values for C and n are 0.1754 and 0.4493, res pectively, with a correlation factor r 2 of 0.98343. Although the correlation can be considered correct, the values of n are lower than expected (literature values are between 1 and 2), probably because the range of density analyzed in this work is narrow [START_REF] Gibson | Cellular solids: structure and properties, Cambridge solid state science series[END_REF]. However, as explained before, the range of density values was directly related to core-back opening course, and this parameter could not be varied during all the tests. The Young's modulus of the whole structural material can be obtained using the next expression [START_REF] Fleck | End compression of sandwich columns[END_REF], considering a mixing law in which the core and solid skin are placed in a series arrangement:

( ) [ ] (4)
, where E e corresponds to the Young's modulus of the structural material (foamed core + solid skin), and E S and E f are the Young's modulus of the solid skin and the foamed core. In addition, t s is the half thickness of the solid skin, and t f is the total thickness of the foamed core. Thicknesses of the solid outer skin and foamed core are obtained directly from the SEM images presented in figure 2.

Table 3 shows the average results obtained in both groups of samples, when applying equation [START_REF] Behravesh | Extrusion of polypropylene foams with hydrocerol and isopentane[END_REF]. As seen in table 3, the results obtained for the Young's modulus of the structural foam are in the range of 1.3 GPa, with no variations between both groups. This value is about 35 % lower than the Young's modulus of the solid material (E s = 2.2 GPa).

Group CBA t s (µm) t f (µm) E f (MPa) E e (MPa)
To summarize, the mechanical properties of the structural foams depend strongly on the density, and as expected, Young's modulus is increased with foamed core density. On the other hand, cell size does not have a great influence in the mechanical parameters measured. In our work, values of the Young's modulus are increased only a 10 % when cell radius is reduced about 100 % (see values of tables 1 and 2). However, it is known that reducing drastically the cell size up to sub-micrometrical or nanometrical sizes has a great influence in the mechanical properties, especially in the Young's modulus [START_REF] Ruiz | Foaming behaviour and compressive properties of microcellular nanostructured polystyrene[END_REF].

The final density of the samples does not vary significantly, indicating that the quantity of gas dissolved during the plasticization process remains constant during all tests. The thickness of the structural foamed samples does not present great variations, because is mainly controlled by the injection parameters that have been unchanged in this work. Only the average cell size is mainly controlled by the type CBA employed. When employing citric acid CBA cell sizes are drastically reduced, but no great impact in the mechanical properties is observed.

-CONCLUSIONS

In this work, the mechanical and morphological aspects of microcellular injected PP plates using Gas Counter Pressure (GCP) and core-back have been analyzed. Two different Chemical Blowing Agents (CBA) were employed, with a fixed proportion of 2% wt. and analyzing the morphological parameters and compression Young's modulus of the obtained samples.

Structural foamed samples obtained presented a homogeneous cellular structure, with similar values of density (about 0.5 g/cm 3 ) and thickness of the solid outer skin, indicating that both aspects are dependent on the injection conditions. However, the average cell size depended strongly on the type of CBA employed, and using a citric acid based CBA leads to a significant reduction in the cell size compared to the cellular structure obtained when using sodium bicarbonate-based CBA.

Concerning the mechanical properties, a set of low velocity plane compression tests were investigated. The values of the Young's modulus of the foamed core were between 200 MPa, about ten times lower than solid polymer, with a good reproducibility of the results. A theoretical approach was also used, which allows obtaining the Young's modulus of the whole structural foamed samples, which values in the range of 1 GPa, about five times higher than the values found for the foamed core.

Figure 1 -

 1 Figure 1 -Core-back expansion molding process.

Figure 2 -

 2 Figure 2 -SEM micrographs showing the morphology of the structural foamed samples (bar scale 1 mm) a) PP + CBA-1 b) PP + CBA-2

Figure 3 -

 3 Figure 3 -Example of a stress-strain curve obtained in a plane compression test.The results for the mechanical testing for both group samples are presented in table 2, in which the average results of the five samples are presented together with the dispersion values.

̅ (µm) Bubble count Solid skin thickness (µm) N c (cm -3 )  f (g/cm 3 )

	1	41  6	314	700  20	1.4410 5	0.56-0.71
	2	74  11	128	715  10	3.2110 4	0.49-0.53

Table 1 -

 1 Morphological parameters of the injected PP samples.

Table 2 -

 2 Mechanical properties of the microcellular injected PP samples.

Table 3 -

 3 Values of the Young's modulus of the PP structural foam samples obtained from equation 4.

	1	1	350  20	1800  40	345  11	1343  48
	2	2	367  10	1734  20	284  8	1321  39

ACKNOWLEDGMENTS

Authors would like to acknowledge SUMIKA Polymer Compounds France for providing the raw materials, and also to Aurore Claverie and Sebastien Huck (MECAPLAST France), for producing the injected plates. Financial support FUI (PLUME Project, France) is also gratefully acknowledged.


L Bossuet
email: lilian.bossuet@univ-st-etienne.fr

V Fischer

P Bayon

Contactless Transmission of Intellectual Property Data to Protect FPGA Designs

Keywords: IP protection, side channel, electromagnetic emanation analysis

Over the past 10 years, the designers of intellectual properties (IP) have faced increasing threats including illegal copy or cloning, counterfeiting, reverse-engineering. This is now a critical issue for the microelectronics industry, mainly for fabless designers and FPGA designers. The design of a secure, efficient, lightweight protection scheme for design data is a serious challenge for the hardware security community. In this context, this paper presents the first ultra-lightweight transmitter using side channel leakage based on electromagnetic emanation to send embedded IP identity discreetly and quickly. In addition, we present our electromagnetic test bench and a coherent demodulation method using slippery window spectral analysis to recover data outside the device. The hardware resources occupied by the transmitter represent less than 0.022% of a 130 nm Microsemi Fusion FPGA. Experimental results show that the demodulation method success to provide IP data with a bit rate equal to 500 Kbps.

I. INTRODUCTION

For digital circuit design the re-use of embedded intellectual properties (IP) is more and more important due to prohibitive cost of ASIC design. Nevertheless the IP business suffers from a lack a security due to the intrinsic form of IPs sales and exchanges. Many dedicated threats target the IP life cycle and result to revenues losses for the IP designers [START_REF] Colombier | Survey of hardware protection of design data for integrated circuits and intellectual properties[END_REF]. The IP threat model includes illegal re-use, illegal sales, cloning (illegal copy) of the IP. The extent of threats targeting IPs is linked to the type of IP: software IPs (typically hardware description language files), firmware IPs (synthesized netlist), and hardware IPs (FPGA bitstream or physical layout).

One of the solutions for the IP designers to protect their intellectual property is to be able to detect the presence of a copy of an IP embedded in a digital device by using IP identification. Works on IP watermarking and IP fingerprinting try to provide the IP identification service. But, most of the time the published solutions are not practical mainly because of the complexity of the watermarking/fingerprinting verification scheme [START_REF] Marchand | IP watermark verification based on power consumption analysis[END_REF]. Efficient IP identification scheme needs to be contactless, rapid and ultra-lightweight. Up to now, these three characteristics are not available in the state-of-the-art. To meet these requirements, in this paper we propose an ultra-lightweight binary frequency shift keying (BFSK) transmitter to forward IP identity (that could generated for example by a feedback shift-register or a physical unclonable function [START_REF] Colombier | Survey of hardware protection of design data for integrated circuits and intellectual properties[END_REF]) discreetly using an electromagnetic channel. Such circuit is usually called "spy circuitry". Using the electromagnetic channel, it is possible to contactless check the presence of an IP inside a digital device.

The rest of this paper is organized as follows. In Section 2, we present related works. In Section 3, we detail the proposed electromagnetic communication of data (i.e. IC/IP information). In Section 4, we present a method to analyze the electromagnetic spectrum to BFSK signal demodulation. In Section 5, we compare published spy circuitries that use a side-channel for IP protection and hardware Trojans. In Section 6, we present two industrial scenarios for the use of the proposed IP protection scheme and in Section 7 we present our conclusions.

II. RELATED WORK

Well-known threat in cryptographic engineering is side channel attacks [START_REF] Kocher | Differential Power Analysis[END_REF], [START_REF] Kamoun | Experimental Implementation of 2ODPA attacks on AES design with flash-based FPGA Technology[END_REF]. Most of the dynamic characteristics of both hardware and software implementations of cryptographic primitives can be used for side channel analysis: computation time, power consumption, electromagnetic radiation, optical radiation, even the sound produced during computation. These side channels can be used as transmission channels to send intellectual property data from a device or an embedded IP. For example in [START_REF] Marsh | Protecting designs with a passive thermal tag[END_REF], the thermal channel representing a contactless communication was used to transfer information from an embedded tag to a remote receiver. However the embedded thermal tag used in this commercial solution requires a relatively large area (255 Spartan-3 slices). In [START_REF] Ziener | Power Signature Watermarking of IP Cores for FPGAs[END_REF], the authors propose using two shift registers to generate a recognizable signature-dependent power consumption pattern to reveal the IP signature. Power consumption was also used in [START_REF] Becker | Side-channel based watermarks for integrated circuits[END_REF] to communicate the IP identity. To reinforce such work, the authors of [START_REF] Kerckhof | Intellectual Property Protection for FPGA Designs with Soft Physical Hash Functions: First Experimental Results[END_REF] propose using the power supply signal of an IP as a physical hash function for fingerprinting.

Related works can be found also in the malicious hardware design such as hardware Trojans. Such systems use side channels to forward secret information such as a symmetric cipher key from cryptographic hardware implementation [START_REF] Lin | Trojan Side-Channels: Lightweight hardware Trojans through Side-Channel Engineering[END_REF], [START_REF] Kutzner | TROJANUS: An Ultra-Lightweight Side-Channel Leakage Generator for FPGAs[END_REF], but also to cause or amplify side-channel leakage of cryptographic hardware [START_REF] Gallais | Hardware Trojans for Inducing or Amplifying Side-Channel Leakage of Cryptographic Software[END_REF].

Expect [START_REF] Marsh | Protecting designs with a passive thermal tag[END_REF], all the related works use power consumption as a communication channel which is not contactless. Unlike the proposed solution, all the related works are not lightweight and rapid as the section 5 of this paper will show.

III. PROPOSED EM COMMUNICATION OF IP DATA

Electromagnetic radiation is an efficient side channel since, unlike measurement of power consumption, electromagnetic radiation can be measured locally. One of the main advantages of this side channel is that it is impossible to hide the leak concerning electromagnetic radiation by using a global countermeasure. Moreover the electromagnetic test bench is not expensive (less than US$ 10K without an oscilloscope, which is the most expensive component). Last but not least, a spectral analysis of the electromagnetic radiation provides information on the oscillating structure such as a ringoscillator [START_REF] Bayon | EM leakage analysis on True Random Number Generator: Frequency and localization retrieval method[END_REF]. For all these reasons, we use the electromagnetic channel for our IP identification scheme. To this end, we designed an ultra-lightweight BFSK transmitter which could be activated outside the device by an ID checker or internally by a specific event (e.g. specific input sequence, internal data value, system state). Note that an enable signal is required to reduce the power consumed by the ring oscillator. Moreover, a permanently activated transmitter could be detected more easily by a spectral analysis of electromagnetic emanations of the device and could also cause local heating and premature aging of the chip.

BFSK is one of the common modulation schemes used in digital communication. The binary data are sent using a sinusoidal carrier at two frequency tones f 0 and f 1 , representing high ('1') and low ('0') logic levels. The binary data arriving at the transmitter input at certain bitrates determine the commutation of the tones at the transmitter output. The proposed BFSK transmitter uses a dedicated configurable ring-oscillator, as shown in Fig. 1. The configurable ringoscillator is designed using one multiplexor, N+K delay elements, and a feedback chain controlled by a NAND gate for activation of transmission to reduce power consumption. Actually, the transmitter is used only during a short time when the enable signal is high, and it consumes power only during this small piece of time. The power consumption of this transmitter is thus completely negligible.

Input data controls the multiplexor, as shown in Fig. 1. When input data is low, the ring oscillator uses N delays and its oscillation frequency is f 0 . When input data is high, the ring oscillator uses N+K delays and its oscillation frequency is f 1 . Since the ring oscillator's oscillation frequency decreases with an increase in the number of delay elements, frequency f 0 is higher than frequency f 1 . These two frequencies have to be selected according to the bandwidth of the electromagnetic analysis platform, which is used to acquire and measure the transmitted signal. The bandwidth of our test bench, which is described in Section 4, was limited to 100 MHz and 1 GHz by the low-noise amplifier.

The proposed BFSK transmitter was implemented in Microsemi FUSION flash based FPGA (130 nm CMOS technology) containing 600K logic gates (M7AFS600). The device contains 13 824 tiles, each tile can be used to implement one D-flip-flop or one configurable multiplexorbased logic block implementing any 3-input logic function.

The configurable number of delays in the ring oscillator of the proposed BFSK transmitter makes it possible to select precisely the two frequencies f 0 and f 1 using parameters N and K. Table I lists the ring oscillator frequencies and the number of Fusion tiles used by the BFSK transmitter for five values of N and K, with N ranging from 0 to 4, and K ranging from 1 to 5. According to Table I, f 0 can be chosen between 119 MHz (N=4) and 385 MHz (N=0) and f 1 can be chosen between 70 MHz (N=4, K=5) and 280 MHz (N=0, K=1). The exact value of f 0 depends on the number of delay elements, but also on the placement and routing of the transmitter. For the values N and K listed in Table I, the frequency variation was less than 1.7% (the maximum frequency deviation in Table I is The number of tiles used by the BFSK transmitter is very low, i.e. from 3 tiles (N=0, K=1) to 11 tiles (N=4, K=5). These values are equivalent to less than 0.022% and less than 0.080% of the total number of tiles included in the targeted 600K-gate FUSION FPGA, respectively. This very small number of tiles is very promising for good dissimulation of the BFSK transmitter inside the sea of gates/tiles. In order to estimate the number of resources needed for implementation with Xilinx SRAM FPGA or Altera SRAM FPGA, Table I gives the number of 4-input look-ups (LUT4) used by the BFSK transmitter with such FPGAs.

To evaluate the logical resources needed by the BFSK transmitter in ASIC implementations, the right hand column in Table I gives the number of equivalent gates (EG) in the transmitter. The gate count was estimated using the Virtual Silicon standard cell library based on the UMC L180 0.18 µm 1P6M Logic process (UMCL18G212T3 [START_REF]0.18 μm VIP Standard Cell Library Tape Out Ready[END_REF]). The delay gates are replaced by more efficient standard NOT gates. The gate count of a standard NOT gate is 0.67 EG, and that of the standard multiplexor, 2.33 EG. The standard NAND gate uses 1 EG. So the number of gates of the whole BFSK transmitter ranges from 4.67 EG (N = 0, K = 1) to 10.03 EG (N = 4, K = 5). Note that one flip-flop requires between 5.33 EG and 12.33 EG to store a single bit [START_REF]0.18 μm VIP Standard Cell Library Tape Out Ready[END_REF]. Such a transmitter is clearly ultra-lightweight in both FPGA and ASIC implementations. The small logical resources requirement of the proposed spy circuitry makes reverse engineering it harder, although not impossible [START_REF] Torrance | The state-of-the-art in semiconductor reverse engineering[END_REF]. Even with recent CMOS technologies, the attacker can reverse engineer ICs using a scanning electron microscope and an automatic tool for circuitry extraction [START_REF] Torrance | The state-of-the-art in semiconductor reverse engineering[END_REF], [START_REF] Subramanyan | Reverse Engineering Digital Circuits Using Structural and Functional Analyses[END_REF]. Nevertheless, the smaller the piece of hardware used for BFSK transmitter the harder it is to detect during reverse engineering. Detection of the transmitter using standard Trojan detection methods [START_REF] Agrawal | Trojan Detection using IC Fingerprinting[END_REF], [START_REF] Jin | Hardware Trojan Detection using Path Delay Fingerprint[END_REF] is not feasible because the transmitter does not change the data path of the circuit and because of the ultra-low signal-to-noise ratio on the electromagnetic channel, as shown in our experimental results below (Section 4).

IV. EXPERIMENTAL RESULTS

The electromagnetic radiation of the device was evaluated using the near-field electromagnetic analysis test bench described in [START_REF] Bayon | EM leakage analysis on True Random Number Generator: Frequency and localization retrieval method[END_REF]. The border between the far field and the near field can be considered to be about 23 mm from the device, depending on the hardware concerned. The most important part of the test bench is the acquisition chain. It determines the signal to noise ratio and measurement precision.

The chain, as presented Fig. 2, is composed of:  A Langer magnetic probe with a frequency range of from 30 MHz to 3 GHz and a spatial resolution of approximately 500 µm.  A Miteq low-noise amplifier with a frequency range of from 100 MHz to 1 GHz.  A LeCroy oscilloscope with a frequency range of up to 3.5 GHz and a sampling rate of up to 40 GS/s.

As presented in Fig. 2, the device to be tested (the board) is fixed to a XYZ table with repeatability of movement of 1 µm. The test bench, including the acquisition chain, XYZ table, FPGA configuration and power supply variations, is controlled by a computer. Electromagnetic analysis of IC is contactless, local, and can be spatial or/and temporal. This last point makes it possible to perform frequency analysis of the electromagnetic emanation. In the your bench the spectral range is limited to 100 MHz and 1 GHz. Standard devices aimed at direct BFSK demodulation cannot be used for these relatively high frequencies. Available integrated BFSK demodulators are limited to a few dozen megahertz. For this reason, we developed a dedicated BFSK demodulation scheme for our needs, in which a spectral analysis of the low noise amplifier output (a component of the test bench) is performed to measure the f 0 and f 1 spectral contribution. The transmitted high (low) level is detected when f 1 spectral contribution is higher (lower) than that of f 0 . Fig. 3 illustrates the spectral analysis of the BFSK transmitter's electromagnetic emission when N = 1 and K = 1, which corresponds to a small transmitter with high frequencies in Table I. For the spectral analysis, a 16 384 points FFT was computed. Fig. 3 presents a zoom (X and Y axis) on the global spectrum of the local EM emanation of the circuit when the BFSK transmitter sends a '0' (in blue) and when the BFSK transmitter sends a '1' (in red). Notice also that we placed a small antenna in the close vicinity of the ring. The amplitude of the spectral rays at 180, 200, 220, 240, 260, 280 MHz are very high compared to the two spectral rays at 189.2 MHz and 272.2 MHz. However, we have cut the upper part of the spectrum in order to see the two interesting rays (at 189.2 MHz and 272.2 MHz). These frequencies correspond to the spectral contribution of the BFSK frequencies f 0 and f 1 . Actually, according to Table I, for Microsemi FUSION FPGA, the frequency f 0 is around 272 MHz and f 1 is around 189 MHz. Slight variations in the frequency values (around the values presented in Table I) are due to environmental variations (mainly temperature). These variations are too small to jeopardize the success of demodulation because the BFSK designer knows the chosen parameters N and K and targeted frequencies f 0 and f 1 for transmission of low and high levels. Without knowledge of the BFSK parameters, the electromagnetic transmission cannot be easily detected because it cannot be distinguished from spectral noise. The signal-to-noise ratio of the BFSK transmission in Fig. 3 is -135 dB for a 1 GHz bandwidth. Such an ultra-low SNR represents efficient protection against unwanted BFSK transmitter detection via a side channel. However, knowing the N and K parameters, the BFSK designer can calibrate the demodulation (determine the two frequencies) by electromagnetic analysis of the ring oscillators based on the differential spectral analysis as described in [START_REF] Bayon | EM leakage analysis on True Random Number Generator: Frequency and localization retrieval method[END_REF].

The spectral contribution of the two BFSK frequencies during transmission (which is limited by the transmitter enable signal) at the low-noise amplifier output is measured to determine the transmitted bit sequence. In order to apply the demodulation technique described here, called coherent demodulation, theoretically the bitstream rate (bitrate) is limited to 0.5 times the frequency difference between f 0 and f 1 . In the case shown in Fig. 3, theoretically the maximum bitrate is 41.5 Mbps (theoretically limited to 53 Mbps when N = 0 and K = 1). The coherent demodulation scheme cannot be used for higher bitrates, instead more sophisticated and more expensive demodulation techniques should be used. To reduce the size of the BFSK transmitter, a time-interleaving transmitter should be used to increase the bitrate. However, the high transmission rate is not so important, because data (e.g. IP identifier or stolen data) can be stored at a higher data rate before being transmitted at a slower rate. The small area and unobtrusiveness of the BSFK transmitter are much more important properties than transmission rate.

For the coherent demodulation of the electromagnetic radiation, we propose a slippery window spectral analysis. Indeed, overall spectral analysis masks the effects of the nonstationarity of the signal and therefore provides no information about its temporal evolution. Slippery window spectral analysis is a three-dimensional representation of the signal: amplitude, frequency, and time. It requires two quantities Fw, the width of the FFT window frame and the difference Δτ between two frames. For our experiment, we chose Fw equal to 16 384 points (2 14 -point FFT) and Δτ equal to 100 points. For each frame, the FFT provides the software demodulator with the amplitude of signals f 0 and f 1 which enables the demodulator to distinguish between a transmitted '1' or '0'.

To illustrate data transmission from the circuit via the EM channel, we used a shift register that stored the following 16-bit sequence: "1011110111011100". The clock frequency of the shift register is 500 kHz. When the enable signal of the transmitter is given, the sequence is sent cyclically to the BFSK transmitter, which transmits it via the electromagnetic channel. Fig. 4 gives the result of the coherent demodulation obtained at a 500 Kbps bit rate, which served as a proof of concept.

16-bit sequence 16-bit sequence

V. COMPARISON WITH STATE OF THE ART SPY CIRCUITRIES USING A SIDE-CHANNEL

Table II compares the implementation of the proposed ultralightweight BFSK transmitter with other recently published state of the art methods. Table II gives the spy circuitry application (App.) for each reference; this may be IP protection (IPP) or hardware Trojan (HT) or both (for the presented work, PW). In addition, Table II gives the year of publication (YoP), the side channel used, the hardware resources required only for the leakage generator (for example we do not take the hardware used for IP watermark generation or the Trojan's payload into account). Unfortunately, the principles compared do not use the same hardware. For the sake of correctness, we give the implementation results as they are presented in the referenced papers. Nevertheless, the implementation bitrate of these previously published works can be roughly compared with our proposed solution. Based on published data, we computed the bitrate of all the proposals by using the number of clock cycles needed to send information via the side channel. For all the references presented in this table, the bitrate was computed using a 500 KHz frequency for data synchronization (same frequency is used for data transmission in Fig. 4).

As can be seen in Table II, the proposed work reaches the highest bitrate. The reason for such a good performance is first that we use a spectral analysis of the local electromagnetic leak based on a simple frequency modulation. Except for [START_REF] Kamoun | Experimental Implementation of 2ODPA attacks on AES design with flash-based FPGA Technology[END_REF], all the other solutions use a global measurement of power consumption, which reduces the signal-to-noise ratio of the information leaked via the side channel. Our proposal is clearly the smallest spy circuitry ever published. Although solutions based on circular shift-registers are well adapted to last generation FPGA families, since the 16-bit shift registers can be designed using only one look-up table, they are not suitable for ASIC technologies. Currently, an ASIC implementation of a 16-bit shift register requires 16 flip-flops whereas the solution we propose occupies an area equivalent to only one D-flip-flop.

In this paper, we present the proposed spy circuitry for IP protection, but it can also be used for hardware Trojan. Most of the other proposals could also be used for both applications. Note that in 2012, Kasper et al. proposed to use the work initiated in [START_REF] Lin | Trojan Side-Channels: Lightweight hardware Trojans through Side-Channel Engineering[END_REF] for hardware Trojan or IP watermarking implementation [START_REF] Kasper | Side Channels as Building Blocks[END_REF]. However, by using electromagnetic emanation and a configurable ring oscillator, the proposed solution is the most convincing for industrial applications (e.g., those aimed at IP protection) because of its very small area and high bitrate. The first scenario is the identification of embedded IP in the supply chain. This identification is used in order to be sure to don't use counterfeiting (fake) devices. In recent years, the issue of counterfeiting of integrated circuits has increased considerably. For example, the number of counterfeit electronic circuits collected by U.S. Customs from 2001 to 2011 has increased by around 700 [START_REF] Gorman | Counterfeit Chips on the Rise[END_REF]. Between 2007 and 2010, U.S. Customs collected 5.6 million counterfeit electronic products [START_REF] Agma | Alliance for Gray Markets and Counterfeit Adatement[END_REF]. Overall, the estimation of counterfeiting is about 7% of the semiconductor market [START_REF] Pecht | Bogus! Electronic manufacturing and consumers confront a rising tide of counterfeit electronics[END_REF] which represented a loss of about $ 22 billion in 2014 for the legal industry.

It is therefore crucial and strategic to be able to detect counterfeit IC as soon as possible in the supply chain (this is particularly crucial for military and space grad devices). Fig. 5 shows a possible framework to manage the device under test (control the enable signal) and check the IP identification by using an EM probe, an amplifier and a dedicated acquisition system including a spectral analysis and the proposed demodulation method. Due to the high bit rate of the proposition solution the identification of the ID requires less than 1ms (with 500Kps bit rate as in Fig 4 .). This counterfeiting detection could be completed by other physical (invasive or not) and optical inspection [START_REF] Tehranipoor | Counterfeit Integrated Circuits -Detection and Avoidance[END_REF]. The second scenario occurs when a IP designers would like to verify the illegal presence of its IP inside a device (ASIC or FPGA). In this case the proposed transmitter provides to the Validity check ok identification scheme a data like a PUF [START_REF] Bossuet | A PUF based on a transient effect ring oscillator and insensitive to locking phenomenon[END_REF] or a watermarking. Watermarking is a technique of steganography which provides the ownership of an IC (or an IP) by checking for the presence of hidden information called the watermark [START_REF] Kocher | Differential Power Analysis[END_REF]. Most of the watermarking methods proposed in the literature need a complex verification scheme. Nevertheless it is possible to use power consumption as proposed in [START_REF] Kerckhof | Intellectual Property Protection for FPGA Designs with Soft Physical Hash Functions: First Experimental Results[END_REF] and [START_REF] Kocher | Differential Power Analysis[END_REF] but it is easy and cheap to use global countermeasure in order to mask the power consumption due to the watermark [START_REF] Kamoun | Correlated Power Noise Generator as a Low Cost DPA Countermeasure to Secure Hardware AES Cipher[END_REF]. Using electromagnetic emanation in this scenario is better because as electromagnetic emanation is local it is really hard to mask it by using a global countermeasure. Moreover, in this paper we have shown that due to the SNR of BFSK signal, it is unrealistic to try to detect it without the precise knowledge of the used frequencies for data transmission.

VII. CONCLUSION

IP protection has become crucial topics for hardware security due to the lack of trust in IP market. In this paper, we have presented an ultra-lightweight transmitter of IP identity using the electromagnetic side channel. Based on a configurable ring oscillator, our solution exploits a BFSK signal to transmit information by way of the electromagnetic channel. By performing a slippery window spectral analysis of the near field electromagnetic emanations captured locally over the BFSK transmitter circuitry, the proposed transmission achieves a high bitrate (experimentally 500 Kbps and theoretically limited to 53 Mbps with a Microsemi Fusion FPGA), which has not been achieved before. Moreover, the transmitter occupies very small area. For the highest frequency and data rate, our solution requires 4.67 equivalent gates, which it is less than the requirement of a small D-flip-flop. Such a small requirement of logical resources makes reverse engineering of the chip very difficult and detection of the transmitter using standard Trojan detection methods is not feasible.

2

 2 MHz when N = 4).

Fig. 1 .

 1 Fig. 1. Architecture of the ultra-lightweight digital BFSK transmitter based on a configurable ring oscillator.

Fig. 2 .

 2 Fig. 2. Near-field Electromagnetic analysis test bench.

Fig. 3 .

 3 Fig. 3. Parts (zoom in x-axis and y-axis) of the electromagnetic emanation spectral analysis of the proposed BFSK transmitter when it sends a high level at f1=189.2MHz (in red) and when it sends a low level at f0=272.2 MHz (in blue).

Fig. 4 .

 4 Fig. 4. Coherent demodulation of the 16-bit code word cyclically sent by the BFSK transmitter at 500 Kbps, the top (bottom) graph corresponds to the spectral contribution of signal at f1 (f0).

  VI. INDUSTRIAL SCENARIOS USING THE PROPOSED IPPROTECTIONAccording to the previous section, in comparison with other works, our propose goes clearly towards using a spy circuit in an industrial context for IP protection. Two industrial scenarios are presented in the following.

Fig. 4 .

 4 Fig. 4. Rapid and contactless IP identification in the supply chain by using EM transmission of IP ID.

TABLE I

 I

			. HARDWARE IMPLEMENTATION RESULTS OF THE
			BFSK TRANSMITTER
	N	K	f0 (MHz)	f1 (MHz)	Fusion Tiles	LUT4	EG
		1	385	280	3	2	4.67
		2	383	210	4	3	5.34
	0	3	384	151	5	4	6.01
		4	385	130	6	5	6.68
		5	381	111	7	6	7.35
		1	272	189	4	3	5.34
		2	272	156	5	4	6.01
	1	3	270	120	6	5	6.68
		4	271	106	7	6	7.35
		5	269	93	8	7	8.02
		1	168	144	5	4	6.01
		2	169	124	6	5	6.68
	2	3	169	100	7	6	7.35
		4	168	91	8	7	8.02
		5	168	79	9	8	8.69
		1	146	128	6	5	6.68
		2	147	112	7	4	7.35
	3	3	146	92	8	5	8.02
		4	145	84	9	6	8.69
		5	144	74	10	7	9.36
		1	123	110	7	6	7.35
		2	121	98	8	7	8.02
	4	3	122	83	9	8	8.69
		4	121	77	10	9	9.36
		5	119	70	11	10	10.03

ACKNOWLEDGMENT

The work presented in this paper was realized in the frame of the SALWARE project number ANR-13-JS03-0003 supported by the French "Agence Nationale de la Recherche"and by the French "Fondation de Recherche pour l'Aéronautique et l'Espace".