Viet-Duc Le
email: viet-duc.le@ensam.eu

Etienne Pessard

Franck Morel

Interpretation of the fatigue anisotropy of additively manufactured TA6V alloys via a fracture mechanics approach

Keywords:

 for modelling small fatigue crack growth. Finally, a simulation of the generalized Kitagawa-Takahashi is presented and good agreement with the experimental data is shown.

Introduction

During recent years, additive manufacturing (AM) technologies have continued to grow steadily due to industrial and academic developments. One of the main advantages of this technology is the capacity to fabricate complex geometries directly from a computer-aided design. Among several metal AM processes, Selective Laser Melting (SLM) is one of the most commonly used. In this process, rather than single points, a line-shaped scanning pattern is used to melt powder within a powder-bed layer thanks to the energy of the laser source. Upon completion of a single layer, a new powder-layer is deposited and subsequent melting is pursued to continue the manufacturing process to build the 3D geometry.

The TA6V (Ti6Al4V) titanium alloy is one of the most common alloys manufactured by the SLM process and has applications in many industries, in particular the aeronautic and the medical domains. The mechanical properties of TA6V alloy obtained by this process have been intensively studied in the literature [START_REF] Leuders | On the fatigue properties of metals manufactured by selective laser meltingthe role of ductility[END_REF][START_REF] Edwards | Fatigue performance evaluation of selective laser melted Ti-6AL-4V[END_REF][START_REF] Xu | Ti-6AL-4V additively manufactured by selective laser melting with superior mechanical properties[END_REF][START_REF] Gong | Influence of defects on mechanical properties of ti 6al 4v components produced by selective laser melting and electron beam melting[END_REF]. For the monotonic tensile behaviour, it has been observed that the tensile strength of AM TA6V is equivalent or even higher than the alloys obtained by conventional processes [START_REF] Bian | Mechanical properties and microstructural features of direct laser-deposited Ti-6AL-4V[END_REF][START_REF] Rafi | A comparison of the tensile, fatigue, and fracture behavior of Ti-6AL-4V and 15-5 PH stainless steel parts made by selective laser melting[END_REF][START_REF] Facchini | Ductility of a Ti-6AL-4V alloy produced by selective laser melting of prealloyed powders[END_REF]. However, in terms of fatigue behaviour, the AM TA6V alloy shows significantly lower fatigue strength [START_REF] Leuders | Structural components manufactured by selective laser melting and investment casting impact of the process route on the damage mechanism under cyclic loading[END_REF]. This fact is due on one hand to the high roughness level of the as-built surface [START_REF] Chan | Fatigue life of titanium alloys fabricated by additive layer manufacturing techniques for dental implants[END_REF][START_REF] Wycisk | Effects of defects in laser additive manufactured Ti-6AL-4V on fatigue properties[END_REF]. On the other hand, the presence of microstructural defects, especially porosity, can also largely reduce the fatigue strength [START_REF] Günther | Fatigue life of additively manufactured Ti-6AL-4V in the very high cycle fatigue regime[END_REF][START_REF] Greitemeier | Fatigue performance of additive manufactured TiAl6v4 using electron and laser beam melting[END_REF].

The principal aim of this paper is to contribute to the understanding of the impact of porosity on the fatigue behaviour of TA6V alloy obtained by the SLM process and to compare the influence of porosity to that of the microstructure. In the literature, the effect of pore size on the fatigue strength of TA6V alloy fabricated by AM processes has been investigated by numerous researchers [START_REF] Wycisk | Effects of defects in laser additive manufactured Ti-6AL-4V on fatigue properties[END_REF][START_REF] Günther | Fatigue life of additively manufactured Ti-6AL-4V in the very high cycle fatigue regime[END_REF][START_REF] Greitemeier | Fatigue performance of additive manufactured TiAl6v4 using electron and laser beam melting[END_REF][START_REF] Fatemi | Multiaxial fatigue behavior of wrought and additive manufactured Ti-6AL-4V including surface finish effect[END_REF][START_REF] Fatemi | Torsional fatigue behavior of wrought and additive manufactured Ti-6AL-4V by powder bed fusion including surface finish effect[END_REF]. However, the results from these studies are quite scattered, principally due to the large variability in the fabrication process. AM processes require the use of many parameters that can affect the material properties, such as building direction, scan speed, laser power, scan pattern, temperatures, the initial powder properties (e.g. granulometry, new/recycled powder ratio) and also the postheat treatments (e.g. annealed, hot isostatic pressing, etc…). Consequently, the variability in the microstructure and porosity (including the porosity density, pore size and pore morphology) is very significant. Even with the same process conditions and same operator, the variability between fabricated batches can also be notable [START_REF] Romano | Fatigue properties of AlSi10mg obtained by additive manufacturing: Defect-based modelling and prediction of fatigue strength[END_REF]. One possible explanation for this variability is that pore formations are complex phenomena and not fully controlled. Even with the same process parameters, the largest pore size that governs the fatigue strength can vary significantly between batches and between specimens. For this reason, the number of specimens must be large enough in order to cover this variability and to reliably characterize the effect of porosity.

A vast experimental fatigue campaign for uniaxial loads with a load ratio of R = 0.1 has been undertaken. Seven different specimen batches have been investigated. They were obtained using different AM building directions (horizontal, vertical and diagonal). All of the specimens were machined to avoid the effect of the as-built surface. The fabricated fatigue specimen size is relatively large (with a diameter of 8 mm and a gauge length of 24 mm) in order to be representative for industrial components. In total, 77 fatigue specimens were tested and all the fatigue failure surfaces were examined by SEM. In addition, the microstructure and the porosity of the investigated alloy was observed by optical microscopy on polished samples taken from failed fatigue specimens.

In the first part of this paper, the characterization by optical microscopy of the microstructure and the porosity is presented as well as the tensile strength properties, to demonstrate the differences between the three building directions. More specifically, the microstructural anisotropy and the pore morphology anisotropy are highlighted. In the second section, the Wöhler curves are shown followed by characterization of the fatigue damage mechanisms. The size of pores at the crack initiation sites is measured and compared with the porosity characterized by the metallographic methodology. A link between the two pore size distributions is established by using the Murakami approach.

In the last section of this paper, the effect of pore on the fatigue behaviour is presented. Several analyses are shown to highlight the effect of the porosity, including the characterization of the critical pore size distribution and generalized Kitagawa-Takahashi maps linking the three variables: stress, number of cycles to failure and pore size. The influence of the pore size on the fatigue life is then modelled by using two fracture mechanics approaches: the Paris law and the Caton approach ( [START_REF] Caton | The effect of solidification rate on the growth of small fatigue cracks in a cast 319-type aluminum alloy[END_REF][START_REF] Caton | The influence of heat treatment and solidification time on the behavior of small-fatigue-cracks in a cast aluminum alloy[END_REF][START_REF] Caton | Demonstration of an endurance limit in cast 319 aluminum[END_REF]). A comparison of these two simulations is also shown. Finally, the generalized Kitagawa-Takahashi map is simulated by using the Caton approach and good correlation is obtained with the experimental data.

Specimen fabrication, metallurgical characterization and tensile strength

Specimen fabrication

The TA6V alloy used in this study is grade 23 ELI dedicated to medical and aeronautical applications. The chemical composition in weight percentage is shown in Table 1.

The TA6V powder, obtained by inert gas atomization, was supplied by AP&C with a 20-63 μm particle range. This particle distribution was confirmed by laser granulometry (Mastersizer Malvern, results: D10 = 30 μm, D50 = 43μm, D90 = 61 μm). This implies that 10% of the volume of the power particles are lower than 30μmand 10% higher than 60 μm.

The particle shape has also been observed by Scanning Electron Microscopy (SEM). Fig. 1) shows that the particles have a mostly spherical shape with no satellites and no agglomeration.

The specimens were fabricated by the Jules Verne Research and Technology Institute in France (IRT Jules Verne), using a SLM280 HL machine with a fabrication chamber volume of 278 × 278 × 325 mm 3 . The fabrication process was optimized thanks to a vast experimental campaign in which several process parameters such as laser power, scan speed, hatch distance and focus were taken into consideration. A porosity density level of approximately 0.01% (determined via optical microscopy observations on polished samples) is obtained.

The process parameter used to fabricate the specimens tested in this work were: a layer thickness of 30 μm, laser power of 175 W, beam diameter of 80 μm, a scanning velocity of V = 775 mm/s. The hatch distance (distance between two laser tracks) was 120 μm. The building platform is heated at 200 °C. The scanning strategy is a stripe pattern type of 10 mm width to be more representative of production job. The traverse (or scan) direction of the laser changes by 67 °between each layer.

The fatigue specimen geometry after machining (shown in Fig. 2) is given as per the ISO 2072:2010 standard and is the same for all batches. The specimen shape produced by SLM is the final geometry with a machining allowance of 1 mm. The SLM contour with a thickness of approximately 200 μm is therefore totally removed by the machining operation.

Table 1

Referred chemical compositions of the investigated TA6V alloy.

In total, seven fatigue specimen batches with three building directions were investigated. Note that the process parameters were the same for all of the specimen batches. The building directions are 0°(horizontal), 45°(diagonal) and 90°(vertical). The reason for fabricating several batches with the same building direction is to investigate the variability due to the process. Table 2 summarizes the fabricated fatigue specimen batches. Fig. 3 shows batches corresponding to the three building directions. It can be seen in these figures that a support structure was used for the whole specimen length for the 0 °specimens while for the 45 °and 90 °specimens, a support structure was used only for one of the specimens heads.

After fabrication, all of the specimens were post-heat treated by a typical heat treatment (annealed at 850 °C for 2 h following by slow cooling within the furnace) to relax the residual stresses. All specimens were machined after the heat treatment. The Ra roughness of the machined surface is lower than 0.8 μm. The roughness measurements were undertaken using a Bruker ContourGT-K 3D Optical profilometer and the ISO 4287 standard.

Microstructural characterization

It is well known in the literature that the microstructure of TA6V alloy, obtained by AM processes, is columnar [START_REF] Thijs | A study of the microstructural evolution during selective laser melting of Ti-6AL-4V[END_REF] with grains growing in the building direction. In the present work, all of the specimens were fabricated and post-heat treated with the same parameters and heat treatment conditions. Therefore, the microstructures of all investigated batches is expected to be the same for all specimen orientations. The microstructure was observed by using an optical microscope on polished and chemically etched samples (by using Kroll solution for 20 s). The observed samples were cut from the fatigue specimens after fatigue testing. The observations were firstly made on planes parallel to the building direction as shown in Fig. 4(a)-(c). It can be seen that the microstructure has the same columnar structure for all of the three building directions. However, it should be noted that the observations on planes perpendicular to the specimen axis are not the same and that these are planes of maximum principal stress, which are often considered as being critical planes in terms of the fatigue behaviour of defect-containing alloys (Fig. 4(a),(d) and (e)).

Porosity characterization

Porosity characterizations were conducted on 3 batches 0°-B1, 45°-B1 and 90°-B1. The analysed pores were observed on the planes perpendicular to the specimen axis, i.e. principal stress planes.

The characterization methodology proposed by Murakami [START_REF] Murakami | Effects of small defects and nonmetallic inclusions[END_REF] was used. As per this method, the size of the biggest pore observed in each standard inspection area of size S 0 is measured. The resulting distribution corresponds to the extreme value distribution related to the area S 0 . In the present work, the inspection area S 0 is equal to the area of the taken images, S 0 = 2.34 mm 2 . The area Fig. 3. Batches corresponding the the three building directions [START_REF] Le | Influence of porosity on the fatigue behaviour of additively fabricated ta6v alloys[END_REF]. and the Maximum Feret diameter of the pores is used to define the pore size. For each batch, 300 images of size S 0 were analysed. Theses images were taken from samples cut from 3 fatigue specimens in the specimen's gauge length after fatigue tests. The three specimens were intentionally chosen so that their fatigue strengths are representative of the fatigue strength of the batch.

Firstly, it was observed an high anisotropy of the geometry of the biggest pores. It can be seen that the largest pores detected by this methodology seem to be formed by the balling (or spatter) effect. This phenomenon has been widely investigated in the welding processes [START_REF] Wu | Understanding of spatter formation in fiber laser welding of 5083 aluminum alloy[END_REF][START_REF] Zhang | Observation of spatter formation mechanisms in high-power fiber laser welding of thick plate[END_REF]. In additive manufacturing, the balling effect has also been studied [START_REF] Gunenthiram | Experimental analysis of spatter generation and melt-pool behavior during the powder bed laser beam melting process[END_REF][START_REF] Koutiri | Influence of SLM process parameters on the surface finish, porosity rate and fatigue behavior of as-built inconel 625 parts[END_REF][START_REF] Darvish | Reducing lack of fusion during selective laser melting of CoCrMo alloy: Effect of laser power on geometrical features of tracks[END_REF][START_REF] Anwar | Study of the spatter distribution on the powder bed during selective laser melting[END_REF][START_REF] Khairallah | Laser powder-bed fusion additive manufacturing: Physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones[END_REF][START_REF] Kasperovich | Correlation between porosity and processing parameters in TiAl6v4 produced by selective laser melting[END_REF]. This effect can be described as small balls of molten metal created near the laser melting zone and blown by gas flow. It is shown in Fig. 5 that the size of "balling" particles can be up to 200 μm, which is much larger than the powder particle size (controlled by sieving with a maximum diameter of 100 μm). Because the laser power was optimized for the powder particles size (by considering several criteria such as porosity density, residual stresses and the distortion of as-fabricated parts), the "balling" particles cannot be totally melted and this fact results in the formation of lack-of-fusion pores.

The geometry of the "balling effect" related pores has been investigated by Kasperovich et al. [START_REF] Kasperovich | Correlation between porosity and processing parameters in TiAl6v4 produced by selective laser melting[END_REF] using micro-tomography observations. These authors showed that the pores relating to the balling phenomenon have a very spread geometry, i.e. a large area but small thickness. The plane on which pores spread out is generally perpendicular to the building direction. A same observation was also showed in the work of Romano et al. [START_REF] Romano | Qualification of {AM} parts: extreme value statistics applied to tomographic measurements[END_REF] for an AlSi10Mg alloy obtained by the SLM process.

From the microstructural and porosity characterizations, Fig. 6 shows a sketch illustrating the microstructural anisotropy and the pore morphology anisotropy for the 3 building directions. This figure is important as it illustrates the relationship between the microstructure, the build direction and the pore shape. Specifically it shows the columnar grains are elongated in the direction parallel to the building direction, and that the pores are large but spread out in planes perpendicular to the building direction. One of the aims of this work is to answer the question as to which is more important to consider in fatigue, the oriented columnar microstructure or the porosity population.

The pore size distributions for the three batches 0°-B1, 45°-B1 and 90°-B1 measured on the planes perpendicular to the specimen axis are shown in Fig. 7. The y-axis is presented by using the reduced variable, Y, defined as a function of the empirical cumulative probability, F, given as the following:

= -- Y ln lnF ( ) (1)
The Gumbel distribution [START_REF] Gumbel | Statistical theory of extreme values and some practical applications: a series of lectures Applied mathematics series[END_REF] is used to fit the experimental data and shown in Fig. 7 as the linear regression lines of the non-zero experimental points. Greater details concerning this methodology can be found in references [START_REF] Le | Simulation of the Kitagawa-Takahashi diagram using a probabilistic approach for cast al-si alloys under different multiaxial loads[END_REF][START_REF] Le | Multiaxial high cycle fatigue damage mechanisms associated with the different microstructural heterogeneities of cast aluminium alloys[END_REF].

It can be observed that in general the pore size of the 0°-B1 specimens is smaller than that of 45°-B1 specimens, both in terms of area and the Feret diameter. The pore size of the 90°-B1 batch terms of area is slightly larger than the 0°-B1 but smaller than the 45°-B1 batch. In terms of Feret diameter, the pore size of 90°-B1 batch is the smallest. This trend is surprising because given the observed anisotropy of the large pore geometry as shown in Fig. 6, it would be expected that the pore size of the 90°batch is the largest. In Section 3.3 presented below, it will be shown that the measured pore size distributions on polished samples are in good agreement with the pore size observed at the crack initiation sites for the 0°-B1 and the 45°-B1 batches. However, for the 90°-B1 batch, the pore size measured on polished samples under-estimates the porosity of the batch.

This disagreement could be due to the small thickness of the pores. At 90°the elongated direction of the pores is parallel to the Fig. 5. Pores formed by the balling effect observed on a 0°(left) and 90°(right) specimen.

Fig. 6. Illustration of the microstrucural anisotropy and the pore morphology anisotropy of the investigated alloy [START_REF] Le | Influence of porosity on the fatigue behaviour of additively fabricated ta6v alloys[END_REF].

observation planes (see Fig. 6). The probability that an observation plane cuts a large pore is therefore small at 90°. The surface size based on 300 images with S 0 = 2.34 mm 2 is probably not enough to determine the largest defect size and correctly characterized the pore distribution of the batch.

In order to verify this hypothesis, a comparison of the pore size distribution between measurements on the perpendicular planes and measurements on the parallel planes for the 90°-B1 specimens is shown in Fig. 8. It can be seen that the pore size measured on the parallel planes is larger than the one measured on the perpendicular planes, especially in terms of Feret diameter.

Tensile behaviour

Monotonic tensile tests were conducted using an INSTRON servo-hydraulic machine on flat specimens with rectangular crosssection of 2 mm × 6 mm and gauge length of 25 mm. The strain was measured thanks to an extensometer with an initial length L 0 =25 mm and a displacement range = ± L Δ 2.5 mm. All of the tests were conducted with a constant strain rate of 0.005 min -1 at ambient temperature.

The comparison of the ultimate tensile strength, the yield strength and the fracture strain of the three building directions is shown in Fig. 9. The values of the yield strength, σ Y , and ultimate strength, σ UTS , were normalized by dividing by the mean value of the yield Fig. 7. Pore size distributions measured on polished samples of the three batches 0°-B1, 45°-B1 and 90°-B1. The observation planes are perpendicular to the specimen axis as shown in Fig. 6.

Fatigue behaviour and damage mechanisms

All of the fatigue tests presented below were carried out at ambient temperature and pressure in laboratory air. The fatigue tests were conducted with a constant stress amplitude, a load ratio R = 0.1 and a frequency of 20 Hz. A maximum fatigue life of × 2 10 6 cycles was used. The stopping criterion was chosen to be the complete rupture of the specimen. The run-out specimens that survived × 2 10 6 cycles were re-tested at a higher load.

Wöhler curves

Fig. 10a-c shows the S-N curves for the 7 investigated batches. The values of the maximum stress, σ max , are normalized by the average value of the ultimate strength, σ UTS , for the three building directions as defined by the following equation. Note that the maximum difference in the ultimate tensile strength between the three building directions is lower than 5%.

= + + °°°σ σ σ σ σ ( ) / 3 max n max UTS UTS UTS , , 0 ,45 ,90 (2)
The regression lines correspond to the Stromeyer theoretical model [START_REF] Kasi | The determination of fatigue limits under alternating stress conditions[END_REF] as given by Eq. (3). The three parameters A B , and σ max n 0,

were identified to best fit the experimental data (see Table 3).

= -× -

Log N A B Log σ σ ( ) ( ) max n max n , 0 , (3)
Firstly, a large degree of variability can be seen between the batches with the same building direction, in spite of the fact that the same process parameters were used. Measurements of the fabrication conditions such as the oxygen level or the plate and chamber temperatures also showed that there was no notable difference between batches. Similar variability was also seen in the work of Romano et al. [START_REF] Romano | Fatigue properties of AlSi10mg obtained by additive manufacturing: Defect-based modelling and prediction of fatigue strength[END_REF] for the SLM process on an aluminium alloy. The scatter in the S-N curves is relatively high compared to wrought TA6V alloy [START_REF] Nalla | Influence of microstructure on high-cycle fatigue of Ti-6AL-4V: Bimodal vs. lamellar structures[END_REF], especially for the 0°-B2 and 45°B2 and B3 batches. It is interesting to note that these batches show a better fatigue strength compared with the other batch with the same building direction. In other words, it seems that the fatigue behaviour of the batch with a better fatigue resistance is more scattered than the one with lower fatigue strength. Fig. 10(d) shows the comparison of the fitted Stromeyer curves of all of the batches as well as the fatigue data bands of the TA6V alloy obtained by the wrought and cast processes documented in [START_REF] Donachie | Titanium: A technical guide[END_REF]. It can be seen in Fig. 10(d) that generally, the fatigue strength of the 90°batches is the lowest. However, no clear tendency is observed between the 0°and 45°specimens. In the literature, the results concerning the effect of the building direction of machined specimens show different trends. In the work of Wycisk et al. [START_REF] Wycisk | Effects of defects in laser additive manufactured Ti-6AL-4V on fatigue properties[END_REF], no effect of the building direction was observed between the 45°and 90°specimens. However, in the work of Edwards and Ramulu [START_REF] Edwards | Fatigue performance evaluation of selective laser melted Ti-6AL-4V[END_REF] and Nicolet to [START_REF] Nicoletto | Anisotropic high cycle fatigue behavior of Ti-6AL-4V obtained by powder bed laser fusion[END_REF], the fatigue strength of the 90°specimens is lower than other directions. Comparing to the wrought and cast TA6V alloys, it can be stated that the fatigue strength of the 0 °and 45 °batches is comparable to the wrought TA6V while for the 90 °batches, the fatigue strength slightly higher than the cast TA6V alloy.

Fatigue crack initiation mechanisms

All of the fatigue failure surfaces were analysed using a conventional ZEISS scanning electron microscope. It was observed that all of the fatigue cracks initiated from only one site. Almost all the fatigue cracks initiate from a pore (76 out of 77 specimens). Only one Fig. 9. Tensile behaviour of the investigated alloy fabricated with the three building directions [START_REF] Le | Influence of porosity on the fatigue behaviour of additively fabricated ta6v alloys[END_REF].

specimen shows crack initiation with no detectable pore.

Most pores at the crack initiation sites are "lack of fusion" pores, located on the specimen surface (66 of 77 specimens) as shown in Fig. 11. Interestingly, it can also be seen in Fig. 11 for 45°and 90°specimens that "balling" particles seem to be the origin of the formation of the critical pore.

An other interesting observation is that the anisotropy of the pore morphology, related to the balling phenomenon (mentioned previously in Section 2.3), can also be seen in the analysis of the failure surface. Fig. 12 shows some pores at the crack initiation sites of 0°(left) and 90°(specimens).

It can be observed that the critical pores for the 0°specimens are generally elongated (i.e. the Minimum Feret diameter is much lower than the Maximum Feret diameter) while for the 90°specimens, the critical pores are more equiaxed.

Pore size distribution used for fatigue strength prediction

It was shown in the characterization of the pore size distribution obtained by optical microscopy on polished samples, taken from three batches 0°-B1, 45°-B1 and 90°-B1 (Section 2.3), that the pore size of the 45°-B1 batch is the largest, while the difference in the pore size between the 0°-B1 batch and the 90°-B1 batch is not clear. However, the order of the fatigue resistance for these three batches is not in agreement with the pore size distributions. In fact, the fatigue strength of the 90°-B1 is the lowest, while the 0°-B1 batch shows the highest fatigue strength (Fig. 10(d)). In order to understand this discrepancy, a comparison of the pore size distributions characterized on polished samples with the distributions of pore size measured at the crack initiation site is shown in Fig. 13.

Firstly, it can be seen that the pore size distributions measured on the fatigue failure surfaces is in a good agreement with the fatigue strength rank of the three batches (i.e. the 90°-B1 batch shows the lowest while the 0°-B1 batch shows the highest fatigue strength).

Secondly, it seems that the pore size distribution of the 90°-B1 batch characterized on polished samples by optical microscopy under-estimates the porosity level of the batch. In fact, the critical pore size distribution measured on the failure surfaces of the 90°-B1 batch is the largest while the distribution from polished samples for this batch shows a relatively low pore size compared to the other batches. One possible explanation is due to the anisotropy of the critical pore morphology (shown in Fig. 12). That is, the possibility that a OM observation plane cuts a large pore is relatively small for the 90°building direction. Consequently, the measured pore size distribution for the 90°batch is not representative of the porosity density of the batch, even though a relatively large sample size was examined (i.e. 300 images of size S 0 = 2.34 mm 2 were analysed for each batch).

It can also be seen in Fig. 13 that the pore size distributions measured on polished samples is smaller than the ones measured at the crack initiation sites on the fatigue failure surfaces. This difference is widely observed in the literature ( [START_REF] Le | Etude de linfluence des htrognits microstructurales sur la tenue en fatigue grand nombre de cycles des alliages daluminium de fonderie[END_REF]) and can be explained by the fact that the material volume characterized by the metallographic method is generally smaller than the gauge volume in the fatigue specimens. Hence, as per the extreme value theory ( [START_REF] Gumbel | Statistical theory of extreme values and some practical applications: a series of lectures Applied mathematics series[END_REF][START_REF] Weibull | A statistical theory of the strength of materials[END_REF]), the pore size extreme distribution is larger when the volume of material characterized is higher.

In order to illustrate this "size effect" and to link the two distributions for each building orientation, it is proposed, as done in previous work [START_REF] Le | Etude de linfluence des htrognits microstructurales sur la tenue en fatigue grand nombre de cycles des alliages daluminium de fonderie[END_REF][START_REF] Osmond | Effect of porosity on the fatigue strength of cast aluminium alloys: from the specimen to the structure[END_REF], to use an extension of the Murakami approach [START_REF] Murakami | Effects of small defects and nonmetallic inclusions[END_REF] proposed by Makkonen et al. [START_REF] Makkonen | Size effect in fatigue based on the extreme value distribution of defects[END_REF]. As per this approach, the cumulative distribution function of the critical pore size at the crack initiation site associated with the loaded volume V is given by:

⎜ ⎟ = ⎡ ⎣ ⎢ - ⎛ ⎝ - -′ ⎞ ⎠ ⎤ ⎦ ⎥ F area V exp exp area μ β ( , ) ( ) (4)
where ′ = + μ μ βln V V ( / ) 0 . μ and β are the location and scale parameter of the Gumbel distribution, metallographically identified by the line of best-fit in Fig. 7.

The reference volume V 0 characterized by the metallographic method can be calculated by = × V max area S ( ) 0 0 where max area ( ) is the maximum pore size observed and S 0 = 2.34 mm 2 is the inspection area. The loaded volume V is the sub-surface hollow cylinder with a height of 24 mm, i.e. the specimen gauge length. The thickness of the hollow cylinder is 250 μm corresponding to the maximum pore size observed on the failure surfaces. The reason for using the sub-surface volume is that it has been shown in the literature ( [START_REF] Le | Multiaxial high cycle fatigue damage mechanisms associated with the different microstructural heterogeneities of cast aluminium alloys[END_REF][START_REF] Dezecot | 3d characterization and modeling of low cycle fatigue damage mechanisms at high temperature in a cast aluminum alloy[END_REF][START_REF] Serrano-Munoz | Influence of surface and internal casting defects on the fatigue behaviour of a357-t6 cast aluminium alloy[END_REF]) for the high cycle fatigue (HCF) regime that surface pores are much more detrimental than internal pores. The dominate fatigue crack often initiates from surface pores, even though bigger pores are present in the bulk. Fig. 14 shows the result of the prediction for the three batches 0°-B1, 45°-B1 and 90°-B1 for which the metallographical pore size characterization was done.

Good agreement between the prediction and the experimental data can be seen for the 0°-B1, 45°-B1 batches. However, for the 90°-B1 batch, the prediction under-estimates the critical pore size. This observation emphasizes both a volume effect and the difficulty to correctly characterize the critical pore size distribution from the metallographic approach for each building direction.

The effect of porosity on the fatigue strength

In this section, the effect of pore size on the fatigue behaviour of all of the batches will be analysed and modelled. As shown in the Fig. 13. Comparison of the pore size distributions obtained from observations on polished samples using optical microscopy (square data points labelled "OM") and on failure surfaces(circular data points labelled "Crack initiations"). Fig. 14. Correlation between the maximum pore size distributions predicted from the metallographical pore size characterization and the pore size distributions measured on the crack initiation sites. previous section, the pore size distribution obtained by the metallographic method should not be assumed representative for all the building directions. For this reason, only the pore size distributions measured on the crack initiation sites will be considered in these analysis.

Correlation between the critical pore size distributions and the S-N curves

The pore size defined as the area measured at the crack initiation sites are plotted in Fig. 15 for all fatigue batches. By comparison with the S-N curves shown in Fig. 10(d), it can be seen that the lowest fatigue strengths, obtained for the 90°b atches are in a good agreement with the largest pore sizes compared with the 0°and 45°batches. For the 0°and 45°batches, no difference in the S-N curves is clearly visible. Furthermore, because of the large scatter in the Wöhler curves, the effect of the pore size on the fatigue strength can not be readily highlighted by correlation with the pore size distributions. In order to better understand the pore size effect, a map that links the applied stress (S), the pore size measured at the crack initiation site ( area ) and the fatigue life (N f ) is presented.

Generalized Kitagawa-Takahashi diagram: stress, pore size, fatigue life maps

The original Kitagawa-Takahashi diagram was proposed to show the relationship between the defect size and the fatigue strength at the infinite fatigue life of the material. The generalized Kitagawa-Takahashi diagram, initially proposed by Ciavarella et al. [START_REF] Ciavarella | On the possible generalizations of the Kitagawa-Takahashi diagram and of the el haddad equation to finite life[END_REF] according to the authors knowledge, is used to show the relationship between three parameters: the fatigue life, the applied stress and the defect size. The authors presented this relationship in the "defect size -fatigue strength" plane using the fatigue life as isocurves.

In the present work, the fatigue tests focused on the limited fatigue life domain with only a few data points were obtained for the "infinite" fatigue life. The run-out number of cycles was defined as × 2 10 6 cycles. Specimens that survived × 2 10 6 cycles were considered to have an "infinite" fatigue life.

As discussed in the previous section, due to the large scatters in both the pore size distributions and the S-N curves, the classical diagrams such as the original Kitagawa-Takahashi or the Gumbel pore size distribution can not clearly highlight the role of the pore size on the fatigue behaviour. The generalized Kitagawa-Takahashi is hence an interesting way to analyse the data obtained at different fatigue lives.

Fig. 16(a) shows the generalized Kitagawa-Takahashi diagram for all pores for the three building directions. A linear interpolation in terms of Log N ( ) f was used to fill the maps. The Y-axis scale is the same for all of the three maps. The run-out specimens are also presented with a fatigue life of × 2 10 6 cycles and with a pore size measured on the failure surface of the same specimen after re-tested at higher load.

Considering this representation of lifetime behavior, two interesting observations can be highlighted. Firstly, the yellow arrows shown in the Fig. 16(a) point to perturbations, or specimens that have a better fatigue strength compared to neighbouring data. Observation of the fatigue failure surfaces, for these specimens showed that these points correspond to crack initiation from an internal pore as shown in Fig. 16(b). This effect of initiation from an internal pore has been widely accepted in the literature ( [START_REF] Serrano-Munoz | Influence of surface and internal casting defects on the fatigue behaviour of a357-t6 cast aluminium alloy[END_REF][START_REF] Murakami | Stress intensity factors handbook[END_REF]).

The second observation concerns the effect of the pore morphology. It was mentioned previously in Section 3.2 that, for some of the 45°-B2 specimens, the crack initiated from spherical gas pores. By examining the map for the 45°batches (see Fig. 17(a)), a slightly better fatigue life can be seen for the specimens in which the final crack initiates from a spherical gas pore (indicated by yellow arrows).

The effect of pore morphology on the fatigue strength has been discussed in the literature ( [START_REF] Dezecot | 3d characterization and modeling of low cycle fatigue damage mechanisms at high temperature in a cast aluminum alloy[END_REF][START_REF] Le | Investigation of the effect of porosity on the high cycle fatigue behaviour of cast al-si alloy by x-ray microtomography[END_REF]) for both experimental and numerical studies.

By removing the data related to the internal crack initiation and the spherical gas pores and by considering only the surface "lackof-fusion" pores, the map shown in Fig. 18 is constructed. Better agreement can be observed between the stress, the pore size and the fatigue life, when compared to the original map. The effect of the oriented microstructure does not seem pronounced. As an example, the two points enclosed by the white ellipse in Fig. 18 correspond to two critical pores with similar size and geometry. The fatigue life of the two specimens are similar (67619 and 69350 cycles) however the first one is a 45°specimen and the second is a 90°specimen.

Modelling the effect of pore size on the fatigue life

As shown in the previous section, the fatigue resistance of the investigated alloy is controlled mainly by the pore size while the effect of the local microstructure does not seem pronounced. In this section, a model to predict the fatigue life, taking into account the pore size will be presented. Note that, in this section, crack initiation from internal pores (observed on 5 specimens out of a total of 77) and from spherical gas pores (observed on 5 specimens out of a total of 77) are not taken into account.

In the literature, several approaches for modelling the fatigue life, taking into account the defect size have been proposed. Among them, the fracture mechanics-based approach using the Paris law ( [START_REF] Paris | A rational analytic theory of fatigue[END_REF][START_REF] Paris | A critical analysis of crack propagation laws[END_REF]) is largely used to model the crack growth rate and predict the fatigue life. Similarly, ) proposed a model for small fatigue crack growth and applied it with success to cast aluminium alloys. In the present study, these two approaches will be applied and compared for the investigated alloy.

• The Paris law:

= da dN C K (Δ ) n (5)
where K Δ is the cyclic stress intensity factor (SIF). For microstructural crack-like pores, this factor can be calculated by the Murakami equation ( [START_REF] Murakami | Stress intensity factors handbook[END_REF]) given by:

= K Y σ π area Δ Δ ( 6
)
where Y is equal to 0.65 for surface pores. By substituting Eq. ( 6) into Eq. ( 5) and given that = -σ Rσ Δ

(1 ) max , the initial crack growth rate from a pore can be determined by the following equation:

= ′ d area dN C σ area ( ) ( ( ) )
max n 1/2 [START_REF] Facchini | Ductility of a Ti-6AL-4V alloy produced by selective laser melting of prealloyed powders[END_REF] where 1)) n n/2 and R is the load ratio.

′ = - C C Y R π ( (
• The Caton approach:

Caton assumed that the small crack growth rate is a function of the macroscopic maximum strain ε max , the stress amplitude σ a and the yield strength of the material σ y , given by the following equation:

⎜ ⎟ = ⎛ ⎝ ⎜ ⎛ ⎝ ⎞ ⎠ ⎞ ⎠ ⎟ da dN D ε σ σ a max a y s t ( 8
)
where a is the defect size.

For the tests undertaken in our work the macroscopic maximum strain ε max can be replaced by σ E / max (where E is the Young's modulus) because the maximum applied fatigue loads always result in a maximum stress which is lower than the macroscopic yield strength of the material. If the Murakami's parameter area is used to define the defect size, the initial crack growth rate from a pore is given by:

= ′ ′ d area dN D σ area ( ) ( ( ) ) max s t 1/2 (9) where ′ = - ( ) D D R Eσ st 1 2 y and ′ = t st 2 .
It can be seen from Eqs. ( 7) and ( 9 1/2 for the Caton approach can be used as variables to estimate the initial crack growth rate from a pore and therefore its final fatigue life. In Fig. 19, the relationship between the fatigue life and these terms for all the specimen batches is shown. In order to compare with the raw data in the S-N curves (shown in Fig. 19(a)), the variable is normalized with reference pore size of = area 150 μm.

The value of = s 1.67 has been determined so that the scatter of the corrected S-N curves is the lowest. One of the most interesting advantages of the "corrected stress versus fatigue life" diagram is that the validity of the criteria under investigation can be estimated with limited influence of the model parameter identification. For the Paris law, no parameters need to be identified. However for the Caton approach, only one parameter s is necessary (instead of 3 parameters ′ D s , and ′ t for the full model). In the literature, similar diagrams have been used to investigate the effect of the pore size and to validate a criterion. For example, in the work of Gunther et al. [START_REF] Günther | Fatigue life of additively manufactured Ti-6AL-4V in the very high cycle fatigue regime[END_REF], the authors plotted the K I max It can be clearly seen in Fig. 19 that the scatter in the corrected S-N diagrams is largely reduced, compared to the raw data. The Caton approach shows smaller scatter compared to the Paris law. Another observation is that no notable difference can be observed between the three building directions. A similar observation was shown in Section 4.2. In other words, the effect of the oriented microstructure, that is the columnar microstructure, is not the key factor controlling the fatigue strength.

To model the generalized Kitagawa-Takahashi map, the Caton approach is chosen. Caton et al. showed in their work ( [START_REF] Caton | Demonstration of an endurance limit in cast 319 aluminum[END_REF]) that the proposed model under-estimates the fatigue strength of the material for the unlimited fatigue life domain, so a combination with the crack growth threshold intensity factor K Δ th seems necessary. A second term representing the asymptotic fatigue limit line is added to Eq. ( 8) as shown in Eq. ( 12). This modification was also used in the work of Romano et al. [START_REF] Romano | Fatigue properties of AlSi10mg obtained by additive manufacturing: Defect-based modelling and prediction of fatigue strength[END_REF].

⎜ ⎟ = ⎛ ⎝ ⎜ ⎛ ⎝ ⎞ ⎠ ⎞ ⎠ ⎟ ⎛ ⎝ - ⎞ ⎠ > da dN D ε σ σ a K K K K 1 Δ Δ ifΔ Δ max a y s t th p th ( 12
)
Four material parameters D s t , , and p, were identified to best-fit the experimental data. The resulting generalized Kitagawa-Takahashi maps is shown in Fig. 20a-c. Fig. 20(a) shows the predictions in log-log scale using the iso-value curves. Fig. 20(b) is nothing but Fig. 18 presented in log-log scale using the iso-value curves. Fig. 20(c) shows the error map in Log-Log scale. The error between the model and the raw data is defined as follow:

= N N Error f Exp f Sim , , (13)
where N f Exp It can be noticed that the discontinuities of the iso-value curves in log N ( ) f represented in Fig. 20(b) are more pronounced in the unlimited fatigue life domain in which the lack of experimental data is probably the cause.

In Fig. 20(a), the utilized model does not permit predicting the fatigue limit in the region without defect. This is simply because the LEFM based criterion used to model the fatigue strength at the unlimited fatigue life uses a constant stress intensity factor, K Δ th ,that is not appropriate to model the fatigue strength of the porosity-free material. In the present study, this fact is not really problematical because all the analysed fatigue specimens showed a crack initiation from pore. It can be seen in the error map that the predictions are quite good in this region with the error lies between 0.8 and 1. ).

Conclusion

This paper deals with the fatigue behaviour of TA6V titanium alloy fabricated by the SLM process. A great deal of experimental data was obtained for 7 specimen batches with 3 building directions. Several analyses have been done in order to finely characterize the effect of pore size on the fatigue strength. A generalized Kitagawa-Takahashi diagram is shown to be an appropriate way to analyse the stress, fatigue life and pore size relationships.

The principal results are as follow:

• It has been highlighted that the oriented microstructure (columnar grains) and the pore geometry anisotropy ("balling effect" related pores) are linked to the building direction.

• The characterization of the porosity by optical observations on polished samples under-estimates the critical pore size for the 90°s pecimens. For the 0°and 45°, the measured pore size distributions by this methodology is in agreement with the order of the fatigue resistance and can be used to predict the critical pore size distributions by the Murakami approach.

• The effect of the pore size on the fatigue behaviour, investigated using the generalized Kitagawa-Takahashi diagrams, showed that he fatigue resistance is mainly controlled by the pore size, characterized by the Murakami's parameter area . The effect of the oriented microstructure is not pronounced. • Two fracture mechanics approaches, the Paris law and the Caton model, were used to predict the fatigue life. By using the "corrected" S-N curves, the Caton approach shows better agreement than the Paris law. The predictions of the fatigue life obtained by using the Caton approach was correlated with the experimental data thanks to the generalized Kitagawa-Takahashi diagram and shows relatively good results.

Fig. 1 .

 1 Fig. 1. SEM analysis of powder particles.

Fig. 2 .

 2 Fig. 2. Final geometry of the fatigue specimens [50].

Fig. 4 .

 4 Fig. 4. (a) to (c): observation planes parallel to the building direction; (d) and(e): observation planes perpendicular to the specimen axis for 45°and 90°specimens.

Fig. 8 .

 8 Fig.8. Comparison of the pore size distributions measured on perpendicular planes and parallel planes for the 90°-B1 specimens.

Fig. 10 .

 10 Fig. 10. Wöhler curves for all of specimen batches: (a) 0°batches, (b) 45°batches, (c) 90°batches and (d) superposition of the Stromeyer theoretical curves of all batches. The arrows indicate the run-out specimens. The fatigue data bands for wrought and cast TA6V is documented in [35].

Fig. 11 .

 11 Fig. 11. Fatigue crack initiation site observed on a 0°(left), 45°(middle) and 90°(right) specimen [50].

Fig. 12 .

 12 Fig. 12. Anisotropy of the pore geometry observed on the fatigue failure surfaces for 0°(left) and 90°(right) specimens.

Fig. 15 .

 15 Fig. 15. Pore size distributions measured on the fatigue failure surfaces.

Fig. 16 .

 16 Fig. 16. (a) Generalized Kitagawa-Takahashi diagram or the Sarea -N f map for all pores. The arrows indicate specimens for which crack initiation occurred from internal pores and (b) SEM images showing crack initiation from an internal pore [50].

Fig. 17 .

 17 Fig. 17. (a) Generalized Kitagawa-Takahashi map for 45°batches showing the effect of spherical gas pores (indicated by arrows) and (b) SEM images of crack initiation from a spherical gas pore.

Fig. 18 .

 18 Fig.18. S-area -N f map for only surface lack-of-fusion pores (i.e. without internal pores and spherical gas pores); circular dots = 0°specimens, square dots = 45°specimens, triangle dots = 90°specimens.

  ) that mathematically, the Caton equation is an extension of the Paris law with the Murakami SIF. When = s 1 the Caton model simplifies to the Paris law. Hence, the terms Yσ area

,

  against the fatigue life while in the work of Masuo et al.[START_REF] Masuo | surface roughness and HIP on fatigue strength of Ti-6AL-4V manufactured by additive manufacturing[END_REF] and Yamashita et al.[START_REF] Yamashita | Defect analysis and fatigue design basis for ni-based superalloy 718 manufactured by selective laser melting[END_REF], a normalized stress σ σ / w based on the Murakami empirical model =

,

  and N f Sim , are respectively the experimental and simulated number of cycles to failure.

Fig. 19 .

 19 Fig. 19. Comparison of the scatter of the S-N curves between (a) the raw data, (b) the Paris law corrected-stress data with the Murakami SIF, (c) the Canton corrected-stress data with = s 1.67. The Y-axis of the three diagrams have the same scale.

Fig. 20 .

 20 Fig. 20. The generalized Kitagawa-Takahashi diagrams iso-values curves in log-log scale, modelled using the Caton approach (a) and built from the experimental data (b); The error map in log-log scale (c).

Table 2

 2 Summary of the fabricated fatigue specimen batches.

	Building Dir.	Batch Ref.	Batch specimen nb.	Total
	0	0°-B1	10	20
		0°-B2	10
	45	45°-B1	14	26
		45°-B2	6
		45°-B3	6
	90	90°-B1	20	31
		90°-B2

Table 3

 3 Parameters of the Stromeyer equation fitted for each batch.

	Batch	A	B	σ max n 0,
	0-B1	4.250	0.742	0.450
	0-B2	4.204	1.022	0.550
	45-B1	3.861	1.150	0.400
	45-B2	4.042	1.049	0.600
	45-B3	3.309	2.106	0.500
	90-B1	3.908	1.238	0.275
	90-B2	3.728	1.193	0.400

  6 (i.e.

	N 0.8 f Sim ,	⩽	N f Exp ,	⩽	N 1.6 f Sim ,

Acknowledgement

This study is part of the FATAL project, managed by IRT Jules Verne (French Institute in Research and Technology in Advanced Manufacturing Technologies for Composite, Metallic and Hybrid Structures). The authors wish to associate the industrial and academic partners of this project; Respectively ACB, Arts et Métiers ParisTech, CNRS, DAHER, Ecole Centrale de Nantes, Europe Technologies, FIVES, Renault and University of Nantes.


S Salhi

Florent Bonnefoy

Stéphane Girard

Maxime Bernier

Nicolas Barbot

Romain Siragusa

Etienne Perret

Frédéric Garet

Enhanced THz tags authentication using multivariate statistical analysis

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

I. INTRODUCTION

Principal Component Analysis (PCA) is mainly used for the classification of materials and authentication in various applications domains [START_REF] Güler | Evaluation of graphical and multivariate statistical methods for classification of water chemistry data[END_REF], [START_REF] Behera | Palm print Authentication Using PCA Technique[END_REF] and has recently attracted research in the THz domain [START_REF] Sommer | Terahertz Spectroscopy of Aged Epoxy Resin Adhesive[END_REF], [START_REF] Xie | Terahertz data combined with principal component analysis applied for visual classification of materials[END_REF]. This method is also commonly used to reduce the dimension of a statistical set of samples [START_REF] Hardle | Applied Multivariate Statistical Analysis[END_REF]. However, more efficient multivariate statistical methods have been developed specifically for classification like the Linear Discriminant Analysis (LDA) [START_REF] Hardle | Applied Multivariate Statistical Analysis[END_REF]. In this study, we propose to use both PCA and LDA to address the authentication of THz tags with the aim to fight counterfeiting of manufactured goods. First, PCA is used to both reduce the dimensions of the statistical measurements set and classify the different tag signatures. Secondly, LDA is used to increase the discriminative rate of these tags. The THz tag structure used in this study is schemed in Fig. 1. A 1D rectangular diffraction grating is engraved on the surface of a plastic card. The "substrate" layer acts as a dielectric waveguide whose spectral response depends on the geometrical parameters of the grating engraved on its top. In the present case these parameters are: period Γ = 800 µm, groove's depth p = 280 µm, card thickness ecard = 775 µm, card refractive index ncard = 1.73 and absorption coefficient α = 17 cm -1 @ 1 THz. We use the presence of guided modes appearing as narrow lines in the transmitted spectra for example, more precisely their frequency positions and width, as a potential unique signature for the tag [START_REF] Bonnefoy | Authentication in the THz domain: a new tool to fight counterfeiting[END_REF], [START_REF] Ali | Potential of chipless authentication based on randomness inherent in fabrication process for RF and THz[END_REF]. In order to prove that such a signature can be used to unitary authenticate a tag, 44 tags were realized with the same pattern and protocol; these tags can be supposed as theoretically identical. However, due to random uncertainties during the fabrication process, each signature presents tiny variations that we take advantage using PCA and LDA.

II. RESULTS

Each tag was repeatedly measured in transmission one hundred times using classical THz-TDS system in the frequency range 200 GHz -800 GHz (spectral range of the maximum dynamics), with a frequency resolution δf = 8.6 GHz and under an incident angle θ = 10°.

Here, we thus considered the 4400 measurements of the twin tags, each containing 70 frequency points, as a statistical dataset thereby composed of 3.08 10 5 samples. We then first calculate the covariance matrix of this dataset [START_REF] Hardle | Applied Multivariate Statistical Analysis[END_REF]. Afterwards we diagonalize this matrix to obtain the associated eigenvalues and eigenvectors. From these latter, it is possible to get the total explained variance and the privileged directions of maximum variance, which are known as the principal components axes in PCA. The hundred extracted spectral responses of each tag were then projected on the two first principal components of PCA, corresponding to the maximal variance (Fig. 2). Each set of hundred measurements leads to a cluster whose probability density can be fitted by a bivariate normal distribution. Note that the cluster size is directly related to the measurement reproducibility. We proceed to a dimension reduction with PCA as a pre-processing step before a Gaussian classification and we evaluate the error rates. For that, a Gaussian classifier is then implemented using this bidimensional model [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. The error probability ε k,t FN to wrongly authenticate the tag k as the tag t (False-Negative case, FN), is directly linked to the clusters size and proximity. ε k,t FN is numerically calculated from the N measurements of the k th tag using: With fk,t(xi): the bivariate normal distribution associated to the tag k and tag t. The probability Ek of confusing the k th tag with any others is therefore given by: = ∑ , , ≠

, = ∑ < , ~ (1)
where M=44 and N=100 in the present case. As expected (cf. Fig. 2), isolated clusters present the lowest error rates (e.g. E4 < 0.02%, green ellipse), while clusters packed each other are more hardly authenticated due to their close vicinity (e.g. E25 < 2.7%, red ellipse). We found a global error rate, calculated by averaging the error rates Ek of 0.5%. This means that PCA proposed in this study can unitary authenticate 1 tag randomly chosen among 44 others, with a success score of 99.5%. In other hand, we show that the 5 first principal axes can explain 99.5% of the total variance. In such a case, the dimensions of the statistical problem could be then reduced by about a factor 15, with an information loss lower than 0.5%. We now use LDA [START_REF] Martinez | PCA versus LDA[END_REF] to try to better separate all the different classes. LDA permits to find the optimal base to project the clusters (blue and red scatter plot in Fig. 3). Let notice that using LDA, the obtained classes, corresponding to such projections can be much better separated than using PCA (see blue and red Gaussian curves). On the other hand, PCA remains necessary to reduce the dimension of the statistical sample before the classification process but not for efficient classification. In Fig. 4, we also plot the rate of correct classification versus the number of axis we take into account for classification. Let notice that this rate can reach 99.994% using only the 3 first axis from LDA.

To evaluate the robustness of the classification, we considered a training database constituted of 70% from the original 4400 measurements while the other 30% will constitute the test database. We then calculate the correct classification rate of each sample in the test database. The whole process is randomly repeated 30 times to test the robustness of the classification process. Thereby we also showed that in these conditions, only 6% of the 4400 measurements, with at least 1 measurement on each tag, is sufficient to obtain the previously mentioned classification rate. Let notice that this result can be different if changing the original database measurements i.e. the set of tags and/or THz reader.

III. CONCLUSION

In this study, we demonstrate the ability to use dielectric planar waveguides associated with a diffraction grating as a tag for unitary authentication in the THz domain. More precisely we take advantage of the very small geometrical differences (λ/10 -λ/100) intrinsically and randomly induced by the fabrication process. To analyze the tags signatures we combined PCA and LDA to respectively reduce the dimension of the statistical sample and to efficiently classify them. A correct classification rate of 99.994% has been obtained.

Fig. 1 :

 1 Fig. 1: a, Picture of the proposed tag and b, scheme of the structure.

Fig. 2 :

 2 Fig.2: Projections on the principal component space of the hundred measurements performed on each tag. For sake of clarity, each set of measurements has a different colour to its closest neighbours letting appear clusters, which are referenced by the tag number. The explained variances are 68% and 22% for the 1 st and 2 nd principal components, respectively.

Fig. 3 :

 3 Fig. 3: Projections of two different classes (blue and red dots) in the principal component space from PCA (1 st and 2 nd axis, left) and on the two first axes from LDA (1 st and 2 nd axis, right).

Fig. 4 :

 4 Fig. 4: Correct classification rate versus the number of LDA axes considered.

ACKNOWLEDGMENT

This work is supported by the ANR AUSTRALE project, grant ANR-18-CE39-0002 of the French Agence Nationale de la Recherche.


Target tracking in the time-frequency domain for a driving aid application Zineddine Bouhroum * † , Thierry Chonavel * † , Ronan Fablet * † and Pascal Deloof ‡ for ACC (Automotive Cruise Control). The associated processing relies on the demodulation of the received signal by the emitted waveform, leading thus to harmonic signals with distinct frequencies during increasing and decreasing frequency periods. These frequencies are linear combinations of the Doppler and the delay of the target, such that the latter can be estimated using a Fourier transform. This approach is however waveform-specific. A more general purpose detection and estimation step consists in performing matched filtering at the outputs of a bank of Doppler demodulators, thus reconstructing cross-ambiguity function between emitted waveform and received signal. Detection threshold can account for ambiguity function side lobes level.

In this paper, we propose a general purpose method for target tracking. This is done by tracking the waveform ambiguity function displacement inside the time-frequency domain, by considering the observed cross-ambiguity function at each observation time. By using a state space model, we are able to perform the target parameters estimation and filtering at the same time. In classical techniques, the observation is restricted to an initial guess for target delay and Doppler, which amounts to focusing on the maximum of the cross-ambiguity function.

Abstract-Driving

aid is an important emerging radar ap-plication. Among associated technological issues, detecting and tracking potential obstacles (cars, pedestrians, static objects . . . ) is of key interest. In this context, we propose an original tracking method, based on particle filtering, where the observation of the state vector is the crossambiguity function of the backscattered signal. Simulations for various waveforms demonstrate the im-provement brought compared to the conventional approach, in which a detection step precedes the tracking. Index Terms-radar; tracking; driving aid; ambiguity func-tion; cross-ambiguity function; particle filtering.

I. INTRODUCTION

About 80% of road accidents are due to human errors. European statistics show that in most cases, the root cause of the accidents is linked to a bad perception of the environment, followed by a wrong action of the driver. These factors make the development of perception and driving assistance an area of very active research. These systems aim at improving the perception of potential obstacles, even in cases of partial occlu-sions or bad weather conditions (fog for example). Detecting and tracking both static and moving, potential obstacles, is among the main challenges to be addressed.

A typical radar system relies on the emission of a waveform and on the analysis of the backscattered signal, with a view to extracting information on the position (related to the time delay), and the velocity (related to the Doppler shift) of the objects. The accuracy of measurement depends both on the selected waveform and the processing method.

For target tracking, the classical method involves two suc-cessive steps: first, the detection of the potential targets and the estimation of position and velocity parameters independently at each time, second, a filtering technique, typically a Kalman filter o r a nother B ayesian fi l ter [START_REF] Ramachandra | Kalman filtering techniques for radar tracking[END_REF].

In practice, particular structures of the waveforms under consideration are generally exploited. For instance, FMCW (Frequency Modulated Continuous-Wave) waveforms are a typical example [START_REF] Griffiths | New ideas in FM radar[END_REF]. Composed of two successive linear frequency-modulated pulses (Chirp) [START_REF] Levanon | Radar Signals[END_REF], with increasing and decreasing frequency respectively, these waveforms are used In our approach, on the contrary, we fully exploit available data in the tracking procedure by considering the whole crossambiguity function. More precisely, the target contribution appears as a weighted and time-frequency shifted version of the ambiguity function in the cross-ambiguity function. By doing so, although the procedure is not waveform specific, specific time-frequency properties of the waveform under consideration are accounted for. Thus, we can expect higher parameter tracking performance with our approach than with the usual estimate then track approach.

Anticipating on the simulation part of the paper, Figure 1 shows how this approach can be fruitful. A comparison of the tracking performance for three standard waveforms (chirp, Barker codes and pulse train waveforms) demonstrates that a better precision is achieved in all cases with the proposed approach. In addition, accounting for the shape of the ambiguity function is a natural way to incorporate information about noisy signal amplitude statistics around the maximum of the cross-ambiguity. Clearly, this is not achieved in the classical approach. It may be point out that associating detection and tracking in a same processing step has already been considered in Track Before Detect (TBD) and Bayesian TBD techniques (see e.g. [START_REF] Ristic | Beyond the Kalman filter, particle filters for tracking applications[END_REF]). In TBD we are however more concerned with detection than with estimation. To our best knowledge, there is no study where joint tracking and estimation is considered and the ambiguity function spreading is large compared to radar image sampling steps.

The paper is organized as follows. In Section II, we recall the notion of ambiguity and cross-ambiguity function. In section III, we address the tracking problem and its general Bayesian solution. In Section IV, we describe the tracking method through our state space formulation of the target tracking problem. In section V, we evaluate the performances of our approach compared to the classical one.

II. AMBIGUITY AND CROSS-AMBIGUITY FUNCTION

A. Ambiguity function

The ambiguity function is a two-dimensional function of time delay and Doppler frequency χ(τ, f ). In a radar application, it is defined as the time response of an adapted filter to emitted signal, when this one is received shifted in time by a delay τ and in frequency by a Doppler shift f relative to the nominal values (zeros) [START_REF] Levanon | Radar Signals[END_REF]. Formally, it is given by:

χ u (τ, f ) = +∞ -∞ u(t)u * (t -τ ) exp(i2πf t)dt ( 1
)
where u is the complex envelope of the signal. A target moving away from the radar implies a positive f. The ambiguity function fully characterizes a given waveform u.

The more localized the ambiguity function in the delay-Doppler space, the more accurate the estimation of the range-Doppler parameters. For a given signal energy, the total volume under the normalized ambiguity function (squared) is constant and independent on the signal waveform. Consequently, a trade-off has to be achieved between the peakness of the ambiguity function and the level of the side-lobes. This results in the proposal of various waveforms.

B. Cross-ambiguity function

The cross-ambiguity function is the time response to the backscattered signal r of an adapted filter to the delay-Doppler shifted emitted signal [START_REF] Levanon | Radar Signals[END_REF]:

χ r,u (τ, f ) = +∞ -∞ r(t)u * (t -τ ) exp(i2πf t)dt (2)
The received signal is the sum of the contributions of the different targets, where each target backscatters the signal with specific delay and Doppler parameters corresponding to its position and speed, and an unknown scaling coefficient corresponding to the attenuation and the phase shift of the signal:

r(t) = n α n u(t -τ n ) exp(i2πf n )
The cross-ambiguity function is then the sum of translated ambiguity functions in the direction [τ n f n ]

T corresponding to the delay and Doppler shift of each target:

χ r,u (τ, f ) = n α n χ u (τ -τ n , f -f n ) (3)

III. NONLINEAR BAYESIAN TRACKING

The tracking problem is largely explained in [START_REF] Ramachandra | Kalman filtering techniques for radar tracking[END_REF]. A target motion model is given by the evolution of its state vector sequence defined by the state equation

X k = f k (X k-1 , V k-1 )
, where f k is a possibly nonlinear function, and V k-1 a process noise sequence. Target tracking resorts to recursively estimate X k from the observation modeled by the equation:

Y k = h k (X k , W k )
, where h k is a possibly nonlinear function, and W k a measurement noise sequence. Within a Bayesian setting, it amounts to stating the probability density function p(X k |Y 1...k . It may be obtained recursively in a two-stage prediction-update scheme :

• Using the state equation, the prediction step determines the prior density function at time k via Chapman-Kolmogorov equation:

p(X k |Y 1...k-1 ) = p(X k |X k-1 )p(X k-1 |Y 1...k-1 )dX k-1 , (4)
• From the new observation , an update step of the estimate is carried out. It exploits the likelihood function defined by the observation model. The updated density function is then given by Bayes rule:

p(X k |Y 1...k ) = p(Y k |X k )p(X k |Y 1...k-1 ) p(Y k |X k )p(X k |Y 1...k-1 )dX k . (5)
These two steps modeled by the recurrence equations ( 4) and ( 5) form the optimal Bayesian solution for target tracking. It can be determined analytically only for particular cases such as the linear Gaussian problem (Kalman filter). Generally, the optimal Bayesian solution is approximated by linearization (Extended Kalman filter), unscented transformation (unscented Kalman filter), or Monte Carlo simulation (particle filter) [START_REF] Arulampalam | A Tutorial on Particle Filters for Online/Non-Gaussian Bayesian Tracking[END_REF].

IV. TRACKING METHOD

A. State Model

We consider a single target moving along a one-dimensional axis corresponding to the direction of the radar. It is often the case in an automotive application as vehicles to be tracked move along the same direction as the vehicle equipped with the radar.

To model the dynamics of this target in the delay-Doppler space, the state vector of the target, at time step k, k = 1, . . . , K, is given by:

τ (k) f d (k) a rel (k)
T , where τ is the delay corresponding to the position of the target (distance between the radar and the target), f d is the Doppler frequency corresponding to the velocity of the target relatively to the velocity of the tracker vehicle, and a rel the relative acceleration of the target. From the formulation of a nearly uniform accelerated motion model, given in [START_REF] Bühren | A global motion model for target tracking in automotive applications[END_REF] among others, where we consider v rel = v obj -v ego , a rel = a obj -a ego , and the linear equations τ = 2d/c, f d = 2v rel /λ 0 , the state equation resorts to:

X k =   τ (k) f d (k) a rel (k)   = F X k-1 + QV k (6)
where:

F =    1 T f0 T 2 c 0 1 2T λ0 0 0 1   , Q =    T 3 3c T 2 λ0 T   .
T is the sampling period, c the celerity of light, f 0 the carrier frequency and λ 0 the corresponding wavelength. V k denotes the noise model of relative accelerations. The prior distribution p(X k |X k-1 ) is issued from the state equation [START_REF] Bühren | A global motion model for target tracking in automotive applications[END_REF]. Considering a white centered Gaussian noise V k , p(X k |X k-1 ) will be a Gaussian distribution, with mean F Xk-1 , where Xk-1 is the state estimator at time step k -1, and a variance related to Q.

B. Observation Model

The observation that we consider consists of a discretized version of the cross ambiguity function :

Y k (a, b) = χ r,u (a.δτ, b.δf ) = n α n (k)χ u,u (a.δτ -τ n (k), b.δf -f n (k)) + W k (a, b) (7)
where W k is the observation noise. a ≤ A, b ≤ B are integers. The observation is then a matrix A × B.

Given a state vector, the associated prediction of the discretized cross ambiguity function χ r,u can be compared to the real noisy observation Y k up to the unknown scaling factor α n (k). The likelihood model can be then be defined from a metric between the observation Y k and the prediction Y X k .

Formally, it is stated as a Gibbs distribution, like in image processing tracking techniques [START_REF] Rathi | Tracking Deforming Objects using Particle Filtering for Geometric Active Contours[END_REF] and is given by:

p(Y k |X k ) ∝ exp(-γ∆(Y k , Y X k ))
where ∆(Y k , Y X k ) is the selected metric between the true observation Y k and the predicted observation Y X k , and γ is a parameter that determines the shape of the distribution.

α n (k) is an unknown propagation coefficient. It possibly changes at each time step. A solution to this problem could be to choose a power invariant metric as can be the mutual information [START_REF] Rusakoff | Image similarity using mutual information of regions[END_REF]. Here, we estimate α n (k) in the mean square sense. Then, for a single target, our definition of ∆ reduces to

∆(Y k , Y X k ) = min αn(k) ||Y k -Y X k | | 2

C. Tracking Algorithm

Given the non-linearity of the observation model, we apply a sampling importance resampling particle filter (SIRPF) [START_REF] Arulampalam | A Tutorial on Particle Filters for Online/Non-Gaussian Bayesian Tracking[END_REF] , which will use the prior density p(X k |X k-1 ) given by the state equation ( 6) as the importance density, and p(Y k |X k ), given by the equation of observation [START_REF] Griffiths | New ideas in FM radar[END_REF], as the likelihood distribution.

As explained in [START_REF] Kyriakides | Target Tracking Using Particle Filtering and CAZAC Sequences[END_REF], every particle is a possible case of the state X k of the target, that we will note by X m k where m = 1 . . . M is the particle index and M is the number of particles used. At time k -1, The state X k-1 is described by M particles X m k-1 and M weights ω m k (m = 1 . . . M ), each weight expresses the probability that have the correspondant particle to be the state X k-1 , given the observation. We independently propagate each particle to time step k using [START_REF] Bühren | A global motion model for target tracking in automotive applications[END_REF], which amounts to sampling from the prior distribution p(X k |X k-1 = X m k-1 ). A predictor of the state X k , before having the observation Y k , could be given by:

Xk -= M m=1 ω m k-1 X m k
While having the observation Y k , we compute recursively the new weigts using the likelihood probability function p(Y k |X k ):

ξ m k = ω m k-1 p(Y k |X k = X m k ) m = 1 . . . M (8) ω m k = ξ m k M m=1 ξ m k (9)
The state estimator is given by: A comparison with the classical approach is carried out. For the classical approach, the detection step consists in the localization of the global maximum of cross-ambiguity function. Regarding the filtering step, a Kalman filter described in [START_REF] Arulampalam | A Tutorial on Particle Filters for Online/Non-Gaussian Bayesian Tracking[END_REF] is considered. For our approach, the particle filter exploits 300 particles. A quantitative evaluation is reported for three different waveforms: a chirp (linearly modulated pulse), a pulse train (6 pulses) and Barker codes waveform ([1 1 1 1 1 -1 -1 1 1 1 -1 1]). Their ambiguity functions are given figure 2. We can see that for the Barker codes waveform and more clearly for the pulse train, the side-lobes are relatively high.

Xk = M m=1 ω m k X m k
Can be observed on figure 3 the improvement of tracking performance with the new method for the three waveforms. For the pulse train one, we can especially observe a failure of the tracking with the classical method. This is probably due to the fact that the pulse train ambiguity function has high side lobes. The classical approach based on estimating then filtering can have for consequence to track the side lobe and not the principle. This problem does not appear in our approach because we consider the entire ambiguity function. Table I gives a quantitative evaluation of tracking performance for the two approaches and for the three waveforms.

VI. CONCLUSION

In this paper, we have developed a new solution for target tracking from cross-ambiguity observation. We have obtained significant performance improvement compared to the classical detection-tracking strategy and improved robustness with respect to waveform diversity.

In future work, we will adapt this approach for scan antennas that are considered in new ACC radars. In this context, target bearing will also be considered and antenna pattern will play a role in bearing estimation analog to that played by ambiguity function in delay-Doppler estimation.

With a view to speeding up the algorithm, in future work, we will also refine cross ambiguity sampling strategy to make the amount of information to be handled as small as possible, while retaining high tracking accuracy. Using unscented Kalman filter or unscented particle filter is also a possible way to achieve reduced processing complexity. An extension to multiple target tracking will also be considered to comply with automotive environment.

Fig. 1 .

 1 Fig. 1. Distance tracking performance of classical (dashed line) and new (solid line) method for three waveforms : 1. Chirp, 2. Barker codes, 3. Pulse train

Fig. 2 .

 2 Fig. 2. Ambiguity function of the three simulated waveforms (modulus)

TABLE I QUANTITATIVE

 I COMPARISON OF TRACKING PERFORMANCES BETWEEN CLASSICAL AND OUR APPROACH

	waveforms	Distance error (m)	Speed error (km/h)
		Classical method New method Classical method New method
	Chirp	0.15	0.09	0.85	0.57
	Barker Codes	0.83	0.09	2.88	0.87
	Pulse train	4.93	0.06	9.21	0.47