K Safa

G Gary

Accounting for the effect of local punching at the bar/ specimen interface in SHPB experiments

The object of this work is to provide a displacement correction for the local punching due to an axial loading at the end of a bar. For this purpose, an analytical calculation of the dynamic 3-D indentation of an elastic isotropic semi-infi nite bar is carried out. It provides a fi rst-order correction for the displacement obtained through the 1-D wave analysis commonly used in SHPB processing. This correction applies to the displacements as they are calculated at both ends of a specimen sandwiched between the two bars of a SHPB device. It improves the results obtained for the dynamic behavior of the specimen, in particular at early instants of loading.

I. INTRODUCTION

The SHPB device, called also Kolsky apparatus, is widely used to measure the stress-strain behavior of materials at high strain rates of loading. This behavior is derived from the measured forces and velocities at the specimen faces. These quantities are obtained through the knowledge of the axial strains that are recorded by strain gauges glued on the lateral faces of the input and output bars (see formulas (I) underneath). These four quantities (forces and speeds) provide redundant measurements allowing for the use of inverse methods to investigate the specimen behavior [l]. A particular and simple use of such methods is to compute (using a 1-D transient calculation) the initially elastic response of the specimen. This method [START_REF] Zhao | On the use ofSHPB technique to determine the dynamic behaviour of the materials in the range of small strains[END_REF] was proposed for an optimized determination of the relative position of the origin of the three waves involved in the SHPB. It has been used for years in our laboratory where we have observed that the best value of the Young's modulus of the specimen was always smaller (in particular with small diameter specimens) than the expected or known value. We have suspected that this systematic error was due to an imperfect measurement of the average strain of the specimen. The bar face in contact with the specimen (when it has a smaller diameter than the bar) does not remain plane as it is non uniformly loaded, i.e. submitted to a local "elastic punching". The corresponding displacement fi eld is restricted to the impact zone and is not recorded by the strain gauges located far from the bar end.

We are interested in the evaluation of the non-planar deformation of the loaded surface induced at the end of a bar by a central axial load. The elementary 1-D theory of wave propagation in bars, and more advanced theories that consider wave dispersion due to lateral inertia effects in the bars [START_REF] Skalak | Longitudinal impact of a semi-infi nite circular elastic bar[END_REF][START_REF] Folk | Elastic strain produced by sudden application of pressure to one end of a cylindrical bar. I. Theory[END_REF][START_REF] Jones | Axially symmetric cross-sectional strain and stress distributions in suddenly loaded cylindrical elastic bars[END_REF], do not tackle this problem. They rather investigate wave propagation at large distances from the bar end [START_REF] Kennedy | Longitudinal wave propagation in a circular bar loaded suddenly by a radially distributed end stress[END_REF]. As suggested by some authors [START_REF] Gama | Hopkinson bar experimental technique: A critical review[END_REF], the determination of the local displacement at a bar end is of practical importance for the use of the SHPB device. This displacement has been determined through a study described in details in [START_REF] Safa | Mise au point d'un essai de compaction dynamique. Application au beton[END_REF]. Because the restricted number of permitted pages, only the basis and the general ideas are recalled in the present paper.

Numerical and experimental studies are also carried out. They provide estimations of the consequence of the proposed displacement correction on the dynamic response of a specimen. The correction is especially significant at the beginning of the loading where it permits an estimation of the Young's modulus of the tested material.

FORMULATION OF THE PROBLEM

Experiment

In a typical SHPB device (for compression testing), three waves are involved during the experiment: an incident compressive wave E;(a, t) generated by the impact of the striker, a reflected (tensile) wave E,(a, t) due to the lower impedance of the specimen (both recorded at gauge A abscissa a) of the input bar), and a transmitted compressive wave E1(b, t) recorded at gauge Babscissa b) of the output bar).

By use of wave theory, the waves measured at A and B are obtained at the bar-specimen interfaces where they are denoted E,(t), E,(t) and E1(t) [START_REF] Kolsky | An investigation of mechanical properties of materials at very high rates of loading[END_REF].

The measured forces and velocities are:

F,(t) = A,E(E;(t) + E,(t)) V; (t) = -Co; ( E; (t) -E, (t))
at the input interface

F0(t) = A0EE1(t) Va(t) = -co, E 1(t)
at the output interface

(1)

Where c0; and c0, are the velocities of the longitudinal waves and Ai and A0 are the cross-sectional areas of the input and output bars, respectively.

In Figure 1, the sample is shown together with the output bar, with a deformed confi guration of the bar end under a distributed force F(t). The bar end displacement derived from relations [START_REF] Rota | Application de methodes inverses au depouillement de l'essai aux barres de Hopkinson[END_REF] is u(t). The additional displacement p(t) resulting from the local elastic 3-D axisymmetric deformation of the bar is referred to as the elastic punching. It appears when the diameter of the specimen is smaller than that of the bar. The determination of the elastic punching p( t) provides a correction for the displacement u(t) given by SHB formulas, allowing for the measurement of the real displacement at the specimen-bar interface: d(t) = u(t) + p(t).

Theoretical problem

Following the usual SHPB practice, the friction between the sample and the bar faces is neglected and the uniaxial stresses within the specimen are assumed to be uniformly distributed through the cross section.

Based on these two assumptions, the problem becomes that of a semi-infinite bar with traction free lateral faces, axi-symmetrically loaded at its end by a uniform distribution of time dependant normal stress. The resulting force F( t) is assumed to be known, as obtained through relations (1) in the case of a SHPB testing.

The problem can be divided in two parts. At early instants, the waves that are produced by the loading are the same as they would be in a semi-infi nite solid, as long as the refl ected waves at the free cylindrical surface of the bar are not superposed with them. The solution of the first part is then valid for a very short time (if compared with the SHPB test duration). An explicit solution is proposed in [START_REF] Safa | Mise au point d'un essai de compaction dynamique. Application au beton[END_REF]. A short time later, the radial oscillations near the bar end tend to produce a quasi static state in the radial direction while the early axial waves produced by the loading are already far from the end of the bar. The solution of the corresponding problem provides a sufficient correction for SHPB applications.

SOLUTION OF THE PROBLEM

Following the notations of fi gure 1, the corrected displacement for any point under the sample is defi ned by the relation d( t) = u( t) + p(t) where p(t) represents the "elastic punching".

A trick used to solve this problem is based on the superposition principle. Two loadings are applied simultaneously at the end of the bar. These loadings are chosen so that they result in opposite forces F( t) and -F( t). (Figure 2). In order to determine the correction p(t), the dynamic problem is transposed into the static configuration of a self-equilibrated end-loading problem of a semi-infinite cylinder. This configuration is able to provide at any time t the mean value of the correction p(t) to be added to the measured displacement u(t). The method and subsequent calculations are detailed in [START_REF] Safa | Mise au point d'un essai de compaction dynamique. Application au beton[END_REF]. The exact analytical solution obtained is not easy to use for SHPB applications.

In a second step, it is then searched (as also described in [START_REF] Safa | Mise au point d'un essai de compaction dynamique. Application au beton[END_REF]) for an approximate expression in a closed-form that will be more convenient in the frame of SHPB practice. This step leads to the following formula.

Where,

4 1-v F 2 2 ( p = ----[2p -(p + 1 )E(p ) -(p -1 )K p )]
(2) 3.n 2 µ a p a a a a a a pis the average displacement under the loaded area, a is the radius of the loaded area, Pa is the ratio of the radius of the loaded area versus that of the bar (r/a) Fis the applied (known) force, µ and v are Lame's coefficients Functions E( pa), K( pa) are well defined in the literature. They have the following expressions:

E( p a ) = (' V 1 -p� sin 2 () d() ; K( p a ) = ['t d() J o J o j 1 -p� sin 2 ()

NUMERICAL SIMULATION

A dynamic numerical simulation is performed with Abaqus explicit code, for a cylinder of diameter 50 mm submitted to uniform distributions of pressure over diameters of 15 and 25 mm, due to the same force magnitude. The evolution of the loading force with time is given in the insert (Figure 3). Figure 3 shows the average displacement as calculated with formula (2) and as obtained numerically, and fi nally the ID displacement at large distance from the impact zone as it is derived from the standard 1-D wave theory: u( t) = -R} r01 F (r ) dr. This simulation confirms that the closed-form solution is almost equal to the one given by the numerical simulation. The test lasts 85 µ S and simulates the behavior of a brittle material. We purposely have chosen a short loading duration to conveniently visualize the obtained solutions. In the present case, it appears that the amplitude of the correction is in the range of the displacement itself up to a force equal to 40 kN (250 or 90 M pa for a specimen of 15 mm or 25 mm respectively). The case of a 25 mm diameter specimen loaded by the same resulting force indicates that the correction for punching becomes rather small in this case (a/R greater than 0.5). The distance between the corrected displacement d( t ) and the 1-D displacement u( t) given by the SHPB, corresponds to the elastic punching p( t). It vanishes when a= R.

PRACTICAL IMPLEMENTATION OF THE DISPLACEMENT CORRECTION

In Figure 1 a scheme of the shape of one bar face is displayed (for any given time). lf we consider the common case of two identical bars and if we denote by 10 the initial length of the specimen, its current length becomes 10 -2[u(t) -p( t)]. The value 10 -2u (t) given by the SHPB device is therefore under-estimating the length of the specimen. Consequently the real strain of the specimen is less than that given by the classical SHPB analysis.

Note that Eq. 2 may be written in the form:

The punching effect is then equivalent to that produced by a spring of stiffness l/Kr (smaller than that of the bars) acting between the bars and the specimen. It explains that SHPB processing without punching corrections leads to an underestimated value of the Young's modulus of the specimen.

Taking into account punching, the classical formulas (1), remain unchanged concerning the forces. The formulas, at the input interface (left) and at the output interface (right), are:

F;(t) = A;E(E;(t) + E,(t)) dF;(t) V;(t) = -c0;(E;(t) -E,(t)) + � Kp; ( 4 )
In the common particular case of identical input and output bars and a cylindrical specimen in equilibrium, formulas giving the stress are unchanged.

Formulas giving the nominal strain are corrected in the following manner:

2Kpa(t)S s E corr. = E sHPB + f s
where Ss and ls are the cross-sectional area and the length of the specimen, respectively.

EXPERIMENT AL ILLUSTRATION AND APPLICATIONS

(5)

In the processing procedure of SHPB, we use [START_REF] Gary | Correction de dispersion pour !'analyse des petites deformations aux barres de Hopkinson[END_REF], the waves are first transported to the bar ends, taking into account the dispersion. The incident wave being known, we assume that the specimen has a purely 1-D elastic behavior. Taking the real specimen dimensions (area, length and density), a 1-D elastic calculation with a chosen Young's modulus, using the incident wave at bar end as input data, builds the reflected and transmitted waves as they would be produced by the elastic specimen. This calculation is very fast (it appears instantaneous). The "best" Young's modulus is then chosen to find the best match in shape between simulated waves and the beginning of (shifted) real ones. This procedure, called "elastic simulation" is described in [START_REF] Zhao | On the use ofSHPB technique to determine the dynamic behaviour of the materials in the range of small strains[END_REF], and is used for a more precise shift for the waves to be processed. This shift correction is always very small (a few microseconds) and has no influence on the dispersion effects. The important point, for the present paper, is that it allows for the measurement of the apparent elastic modulus. When used without taking account of the punching correction, it leads to an underestimated value of Young's modulus. This is illustrated in Figures 4 to 6 below.

It is observed that the expected real modulus of steel corresponds to waves that are different of real ones and that a (two times) lower modulus produces waves in good agreement with the real ones. This is due to the punching effect. When a test is processed with a specimen that has the same diameter than the bar, there is no punching effect and one recovers the expected modulus.

Figure 6 shows that the punching correction is worth being done. It particularly applies when the specimen diameter is smaller than that of the bars and when the total strain of the test is in a reduced range, as it is the case when low strain-rates are required. In the present test the average strain-rate was around 100/s.

CONCLUSION

The problem of the local punching of a bar axially loaded at one end has been investigated. When a known axial stress is dynamically applied on a circular central part of a bar end, it induces an axial displacement which varies along the radius of the bar. This displacement, depending on the mechanical parameters of the bar and on both parameters of the loading (stress and diameter), is closely approximated by a simplified formula. When applied to SHPB testing, it allows us to calculate the local punching induced by a specimen smaller than the bars. Consequently, it improves the measurement of material properties at low strains, and in particular provides an evaluation of the Young's modulus of elastic materials.

slrain edwi s 0 .2e-3 ----• • • • • ---------• • -------• • • • -------• • • •-------• • • • ------• • • • -0.2e-3 -----•• ••------••--; -----•• liiiSirilllid• : ------• •••••-----••••

Figure 1 .

 1 Figure 1. Schematic view of displacements at the output bar end during a SHPB experiment.

  )=�(!)/JT(I�/ / , : q(t)=-F(t)/7TR2 /,,,,./ u(t)=O

Figure 2 .

 2 Figure 2. Determination of the local displacement p,,,.11 (t) by superposition of states of stress.

Figure 3 .

 3 Figure 3. Numerical, analytical and lD average displacements under disks of pressure of diameters 15 and 25 mm, applied on a 50 mm diameter cylinder, according to the same time-dependant function.

Figure 4 .

 4 Figure 4. Basic recorded waves of a real test.Figure 5. Waves at bar ends.

Figure 5 . 3 Figure 6 .

 536 Figure 4. Basic recorded waves of a real test.Figure 5. Waves at bar ends. (Steel bars, diameter 20 mm, striker speed 2.8 m/s, steel specimen: length 6.5 mm, diameter 5.1 mm) "incident", "reflected", "transmitted" (thin lines) show real waves at bar ends. "tr-sim l" and "ref-siml" show simulated transmitted and reflected waves produced by a purely elastic specimen with Young's modulus= 1.10 11 Pa "tr-sim2" and "ref-sim2" show simulated transmitted and refl ected waves produced by a purely elastic specimen with Young's modulus= 2. 10 " Pa (as expected)

1-v [ J

-2 pa -(p� + l)E(pa) -(p� -l) K (pa) F(t) = K pF(t) n a p � (3)Where K P is function of the specimen, elastic properties and diameter of the bar Therefore K P is constant all along the experiment and from (3) we can deduce that l /Kp is homogeneous to a spring stiffness.


Manuel Trotta
email: manuel.trotta@uca.fr

Dominique Feillet
email: dominique.feillet@mines-stetienne.fr

Alain Quilliot
email: alain.quilliot@isima.fr

Pickup and delivery problems with autonomous vehicles on a ring

Keywords: vehicle routing, mobility, autonomous vehicles

Introduction

It is a matter of a fact that urban mobility (by means of public transport or personal cars) has become a key factor in people's everyday life, making it easier. But urban mobility is also the source of some issues, some of which are traffic congestion, air pollution and noise. Some of these problems may be solved by replacing normal vehicles with autonomous cars, minibuses, and shuttles. However, the use of autonomous vehicles raises legal, ethical, economic and safety issues. Due to these problems it is not likely that autonomous vehicles will totally replace normal vehicles very soon. However, they will probably first be authorized for collective transportation.

There already exist few cases where fully autonomous vehicles are used in public transportation. A collective transportation service is active in the Saclay Plateau neighborhoods. An autonomous Transdev-Lohr i-Cristal shuttle serves four stops, providing users with transportation between the Massy station and the Saclay Plateau to supplement the existing train, regional express, coach and bus services.

On 17 December 2018, Keolis and the European Metropole of Lille launched an electric autonomous shuttle service at the university of Lille in Villeneuve d'Ascq. The service employed two NAVYA electric autonomous shuttles for one year serving four dedicated stops on a 1.4 km route and providing connections to two metro stations. For a recent benchmark study on experimentations with Autonomous Shuttles for Collective Transport, the reader is referred to Antonialli [START_REF] Antonialli | International benchmark on experimentations with autonomous shuttles for collective transport[END_REF].

In these situations there is often a need to schedule the transportation requests that the autonomous vehicles must serve while moving on circular routes (we will call them rings). Pickup and Delivery and Dial-a-Ride problems on circular tracks have already been studied in a certain number of papers, such as [START_REF] Gendreau | Heuristics for the traveling salesman problem with pickup and delivery[END_REF][START_REF] Ilani | A fixed route dial-a-ride problem[END_REF][START_REF] Pimenta | Models and algorithms for reliabilityoriented dial-a-ride with autonomous electric vehicles[END_REF][START_REF] Tzoreff | The vehicle routing problem with pickups and deliveries on some special graphs[END_REF].

The problems we want to study in this paper have the following structure. There are m stations arranged in a circle and served by some capacitated vehicles that can travel in either one or both directions. The vehicles must serve some transportation requests. A transportation request (i, j) (i, j = 1, ..., m, i = j) requires a vehicle to reach station i (if it is not already there), pick up a certain number of people or objects and then move them to station j. The number of requests they can serve at the same time depends on their capacity. An additional station (0) is added, that can be thought of as the depot, i.e. a station where the vehicles stay before starting the service and after serving all the requests. Due to the circular layout, to go from station i to station j the vehicles must pass through all intermediate stations between i and j. There can be identical requests. From this common structure, several different problems can be defined.

The objective of all problems is to minimize either the maximum completion time or the total completion time subject to different constraints, such as the the vehicle capacity, release dates or due dates.

Our contributions

We propose a classification scheme for these problems, insipired by the three-field notation for scheduling problems introduced by Graham et al. [START_REF] Graham | Optimization and approximation in deterministic sequencing and scheduling: a survey[END_REF]. It is made up of three fields separated by a vertical bar. Each field may be a comma separated list of symbols that describe one or more features of the problem.

We investigate the complexity of different variants, presenting some complexity results. We prove the complexity of some them by showing reductions from known problems, while for some others we describe new complexity proofs.

We show the equivalence of some variants, in the sense that each can be polynomially reduced to the other and thus solved with the same algorithm. We also show the relations between some problem variants and arc routing, graph coloring and scheduling problems.

For polynomially solvable problems we present some efficient exact algorithms, while for some NP-hards problems we propose mathematical formulations and greedy heuristics. We also compare these heuristics with integer programming models solved with CPLEX and we show the experimental results.


Maelig Ollivier
email: maelig.ollivier@cea.fr

Laurence Latu-Romain
email: blaurence.latu-romain@cea.fr

A Mantoux

M Martin

T Baron

V Soulière

G Ferro

E Bano

Elaboration of core Si / shell SiC nanowires

Keywords: silicon carbide, carburization, FIB-SEM, nanowires, core/shell, epitaxy, biocompatibility

Silicon nanowires obtained by a top-down approach have been carburized at high temperature and atmospheric pressure with two different gaseous precursors: CH 4 and C 3 H 8 . These processes reveal core Si / shell 3C-SiC nanowires. After being characterized by SEM, FIB-SEM and TEM microscopies, the 3C-SiC layer has been used as seed layer for the growth of epitaxial 3C-SiC on the nanowires. Preferential growth of 3C-SiC on the sidewalls of nanowires has been observed. Thanks to the biocompatibility of SiC compared to Si, this layer could act as a protective shell for biosensors based on Si nanowires transistor.

Introduction

3C-SiC is an important material for the fabrication of electronic devices because of its high physical properties: it has a high thermal conductivity, a high mechanical strength, a high breakdown electric field, and a high electronic mobility [START_REF] Morkoc | Large-band-gap SiC, III-V nitride, and II-VI ZnSe-based semiconductor-device technologies[END_REF]. Therefore it can operate at high power, high temperature and in harsh conditions. In the meantime, the reduction of the size of devices, and particularly one-dimensional (1D) nano-objects, has an increase interest for industrial manufacture and academic research [START_REF] Lu | Nanoelectronics from the bottom up[END_REF]. In the case of 3C-SiC, several groups have successfully synthesized SiC nanostructures such as SiC nanowires (NWs) [START_REF] Attolini | Synthesis and characterization of 3C-SiC nanowires[END_REF], core SiC / shell C NWs [START_REF] Hyun | Growth of 3C-SiC nanowires on nickel coated Si(100) substrate using dichloromethylvinylsilane and diethylmethylsilane by MOCVD method[END_REF], SiC nanotubes (NTs) [START_REF] Latu-Romain | From Si nanowire to SiC nanotubes[END_REF][START_REF] Ollivier | Carburization of Si microwires by chemical vapour deposition[END_REF]. However, there are still some challenges before the elaboration of SiC nanostructures suitable for making good components [START_REF] Zekentes | SiC nanowires: material and devices[END_REF]. Indeed all the above methods lead to the same result: 3C-SiC is mainly crystalline and has a high density of twins and stacking faults.

In this study, core Si / shell 3C-SiC NWs have been elaborated by carburization of Si NWs. These as-obtained core Si / shell SiC NWs have been used to grow by epitaxy a thick layer of 3C-SiC. Each step of the synthesis has been characterized either by SEM or by TEM.

Experimental details

For this study, Si NWs are obtained by a SF 6 /O 2 plasma etching of a (100) Si substrate. They are 3.4 µm long with a diameter of 200 nm. Fig. 1 shows SEM image of etched Si NWs before carburization. It can be noticed that the sidewalls display a non-negligible roughness, compared to Si NWs obtained by VLS method, due to the etching process itself. However, these etched Si NWs present major advantages such as high reproducibility of the fabrication process, single crystallinity and absence of catalyst. The carburization is performed in a CVD furnace at 1100°C and atmospheric pressure [START_REF] Latu-Romain | From Si nanowire to SiC nanotubes[END_REF][START_REF] Soueidan | Effect of growth parameters on the heteroepitaxy of 3C-SiC on 6H-SiC substrate by chemical vapour deposition[END_REF]. The main flow is composed of H 2 or H 2 -diluted-in-Ar. The H 2 gas leads to a reductive atmosphere useful to remove native oxide at high temperature. CH 4 or C 3 H 8 highly diluted in the main flow can be used as carbon precursors and are sent at 800°C. This high dilution allows limiting the carbon deposit in favor of the carburization of Si. After a plateau which lasts from 1 to 60 minutes, the precursor gas flow is switched off and the system is naturally cooled down.

For the 3C-SiC epitaxial growth step, C 3 H 8 and SiH 4 are introduced into the reactor and the system is heated up to 1300°C. After 2 min of growth, precursors are removed and the system is naturally cooled down under H 2 until ambient temperature [START_REF] Soueidan | Effect of growth parameters on the heteroepitaxy of 3C-SiC on 6H-SiC substrate by chemical vapour deposition[END_REF].

The obtained NWs are morphologically and structurally characterized with SEM (Zeiss Ultra Plus) and TEM (FEI Tecnaï). A dual focused ion beam -scanning electron microscope (FIB-SEM, FEI Helios 450S) is used to prepare the thin lamella for TEM characterizations of the as-grown SiC layer on carburized Si NWs.

Results and discussions

Carburized Si NWs have the same roughness than before carburization (Fig. 2a), and the overall morphology is kept, such as length and diameter.

The thickness of the 3C-SiC layer, measured by TEM, is about ε ≈ 3 nm and the 3C-SiC layer seals entirely the surface (Fig. 2b). This thickness is in accordance with the thickness measured on Si (100) surface. Si core / SiC shell NWs have been elaborated.

Then, the role of the carburization time is studied on the thickness ε of the 3C-SiC layer. In Fig. 2c, the thickness of the 3C-SiC layer as a function of the squared root of the carburization time is reported. It appears a linear relationship between the thickness and the square root time which is typical of a limited-diffusion process. In the case of the carburization of Si, it has been shown that Si is the diffusing specie through the SiC layer. In a previous study, this phenomenon has been used to produce SiC nanotubes from carburization of Si NWs [START_REF] Latu-Romain | From Si nanowire to SiC nanotubes[END_REF].

According to the second Fick's law, and the Einstein-Smoluchowski's relation, the depth ε of the Si surface affected by the carburization is given by Eq. 1:

ε² = 2 • D Si • t ( 1
)
Where t is the time and D Si is the apparent diffusion coefficient of Si through the SiC layer. Applying 1, the calculation of the diffusion coefficient of Si into SiC at 1100°C under CH 4 is D Si = 1.15×10 -17 cm².sec -1 .This value is in accordance with previous studies conducted in similar experimental conditions [START_REF] Ferro | Infrared kinetic study of ultrathin SiC buffer layers grown on Si( 100) by reactive chemical vapour deposition[END_REF][START_REF] Cimalla | The diffusion coefficient of silicon in thin SiC layers as a criterion for the quality of the grown layers[END_REF]. Interestingly, it has to be noticed that with an important gas flow of H 2 during the rise of temperature, smoothing of sidewalls is observed, accompanied by a slight diameter reduction (~ 11 nm). A closer inspection shows that these desoxidized Si NWs display facetted sidewalls with {220} and {400} planes with an irregular octagonal section (see Fig. 3a). This is probably due to the minimization of surface energies.

After carburizing these smooth NWs, no morphological change is visible (Fig. 3a). In order to check if any SiC material was formed and to study the structure of carburized Si NWs after carburization, a thin NWs cross-section lamella was prepared thanks to the FIB-SEM. The result is shown in Fig. 3b. As mentioned earlier, the octagonal shape of the section is clearly visible. Dotted lines represent the {220} planes and the continuous lines represent the {400} planes of Si. Interestingly, the {220} planes are smaller than {400} planes. A high magnification image of the periphery is shown Fig. 3c. One can see that the carburization step leads to a single crystalline, ~ 2 nm thick 3C-SiC layer which covers perfectly the Si NWs. The diffraction pattern inserted on Fig. 3c points out the cubic polytype of the buffer layer in epitaxial relationship with the Si NW. These carburized smooth nanostructures have been then used as template for the epitaxial growth of 3C-SiC. As it can be seen on Fig. 4a, this step leads to rough sidewalls of NWs but it also enhances the faceting of the NWs leading apparently to a square section. The occurrence of this faceting suggests that the deposit is crystalline and in epitaxial relationship with the NWs. Moreover, this deposit covers completely the NWs, from the top to the bottom. A FIB-cut of the head of a NW after SiC epitaxy confirms the square shape after growth (Fig. 4b). On this image is reported the {220} planes and the {400} planes, according to those identified after the carburization step done on smoothed NWs thanks to the difference of length between {400} and {220} planes. The anisotropic growth of 3C-SiC on carburized Si NWs is therefore favored towards {400} directions at the expense of the {220} ones. The thickness of the 3C-SiC layer is around 60 nm on the {400} planes and 10 nm on the {220} planes.

Conclusion

Core Si / shell SiC NWs have been successfully elaborated from etched Si NWs. Firstly etched Si NWs are carburized under either CH 4 or C 3 H 8 , both diluted in H 2 , at atmospheric pressure. Monitoring the H 2 flow allows us to smooth the surface of Si NWs before the carburization step but also to highlight the {400} and {220} planes because of reduction of surface energies. The as-grown 3C-SiC layer obtained by carburization is single crystalline, ~3 nm thick and covers entirely the surface. This layer is then used as seed for the epitaxial growth of 3C-SiC. This growth occurs preferentially on the {400} planes. Thanks to the biocompatibility of the SiC shell associated with a Si core in which electronic transport is excellent, this structure has high potential for bio-sensors.

Fig. 1 :

 1 Fig. 1: SEM image of Si NWs obtained by a SF 6 /0 2 plasma etching of Si(100).

Fig. 2 :

 2 Fig. 2: a) SEM image of Si NWs carburized at 1100°C under CH 4 at atmospheric pressure. b) TEM image of the cross-section of a carburized Si NW. The SiC layer is ~ 3 nm thick, as shown in insert. A thicker layer of about 20 nm appears on the image which is only due to the original roughness of the Si NWs. c) 3C-SiC thickness as a function of the squared root of the carburization time.

Fig. 3 :

 3 Fig. 3: a) SEM image of Si NWs carburized under C 3 H 8 with an important flow of H 2 . Sidewalls are smoothed and NWs exhibit an octagonal irregular cross section. b) TEM crosssection of a carburized Si NW with the octagonal shape. {220} planes (dotted lines) and {400} planes (continuous lines) are highlighted. c) Higher magnification on the edge of the NW showing the 3C-SiC carburization layer, diffraction pattern in insert.

Fig. 4a )

 4a Fig. 4a) SEM image after the epitaxial growth of 3C-SiC on the smoothed and carburized Si NWs. Sidewalls are rough and the epitaxial layer is homogeneous down to the NWs feet. b) SEM image of the square section of one NW after a FIB-cut. The 3C-SiC epitaxial layer has grown preferentially on the <400> directions (continuous lines) at the expense of <220> ones (dotted lines).


M Zhe

O Dezellus

G Parry

M Braccini
email: muriel.braccini@simap.grenoble-inp.fr

J C Viala

Modified 4-Point Bending Test for Adhesion Measurement at the Interface of Iron Coated with Aluminum Casting Alloy

Keywords: Casting, fracture energy, aluminum alloys, steel, interface

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The development of multi-material systems enhances the research and understanding of the physics, chemistry and mechanics of interfaces. In particular, it is often necessary to evaluate the practical adhesion at the interfaces. Indeed, a good adhesion is a necessary condition for the interface to fulfill its function, even if it is a physical one (electrical conductivity, thermal conductivity, etc.). Various methods have been developed to mechanically characterize an interface: they either measure the strength of the interface which is the maximum load (or stress) that causes a fracture; or the interface toughness, which is the mechanical energy per unit area at which the interface separation occurs. The interface strength, measured by a tensile pull test [1], for example, is strongly dependant on crack initiation conditions. The interface toughness is measured by a steady state propagation of a crack in a fracture mechanics specimen [2][3][4] and provides a direct measurement of the interface adhesion energy. These kinds of methods are less sensitive to the specimen details and leads to fewer discrepancies.

The 4-point bending test is increasingly used for quantitative characterization of interface toughness [5][6][7][8] because it can be applied either to the interface between bulk materials (thick coating for example) or to thin films using a modified test where a counter-part is fixed onto the film to increase its stiffness [7]. Moreover, in the case of pure elastic deformations, the analysis is quite simple: during the crack propagation, the load is constant and depends directly on the interface adhesion energy [4]. Nevertheless, the sample preparation can sometimes be difficult, particularly when a counter-part has to be fixed because the film/counter-part interface must be stronger than the interface of interest. Moreover, pure elastic deformation hypothesis is not always fulfilled. Particularly in studies on interface between metals there is always some plastic deformation, and depending on the level of the adhesion energy and on the ductility of the metals, the plasticity will be limited to the crack tip process zone or will extend into the material. When the plastic deformation is limited to the crack tip area, the Griffith approach to brittle fracture can still be applied and the analytical relation between the toughness and the constant load is valid. In the case of an extended plasticity, the energy dissipated in the plastic deformation cannot be neglected [9]. At this point numerical methods are often used: they are based on a finite elements method in which the plastic behavior of the metals can be introduced.

One way to study interface crack propagation in finite elements methods is to introduce cohesive elements at the place of the interface. The cohesive elements are used to model cracks in homogeneous or heterogeneous materials. In the case of interface separation, this kind of approach is easier because the crack path is known and is located at the interface. As a result cohesive elements have been used in adhesive joints [10], brazing joints [START_REF] Braccini | Proceedings of the International Conference on Fracture XII[END_REF] or thin film delaminations [START_REF] Shirani | [END_REF]. Modeling the 4-point bending test, the cohesive elements parameters can be evaluated using an inverse method: results from the modeling are compared to experimental data. In this work this method, combining 4-point bending with cohesive zone modeling, is applied to measure the interface adhesion of a Fe/Al-Si interface. This system is of practical interest because Al-Si cast alloys are intensively used for automotive applications and their local reinforcement by ferrous inserts is a clas-sical method to enhance their ability to sustain high load levels [START_REF] Clyne | An Introduction to Metal Matrix Composites[END_REF][START_REF] Durrant | [END_REF]. In this case, a mechanically strong interface, ensured by a true metallurgical bond, is a key point for load transfer. Such interfaces can be obtained by means of the well-known Al-fin process which consists in two steps: aluminizing of the ferrous insert immediately followed by its insertion in a casting mould [15,16]. This process leads to the formation of a complex interfacial reaction zone with different intermetallic phases of the ternary Al-Fe-Si system. The advanced knowledge already acquired by some of the present authors on the thermodynamic and kinetic aspects of interfacial reactivity in this ternary system [17][18][START_REF] Pierre | Diffusion in Materials[END_REF] not only allows us to produce 4-point bending specimens with a rigorously controlled interfacial chemistry but also to further understand the relationship between the interfacial chemistry and mechanical properties of bimetallic systems.

Experimental

Materials and Procedure

Bimetallic samples were processed from an Al-7 wt% Si-0.3 wt% Mg casting alloy suppliedbyAlcan ® and referred to as A356.2 according to ASTM standards. The composition of this alloy which will be referred to as Al-Si alloy in the present paper is given in Table 1.

Iron of commercial purity (99.5%) was provided by Goodfellow ® as 1 mm thick sheets (typical composition in wt ppm: Mn 3000, Si 1000, C < 800, P < 400, S < 500, Fe balance). These sheets were cut by Electrical Discharge Machining (EDM) into plates of the desired dimension for 4-point bending test specimens: 10 mm × 60 mm. All the plates received the same surface preparation of mechanical abrasion using SiC paper under water to a mean surface roughness of about 3 µm, followed by rinsing and drying. The plates were kept in dichloromethane to prevent air contamination. A key feature to measure the toughness of an interface by 4-point bending test is to ensure that the crack propagates at the interface. In order to facilitate this interfacial propagation, some special treatment must be applied to produce a weak zone, without any metallurgical bond, at the center of the specimens. In this work this was done by depositing, just before immersion in the Al-Si alloy-melt, a thin 2 mm wide layer of Stopyt ® 62A™ from Wesgo ® Metals (Yttrium oxide suspension in an organic liquid -a procedure often used in brazing applications to hinder the spreading of liquid metals). The Al-Si alloy was cut into small pieces that were degreased in an ultrasonic bath of dichloromethane, placed in a ceramic oxide crucible and heated above their melting temperature by direct radio-frequency coupling as shown in Fig. 1. The establishment of a strong metallurgical bond across the liquid solid interface can be obtained only after the removal of the native oxide film layer on the Fe [20]. The removal was promoted by the dissolution of Fe in the liquid alloy and is therefore favored in the case of melts with a low Fe content [18]. However, in this case, the growth of reaction layers was limited by their dissolution in the bath until the solubility limit of Fe was reached (6-8 wt% according to [18,21]). In the present study we have chosen to add 2 wt% Fe to the melt in order to obtain a reasonable compromise between the dissolution of Fe and the growth of the reaction layers.

Before immersion of the Fe plates the bath was regularly stirred and the surface oxide skin removed. The bath temperature was measured with a precision of better than ±0.2 • C by plunging in the melt a K-type (Ni/Cr) thermocouple. Note, that with a bath heated at 805 • C, immersion of a Fe plate led to a temperature decrease of about 10 • C. Therefore, interface reactions were considered to develop at a temperature of 795 ± 5 • C. After a 2 min of immersion, the Fe plates were gently removed from the bath and solidification of the aluminizing layer occurred in air in a few seconds.

Mechanical Testing

The adhesion energy of the bimetallic Al-Si alloy/Fe specimens, G c , was measured using a 4-point bending test [4]. The test was adapted to the system of interest, which consists of a thin and brittle intermetallic layer that was enclosed between a thin aluminizing layer and the Fe sheet [7]. A metallic stiffener (Al 2024-T3) was then bonded with Araldite ® 2015 on the top of the aluminized Fe plate.

During the 4-point bending test, a symmetric crack propagates along the weakest interface of the multilayer specimen. To ensure the localization of the crack at the Al-Si alloy/Fe interface, optimization of the strength of other interfaces, Al-Si alloy/Araldite ® 2015 and Araldite ® 2015/Al 2024, was necessary. This was done by applying a surface treatment prior to joining with the adhesive: this treatment consists of different steps including an anodization treatment step [22,[START_REF] Wegmen | Surface Preparation Techniques for Adhesive Bonding[END_REF]. Next, the aluminized Fe plate was glued to the stiffener at room temperature for 12 h in a Teflon mould in order to ensure the alignment and constant thickness of the Araldite ® layer along the 4-point bending specimen. Finally, a pre-notch was made on the stiffener surface to the Fe plate using a diamond wire saw. As shown in Fig. 2, this pre-notch was located in the middle of the specimen, where a Stopyt ® 62A™ layer had been deposited to facilitate crack initiation in the Fe/Al-Si alloy interfacial zone.

Mechanical testing was performed in two steps: (i) initiation of a crack by 3-point bending test at a speed of 1.2 mm/min up to 1 mm displacement, (ii) followed by crack propagation by means of 4-point bending test geometry at a typical speed of 1.2 mm/min. During mechanical testing two sets of data were recorded: the classical evolution of load vs. displacement but also side view pictures of the sample were taken in order to follow the crack propagation (periodicity of 6 s).

Characterization

After solidification, some bimetallic specimens were cut with a diamond wire saw and diamond polished to a finish better than 1 µm for examination by optical microscopy (OM), scanning electron microscopy (SEM) and electron probe microanalysis (EPMA). After failure in the 4-point bending test, the fracture surfaces were characterized to determine the precise path of the crack. Both sides of the fracture surface (Fe side and Al 2024 side) were observed by OM, SEM and analyzed by X-ray diffraction and EPMA. Observations were also made after failure on sections of the Fe part cut perpendicular to the fracture surface. In that case, before mounting the sample in resin, an electroless hard Ni coating (Nicklad™ ELV 804 from MacDermid) was applied on the fracture surface in order to protect it from detrimental polishing effects such as brittle matter loss or edge rounding off.

The XRD spectra were recorded using standard diffraction equipment (Panalytical MPD-Pro diffractometer equipped with a back monochromator and a X'celerator detector, Cu K α radiation). The OM and SEM observations were made on diamond polished sections. The SEM observations and EPMA analyses were carried out using a Camebax apparatus (Cameca) equipped with an energy dispersive analyzer. The accelerating voltage was 10 kV and the beam current was 7.8 nA. After background subtraction, the counting rates obtained for Al, Fe and Si in at least eight different points were averaged and referred to the counting rates recorded under the same conditions on pure and freshly polished element standards. After corrections for atomic number, absorption and fluorescence, the atomic contents of Al, Fe and Si in the different phases deriving from the Al-Fe-Si ternary system were obtained with accuracy of better than ±0.5 at%.

Results and Discussion

Interface Chemistry

Examination by SEM of the bimetallic sample end sections indicates that the reaction zone formed at the Fe/Al-Si alloy interface, during the immersion process, consists of different intermetallic phases (see Fig. 3). More especially, the following interface reaction layer sequence was found by EPMA from Fe to Al-Si alloy: is θ -Al 13 Fe 4 (Si) [17,21]. Therefore, one can conclude that the first two layers of the interfacial sequence (ηAl 5 Fe 2 (Si) + τ 1-9 -Al 3 Fe 3 Si 2 /θ -Al 13 Fe 4 (Si)) formed during immersion at 795 ± 5 • C and that this sequence corresponds to a diffusion path in the isothermal section of the ternary system at 800 • C [18,[START_REF] Pierre | Diffusion in Materials[END_REF]. As to the outer layer of τ 5 -Al 7.4 Fe 2 Si, it appears to have formed upon cooling as a result of Fe saturation of the liquid alloy film covering the Fe plate just removed from the bath.

4-Point Bending Test

Figure 4 reports a typical load-displacement curve obtained by 4-point bending testing. Note that reproducibility has been checked by testing 7 different specimens.

Initially the load increases linearly as the sample deformed elastically. Next, a plateau was observed at a constant load level. By combining a load-displacement curve and images recording during the test, it can be concluded that during this plateau the crack propagates along the interface until it reaches the inner loading points. Finally, plasticity of the free Fe plate intervenes leading to a slow increase of the load and an important bending.

Because crack propagation occurs during elastic loading, the fracture energy of the interface G c is equal to the strain energy release rate G, which has been establishedby [4]:

G = 3 2 F 2 l 2 b 2 h 3 (1 -ν 2 2 ) E 2 × 1 (h 2 /h) 3 - λ (h 1 /h) 3 + λ(h 2 /h) 3 + 3λ((h 1 h 2 )/h 2 )(h 1 /h+ λ(h 2 / h)) -1 , (1) where λ = E 2 (1-ν 2 2 ) E 1 (1-ν 1 2 ) .
Subscript 1 denotes the aluminum stiffener and 2 the iron plate. Geometric parameters of these equations are the width of the specimen b, the heights of each material composing the beam h 1 and h 2 and the total height of the specimen h = h 1 + h 2 ; l is the spacing between inner and outer load lines: l = 12 mm in this study. Mechanical parameters are the elastic properties of the materials, the Young's modulus E and the Poisson's ratio ν, and the total load F . Numerical values of those parameters are given in Table 2.

Since cracks propagate for loads between 47 N and 53 N, equation ( 1) gives an adhesion energy of 23 ± 3J/m 2 .

Crack Path

Results in terms of phase composition on both sides of the fracture path are illustrated by the XRD patterns presented in Fig. 5. On the Al side the diffraction lines of η-Al 5 Fe 2 ,θ-Al 13 Fe 4 and Al are present. Note that other phases (τ 5 -Al 7.4 Fe 2 Si and τ 1-9 -Al 3 Fe 3 Si 2 ) cannot be clearly identified on this diffraction pattern because of the low intensity of the residual peaks.

On the Fe side only η-Al 5 Fe 2 and Fe are characterized: η-Al 5 Fe 2 is thus the only phase present on both sides of the crack path.

Observations on a polished section perpendicular to the Fe side complete those analyses (Fig. 6). It clearly appears in Fig. 6 that the crack has propagated in the interface reaction zone, at a distance from the η-Al 5 Fe 2 /Fe interface that is typically lower than 1 µm. EPMA analyses performed at several points of the intermetallic layer remaining on Fe substrate also confirm the presence of η-Al 5 Fe 2 as the major phase with some Si in excess attributable to τ 1-9 -Al 3 Fe 3 Si 2 inclusions.

Combining XRD results with metallographic examination and EPMA characterization leads to the conclusion that the crack propagates in the η-Al 5 Fe 2 intermetallic reaction layer, most often at a distance less than 1 µm from the Fe underlying substrate.

Description of the Inverse Method. Application to the 4-Point Bending Test

This section describes the numerical method that has been used in order to obtain a precise value for the toughness of the Fe/Al-Si alloy intermetallic interface, prepared as described above. Indeed the analytical approach leading to equation ( 1) is a simple energy balance based on elasticity calculations on the overall system. Moreover this approach omits the presence of various thin films between the two plates: reaction layers and adhesive bonds that are difficult to account for in the analytical calculation, but are readily included in the numerical calculation. Hence it is only valid for a stationary state propagation of the crack, no crack initiation or transitory crack growth can be described in this approach.

The numerical method developed here is based on cohesive zone modeling. The basic idea was to carry out a numerical simulation of the 4-point bending test, using the finite elements method involving cohesive zone elements. In these simulations, the adhesion energy is described as a local property of the interface even if it is in practice -the work of adhesion. It is possible to fit this parameter, as well as the value of the critical stress leading to the failure of the interface, by using the data from the 4-point bending experiment. If the response of the system can be numerically reproduced with enough precision on a point by point basis, then a good estimate of the interface properties can be expected.

Choice of the Cohesive Zone Model

The numerical calculations have been carried out using the finite element software ABAQUS © , in which the cohesive zone model is implemented through the use of specific elements called cohesive elements.

Cohesive zone models regard a fracture as a gradual phenomenon in which separation takes place across an extended zone in the neighborhood of the crack tip, or 'cohesive zone', with cohesive tractions resisting the separation. The concept of the cohesive zone used in modeling fractures, was introduced by Barenblatt [START_REF] Barenblatt | [END_REF]. Cohesive zone models became popular in the field of material sciences after the introduction of the Xu and Needleman model [25] based on an exponential cohesive law that could be conveniently introduced in finite element codes.

Cohesive zone elements implemented in finite element codes do not represent any physical material, but describe the cohesive forces which occur when material elements (such as grains) are being pulled apart: they are placed between continuum (bulk) elements, as shown in Fig. 7. The traction T = T t t + T n n is the interface force per unit area of the interface, while the separation δ = u +u -= δ t t + δ n n is defined as the difference between the displacements u + and u -, respectively on the upper and the lower faces. As the calculations were carried out in the framework of plane strain, only tractions in the normal and the shear direction were considered.

When damage growth occurs these cohesive zone elements open in order to simulate crack initiation or crack growth. The description of the failure behavior is defined by traction-separation laws. There is a great variety in traction-separation laws [26] which exhibit the same global behavior; we chose a simple bilinear law as depicted in Fig. 8 and detailed in the Appendix. As the cohesive surfaces separate, the traction first increases until a maximum is reached, and subsequently decreases to zero. The energy dissipated at failure at a material point of the interface, G c ,is given as a parameter.

Application to the 4-Point Bending Test

The mechanical model used to simulate the 4-point bending test is depicted in Fig. 9. A state of plane strain is assumed (i.e., the geometry is supposed to be invariant in direction 3). The different materials composing this specimen are supposed to have an elastic behavior, except for the lower part of Fe where the damaging process is taking place: the intermetallic layers were modeled with cohesive zone The different layers that were taken into account are, from the top to the bottom of the multilayer: the elastic Fe layer, the cohesive zone at the bottom of Fe layer, the glue and the Al 2024 stiffener. The summary of the various layers properties can be found in Table 2. The properties chosen for the cohesive zone are discussed in the Appendix.

The aim of the numerical calculation was: (a) to check consistency with the simple analytical approach and determine the maximal traction parameters T 0 n and T 0 t by an inverse method for a given value of G c .

(b) To compute precisely the stresses components in the iron plate during loading in order to check the validity of the linear elasticity assumption.

The set of data generated by the numerical calculation was fitted to the corresponding set of data from the experiment in order to identify some of the parameters of the cohesive zone. The set of data was the reaction force as a function of the applied displacement of the upper loading point (F(ū)). The experimental data show that the curve F(ū) was at first linear up to a given value of the applied displacement ūc , for which the crack was beginning to propagate. Then the curve reached a short plateau F = F p . It is worth noting that in the present case F p >F(ū c ), which means that the force transmitted by the loading point has to increase in order to propagate the crack. The experimental value of the plateau force F p provides a criterion for fitting the cohesive zone parameters using the numerical results.

In order to simplify the approach, the same value was assumed for the maximal normal and shear traction components: T 0 n = T 0 t = T 0 .Thev alueG c was set to 23 J/m 2 , resulting from the analytical calculation carried out in Section 3. Numerical calculations have been systematically carried out with increasing values of T 0 , and the results have been compared to the experimental data. The results for T 0 ={30, 40, 50, 60} MPa are reported in Fig. 10, together with the experimental data. The load F is plotted versus the applied displacement ū. Figure 10(a) shows a global view of the loading curve, whereas Fig. 10(b) focuses on the plateau. The numerical calculations produce a loading curve which is in good agreement with the experimental curve from the initial elastic regime up to the plateau. The critical displacement value ūc , for which the response of the multilayer system becomes nonlinear, does not seem to be affected by the value of T 0 . The best agreement is found for a maximal traction T 0 between 40 and 50 MPa. It is worth noting that this value for the maximal stress T 0 is quite close to some estimated in previous experimental studies on similar systems but with different testing geometries [20,27]. However, this experimental value was indeed a value for the normal traction T 0 n .An experimental value for the shear traction would also be necessary in order to check the validity of our assumption T 0 n = T 0 t = T 0 . It can be concluded that the results from the finite element calculations using a cohesive zone model are quite consistent with the analytical results. The inverse method based on the comparison of the numerical and experimental output allows for the determination of a value for the maximal traction T 0 . The shapes of the numerical curve and the experimental curve are quite similar, up to the end of the plateau. However, after the plateau the two evolutions are different. This discrepancy originates in large plastic deformations occurring in the real sample in the final stage of the testing, which are not taken into account in the numerical model as it is discussed in next section.

Discussion

One central assumption of both the analytical and the numerical model is that the deformations occurring in the structure are essentially elastic. It is possible to check the validity of these assumptions by analyzing the values of the stress tensor components that have been computed numerically during the loading. The evolution of the Mises stress invariant during loading has been reported in Fig. 11 for two points of the iron plate: one near the top surface and one near the bottom surface (above the cohesive zone). The experimental value of the yield stress obtained by a tensile test is also reported. These results suggest that there is no plastic deformation in the iron during the plateau, which justifies a posteriori the hypothesis of elastic deformations. The calculations also indicates that the yield stress is reached in iron for applied displacements above ū = 0.45, starting at the lower and upper surfaces of the iron layer and progressively extending to the center of the layer. This results explains why the behavior of the numerical model is so different from the one observed experimentally for larger applied displacements.

Finally, it is worth noting that in the configurations where some plasticity appears for the lower loading, no plateau was found during the experiment and the numerical approach. This was the only one to retrieve the value of G c from experimental data. Accounting for plastic deformations in the calculation is not difficult if the plastic behavior of the material is determined by a tensile test for example (yield stress, hardening curve).

Conclusion

A modified 4-point bending test has been used to measure the fracture resistance of an Fe/Al-Si alloy interface. This interface was obtained by the immersion of Fe plates in an Al-7 wt% Si melt leading to the formation of mainly three phases in the interfacial reaction zone: a 3 µm thick η-Al 5 Fe 2 (Si) layer in contact with Fe; a 2 µm thick θ -Al 13 Fe 4 (Si) layer and a 6 µm thick τ 5 -Al 7.4 Fe 2 Si layer in contact with the Al-Si alloy film. Accurate characterization of the crack path by different techniques has clearly shown that crack propagates only in the ηAl 5 Fe 2 (Si) layer. Analytical calculation based on beam bending theory leads to a value for the fracture energy, associated with this crack growth, equal to 23 J/m 2 . Thus, it appears that the ηAl 5 Fe 2 (Si) phase is the weak layer in the interface reaction zone and should be maintained as thin as possible by controlling the solidification rate as well as during post thermal treatment.

To complete the analytical approach for 4-point bending data analysis a FEM numerical method was developed using cohesive zone modeling. An inverse method based on the comparison of numerical results to the experimental loaddisplacement curves allows one to calculate the interfacial properties. The cohesive zone model presented in this paper could be used for the evaluation of the interface energy in the case where the interface does not show a perfectly brittle behavior, as it is assumed in the analytical approach. One perspective would be to extend it to samples with an extended plastic deformation, i.e., when the metal surrounding the interface undergoes large plastic deformations.

Appendix: Boundaries Conditions and Some Details of the Numerical Model

• In the numerical approach, only half of the system was modeled, for symmetry reasons. The symmetry plane cuts the picture along the AC direction (Fig. 9). The portion BC must be set free of motion in the two directions, since a cut was made along this segment in the experimental conditions. Symmetry conditions were prescribed along AB (u 1 = 0a n dθ 3 = 0, where u 1 is the displacement in direction 1 and θ 3 is the rotation around direction 3). A cohesive zone thickness of 15 µm was chosen in accordance with experimental observations. • The two alumina rods where loading was applied were modeled as rigid cylinders. Contact conditions without friction were set up between the alumina rods and the upper and lower face of the multilayer. The lower loading point was fixed, whereas the upper one could move vertically. Hence, the experimental loading conditions were accurately reproduced. The vertical displacement of the alumina rod was prescribed during loading, and the reaction force on it could be evaluated at each step. • Large deformations were taken into account in the simulation. The elements chosen for the glue and the Al 2024 layers were standard four nodes plane strain elements. For the Fe layer, four nodes incompatible modes plane strain elements were chosen in order to give a good representation of bending. The analysis was quasi-static. • The traction-separation law describing the behavior of the cohesive elements at the interface between Fe and Al 2024 is depicted in Fig. 8 where δ m is an effective displacement define in ABAQUS © as δ m = δ 2 n + δ 2 t . It is a bilinear law. In this kind of model, the cohesive zone has a stiffness described by a stiffness matrix K, corresponding to the slope of the first linear (increasing) part of the traction/separation law before the onset of damage:

T = T n T t = K nn K nt K nt K tt δ n δ t .
This stiffness is more to be seen as a numerical artifice than a real property of the interface. A quasi-rigid representation has been chosen here. The total separation work is δ f m T 0 2 , δ f m being the final separation when failure of the element is completed. The stiffness components K nn and K tt are chosen such that the pseudo elastic energy (area δ 0 m T 0 2 under the increasing part of the curve, where δ 0 m is the separation corresponding to the maximum traction T 0 ) is limited to 1% of the total separation work. In this condition, the total separation work can be assimilated into G c without ambiguity.

This leads to the following values of the stiffness components:

K nn = 100 (T 0 n ) 2 2G c ,K tt = 100 (T 0 t ) 2 2G c ,K nt = 0,
where T 0 n and T 0 t are the maximum normal (respectively shear) traction that the interface can transmit without being damaged.

The onset of damage for the interface is assumed to happen when the following equality is fulfilled:

T n T 0 n 2 + T t T 0 t 2 = 1.
From this point, damage growth in the element leads to a decrease of the stress components T n and T t . We chose a linear decrease according to an effective displacement define in ABAQUS © as δ m = δ 2 n + δ 2 t . • During the specimen preparation, a weak zone was realized by the deposition of a Stopyt ® 62A™ layer at the center of the specimen (segment BD in Fig. 9) that opened during the preliminary 3-point bending test. The resulting crack was estimated to be 4 mm on each side. This has been taken into account in the simulations, by changing the material characteristics of the cohesive zone elements lying along this preliminary crack region of 4 mm. Very low T 0 n , T 0 t and G c values have been chosen for this portion of the interface. The elements along this weak line were completely damaged in the very first steps of the loading, thus initiating a 4 mm crack at the beginning of the four points bending test.

• In order to give numerical stability to the calculations, a very small viscosity coefficient has been added to the cohesive zone elements. This is a classical method that is readily implemented in ABAQUS © . Caution was taken to check that the influence of this parameter on the results was very small.

Figure 1 .

 1 Figure 1. Manufacturing process of bimetallic sample by immersion of a Fe plate in a molten bath of Al-Si alloy.

Figure 2 .

 2 Figure 2. Schematic of 4-point bending test specimen with pre-notch and crack initiation.

Figure 3 .

 3 Figure 3. (a) Typical morphology of the interface between Fe and Al-Si alloy after 2 min of contact at 800 • C. (b) EPMA concentration profiles of the Fe/Al-Si transition zone shown in (a).

Figure 4 .

 4 Figure 4. Typical load-displacement curve obtained during 4-point bending test of Al-Si alloy/Fe bimetallic specimen. Crack propagation starts at point A and ends in point B when it reaches the inner loading points (see Fig. 2).

Figure 5 .

 5 Figure 5. XRD pattern recorded of the fracture surfaces on both sides after 4-point bending test.

Figure 6 .

 6 Figure 6. Transverse section of the fracture surface at Fe side after 4-point bending test and coating by Ni electroless plating (SEM). Phase identification has been performed by EPMA.

Figure 7 .

 7 Figure 7. Example of a structure made of a bimaterial submitted to a remote traction. The nominal traction stress vector T , representing the interface loading. The corresponding separation is denoted δ.

Figure 8 .

 8 Figure 8. Constitutive traction separation law for the cohesive zone. The damage initiation is occurring for a given value δ f m of the equivalent separation. δ f m is the value of the final separation. G c as defined in the code is the energy dissipated between δ 0 m and δ f m .

Figure 9 .

 9 Figure 9. Mechanical model for the 4-point bending test of the Fe/Al 2024 specimen. Half of the model is depicted, due to symmetry. Materials and dimensions are described as well as boundary conditions.

Figure 10 .

 10 Figure 10. Reaction force versus displacement of the upper indenter, for G c = 23 J/m 2 and various values of T 0 . The evolution of the loading obtained experimentally is reported as well.

Figure 11 .

 11 Figure 11. Mises equivalent stress for two points of the iron plate (one near the upper surface and one near the bottom surface) and yield stress experimentally determined by tensile test, evidencing that the deformations remain elastic up to the plateau.

Table 1 .

 1 Composition

		(wt%) of the AlSi7Mg0.3 alloy provided by Alcan ®
	Element Si	Fe	Cu	Mn	Mg	Ni	Zn	Pb + Sn Sb	Ti	Others
	Min.	6.7				0.3
	Max.	7.3 0.14 0.02 0.04 0.4	0.02 0.04 0.02	0.16 0.15 0.1
	±	0.2 0.01 0.01 0.01 0.02 0.01 0.01 0.01	0.02 0.01

Table 2 .

 2 Numerical values of geometric and mechanical parameters for the 4-point bending tests The Young's modulus, the Poisson's ratio and the thickness of the different materials are denoted by E, ν and h, respectively, with an appropriate subscript.

	Layers	Thickness h (mm)	Young's modulus E (GPa)	Poisson's ratio ν
	Iron (subscript 2)	1	200	0.29
	Al 2024 (subscript 1)	2	73	0.33
	Glue	0.4	2	0.4
	elements.

Acknowledgements

The authors like to thanks the 'cluster MACODEV de la région Rhône-Alpes' for funding part of this work.


Thierry Duffar
email: thierry.duffar@grenoble-inp.fr

La ligne triple grain-grain-liquide lors de la solidification du silicium multi-cristallin

The grain-grain-liquid triple phase line during solidification of multi-crystalline silicon

Introduction

The grain structure of multi-crystalline Si ingots impacts the conversion efficiency of photovoltaic cells, because grain boundaries are location of lattice defects that interact with minority carriers [START_REF] Ciszek | Proc. 23rd PV Specialist Conf[END_REF][START_REF] Nelson | The Physics of Solar Cells[END_REF]. Furthermore, grain boundaries act as getters that attract impurities, which has the drawback of increasing the grain boundary electrical effect but the advantage of cleaning the surrounding grains. These effects have gained importance as industry uses less pure silicon raw material [START_REF] Dubois | [END_REF].

It is therefore important to know how the grain structure develops in Si multi-crystalline ingots. This topic has been studied for a long time (see a review in [START_REF] Duffar | [END_REF]) and the multi-silicon structure is typically composed of three types of grains: [19] Grain boundary twisted around 110 or 100 axis From 0.6 to 1.3 J m -2 (except for Σ3 twin) See plots in [21,22] a Ref. [23] r e p o r t st h a ta n i s o t r o p yi sl e s st h a n1 5%i nt h eo r d e r{111} < {110} < {100}. This is confirmed by numerical simulation [24].

-Long columnar grains with a diameter of the centimeter order, which develop along the growth direction. This is the preferred structure because grains are large and do not impact very much on the properties. -Small equiaxed grains of millimeter size. They should be avoided as they decrease the photovoltaic efficiency and are associated with SiC precipitates that induce severe problems during subsequent wafering. They have been shown to be related to the nucleation of Si on SiC precipitates [5], when C pollution occurs, but this needs a significant undercooling of the solid-liquid interface during growth, at least 0.1 K [6].

-Σ3 twins that are not electrically active by themselves but, when crossing or multiplying, produce Σ9, Σ27, etc. boundaries. It has been shown that four levels of successive twinning result in a random structure of incoherent grains [7].

The same authors have measured that two cell wafers cut from a directionally solidified ingot are composed of only 50 to 60% of random grains, the remaining grain boundaries being of the type Σ3 n . This result has been confirmed in the case of a Si ingot grown by cold crucible continuous casting, under much higher growth rate and thermal gradient [8], demonstrating the importance of twinning in photovoltaic Si solidification. A mechanism of twinning in multi-crystalline Si has been proposed recently [9], based on the existence of undercooled facets at the interface.

It follows that the grain structure in multi-crystalline silicon largely depends on a series of phenomena related to faceting/undercooling of the solid-liquid interface during growth. This paper discusses the solid-liquid interface structure based on the morphology of the grain-grain-liquid triple lines. Very classical kinetic and energetic considerations are used. The most important is that, according to the Jackson criterion [START_REF] Jackson | Liquid Metals and Solidification[END_REF] and as demonstrated in many occasions, facets in Si are only of the {111} family. Ref. [START_REF] Rudolph | [END_REF] reports possible {110} facets, but this is still to be confirmed.

Growth kinetics in Si are well known:

-In the case of a rough interface, the relationship between the interface growth velocity, V I , and the undercooling, T , is linear; for the 100 direction, [12,13] give:

V I{100} ≈ 0.12 T in m s -1 and K

which gives an undercooling of 10 -4 K for a typical solidification rate of 10 µm s -1 .

-In the case of a faceted {111} interface where two-dimensional nucleation occurs, [14] has reviewed the existing litera- ture and proposes the relation:

V F = 0.04956e -46.67 T (0.3 T ) 2/3 in m s -1 and K ( 2
)
which gives an undercooling of 5 K for a typical solidification rate of 10 µm s -1 .

-In the case of a faceted {111} interface with unsaturable steps (for example due to dislocations), [15] proposed the relation:

V F ≈ 0.3 × 10 -3 T 2 in m s -1 and K (3)
which gives an undercooling of 0.2 K for a typical solidification rate of 10 µm s -1 .

From the energetic point of view, Table 1 gives a selection of interfacial energies in Si, relevant to our considerations. It can be seen that, except for the Σ3 twin, the grain boundaries have similar energies, whatever their relative mis-orientation.

Bristowe [START_REF] Bristowe | Properties of Crystalline Si[END_REF], in a review of molecular dynamics and quantum mechanics simulations, shows that an amorphous layer, 0.25 nm thick, exists between twist grains. The energy of tilt grains is lower and it is not clear whether pre-melting exists at grain boundaries close to the melting point in this case. Following [START_REF] Priester | Les joints de grains, de la théorie à l'ingénierie[END_REF], it is likely that only the Σ3 boundary remains crystalline close to the melting point. In the following, it will be considered that, for all other boundaries, the phenomena occurring on one side of the groove are not influenced by what happens on the other side.

In the second part, the morphology of the boundary surrounding a single grain will be discussed. This will allow, in the third part, a simple computation of the faceted or rough structure of the solid-liquid interface as a function of the growth parameters. Several types of Grain-Grain-Liquid Triple Phase Lines (GGL-TPL) will be presented in the fourth part and it will be shown that each of them induces a particular orientation of the grain boundary toward the solid-liquid interface.

Conditions for grain boundary faceting

The morphology of the GGL-TPL will depend on whether the interface is rough or faceted at the TPL. In order for the boundary of one grain to be faceted, there is first a crystallographic criterion: as facets are necessarily {111} planes, they can only exist on some parts of the grain boundary around one grain. Fig. 1 shows the case of a grain which is growing in the 100 direction, for the sake of simplicity. It can be seen that the four possible {111} facets can be found in the [011] directions only. In this case the facets make an angle of 54.7 • with the (100) plane, i.e. with the solid-liquid isotherm, T fus . More generally, the locations of facets around a grain depend on the orientation of the grain relatively to the normal to the solid-liquid isotherm: a 100 -orientated grain will show 4 facets, a 110 -orientated grain will show two facets (at 35.3 • from the solidification isotherm), a 111 -orientated grain will show one facet parallel to the solidification isotherm and three facets at the grain boundary (at 70.5 • from the solidification isotherm), and so on.

There is also an energetic criterion: following [START_REF] Voronkov | [END_REF] a facet is possible only if the energy of the system decreases when a growing step reaches the facet boundary. This gives a condition on the angle, ν, between the facet and the grain boundary (see Fig. 1-b). As shown in [9], this angle should be in the range [93 • -176 • ], which range depends on the values chosen for the step, facet-liquid and solid-liquid surface energies (see Table 1).

Proportion of facets as a function of growth conditions

In order to grow upwards with the same velocity, V I , as the displacement rate of the melting isotherm, T m ,t h ef a c e t velocity, V F ,s h o u l db e( s e eF i g .2):

V F = V I cos θ (4)
For a non-dislocated grain, the facet undercooling, T , depends on the facet velocity, V F ,t h r o u g hE q .( 2) (for a dislocated grain, Eq. ( 3) should be used). This undercooling gives the temperature of the coldest part of the facet, which means that the facet size, l, can be calculated as a function of the temperature gradient at the solid-liquid interface, G th :

l = h tan θ = T (V F ) G th tan θ (5)
If the facet size is larger than the half width of the grain, the grain will be totally faceted. Consequently, for a given grain size, there are growth conditions, V I and G th , under which the solid-liquid interface is totally faceted. Fig. 3 is a velocity- 2), ( 4) and ( 5) with a mean angle θ = 54.7 • , which gives the transition between a fully faceted and a partially faceted interface, for two grain sizes: 2 mm, typical of the structure obtained by Cold Crucible Continuous Casting (4C), and 1 cm, as obtained by Directional Solidification (DS).

The horizontal arrow, for a gradient of 20 K m -1 , represents Miyahara's experiments [27], in which the temperature of the solid-liquid interface has been measured with a thermocouple, as a function of the growth rate. For a low growth rate the interface undercooling was negligible, showing a mainly rough interface, and it increased up to 3 K for large velocities, showing a mainly faceted interface. Their grain size was about 2 mm and the agreement with Fig. 3 is surprisingly good. This supposes that their sample did not show dislocation, which is supported by the large measured undercooling.

It can be seen that the two growth techniques, DS and 4C, whose typical growth conditions are represented by stars in the figure, fall into the faceted region: this explains why so much twins can be seen and measured in these wafers [7,8],a s twinning can only occur in the presence of facets. This also explains why equiaxed grains occur in these processes as soon as significant carbon allows SiC precipitation, because nucleation of Si on the SiC necessitates undercooling which is provided b yt h ef a c e t s .T h ed i a g r a mp r e s e n t e di nF i g .3 is quite rough because it does not take into account crystallographic details or the real kinetic of a faceted GGL-TPL. In case of a preferred orientation of the grains, the angle θ should be adjusted accordingly and, for a given process, it will be possible to include more detailed kinetic considerations (see Section 4).

The four possible configurations of the GGL-TPL

As seen in Section 2, the boundary around one grain can be faceted or rough, depending on the crystallographic orientation and conditions. Consequently, when two grains are facing in order to build a grain boundary, three configurations are possible.

(1) When two rough boundaries are facing, it gives a Rough-TPL, see Fig. 4, bottom, and Fig. 5. This situation is very well known in the solidification of metallic alloys. The model developed by Bolling [28] has been applied to the case of Si [9], and has given a groove height of 17 micrometers which gives an undercooling T = 0.017 K (for G th = 1000 K m -1 ) at the groove bottom: clearly these values are negligible. As in the case of metals, the grain boundary angle to the growth direction, φ, is dictated by the kinetic coefficients k 1 and k 2 of the two grains, which depend on their orientations θ 1 and θ 2 toward the growth direction (Fig. 5):

ϕ = π 2 -arcsin k 1 (θ 1 ) k 2 (θ 2 ) (6)
Some data are available in the case of Si, Atwater [29] gives

k 100 k 110 = 1.
2, and Steinbach [30] uses (2) When two facets are facing, this gives a Faceted-TPL, see Fig. 4,t o p ,a n dF i g .6. This configuration has been observed at the surface of EFG ribbons, where marking of the solid-liquid interface occurs because of pulling shaft vibrations [START_REF] Surek | [END_REF]32]. Interface marking during Bridgman growth of multi-crystalline GaSb has also shown Faceted-TPL [33]. It is clear that the undercooling at the bottom of the groove is maximum and the same for both facets. Therefore they should have the same facet velocity, V F , and the grain boundary angle ϕ toward the growth direction n I should follow the bisector of the angle between the two facets, θ f : It depends only on the crystallographic orientation of the two grains. From Fig. 6, considering that the groove bottom moves up at the velocity V I , it comes:

V F = V I sin(θ/2) cos(ϕ) (7)
This allows, through Eq. ( 2)o rE q .( 3), computing the undercooling and then the groove height. In the case where there is no step source on the facet, the undercooling is of the order of 4 K under typical growth conditions, and the groove can be one centimeter deep or even more. When Eq. ( 3) is used, the groove depth is of the millimeter order, which is in agreement with experimental results [34].

(3) When a facet is facing a rough boundary, it gives a Rough-Faceted-TPL, Fig. 7. The rough side fixes a low undercooling at the groove bottom, always smaller than 0.1 K, as computed from [28], which makes nucleation of 2D seeds on the faceted side very rare. Therefore the facet cannot move significantly and the grain boundary simply follows the {111} facet. This configuration has been observed experimentally by Börzönyi [35] in the case of the nematic-smectic faceted interface in C 4 H 9 -(C 6 H 10 ) 2 CN: the facet moved very slowly compared to the growth rate.

(4) A fourth configuration of TPL concerns the intersection of a Σ3 twin boundary with the solid-liquid interface. Crystallographic relationships show that it belongs to the Faceted-TPL but with the particularity that this twin boundary is the only one that does not present an amorphous layer between the solid grains. Its surface energy is much lower than any other solid-solid boundary, so that it necessarily grows always in the same direction, giving the well-known straight boundaries that can be observed on any multi-crystalline Si solar cell. In spite of the fact that it shows a Faceted-TPL, there is no chance that twinning occurs on these facets, because this would create another Σ3 twin boundary plus a Σ9 boundary, which is totally unlikely from the energetic point of view.

Conclusion

The morphology of the grain-grain-liquid triple phase line in multi-crystalline Si has been discussed on the basis of classical kinetic, energetic and crystallographic considerations.

It has been shown that the TPL can be faceted under given conditions and that the facets are likely to overrun the whole solid-liquid interface, especially under relatively high growth rates and low temperature gradients. This enhances the occurrence of twinning and equiaxed growth and therefore importantly controls the grain structure of the ingot.

When two grains are facing, the TPL can show faceted or rough sides, which gives three different configurations, each of them with its particular grain boundary orientation law. The Σ3 twin boundary is a particular case also taken into account.

In summary:

-When the TPL is rough on both sides, the grain boundary orientation results from the kinetics of the two grains (as in classical metallurgy); -When the TPL is faceted on both sides, the grain boundary orientation is the bisector of the two facets; -When the TPL is rough on one side and faceted on the other, the grain boundary follows the facet direction; -When the TPL results from a Σ3 twin boundary, the grain boundary follows the twin plane. This opens the possibility to predict the grain structure in the Si ingot, in the case of pure material. The next step of this research will be the development of a two-dimensional simulation software in order to model the structure of a Si ingot as a function of the growth conditions. Further development could include 3D simulation but a severe limitation is that it should take into account not only the direction of the grain boundaries, but also the behavior of the grain-grain-grain triple lines during growth, which is totally unknown at the moment and should deserve further research.

On the experimental point of view, there is a need to check the validity of the proposed model and to feed it with appropriate physical properties, especially kinetic coefficients as a function of orientation. In situ observation of the solid-liquid interface morphology by synchrotron X-ray densitometry is aiming to check the model and has already confirmed the existence of the Faceted-TPL during Si growth, see the paper published in the same issue of this journal [34].

It is worth to notice that the structure of a given TPL is likely to evolve during the growth process, because the growth conditions, or the orientation of the solid-liquid interface, are changing. Also, multiplication of dislocations with time can have a very strong effect on the facet sizes and then on the evolution of TPL structures. Therefore a Faceted-TPL may become a Rough-TPL and vice versa.

Considering further improvements of the model, impurity effect should be added. This could be through modifications of surface energies, but also by taking into account interface destabilization. Finally the reader should keep in mind that grain boundaries can also be generated in the solid state, through dislocation aggregation or mechanical twinning, phenomena which have not been taken into account here.

Fig. 1 .

 1 Fig. 1. A 100 -orientated grain has four facets all around its grain boundary, each of them corresponding to a [011] direction. (a) Seen from the melt; (b) cut along the 0 11 direction.

Fig. 1 .

 1 Fig. 1. Un grain orienté 100 présente 4 facettes sur sa frontière, chacune dans une des directions 011 . (a) Vu depuis le liquide ; (b) vu en coupe suivant la direction 0 11 .

Fig. 2 .

 2 Fig. 2. Calculation of the facet size.

Fig. 2 .

 2 Fig. 2. Calcul de la taille d'une facette.

Fig. 3 .

 3 Fig. 3. Diagram showing the growth conditions under which the solid-liquid interface is fully faceted. Temperature measurements by [27] are shown on the right and have been performed along the horizontal arrow for 2 mm grains. Stars show the operating conditions of the DS and 4C processes.

Fig. 3 .

 3 Fig. 3. Carte montrant les conditions de croissance sous lesquelles l'interface est facettée. Les mesures de températures de [27] sont montrées à droite et ont été réalisées suivant la flèche horizontale. Les étoiles montrent les conditions opératoires des procédés DS et 4C.

k 100 k

 100 111 = 1.5, with k 111 the kinetic coefficient of surfaces close to {111}.

Fig. 4 .

 4 Fig. 4. Three grains (all 100 -orientated) seen from the liquid, giving a Rough-TPL and a Faceted-TPL.

Fig. 4 .

 4 Fig. 4. Trois grains (tous supposés d'orientation 100 ) vus depuis le liquide, montrant des lignes triples rugueuse et facettée.

Fig. 5 .

 5 Fig. 5. The growth configuration of the Rough-TPL.

Fig. 5 .

 5 Fig. 5. Configuration de croissance d'une ligne triple rugueuse.

Fig. 6 .

 6 Fig. 6. The growth configuration of the Faceted-TPL.

Fig. 6 .

 6 Fig. 6. Configuration de croissance d'une ligne triple facettée.

Fig. 7 .

 7 Fig. 7. The growth configuration of the Rough-Faceted-TPL.

Fig. 7 .

 7 Fig. 7. Configuration de croissance d'une ligne triple rugueuse-facettée.

Table 1

 1 Selection of surface energy values for Si interfaces. Une sélection de valeurs des énergies de surface dans le cas du silicium.

	Tableau 1
	Interface type	Surface energy	References
	Solid-liquid a	0.4 to 0.68 J m -2	[16,17]
	Solid-liquid step energy	0.087 J m -1	[18]
	Solid-solid Σ3 twin	Theory: 0.026 to 0.077 J m -2	[19]
		Experiments: 0.060 J m -2
	Solid-solid Σ9t w i n	0 . 6 4Jm -2	[20]
	Random grain boundary	1.3 J m -2

Acknowledgements

This research has been supported by the ANR-HABISOL-SiX and the OSEO-SNC research projects and through a PhD grant of the MESR attributed in the frame of Grenoble-INP BQR. Authors acknowledge the constructive help and discussions of Prof. D. Bellet, their colleagues of IM2NP-Marseille and the whole Si-X team.