Agnieszka Klimek-Kopyra

Denis Vile
email: denis.vile@inrae.fr

Julie Guiguitant

Hélène Marrou

Trait Diversity of Pulse Species Predicts Agroecosystem Properties Trade-Offs

Keywords: agroecosystem services, functional diversity, leaf nitrogen content, legumes, pulse crops, trait-based agroecology, traits-services trade-offs

Crop diversity management in agriculture is a fundamental principle of agroecology and a powerful way to promote resilient and sustainable production systems. Pulses are especially relevant for diversification issues. Yet, the specific diversity of legumes is poorly represented in most cropping systems. We used the trait-based approach to quantify the functional diversity of 30 pulses varieties, belonging to 10 species, grown under common field conditions. Our aim was to test relationships between traits, yield, and supporting agroecosystem properties. Our experimental results highlighted trade-offs between agroecosystem properties supported by different combinations of traits. Also, results demonstrated the relevance of leaf nitrogen content (LNC), leaf area ratio (LAR), and reproductive phenology to predict most of the trade-offs observed between agroecosystem properties. A comparison with a previous analysis based on literature data collected in diverse agronomic situations suggested that some traits are more plastic than others and therefore contribute differently to frame legumes diversity depending on the conditions of observation. Present results suggested that the implementation of such trait-based approach would rapidly benefit the selection of species/varieties for specific targeted agroecosystem services provisioning under specific (environmental or management) conditions.

INTRODUCTION

Crop diversity management in agriculture is a fundamental principle of agroecology and a powerful way to promote resilient and sustainable production systems [START_REF] Gaba | Multiple cropping systems as drivers for providing multiple ecosystem services: from concepts to design[END_REF]. Agroecosystem management offers multiple alternatives to monoculture, such as diversification with pulses, through intercropping or rotation. Their well-known ability to fix atmospheric nitrogen reduces fossil energy consumption and makes them particularly suitable for low-input systems. They are also particularly important in human nutrition as a source of proteins to complement cereal-based diets [START_REF] Singh | Improving the production and utilization of cowpea as food and fodder[END_REF][START_REF] Tharanathan | Grain legumes -a boon to human nutrition[END_REF]. Yet, the specific diversity of legumes is poorly represented in most cropping systems. While more than 80 pulses species are known to contribute to human diet, the FAO database includes only 11 of them [START_REF] Tiwari | Pulse Foods: Processing, Quality and Nutraceutical Applications[END_REF]. This lack of attention creates a risk of continued erosion of knowledge and genetic resources that could further hamper the integration of these crops into more diversified cropping systems [START_REF] Dansi | Diversity of the neglected and underutilized crop species of importance in benin[END_REF].

To increase diversity within cropping systems, choosing appropriate species is an important prerequisite [START_REF] Brooker | Improving intercropping?: a synthesis of research in agronomy, plant physiology and ecology[END_REF]. The trait-based approach is well suited for designing diversified cropping systems capable of providing specific sets of services. In ecology, the approach has been increasingly employed to link biodiversity to ecosystem functioning and services [START_REF] Lavorel | Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail[END_REF][START_REF] Violle | Let the concept of trait be functional![END_REF][START_REF] Garnier | A trait-based approach to comparative functional plant ecology: concepts, methods and applications for agroecology[END_REF]; a methodology that could be applied to cropping systems [START_REF] Litrico | Diversity in plant breeding: a new conceptual framework[END_REF][START_REF] Martin | Plant functional traits in agroecosystems: a blueprint for research[END_REF][START_REF] Wood | Functional traits in agriculture: agrobiodiversity and ecosystem services[END_REF][START_REF] Damour | A revised trait-based framework for agroecosystems including decision rules[END_REF]. However, trait-based approach requires extensive data collection to develop trait databases for crop species or varieties and document relationships between traits and services [START_REF] Barot | Designing mixtures of varieties for multifunctional agriculture with the help of ecology. A review[END_REF]. As a neglected group of species, pulses are particularly difficult to document [START_REF] Barot | Designing mixtures of varieties for multifunctional agriculture with the help of ecology. A review[END_REF], for implementing a trait-based approach.

In a previous study, [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] set up a database suitable for functional trait approach on pulses. The database gathered disparate data from multiple published and/or public sources and contained data on plant functional traits for 43 pulse species, and agroecosystem properties recorded at canopy level for a subset of these species. The analysis of the functional trait space highlighted three ecological strategies among cultivated pulses that were similar to those already described in spontaneous species in the leaf-seed-height (LSH) scheme [START_REF] Westoby | A leaf-height-seed (LHS) plant ecology strategy scheme[END_REF]. However, reduced data on agroecosystem properties and data disparity did not allow to identify clear link between traits and agroecosystem properties. In an applied perspective of supporting the integration of pulses in cropping systems, a consolidated analysis of trait-properties relationships is still needed to identify which traits would make legume species more productive, competitive, or adapted to low nitrogen or water conditions and favor production services. Trade-offs between plant traits and how they determine the different strategies also require more attention.

The aim of the present study was to identify pulses traits and/or trait combinations that can be used as predictors of yield and supporting agroecosystem services, further accounting for the relationships between traits and between agroecosystem properties. Original data were collected on 30 pulse varieties from 10 different species monitored in similar environmental conditions. Agroecosystem properties were measured at the canopy level, whereas functional traits, i.e., morphological, physiological, or phenological traits supposed to have direct or indirect effects on agroecosystem properties, were measured at the individual or organ level. We combined these original data collected within uniform environmental conditions with literature data collected in diverse agronomic situations (see [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] to assess the effect of trait measurement conditions on trait values and whether it could distort the relationships between traits and agroecosystem properties.

MATERIALS AND METHODS

Plant Material and Growing Conditions

A field experiment was conducted at UCOSEM's experimental farm in Lectoure (lat./long 43.911158/0.666182) on a deep claylimestone soil. Three varieties for each one of the 10 pulses species (Table 1) were planted in two sequential batches/subexperiments on contiguous fields of approximately 1 ha each. According to their temperature requirements: Cool season legumes were sown after the last negative temperatures, on March 14, 2018, and warm season legumes were sown on June 7, 2018, when daily temperature was over 15 • C. Each batch/subexperiment was laid out in split plot design with four replicates, varieties being subplots and species the main plot. Plots were 6 m long and counted four rows. Preceding crop was sunflower, and fields were fertilized with 0/22/8 NPK (350 kg ha -1 ) fertilizer before sowing. Each species was manually inoculated before sowing with an adequate strain of rhizobia (provider: UMR LSTM, Montpellier, France). Sowing density was chosen according to agronomical recommendations for each species and ranged from 10 to 150 plants m -2 . Manual weeding and drip irrigation were performed on a regular basis from sowing. Soil moisture was recorded weekly in each block and for each species at three depths (0-20, 20-40, and 40-60 cm) using TDR probes installed vertically at sowing (mini TRASE TM system with buriable probes, Soilmoisture Equipment Corp., CA) and confirmed that the plants did not experiment any water deficit in none of the two sub-experiments.

Measurements of Individual Plant Traits

Occurrence of phenological stages (date of 50% flowering, date of physiological maturity, i.e., when more than 50% of the pods were yellow-brown) was recorded in calendar days and then converted into thermal units using a base temperature of 0 • C for Cicer arietinum and Vicia faba, 2 • C for Lens culinaris, 3.5 • C for Lathyrus sativus, 4.1 • C for Trigonella foenum-graecum [START_REF] Covell | The influence of temperature on seed germination rate in grain legumes. I. A comparison of chickpea, lentil, soyabean and cowpea at constant temperatures[END_REF][START_REF] Lamichhane | Analysis of soybean germination, emergence, and prediction of a possible northward establishment of the crop under climate change[END_REF], 4.7 • C for Lotus tetragonolobus [START_REF] Moot | Base temperature and thermal time requirements for germination and emergence of temperate pasture species[END_REF], 6 • C for Phaseolus vulgaris [START_REF] López | Effects of heating strategies on earliness and yield of snap beans (Phaseolus vulgaris l.) grown under[END_REF], and 8.5 • C for Vigna sp. [START_REF] Craufurd | Development in cowpea (Vigna unguiculata). I. The influence of temperature on seed germination and seedling emergence[END_REF]. This allowed to compute degree days from sowing to flowering (FLO) and degree days from sowing to maturity (MAT).

Biomass was harvested at flowering stage and maturity. At each sampling date, plants were collected on 0.5 m 2 quadrats, in the two central rows to avoid border effect. Plant number was recorded, and then leaves, stems, flowers, and pods were detached and then dried separately at 60 • C until constant weight. Leaf and stem dry mass ratio over total aboveground plant mass (LMR and SMR, respectively) was calculated at flowering, and harvest index (HI) was calculated at maturity. Leaf area ratio (LAR; calculated as total plant leaf area (LA) over total aboveground mass), total LA (mm 2 plant -1 ), specific LA (SLA, cm 2 g -1 ), and leaflet length (LL, mm) as well as leaflet width (LW, mm) were measured at flowering (mean of the measurement of every leaf of a plant sub-sample), whereas thousand seed weight (TSW, g) was determined at harvest. Total nitrogen content was determined (Vario PyroCube autoanalyzer and Precision mass spectrometer) on leaves, stems, and grain samples collected at maturity and at flowering. Total plant nitrogen content was computed both at flowering and maturity. Leaf nitrogen content (LNC, mg g -1 ), and seed nitrogen content (SNC, mg g -1 ) were determined at flowering and at maturity, respectively.

Phyllochron was calculated from observations of three tagged plants on each plot for which emerged leaves were counted weekly (cool season legumes) or twice a week (warm season legumes).

The rate of soil coverage was estimated by measuring soil cover weekly, using nadiral photos taken at the same location within each plot, with a large angle camera at chest height, and then analyzed with ImageJ software [START_REF] Schneider | NIH image to ImageJ: 25 years of image analysis[END_REF]. Soil cover was modeled as a logistic function of thermal time. The maximal rate of soil coverage was defined as the slope at the inflexion point of the logistic curve.

Measurement of Agroecosystem Properties

Agroecosystem properties were measured on the same experimental plots (Table 2). Grain yield (GY, t ha -1 ), biomass yield (BY, t ha -1 ), and biomass measured at flowering (t ha -1 ) were determined from plant sampling on 0.5 m 2 quadrats. Water use efficiency (WUE) was computed as grain mass produced per mm of water input (rainfall and irrigation). Maximum soil cover (%), LA index (LAI), and the duration (in thermal units) of the phase when soil coverage is above 80% of the maximum (soil cover duration) were recorded. Finally, the capacity to fix nitrogen from the atmosphere was estimated through the percent of nitrogen derived from biological fixation at flowering and at maturity, measured with the natural isotopic dilution method, using barley and maize grown in the experiment outer borders as control plants for cool season and warm season legumes, respectively. From the 15 N abundance (d 15 N) of the samples, the proportion of plant N derived from atmosphere (%Ndfa) was calculated from Eq. 1 [START_REF] Shearer | N2-fixation in field settings: estimations based on natural 15 N abundance[END_REF]: Where β is the 15 N natural abundance of the N derived from biological nitrogen fixation (BNF) in the aerial tissue of the legume, δ15N ref is the 15 N abundance in the control plant (maize or barley) tissues, and δ15N legume is the 15 N abundance in legume tissues. The β values used are detailed in Table 1 [START_REF] Unkovich | Measuring Plant-Associated Nitrogen Fixation in Agricultural Systems[END_REF].

%Ndfa = 100 δ15N ref -δ15N legume δ15N -β ref (1)

Data Analysis

Data were first analyzed through the prism of the framework presented in [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] in order to compare the variability observed with the global functional diversity of pulses species. Data of the present field measurements were projected as supplementary individuals in a probabilistic principal component analysis (PPCA) [START_REF] Tipping | Probabilistic principal component analysis[END_REF] computed with nine common functional traits (DM, DF, HI, LL, LW, LA, SLA, LNC, and TSW) values taken from [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] database. We used Euclidian distance to quantify the distance between the scores of barycenters of each species as calculated from field data and previously published in [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF]. The consistency of species ranking in trait values across literature and experimental data was tested with Spearman's rank correlation analysis.

Second, a principal component analysis (PCA) was performed on experimental data including all traits measured. Agroecosystem properties were added as supplementary variables, in order to investigate how they fit into the phenotypic space defined by the traits. The appropriate number of PCs was determined with Kaiser' criterion which selects components that correspond to eigenvalues larger than 1. A between and within-class correspondence analysis was carried out to evaluate the percentage of variance captured by the intra and inter-specific diversity [START_REF] Doledec | Rythmes Saisonniers et Composantes Stationnelles en Milieu Aquatique. I -Description d'un Plan D'observations Complet par Projection de Variables[END_REF].

Finally, trade-offs among agroecosystem properties were analyzed using a hierarchical clustering analysis (HCA). A classification and regression tree (CART) analysis was done on raw data (n = 120) to identify how functional traits could allow to predict the clusters of varieties obtained from the HCA. Regression trees are prediction models obtained through machine-learning algorithms that recursively partition the data space in order to fit the simplest prediction model within each partition. The resulting partitioning can be represented graphically as a decision tree. The random forest aggregation method, which consists in adding a random component to the choice of variables in the model, was also used. This method does not allow to generate a visual decision tree, but it provides indices proportional to the importance of each variable in the aggregated model and thus its participation to the discrimination between groups of individuals (here, varieties). Mean decrease in Gini impurity (MDG) can be defined as the total decrease in node impurity (weighted by the proportion of samples reaching a given node) averaged across all of the trees that make up the forest. The most important variables to the model will have the largest MDG values; conversely, the least important variable will have the smallest MDG values.

All statistical analyses were performed in the computing environment R 3.5.2 (R Core Team, 2018) using pcaMethods package [START_REF] Stacklies | pcaMethods a bioconductor package providing PCA methods for incomplete data[END_REF], rpart package [START_REF] Therneau | rpart: Recursive Partitioning and Regression Trees[END_REF], ade4 package [START_REF] Dray | The ade4 package: implementing the duality diagram for ecologists[END_REF], and randomForest package [START_REF] Breiman | Random forests[END_REF].

RESULTS

Consistency of Pulses Functional Trait Space

A PPCA was performed on the 43-species-database from [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] using only traits that were also measured in the present field experiment. The first three PC axes of this PPCA explained 70% of total variance observed in the literature (Figure 1). PC1 (34%) was associated with morphological traits, especially leaf traits (LL, LW, LA, and to a lesser extend TSW). PC2 (24%) was positively associated with late flowering (DF) and maturity (DM) and to a lesser extend HI. PC3 (12%) was strongly associated to SLA and LNC. Experimental data variability on PC axes was in the same range as the literature variability. Experimental variance was close to 100 and 60% of literature variance on PC1 and PC2, respectively. On the contrary, the variance of supplementary individuals on the third axis (driven by SLA and LNC) was almost three times higher than the variance of active data from the literature (Appendix 1). For all other traits, field measurements range was always narrower than the one observed in [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF], for a larger set of species (Appendix 1).

Distance between the barycenters of each species as tested in the present experiment and from the literature database was short in the first plane of the PPCA (Figure 1). Species ranking in literature versus present data remained conserved on PC1 (ρ = 0.85, p = 0.003) and to a lesser extent on PC2 (ρ = 0.53, p = 0.10). Species rankings were rather inconsistent on PC3 (ρ = 0.24, p = 0.12). Cool season species had barycenters closely positioned along the three dimensions of the PPCA (Euclidian distance lower than 2; Table 3). Among them, C. arietinum and T. foenum-graecum had the highest Euclidian distance mostly because of the divergence with literature on the third axis, while V. faba diverged mostly on the second axis. For warm season species, distances between barycenters from the literature and from the present experiment were larger, on PC3 but also on PC2 for Vigna sp.

Characterization of Functional Trait Diversity

The first four axes of the PCA performed on 30 individuals and 16 traits from the present experiment (Table 2) captured 80% of the total variance (Figure 2). Interspecies variance explained 81% of the total variance, while intra-species explained 19%.

Leaf attributes were strongly associated with PC1 (40%). On the first principal component (PC1), LA, LL, and plant N content at flowering were negatively correlated with SLA, LNC, SNC, and LAR. This axis discriminated cool season species which were characterized by small leaflets and small LA with high leaf and seed nitrogen content and a large LAR, from warm season species characterized by large leaflet, high LA, and high plant nitrogen content. PC2 (20%) was positively correlated to LMR and negatively correlated to SMR and cover rate. The cover rate was significantly negatively correlated to MAT (ρ = -0.72, p < 0.001). SMR was also significantly negatively correlated to the phyllochron (ρ = -0.60, p < 0.001). L. tetragonolobus and 2 for abbreviations of traits.

P. vulgaris were apart from all other species on PC2, and both expressed a high LMR at flowering. PC3 (11%) and PC4 (9%) were associated to FLO and TSW, respectively. PC3 separated V. faba from L. culinaris and L. tetragonolobus, whereas PC4 discriminated C. arietinum and P. vulgaris.

All agroecosystem properties, except LAI, biomass at flowering and Ndfa at maturity were not well represented in the factor planes formed by the fourth PCs constructed with traits (Figure 2). LAI and biomass at flowering were relatively well represented in the first plane and were significantly positively correlated to LL (ρ = 0.46, p < 0.001; ρ = 0.63, p = 0.01, respectively) and LA (ρ = 0.67, p < 0.001; ρ = 0.83, p < 0.001, respectively) as well as cover rate (ρ = 0.59, p < 0.001; ρ = 0.60, p < 0.001, respectively). Although Ndfa at flowering was poorly represented on the first PC, this property was significantly positively correlated to plant N content at flowering (ρ = 0.41, p = 0.02). Ndfa at maturity was not correlated to Ndfa at flowering (ρ = 0.04, p = 0.83) but was significantly positively correlated to LNC (r = 0.49, p = 0.005) and SLA (ρ = 0.57, p = 0.001) as well as with SMR (ρ = 0.43, p = 0.02). The soil cover duration was inversely correlated to LAI (ρ = -0.21, p = 0.0428) and biomass at flowering (ρ = -0.22, p = 0.033) and significantly positively correlated with FLO (ρ = 0.47, p = 0.009) and to a lesser extent with MAT (ρ = 0.33, p = 0.08). Maximum soil cover was not correlated to soil cover duration (ρ = -0.03, p = 0.86) and was slightly associated with SMR (ρ = 0.35, p = 0.05). Finally, GY and WUE were positively correlated to PC3 which was mostly determined by FLO, whereas BY, and to a lesser extent GY and WUE, were positively associated to PC4 alongside TSW.

Trade-Offs Between Agroecosystem Properties

The classification analysis revealed seven clusters of species and varieties based on their agroecosystem properties (Figure 3). All varieties of a species generally grouped into the same cluster.

V. aconitifolia varieties were all in cluster 3 together with one variety of Vigna unguiculata. This cluster is apart from a group of close clusters composed essentially by cool season species; L. sativus varieties were grouped with L. tetragonolobus (cluster 4), L. culinaris was alone in cluster 5, cluster 6 was dominantly occupied by T. foenum-graecum, and cluster 7 grouped both V. faba and Vigna mungo. These groups were distant from cluster 1 essentially represented by C. arietinum and cluster 2 that contained varieties of P. vulgaris, and of V. unguiculata.

Each cluster was characterized by a combination of values of agroecosystem properties. Cluster 1 grouped species with high production (GY and BY) as well as a high maximum soil cover and soil cover duration, but a low Ndfa and a low biomass and LAI at flowering. Cluster 2 was characterized by a low biomass at flowering as well as a low LAI at flowering and maximum soil cover but a high Ndfa at flowering and a long soil cover duration. Cluster 3 varieties expressed the lowest yield as well as the lowest WUE but high Ndfa at maturity, a high LAI at flowering, and a long cover duration. Cluster 4 varieties expressed the highest WUE concomitantly with a high GY and a low Ndfa at both stages. Varieties of cluster 5 expressed high Ndfa both at flowering and maturity and high maximum soil cover. Yet, they showed a low maximum soil cover and soil cover duration. Cluster 6 was characterized by a low BY and GY, low Ndfa both at flowering and maturity, and poor soil covering capacity, yet LAI and soil cover at flowering were high. Finally, cluster 7 was composed of varieties producing the highest biomass at flowering, along with the highest maximum soil cover, LAI, and Ndfa at flowering.

Relevant Traits for Prediction of Agroecosystem Properties and Structuration of Functional Diversity

Classification and regression tree analysis showed very diverse trait combinations to predict the cluster affiliation (Figure 4A). LNC was identified as the most important trait to discriminate the clusters as it had the highest MDG value, which means that it allowed to gain purity in the process of classification for most of the nodes in the CART it was associated with (Figure 4B). LNC was assigned to the first division node and allowed to sort cluster 1 varieties apart from the remaining ones. Cluster 1 varieties have an LNC between 2.7 and 3 mg N g -1 . Varieties of cluster 4 were separated. Half of them fell on the same side of the tree as cluster 1 but exhibited lower LNC value (below 2.7 mg N g -1 ), while the other half had an average LNC above 3 mg N g -1 . Varieties in this second cluster 4 half also had an LMR above 43% concomitantly with low plant N content at flowering and at maturity and a longer cycle. LMR, DDR, and plant N content at maturity had strong MDG unlike the plant N content at flowering. Varieties with an LMR under 0.43 and having high MAT were assigned either to cluster 3 (characterized by a high plant N content but only at flowering), or cluster 2 or 5 (both characterized by high plant N content at flowering and maturity). Varieties in cluster 2 flowered later than those in cluster 5. Short cycle varieties were dominantly associated to cluster 6 but were also found in cluster 5 if they also expressed a high LA. Finally, cluster 7 was characterized by an LMR under 43%.

Relying on MDG, plant N content at maturity was the second most important trait to discriminate the clusters. Phyllochron, SNC, SLA, and LL were not selected in the tree and had very low mean decrease in Gini value. LAR was not used in the trees, but had a higher mean decrease in Gini value than LA and LMR. Yet, LA and LMR captured more variability in the PCA than LAR. Finally, cover rate did not appear on the tree despite its high mean decrease in Gini value, probably due to its negative correlation with DDR.

DISCUSSION

Traits Related to Resource Acquisition Are More Plastic Across Growing Conditions and Varieties Than Architectural Traits

Through functional traits, researchers usually seek to make generalizable predictions across organizational and spatial scales [START_REF] Adler | Traitbased tests of coexistence mechanisms[END_REF] and therefore usually use data collected in diverse environment and management contexts. On the contrary, field measurements sharing similar environment as in the present study deliver a local picture of functional diversity which variations can be solely attributed to genotypic determinants. Several studies have demonstrated that despite the significant plasticity of most plant traits across environments, the classification of species for a given trait remains fairly consistent in space and time [START_REF] Garnier | Consistency of species ranking based on functional leaf traits[END_REF][START_REF] Kazakou | Are trait-based species rankings consistent across data sets and spatial scales?[END_REF][START_REF] Tribouillois | A functional characterisation of a wide range of cover crop species: growth and nitrogen acquisition rates, leaf traits and ecological strategies[END_REF]. Yet, other studies have reported a significant effect of plasticity on species ranking across environments, particularly with respect to traits associated with light acquisition and nitrogen nutrition [START_REF] Lipowsky | Plasticity of functional traits of forb species in response to biodiversity[END_REF][START_REF] Siebenkäs | Phenotypic plasticity to light and nutrient availability alters functional trait ranking across eight perennial grassland species[END_REF][START_REF] Keenan | Global leaf trait estimates biased due to plasticity in the shade[END_REF]. These inconsistencies may impair the use of simple aggregation functions (means, medians, etc.) for generalized trait-based approaches at the species level. 2 for abbreviations of traits and properties. Soil cover stands for maximal soil cover.

In the present study, we compared common field measurements to the meta-analysis by [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] to analyze how environmental and management conditions could alter trait values and distort their relationships. When using the experimental dataset, where all data were acquired in a common environment, species rankings on the PC axis summarizing the variation of architectural traits were fairly consistent with literature ranking. However, rankings on the axis summarizing SLA and LNC variation diverged between the two datasets. This contrast mainly concerned warm-season species, particularly Vigna sp. for which divergence between the two datasets was also noted in terms of phenology. Number of days required to reach maturity was generally lower in literature than in experimental data. This comparison between local measurements and literature averages suggests that some traits are more plastic than others and therefore contribute differently to frame legumes diversity depending on the conditions of observation.

Intraspecific Variation Contributed Only Marginally to Global Variability

Results from the PCA showed that the total trait variance was mainly explained by the variance between species. Variation in intraspecific traits is generally considered of marginal importance compared to differences in interspecific traits [START_REF] Mcgill | Rebuilding community ecology from functional traits[END_REF][START_REF] Hulshof | Variation in leaf functional trait values within and across individuals and species: an example from a Costa Rican dry forest[END_REF]. However, several recent studies have shown that intraspecific variation in traits accounts for a substantial proportion of overall plant community variability [START_REF] Jung | Intraspecific variability and trait-based community assembly[END_REF][START_REF] Albert | When and how should intraspecific variability be considered in trait-based plant ecology?[END_REF][START_REF] Roscher | Contrasting effects of intraspecific trait variation on trait-based niches and performance of legumes in plant mixtures[END_REF]FIGURE 3 | Hierarchical clustering analysis (HCA) of the 30 pulse varieties based on mean value of agroecosystem properties (n = 3). The names of varieties belonging to the same species are written in the same color. Warm and cool season species are depicted with yellow and green bars, respectively. Radars represent the mean (standardized to maximum) value of agroecosystem properties of each cluster. Colors in radar plots correspond to branch color in the HCA. [START_REF] Siefert | A global meta-analysis of the relative extent of intraspecific trait variation in plant communities[END_REF]. The small number of varieties considered in this study may have not been able to highlight such diversity.

Leaf Morphological Traits Confirmed as Variability Drivers Among Pulses [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF] found that seed and leaf traits exhibited the greatest variability among 43 pulse species. In the present study, leaf traits (LA and LL) remained an important variability driver that was correlated to the first PCA axis, whereas seed size, associated to the first axis in the literature analysis, poorly contributed to the total variance across varieties (associated only with the fourth PCA axis) and was decoupled from leaf size, which is consistent with LHS scheme [START_REF] Westoby | A leaf-height-seed (LHS) plant ecology strategy scheme[END_REF].

The experimental data confirmed that LL is an important determinant of legume diversity. This trait was strongly related to LNC and to a lesser extent SLA. This finding can be related to the well know leaf economic spectrum (LES; [START_REF] Wright | The worldwide leaf economics spectrum[END_REF]. The LES reflects a trade-off between a leaf 's lifespan and its maximum net photosynthetic rate. Yet, relations between LNC and photosynthesis in grain legumes have been questioned [START_REF] Adams | Legumes are different: leaf nitrogen, photosynthesis, and water use efficiency[END_REF]. Thus, the fact that this axis reflects a gradient in leaf photosynthetic capacity is uncertain. However, a high LNC could have other biological implications such as remobilization to the reproductive organs [START_REF] Davies | Remobilisation of carbon and nitrogen supports seed filling in chickpea subjected to water deficit[END_REF][START_REF] Schiltz | Dynamics of exogenous nitrogen partitioning and nitrogen remobilization from vegetative organs in pea revealed by 15 N in vivo labeling throughout seed filling[END_REF]. Indeed, SNC was correlated with LNC in our experiment. Small leaflet associated with a high LNC could support the growth of N-rich legume seeds through remobilization. 3). Root nodes represent single input variables (traits) and related split point used to make the prediction; the best discriminant is chosen at each step of the classification procedure. Leaf nodes contain mean, number, and percentage of observations of the predicted variable (cluster membership). (B) Gini decrease value for each trait obtained through a random forest procedure. A high Gini decrease value means the variable considered allowed to gain purity in the process of classification for most of the nodes in the CART it was associated with during the random forest procedure.

Yield Is Inversely Related With Ndfa and Soil Cover

Nitrogen fixation by legumes is generally expected to contribute to improving soil fertility through the mineralization of crop residues such as senescent leaves or stems, or roots. In the current study, Ndfa at flowering and Ndfa at harvest in grain tissue varied almost independently. Clustering reveled three types of "high fixing species": cluster 7, cluster 3, and cluster 5. On the other hand, cluster 1 or 6 expressed a low N fixation capacity. "High yielding" species were clustered in cluster 1 and, to a lesser extent, cluster 4 both concerned by low fixation. Cluster 5 fixed a lot of N during flowering and produced a decent yield. However, yield in cluster 5 was lower than that of clusters 1 or 4. Conversely, the very low yielding species (cluster 3) expressed a significant increase in relative Ndfa between flowering and maturity. Other multiple species comparisons confirmed the negative relation between BNF and seed yield [START_REF] Liu | A side-by-side comparison of biological nitrogen fixation and yield of four legume crops[END_REF]. Some studies reported BNF over a certain level results in a reduction of yield [START_REF] Mesfin | Mineral fertilizer demand for optimum biological nitrogen fixation and yield potentials of legumes in Northern Ethiopia[END_REF].

The multi-faceted benefits of cover crops are well known [START_REF] Blanco-Canqui | Cover crops and ecosystem services: insights from studies in temperate soils[END_REF], but land cover is often not considered for cash crops and the interaction between cover capacity and yield is rarely discussed. Maximum cover and cover duration appeared to be related to biomass at flowering and LAI as well as to yield properties. Cluster 7 expressed an important cover at flowering as a result of an important biomass and a large LA. However, cluster 1 and 3 had a low biomass at flowering and a relatively high maximum cover that lasts long, probably because of a delayed senescence. Except for cluster 1, most of those species had low yield. Delayed senescence, by prolonging photosynthesis, is usually expected to contribute to increasing the yield of highcarbon crops such as cereals [START_REF] Gregersen | Plant senescence and crop productivity[END_REF]. However, studies suggest that the stay-green trait may have limited value for high N species [START_REF] Ismail | Delayed-leaf-senescence and heat-tolerance traits mainly are independently expressed in cowpea[END_REF][START_REF] Kumudini | Trials and tribulations: a review of the role of assimilate supply in soybean genetic yield improvement[END_REF] as its value lies in an improvement of remobilization capacity. This corroborates the absence of relation we found between yield and duration of ground cover.

Trait-Based Predictions of Agroecosystem Properties

Understanding the links between functional traits and agroecosystem functioning is fundamental for optimal management of agroecosystem services [START_REF] Zhang | Ecosystem services and dis-services to agriculture[END_REF][START_REF] Litrico | Diversity in plant breeding: a new conceptual framework[END_REF][START_REF] Wood | Functional traits in agriculture: agrobiodiversity and ecosystem services[END_REF]. As noted earlier, leaf characteristics and especially LNC were crucial in describing pulse diversity. Other traits such as LAR, total nitrogen, seed size, and cycle length are also important, which is consistent with the study of [START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF]. The decision tree did not highlight all these traits but discriminated specific functional characteristics for each cluster. Clusters 1 and 4 (high yielding species) were characterized by low LNC values which is consistent with their low N fixation. Cluster 7 (high biomass at flowering) was characterized by a large LA but a low biomass investment in leaves. Clusters 6 and 5 (medium maximum soil cover and low value for every other properties) were characterized by a short cycle, with cluster 5 (high nitrogen fixation) having a higher LA. Clusters 2 (medium value for every properties) and 3 (high N fixation at maturity and long cover duration) were characterized by a long cycle and high nitrogen content at maturity or flowering, respectively.

To conclude on the highlighted linkages, a high LNC value does not translate into higher yields but is generally associated with a high N fixation. Investment in LA is important for N 2 fixation at flowering but is associated with rapid senescence and low yield. Biomass investment in leaves and a short cycle coincide with low values for almost every property except maximum soil cover. Finally, long cycle and high N content at maturity are related to delayed senescence and maintenance of N 2 fixation along pod filling.

CONCLUSION

Integrating pulses in resilient and sustainable production systems requires a better understanding of their functional diversity and how this functional diversity scales up to provide specific sets of agroecosystem services. Previous studies demonstrated the usefulness of trait value aggregation at the species-level to implement trait-based approaches for a better knowledge of the plant-agroecosystem functioning relationships and increasing our predictive capacity [START_REF] Lavorel | Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail[END_REF][START_REF] Violle | Let the concept of trait be functional![END_REF][START_REF] Garnier | A trait-based approach to comparative functional plant ecology: concepts, methods and applications for agroecology[END_REF][START_REF] Kazakou | Are trait-based species rankings consistent across data sets and spatial scales?[END_REF][START_REF] Martin | Plant functional traits in agroecosystems: a blueprint for research[END_REF][START_REF] Wood | Functional traits in agriculture: agrobiodiversity and ecosystem services[END_REF][START_REF] Guiguitant | Evaluation of pulse crops' functional diversity supporting food production[END_REF]. Our experimental results highlighted trade-offs and synergies between agroecosystem properties supported by different pulses species. We were able to demonstrate the relevance of LNC, LAR, and phenology to predict most of the trade-offs observed for agroecosystem properties in those conditions.

However, the present study advocates further for a systematic characterization of the functional diversity (and its plasticity), including intra-specific variability. We argue that the trait-based approach may be easier to implement in crop species where varieties origin and genetics are often better described and accessible, compared with the huge genetic variation encountered in populations of spontaneous species. Application of such trait-based approach would rapidly benefit the selection of species/varieties for specific targeted agroecosystem services provisioning under specific (environmental or management) conditions. of the manuscript. DV and HM contributed substantially to the manuscript through extensive revision. HM designed the experiment. DV participated to data analyses.

FIGURE 1 |

 1 FIGURE 1 | Principal component analysis (PCA) performed on nine functional traits collected on 43 pulses species. (A,C) Individual species projection on the first three axes of the PCA; light gray points represent the 43 species from literature (Guiguitant et al., 2020) and colored points are from the current experimental study projected as supplementary individuals. Transparent points are the data of the four plots for the three varieties of the 10 species, the larger opaque points represent the barycenter of each species, and the horizontal and vertical bars represent the standard deviation along axes. Visualization of the variables and correlation circles on (B) PC1-PC2 and (D) PC2-PC3 planes. See Table2for abbreviations of traits.

FIGURE 2 |

 2 FIGURE 2 | Principal component analysis (PCA) performed on traits measured on three varieties of 10 species. Visualization of the variables and correlation circles on (A,C) PC1-PC2 and (B,D) PC3-PC4 planes. Gray scale represents the contribution of each variable to axes construction. Orange variables are agroecosystem properties, used as supplementary variables. See Table2for abbreviations of traits and properties. Soil cover stands for maximal soil cover.

FIGURE 4 |

 4 FIGURE 4 | (A)Classification and regression tree (CART) regression for prediction of membership to one of the seven clusters based a hierarchical clustering analysis of the 30 pulse varieties based on agroecosystem properties (see Figure3). Root nodes represent single input variables (traits) and related split point used to make the prediction; the best discriminant is chosen at each step of the classification procedure. Leaf nodes contain mean, number, and percentage of observations of the predicted variable (cluster membership). (B) Gini decrease value for each trait obtained through a random forest procedure. A high Gini decrease value means the variable considered allowed to gain purity in the process of classification for most of the nodes in the CART it was associated with during the random forest procedure.

TABLE 1 |

 1 List of the 10 species included in the experiment with a description of the three varieties associated to each species, the abbreviation (Abb.) used in the study, sources of the seeds, and the ß value used for %Ndfa computation.

	Species	Varieties	Abb.	Source	β ( )
	Vigna mungo	ISRA 58-77 (V1), ISRA 66-41 (V2), ISRA 67-30 (V3)	Vmu	ISRA, Sénégal	-1,75
	Vigna unguiculata	YACINE (V1), SAM (V2), IT98K-1092-1-1 (V3)	Vun	ISRA, Sénégal	-1,61
	Vigna aconitifolia	SH_VM (V1), SH_VA (V2), Tvu 45 35 (V2)	Vac	ISRA, Sénégal	-0,91
	Lotus tetragonolobus	Landrace	Lte	Palmbeach medicinal herb, rareexoticseed, hobbyseed	-0,12
	Phaseolus vulgaris	Red kidney (V1), zorro (V2), michelet (V3)	Pvu	Epi de Gascogne, France	-2,16
	Cicer arietinum	Elmo (V1), billy bean (V2), orion (V3)	Car	Epi de Gascogne, France; USDA Pullman, WS	-1,75
	Lens culinaris	Belezana (V1), Richela (V2), Anicia (V3)	Lcu	Epi de Gascogne, USDA Pullman, WS	-0,56
	Trigonnella foenum-graecum	FENUSOL, landrace	Tfo	Epi de Gascogne, France; ICARDA, Morroco	-0,42
	Lathyrus sativus	Wassie (V1), Almenaza (V2), Chicharo (V3)	Lsa	Epi de Gascogne, France	-0,38
	Vicia faba	Aguadulce (V1), alfia (V2), lobaba (V3)	Vfa	ICARDA, Morroco	-0,5

TABLE 2 |

 2 List of agroecosystem properties and traits measured in the experiment with their respective abbreviation and unit.

	Variable	Abbreviation	Unit
	Agroecosystem property
	Grain yield	GY	t ha -1
	Biomass yield	BY	t ha -1
	Leaf area index	LAI
	% Nitrogen derived from	Ndfa_MAT	%
	atmosphere at maturity
	% Nitrogen derived from	Ndfa_FLO	%
	atmosphere at flowering
	Water use efficiency	WUE	kg ha -1 mm -1
	Biomass at flowering	B_FLO	t ha -1
	Maximal soil cover	Maximum soil cover	%
	Soil cover duration	Soil cover duration	Thermal unit
	Traits
	Leaf area	LA	mm 2
	Specific leaf area	SLA	cm 2 g -1
	Leaflet length	LL	mm
	Leaflet width	LW	mm
	Leaf nitrogen content	LNC	mg g -1
	Plant N content at flowering	Nplant_FLO	g plant -1
	Plant N content at maturity	Nplant_MAT	g plant -1
	Phyllochron	Phyllochron	Thermal unit
	Soil cover rate	Cover rate
	Degree days to maturity	MAT	Thermal unit
	Degree days to flowering	FLO	Thermal unit
	Thousand seed weighed	TSW	G
	Leaf area ratio	LAR	-
	Leaf mass ratio	LMR	-
	Stem mass ratio	SMR	-
	Harvest index	HI	-
	Seed nitrogen content	SNC	mg g -1

TABLE 3 |

 3 Euclidian distance between literature and current study coordinates of the barycenter of each species on the firsts three PCA axes.

	Species	PC1	PC2	PC3	Three PCs
	Cicer arietinum	1.1	0.39	1.25	1.68
	Lathyrus sativus	1.0	0.18	0.11	1.03
	Lens culinaris	0.4	0.25	1.29	1.39
	Lotus tetragonolobus	0.9	1.41	1.74	2.40
	Phaseolus vulgaris	0.4	0.92	3.74	3.87
	Trigonella foenum-graecum	0.8	0.44	1.73	1.97
	Vicia faba	1.1	1.20	0.22	1.64
	Vigna aconitifolia	1.0	2.61	2.96	4.08
	Vigna mungo	0.6	1.21	1.76	2.23
	Vigna unguiculata	0.7	1.23	2.39	2.77

Frontiers in Plant Science | www.frontiersin.org

March 2021 | Volume 12 | Article 636915

ACKNOWLEDGMENTS

The authors thank Laurent Palau and UCOSEM team in Lectoure for having the field experiment on their farm and ensuring proper crop husbandry, and for providing genetic resources. The authors are also grateful to the LSTM lab (B. Brunel) and the USDA ARS winter legume breeding program at Pullman (WS) (R. McGuee) for providing, respectively, rhizobial inoculum as well as plant genetic resources. The authors also thank ISRA Senegal for sharing Vigna accessions.

FUNDING

This work was funded by the European Programme ARIMNet/ERA-NET, FP7-KBBE, session 2017-2020, SEMIARID, grant agreement no. 618127.

DATA AVAILABILITY STATEMENT

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS

All authors conceived the ideas and approach together. JG conducted the experiments and data analyses and led the writing

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Yingning He

Douglas Silva De Vasconcellos

Véronique Bardinal

David Bourrier

Eric Imbernon

Ludovic Salvagnac

Adrian Laborde

Xavier Dollat

Thierry Leichlé
email: tleichle@laas.fr

Lateral Porous Silicon Interferometric Transducer for Sensing Applications

Keywords: porous silicon, interferometry, lab on a chip, biosensor I

In this work, we demonstrate the transducing ability of lateral porous silicon membranes (LPSi) using optical interferometry. To this aim, we carry out Fourier-transform infrared spectroscopy measurements using a microscope equipped with an appropriate objective in order to overcome the difficulty to obtain interference signal from LPSi membranes with small dimensions. Reflectance spectra are recorded while filling the membrane with various solvents and their analysis provide estimations of the effective optical thicknesses and the resulting index of refraction of the fillers in the NIR range. The results show that the various solvents can appropriately be identified and discriminated through the derived refractive indices, thus providing a proof-of-concept of the LPSi interferometric transducer.

INTRODUCTION

Sample analysis usually involves the detection of a specific biochemical compound. This is achieved by transducing a biological recognition event into a measurable and quantifiable signal. Common signal transducers include electrical, electrochemical, optical, mechanical, and thermal sensors. The optical sensing scheme enables to probe surfaces and films in a nondestructive manner; it offers advantages in speed, sensitivity, and robustness, as well as permitting in situ sensing and real-time measurements. Out of the various materials available to constitute optical-based biosensors, porous silicon has high surface area, convenient surface chemistry, and has been widely studied in interferometric sensing configuration [START_REF] Harraz | Porous silicon chemical sensors and biosensors: A review[END_REF].

Porous silicon-based interferometric sensors and biosensors consist of vertical porous silicon layers created into a silicon wafer where the pores are oriented perpendicular to the wafer surface [START_REF] Dhanekar | Porous silicon biosensor: current status[END_REF]. Vertical porous silicon is most often created by silicon anodization. For biosensing purpose, the surface of the porous silicon is modified with biomolecules. Reflection of white light at the top and bottom of the porous silicon layer results in an interference pattern (Fabry-Pérot fringes) that depends on the refractive index of the porous silicon matrix. Interactions of molecular species with their recognition partners immobilized onto the porous silicon surface induce a change in the refractive index of the porous layer, giving rise to a wavelength shift in the fringe pattern.

It is possible to carry out biosensing using porous silicon via optical measurement through a fluidic cell, as demonstrated by works where vertical porous silicon is integrated into microfluidics for real-time detection [START_REF] Janshoff | Macroporous p-type silicon Fabry-Perot layers. Fabrication, characterization, and applications in biosensing[END_REF]. These structures have been vastly studied for biosensing of various elements, e.g. DNA, protein, enzyme activity and bacteria [START_REF] Canham | Porous Silicon Optical Biosensors[END_REF].

In a previous work [START_REF] He | Fabrication of lateral porous silicon membranes for planar microfluidics by means of ion implantation[END_REF], we have developed a unique process to fabricate porous silicon membranes with lateral pores, i.e. with horizontal as opposed to above described vertical pores, and integrate them into planar fluidic channels. While the use of porous layers with vertical pores for biosensing in two-dimensional microfluidics results in flow-over (FO) sensing configurations, porous silicon membranes with lateral pores enable liquid to flow through (FT) (Figure 1). Theoretical works have shown that, with similar sensor footprints, flow-through (FT) sensor offers up to 20-fold improvement in response time [START_REF] Escobedo | Flow-through vs flow-over: analysis of transport and binding in nanohole array plasmonic biosensors[END_REF]. Besides, the experimental comparison between FO and FT sensors, later investigated using porous silicon microcavities and freestanding membranes [START_REF] Kumar | Investigation of non-specific signals in nanoporous flow-through and flow-over based sensors[END_REF], revealed some issues encountered with FO sensors where analytes or contaminants can be trapped within the closed-ended pores. Thus, following our previous results on fabricated lateral porous silicon (LPSi) membranes, the aim of the present work is to demonstrate the sensing ability of lateral porous silicon via interferometric measurements to offer solutions for future on-chip flowthrough biosensors. Optical measurements were performed on a VERTEX 70 Fourier Transform Infra-Red spectrometer (FTIR) (Bruker Optics) equipped with a tungsten light source, a Quartz beam splitter (T502) and a Si-diode detector (SiD 510), which were used to cover a spectral range from 600 to 1100 nm. The spectrometer was connected to a HYPERION microscope stage holding various objectives and apertures. A 36× objective was used for localized measurements.

Microfluidic chips bearing lateral porous silicon membranes were fabricated following the process described in [START_REF] He | Fabrication of lateral porous silicon membranes for planar microfluidics by means of ion implantation[END_REF]. Each 1.6 cm×1.6 cm microfluidic chip integrated two 5 µm deep microchannels separated by a 10 µm thick and ~4 µm high porous silicon membrane with average porous diameter of ~15 nm. The porous membrane was obtained through silicon anodization at 200 mA cm -2 during 1 min. The porous silicon chip was loaded into a sample holder which was subsequently fixed on the microscope stage (Figure 2). The sample holder was connected to a 4-channel reservoir (FLUIWELLTM, Fluigent) with pressure sources (MFCS-8C, Fluigent) to control the fluid flow inside the chip.

Various solvents loaded in the reservoir were pumped into the microchannels and subsequently pushed through the LPSi membrane via positive pressure to ensure the full infiltration of the analyte into the membrane. The observation window was downsized to a rectangular area with dimensions of 10 × 100 µm 2 using adjustable x-and yapertures. This window size was chosen so as to only observe the 10 × 250 µm 2 LPSi membrane. Reflectance data acquisition was carried out with a spectral resolution of 4 cm -1 . In order to test the interferometric transduction capability of the LPSi membranes, three kinds of solvents (water, ethanol, acetone) were tested. It is important to note that, prior to the measurements, a reference spectrum was taken onto a gold mirror using the same experimental conditions (i.e. sample stage height and observation window size). Besides, before injecting test solutions, the chip was exposed to oxygen plasma in order to increase the wettability of the microchannels and the nanopores.

III. RESULTS AND DISCUSSIONS

A. Validation of Optical Transduction

The relationship between the optical path length 2nLcosθ, also referred to as the effective optical thickness (EOT), and the Fabry-Pérot fringes [START_REF] Iizuka | Elements of Photonics[END_REF] is:

(1)
Where, n and L are the refractive index and the thickness of the PSi layer, θ is the incident angle of the light at the PSi/Si interface, m is an integer corresponding to the spectral order of the fringe and λ max is the wavelength of the fringe maximum.

Since the refractive index of porous silicon membranes depends on the pore-filling medium, we have tested the ability of our platform to detect a change of refractive index through the effective optical thickness by loading the membranes with liquids of various optical properties. Evaporable solvents (water, ethanol and acetone) were tested according to the procedure described in section II. Figure 3 presents the resulting reflectance spectra. Compared to empty and dried porous silicon, the presence of all three liquids clearly led to a red shift of the spectrum corresponding to an increase in refractive index. The wavelength shifts from bare porous silicon are in the following order: water (n~1.328) < acetone ≈ ethanol (n~1.35).

B. Analysis of Reflectance Spectra

We have analyzed the Fabry-Pérot fringes by means of two methods in order to derive the measured effective optical thickness.

First, we have plotted the spectral order of each fringe (peak number) as a function of frequency, as shown in Figure 4, where the slope of the linear fit equals the effective optical thickness 2nLcosθ [START_REF] Sailor | Porous silicon in practice: preparation, characterization and applications[END_REF]. Since the membrane thickness (L) can be measured by SEM and θ is given by the numerical aperture of the objective, we can derive an experimental value for the refractive index. The spectrum analysis was also carried out by Reflectometric Interference Fourier Transform Spectroscopy (RIFTS) [START_REF] Sailor | Porous silicon in practice: preparation, characterization and applications[END_REF]. This method computes the frequency spectrum of an input waveform by fast Fourier transform, yielding a peak whose position along the x-axis corresponds to the effective optical thickness.

The theoretical refractive index of the porous silicon layer can be estimated by Bruggeman Effective Medium Model:

(2

)
Where P is the porosity of porous silicon (40% in the case of the fabricated LPSi membranes), n fill is the refractive index of the medium filling the pores, n skeleton is the refractive index of the porous material (i.e. the refractive index of silicon, 3.673 [START_REF] Aspnes | Dielectric functions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV[END_REF]) and n layer is the refractive index of the composite porous silicon. 1 provides the refractive indices derived from the experimental methods and from theory. Since the dispersion of the index of refraction is fairly small in the wavelength range considered in this study [START_REF] Sailor | Porous silicon in practice: preparation, characterization and applications[END_REF], it was assumed constant for all solvents and the index of porous silicon as a function of filling solvent was obtained using values at 830 nm and 20 °C [START_REF] Rheims | Refractive-index measurements in the near-IR using an Abbe refractometer[END_REF]. Data show that values obtained from the experimental methods are in accordance with the theoretical ones. The composite refractive index (n layer ) is in order of air < water < acetone ≈ ethanol, which, indeed, corresponds to the shifts of measured interference patterns, thus demonstrating the possibility to discriminate the various solvents with lateral porous silicon.

IV. CONCLUSION

In this work, we have demonstrated the Fabry-Pérot interferometric transduction capability of lateral porous silicon membranes bridging microfluidic channels. The specificity and difficulty of our study is that the LPSi membranes exhibit small dimensions, which increases the difficulty to obtain exploitable interference signals. Nevertheless, this issue was overcome by optimizing the experimental setup and the optical measurement conditions. We were able to identify the presence of solvents with different refractive index, thus providing a proof-of-concept of the LPSi interferometric transducer.

Fig. 1 .

 1 Fig. 1. Two interferometric sensing configurations with porous silicon: (a) a typical flow-over configuration (FO) with a surface-based sensor within a microfluidic channel, (b) a flow-through configuration (FT) with an array of lateral nanopore tubes adequately connecting both microchannels.

Fig. 2 .

 2 Fig. 2. Experimental setup for microscale analysis and fluid-managing system used for optical detection.

Fig. 3 .

 3 Fig. 3. Experimental reflectance spectra of a LPSi membrane filled with different solvents: water, acetone, and ethanol. The spectrum marked as air means that no liquid fills the pores.

Fig. 4 .

 4 Fig. 4. Fringe number as a function of corresponding spatial frequency. The inset shows the spectrum replotted in frequency unit and the peaks of the interference pattern corresponding to successive peak number.

TABLE I .

 I EXPERIMENTAL AND THEORETICAL INDICES OF REFRACTION OF LATERAL POROUS SILICON MEMBRANES FILLED WITH VARIOUS SOVENTS

Table

ACKNOWLEDGMENT

The authors acknowledge the Agence Nationale de la Recherche (ANR-17-CE09-0024-01). This work was partly supported by the French RENATECH network.


Anna Carla Araujo

B G Jun

R E Devor. S G Kapoor

R E Devor

M B G Jun

Experimental Investigation of a Combined Drilling and Thread Milling Process: Thrilling

Keywords: Drilling, thread milling, thrilling, cutting force, thread quality

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

INTRODUCTION

Threads on a workpiece can be produced in a variety of ways based on two basic principles: plastic working and metal cutting. Threads made by plastic working are stronger but lack the accuracy required for many applications. Threads made by cutting on the other hand provide better accuracy and finish, and threads can be made even with brittle materials. There are two common cutting processes for internal thread generation, which are cut tapping and thread milling [START_REF] Stephenson | Metal cutting theory and practice[END_REF]. Recently, a thread cutting tool called thriller has been developed by Emuge Corp., which can perform both drilling and thread cutting operations. Thus, after drilling a hole, no tool change is required to make threads, leading rapid thread generation. This process needs to be understood for effective application of the process.

The thrilling process combines the drilling and thread milling processes. Thus, the bottom of a thriller is similar to a drill, and along the helix of the tool has the tool geometry of a thread milling tool. Since one tool performs two different processes consecutively, the quality of the threads generated would depend on both drilling and thread milling operations. However, though thread making processes both based on plastic working [START_REF] Cao | Investigation of thread tapping load characteristics through mechanistics modeling and experimentation[END_REF][START_REF] Ivanina | Influence of parameters of the cutting part of taps on threading accuracy[END_REF][START_REF] Warrington | Experimental investigation of thread formation in form tapping[END_REF] and thread cutting [START_REF] Araujo | A model for thread milling cutting forces[END_REF][START_REF] Smith | Advanced machining: The handbook of cutting technology[END_REF]] have been extensively studied by many researchers, the thrilling has not been investigated at all. Also, very little is known about cutting characteristics of the thrilling process.

The objective this paper is to investigate the thrilling process. Cutting forces during thrilling are measured and related to the cutting characteristics of the process. The quality of the threads generated at different drilling and thread milling conditions is examined to understand the relationship between drilling and thread milling operations as one affects the other. Two different workpiece materials are considered to investigate the effect of workpiece materials.

The paper is organized as follows. The unique characteristics of the thrilling process are first described, followed by examination of the tool geometry. Then, the cutting characteristics during both drilling and thread milling are presented, followed by study on the quality of the threads generated.

DESCRIPTION OF A THRILLING PROCESS

The machining process of thrilling can be divided into seven stages as shown in Fig. 1. The first stage involves positioning the tool above the workpiece where the hole with threads needs to be made. Once the tool is positioned, it drills down into the workpiece to the desired depth of the hole. Then, the tool retracts up one third of the pitch. In order to make the thread, while rotating at the prescribed rotational speed and retracting the tool half the pitch, the tool moves radially into the workpiece to a proper radial depth of cut through a helical "run-in loop". Once the cutting edge is engaged into the workpiece at the desired radial depth of cut, the tool is moved up helically for one complete revolution. The number of threads generated in this one revolution corresponds to the number threads engaged in the workpiece. If generation of the threads for the full depth of the hole is achieved in several steps, the stages 4-6 are repeated until all the threads are generated. Then, the tool moves back to the center of the hole, and, finally, the tool is retracted out of the hole.

During the seven stages, major cutting occurs at stages 2 (drilling operation) and 6 (thread milling operation). The drilling operation affects the quality of the hole, which in turn affects the thread milling process. Since both of these operations influence the whole process, their interaction may play a role in determining the final thread quality.

THRILLING TOOL GEOMETRY

The geometry of a thriller is shown in Fig. 2. Each flute has a chisel edge, drilling cutting lip, milling cutting edge, and thread milling cutting edge. The chisel edge does not cut but only spread the material sideways by an indentation mechanism [START_REF] Altintas | Manufacturing automation: Metal cutting mechanics, machine tool vibrations, and CNC design[END_REF]. It has a width of 2w and an edge angle of ψ c . The cutting lips have an offset from the center of the tool due to the chisel edge. The cutting lips remove the material with a constant chip thickness and the tool is fed into the workpiece at a feedrate of f d . The chisel edge and cutting lips have the shape of a drill. The cutting lips at the point angle κ and with the height h d meet with the milling cutting edge. The height of the milling cutting edge is denoted as h e . The cutting tips, though, do not interfere with the thread milling operation. The milling cutting edge separates the drilling part from the thread milling part. In the drilling stage, it does not cut but only acts as a helical flute as in drilling during the drilling operation. Then, the milling cutting edge cuts as a milling tool so that the cutting lips do not interfere with the thread milling operation and allow the thread milling operation. During the thread milling operation, the cutting forces generated include the forces due to thread making and also cutting by this milling cutting edge.

The external diameter of the thread milling cutting edge is smaller than the milling cutting edge diameter.

The threads on a thriller have the maximum diameter on the outside d o and minimum diameter d i , which can be obtained as

d i = d o - 2h t
, where h t is the thread height. The thread pitch is denoted as p. Figure 3 shows an example of the thread geometry of an M10 thrilling tool. The dashed line shows the height of a M10 thread. Since the drilled hole needs to be greater than the thriller thread diameter, the milling cutting edge height or radius is greater than the thread tooth height or radius. To make an M10 thread in the workpiece, the distance to be moved into the workpiece during the thread milling operation is equal to the thread height plus the difference between heights of the milling cutting edge and thread tooth. For the thriller in Fig. 3, this distance is equal to 0.9 mm. The experiments were conducted with an Emuge thriller of size M10x1.5 standard threads. Two different workpieces (Brass and Al6061) were used for the experiments. Castrol 6519 at 5% dilution was used as cutting fluid. The quality of the threads was analyzed using Mitutoyo Contracer CBH-400. The helix angle was measured using an optical microscope (Hirox CX-10C) and TechDig software (Fig. 5). Table 1 shows the thriller geometry.

THRILLING FORCE EXPERIMENTS Linear Thread Milling Experiments

In order to study the force characteristics of the thrilling process, experiments were first conducted by simplifying the thread milling operation as shown in Fig. 6. By rolling out straight the curvilinear tool-workpiece path in Fig. 6(a) and neglecting the motion in the z-direction, the thread milling operation can be simplified as in Fig. 6(b) such that the tool moves linearly in the x-direction similar to a milling process. This simplification is a good approximation of the process as the z-direction velocity is very small compared to the cutting velocity and the radius of the curvature of the hole is large compared to the chip thickness. This simplification allows more thorough examination of the force profiles and the effects of the number of pitches engaged. The only major difference with this simplification is that the cutting flue exits the workpiece sooner than in actual thread milling operation [START_REF] Araujo | A model for thread milling cutting forces[END_REF]. In order to see the effect of number of threads engaged, linear thread milling experiments were conducted at different axial depth of cuts (2.9 -10.4 mm) with brass workpiece. Initially, experiments were conducted with only the milling cutting edge and no thread engaged. Then, experiments were conducted with increasing number of threads engaged. Figure 7 shows the peak and average forces averaged over 100 revolutions during linear thread milling experiments, plotted against the number of threads engaged. "Zero" thread engaged means that only the milling cutting edge was engaged. It shows that the force in the normal-to-feed direction is the highest. Also, as the number of engaged threads increases, the force in the normal-to-feed direction increases much more than in the feed direction. In actual thrilling operations, the feed direction represents the tangential direction. Thus, increase in the number of engaged threads would result in small increase in the torque, leading to less torque requirements. However, the radial forces would be high, which may result in more tool deflections. Because of the helix angle and small radial depth of cut, though the average force increases linearly with the number of threads engaged, the peak force does not. Figure 8 shows the force profiles with (a) 0, (b) 1, (c) 3, and (d) 5 threads engaged. As shown, though peak force does not increase as much, the area under the profile is larger with increased number of threads engaged, resulting in increased average forces. Also, similar to the results in [START_REF] Araujo | A model for thread milling cutting forces[END_REF], the number of engaged threads cause the force profiles to have "local" peaks. In Fig. 8(d), due to increased number of engaged threads, there is an increased number of these local peaks on the force signature. These local peaks lead to forces with higher exciting frequencies than the tooth passing frequency, which can have a dynamic effect on tool vibrations.

Thrilling Experiments

In order to study the force characteristics during actual thrilling operations, experiments were conducted for two different workpiece materials (brass and Al6061). Table 2 shows the cutting conditions. In order to evaluate the effectiveness of the thrilling process, the thread milling forces were compared to the tapping results given in [START_REF] Chowdhary | Modeling forces including elastic recovery for internal thread forming[END_REF] for the same workpiece materials. Figure 9 shows the torque and thrust during the thread milling operation at the spindle speed of 1000 rpm and the feedrate of 0.02 mm/rev. The forces during the run-in loop, thread milling, and run-out loop are also shown in Fig. 9. An example of the tapping torque and thrust at the spindle speed of 1000 rpm for the same size tool is shown in Fig. 10. Unlike the tapping torque and thrust, thread milling torque and thrust have periodic components, which could cause torsional and axial tool vibrations. Since radial forces are higher than tangential force in thread milling, the forces in the x-and y-directions are quite large also with the periodic components as shown in Fig. 11, which could cause lateral deflections. However, if maximum values of torque and thrust are taken during the thread milling operation, however, comparison of the magnitudes shown in Table 3 show that thread milling torque and thrust are, for most cases, more than ten times smaller than tapping results given in [START_REF] Chowdhary | Modeling forces including elastic recovery for internal thread forming[END_REF]. In tapping, the feedrate is determined by the spindle speed and the thread pitch. However, in thread milling, the feedrate can be any value, which is the reason why the thread milling forces in Table 3 are much lower than tapping forces and can be lower if desired by using lower feedrates. This freedom to select a desired feedrate allows the thrilling or thread milling process more versatile.

THREAD QUALITY EXPERIMENTS

In order to investigate the effect of both drilling and thread milling and their interaction on the thread quality, two levels of spindle speed (1000 and 2000 rpm) and feedrate (0.05 and 0.15 mm/rev for drilling and 0.01 and 0.02 mm/rev for thread milling) for both drilling and thread milling operations are considered. Two different workpiece materials (brass and Al6061) were considered. Table 4 shows the experimental design. The thread quality was evaluated in terms of the thread geometry (pitch and height). In order to measure the thread pitch and height, the block with the threaded hole was cut in half as shown in Fig. 12(a). Then a profilometer (Mitutoyo Contracer CBH-1400) was used to scan along the direction as shown in Fig. 12(b).

The measurements were taken by moving the stylus of the profilometer at a rate of 0.5 mm/s. An example of the profilometer measurement is shown in Fig. 13. The thread pitch and height were measured to be the distances as shown in Fig. 13. A total of four scans were taken for each block to measure the thread pitch and height. Two scans were first taken on one half of the block, and two other scans were then performed on the other half of the block. For each scan, five threads were scanned starting from the bottom of the hole, and the heights and pitches of these five threads were averaged. Then, the average of the four scan results were recorded as the thread geometry for each experiment. Table 5 shows the results of the thread geometry measurements. Table 5 shows that though the thread pitch is quite close to the nominal value (1.5 mm), the thread height is generally smaller than the nominal value (0.75 mm). This may be due to tool deflection caused by the radial or lateral forces during the thread milling operations. In fact, depending on the location of the thread in the hole, the thread height varies. Figure 14 shows the thread heights measured for one of the scans for the Exp. 1 hole. Low thread number on the abscissa indicates the thread near the bottom of the hole. As seen, the thread height increases close to the nominal value as the thread number increases. This indicates the presence of the lateral tool deflection due to radial forces. Note that though forces for Al6061 were generally higher as was shown in Table 3, the thread heights are lower for brass in Table 5. This may be due to the discontinuous nature of the brass chips. It was observed during the experiments that because brass chips are discontinuous, chip evacuation was difficult during the drilling operation and many segments of chips were left in the hole, which may adversely affect the thread milling operation. In order to study the effects of both drilling and thread milling conditions, four effects considered are listed in Table 6 along with the low and high values. The four main effects and their interactions are analyzed. Figure 15 shows the normal probability plot for thread height of brass material. As shown, the effect of the thread milling feedrate on thread height is surely significant. The similar result is observed with Al6061. This is somewhat expected since increase in the thread milling feedrate directly increases the cutting forces. Figure 16 shows the main effects on thread height for brass. It shows clear effects of thread milling feedrate, and increasing the thread milling feedrate decreases the thread height away from the nominal value. For the thread pitch, it seems that the interactions between drilling spindle speed and the thread milling conditions have significant effects (Fig. 17). The fact that the drilling spindle speed has an interaction effect indicates that the pitch geometry depends on the drilled hole conditions. Figure 18 shows the main effects on thread pitch for brass. It shows that the effect of drilling spindle speed is significant and increasing the drilling spindle speed improves the thread pitch dimension. The thread milling spindle speed have a similar effect. This may be softening of the materials during drilling resulting in lower forces and better hole quality.

CONCLUSIONS

The following conclusions can be drawn from the experiments conducted:

1. As the number of engaged threads increases, the force in the feed (tangential) direction increases much less than in the normal-to-feed (radial) direction, leading to less torque requirement. 2. Due to helix angle, though the peak force does not increase as much, the area under the profile is larger with increased number of threads engaged, resulting in increased average forces.

3. Compared to tapping torque and thrust for the same size threads, thread milling torque and thrust are much lower (more than ten times). 4. In tapping, the feedrate is fixed by the thread pitch and the spindle speed. However, in thread milling, any feedrate can be selected, allowing either higher productivity with high feedrates or lower forces with low feedrates. 5. The thread milling feedrate has a significant effect on the thread height due to tool deflections. 6. Drilling spindle speed has a positive effect on the thread pitch. This may be due to softening of materials resulting in lower forces and better hole quality. This indicates that the pitch dimension depends more on the drilled hole conditions.

FIGURE

  FIGURE 1. A THRILLING PROCESS

FIGUREFigure 4

 4 FIGURE 3. THRILLER THREAD GEOMETRY EXPERIMENTAL SETUP Figure 4 shows the experimental setup for thrilling operations. Experiments were conducted on Mori Seiki TV30 CNC milling center. A Kistler

FIGURE

  FIGURE 4. EXPERIMENTAL SETUP

FIGURE 7 .

 7 FIGURE 7. MEASURED FORCES DURING LINEAR THREAD MILLING

FIGURE 8 .

 8 FIGURE 8. FORCE PROFILES DURING LINEAR THREAD MILLING

FIGURE 9 .

 9 FIGURE 9. TYPICAL THREAD MILLING TORQUE AND THRUST

  FIGURE 13. PROFILOMETER MEASUREMENT

FIGURE

  FIGURE 14. PROFILOMETER MEASUREMENT

FIGURE 17 .

 17 FIGURE 17. NORMAL PROBABILITY PLOT OF EFFECTS ON THREAD PITCH FOR BRASS

TABLE 1 .

 1

	THRILLER GEOMETRY
	Geometry	Value
	Point angle κ	140 o
	Cutting lips height h d	1.64 mm
	Milling cutting edge height h e	1.35 mm
	Milling cutting edge diameter d d 8.5 mm
	Helix angle α h	28 o
	Thread pitch p	1.5 mm
	Thread height h t	0.75 mm
	Thread outside diameter d o	8.2 mm
	Thread inside diameter d i	5.2 mm
	FIGURE 5. HELIX ANGLE MEASUREMENT

TABLE 2

 2

		. CUTTING CONDITIONS FOR FORCE
		EXPERIMENTS
		Drilling	Thread Milling
	Exp. #	Speed (rpm)	Feed (mm/rev) (rpm) (mm/rev) Speed Feed
	1	1000	0.05	1000	0.01
	2	1000	0.15	1000	0.02
	3	2000	0.05	2000	0.01
	4	2000	0.15	2000	0.02

TABLE 3 .

 3 COMPARISON OF THREAD MILLING AND TAPPING TORQUE AND THRUST

		Thread Milling
		Brass	Al6061
		Torque Thrust Torque Thrust
	Exp. #	[Ncm]	[N]	[Ncm]	[N]
	1	44.4	60.7	45.4	59.4
	2	52.9	70.3	65.2	86.2
	3	43.3	58.9	42.3	51.1
	4	54.9	69.7	65.4	99.7
			Tapping
	Plug	925	760	1049	766
	Bottom	832	984	913	857

TABLE 4 .

 4

		EXPERIMENTAL DESIGN
		Drilling	Thread Milling
	Exp. #	Speed	Feed Speed	Feed
		(rpm) (mm/rev) (rpm) (mm/rev)
	1	2000	0.15	2000	0.02
	2	1000	0.15	2000	0.02
	3	2000	0.05	2000	0.02
	4	1000	0.05	2000	0.02
	5	2000	0.15	1000	0.02
	6	1000	0.15	1000	0.02
	7	2000	0.05	1000	0.02
	8	1000	0.05	1000	0.02
	9	2000	0.15	2000	0.01
	10	1000	0.15	2000	0.01
	11	2000	0.05	2000	0.01
	12	1000	0.05	2000	0.01
	13	2000	0.15	1000	0.01
	14	1000	0.15	1000	0.01
	15	2000	0.05	1000	0.01
	16	1000	0.05	1000	0.01

TABLE 5

 5

	. THREAD MEASUREMENT RESULTS
		Brass	Al6061
	Exp. #	Height	Pitch Height	Pitch
		(mm)	(mm) (mm)	(mm)
	1	0.691	1.494 0.707	1.511
	2	0.674	1.509 0.738	1.506
	3	0.714	1.509 0.751	1.503
	4	0.736	1.507 0.738	1.504
	5	0.719	1.496 0.750	1.506
	6	0.746	1.525 0.744	1.507
	7	0.679	1.494 0.744	1.507
	8	0.741	1.503 0.734	1.550
	9	0.722	1.508 0.754	1.511

TABLE 6 .

 6

		CONSIDERED EFFECTS
		Variable	Low (1) High (2)
	A	Drill speed	1000	2000
	B	Drill feed	0.05	0.15
	C Thread milling speed 1000	2000
	D	Thread milling feed	0.01	0.02
	FIGURE 15. NORMAL PROBABILITY PLOT OF
	EFFECTS ON THREAD HEIGHT FOR BRASS
	FIGURE 16. MAIN EFFECTS ON THREAD HEIGHT
		FOR BRASS


Anna Carla Araujo

Pedro Manuel Lopes

Calas Pacheco

Marcelo Amorim

DYNAMICAL ANALYSIS OF AN END MILLING PROCESS

Keywords: End milling, force prediction, nonlinear dynamics, nonsmooth systems, chatter, machining

The cutting processes include different manufacturing procedures as milling, turning and drilling. In all of these processes, the machine tool vibration plays an important role concerning the cutting characteristics. Therefore, its correct understanding is essential in order to improve the workpiece surface quality and avoid tool breakage. Chatter and squeal are some undesirable phenomena related to an improper functioning. This article analyses the nonlinear dynamics associated with the end milling tool and its vibration relative to the work piece considering a nonsmooth two-degree of freedom system. End milling process is of concern assuming different machining characteristics. In this work, the cutting force is considered as a composition of contact/non-contact of the tool and the workpiece with stick-slip behavior guided by friction force and the prescribed velocity of the tool holder. Numerical simulations are carried out showing some situations related to proper and improper functioning during the cutting process.

INTRODUCTION

Machining process is associated with a complex dynamics that involves the coupling of different phenomena. During cutting processes, for example, it is common to have temperature variations that induce dramatic changes in expected behaviors. Although all these complexity, some simple models may be used in order to obtain useful information concerning tool and workpiece behaviors as dynamic tool prediction.

The machine tool vibration during the cutting process plays an important role concerning the workpiece surface quality and also the tool durability. Chatter and squeal are some undesirable phenomena related to an improper functioning. During certain cutting conditions, the motions of the workpiece-tool system are characterized by large amplitudes, which are not desirable for obtaining a good surface finish. The undesirable motions, which are often referred to as chatter, can result in wavy surfaces on the workpiece, inaccurate dimensions, and excessive tool wear.

The analysis of machining using a dynamical approach has been done by different research efforts. The first works in this way is due to [START_REF] Tobias | Machine tool vibration[END_REF] , [START_REF] Merritt | Theory of self-excited machine tool chatter[END_REF] and Tlusty [START_REF] Altintas | Manufacturing automation[END_REF] that investigate the chatter on the process. Their research provided instability diagrams to inform the safe parameters to avoid chatter. After them, a new approach has been developed using nonlinear dynamic tools as bifurcation analysis to provide better comprehension of the process (Moon, 1978) .

Another references treated machining operations by dynamical approach. For turning operations, [START_REF] Chandiramani | Dynamics of 2 dof regenerative chatter during turning[END_REF] analyzed the dynamics of cutting considering a two-degree of freedom system (2-dof) and orthogonal cutting modeling to predict chatter. [START_REF] Pratt | Design and modelling for chatter control[END_REF] studied the boring bars for turning, which commonly presents chatter problems and experimentally determined modal properties of the tool. Also dedicated to turning, Kalmar-Nagy, Sthepan and [START_REF] Kalmar-Nagy | Subcritical Hoft bifurcation in the delay equation model for machine tool vibrations[END_REF] showed the existence of a subcritical Hopf bifurcation in the delay-differential equation used to describe the machining equations of motion. The stability of the milling process was investigated by [START_REF] Zhao | Dynamics and stability of milling process[END_REF] , Insperger et al (2003) and [START_REF] Gradisek | On stability prediction for milling[END_REF] by considering single and two degree of freedom systems for different experimental conditions. The mathematical model is represented by delay-differential equations.

Milling process has the peculiarity to have a contact/non-contact behavior due to its geometry and due to tool vibration. This kind of behavior is related to nonsmooth systems that are usually associated with either the friction phenomenon or the discontinuous characteristics as intermittent contacts [START_REF] Savi | Numerical and experimental investigations of the nonlinear dynamics and chaos in non-smooth systems[END_REF] .

This article investigates the milling process representing the end milling tool by a mass-spring-dashpot system and the contact with the workpiece by a nonsmooth contact/non-contact system. The process of cutting is related to a stick-slip behavior that defines whether the chip is being removed from the workpiece. The tool holder displacement is prescribed and the local force that allows the slip motion is related to the workpiece shear stress. Under these assumptions, the equation of motion is represented by a differential equation that is solved employing the Runge-Kutta method. Numerical simulations are carried out showing some situations related to proper and improper functioning during the cutting process.

DYNAMICAL MODEL

The milling process is of concern by assuming that it is a full immersed milling in the x direction. Figure 1(a) shows the general view of the tool and the workpiece, where it is identified the prescribed displacement, X h and the tool tip displacement, X t . This process is modeled by assuming a nonsmooth system that is composed by a primary system that represents the tool and a secondary system, representing the workpiece. The primary system consists of a linear springdashpot-mass oscillator with parameters m, k, c and displacement X t , and presents a gap that separates itself from the secondary system that represents the work-piece.

The workpiece system has weightless slider connected to a linear spring-dashpot system with parameters k s , c s that can present progressive motion when the acting force exceeds F t , here represented by a dry friction, and displacement X p , the tool-work piece contact position.

Similarly to the stick-slip phenomena reported by Marian ( 2001) , the progressive motion X c occurs when the force acting on the slider exceeds the threshold of the dry friction force, considering that there is no cutting fluid. Figure 1(b) presents a schematic picture of the cutting process system related to the milling. As the system may operate in stick-slip phases and contact-non contact situations, its dynamic behavior is nonlinear and its equations should consider each case. The system dynamics may be understood as a three stage motion as represented in Figure 2: non-contact, contact in stick and contact in slip. The non-contact stage may be understood as a two independent systems representing the situation where the tool is not in contact with the workpiece. Figure 2(a) shows that the displacements X t and X h are not connected with X c and X p by the mass.

An auxiliary displacement g is used to locate the work-piece position related to the tool position, this variable remains equal to X t (as in Figure 2a) in non-contact case and in stick phase, as it will be shown on Figure 2. The displacement g changes only when the slider slips simulating the chip removal.

As long as g is greater than X t , the system remains separate. When there is no more gap between X t and X p , the systems are in contact and the dynamic changes, as shown in Figure 2 In stick stage the slider displacement does not change and the system is restricted to a fixed boundary. In slip stage this boundary is free and X c displacement changes. The equations of motion are formulated by considering each kind of response separately.

Non-contact Stage

In non-contact stage, when X t < g, the secondary system has no movement and its displacements does not change Ẋp = Ẋc and the primary system is described by the dynamic equation of the single degree of freedom, as follows:

m( Ẍt -Ẍh ) + c( Ẋt -Ẋh ) + k(X t -X h ) = 0 (1)
The force acting on the tool tip is calculated by:

F t = -k(X t -X h ) -c( Ẋt -Ẋh ) (2)
In order to rewrite the equations 2 for either contact and non-contact situations, the Heaviside function H(.) is included. Note that the damping coefficient when the tool retracts from the workpiece ( Ẋt < 0) is not considered.

F t = -k(X t -X h ) -c( Ẋt -Ẋh ) -H[X t -g] k s (X t -X c ) + H[ Ẋt ]c s Ẋt (3)
The difference between chip displacement X c and work-piece displacement X p when the secondary system is free is constant and it is defined as δ, used to calculate the chip movement X p on slip phase.

Contact -Stick Stage

The chip, represented by the slider, is not removed in stick stage. This case occurs when the friction force F f is less than a maximum friction force F f max supported by the workpiece material, which can be considered as a function of the shear area and the maximum shear stress.

F f max = A s τ s (4)
The friction force is calculated, in stick phase, by:

F f = H[X t -g] k s (X t -X c ) + H[ Ẋt ]c s Ẋt ) (5)

Contact -Slip Stage

During slip stage, when the chip is being removed, the maximum friction force is achieved. The chip position X c follows the workpiece displacement, X c = X p + δ, and g assumes the work-piece position g = X p . The friction force is, in slip phase, the value of maximum friction force F f = F f max .

The force applied in the tool tip in slip phase cannot be calculated by equation ( 4) due to the difference between k s (X t -X c ) + H[ Ẋt ]c s Ẋt ) and F f max . The force F t in slip phase is:

F t = -k(X t -X h ) -c( Ẋt -Ẋh ) -F f max (6)

Equations of motion

Under these assumptions, the system dynamics may be represented by a nonsmooth system as follows:

m( Ẍt -Ẍh ) + k(X t -X h ) + c( Ẋt -Ẋh ) = F (7)
where

F = F f if F f < F f max (8) F f max elseif (9)
and Numerical simulations are carried out in order to evaluate the model capability to describe milling process. The governing equations are solved by considering the forth order Runge-Kutta method with time steps less than 5.10 -4 s. The end milling tool parameters as diameter, run-out, damping [START_REF] Mann | Milling bifurcation from structural asymmetry and nonlinear regeneration[END_REF] ) and system parameters are presented in Table 1. The depth of cut was 10mm and the distance between the tool holder and the workpiece in t = 0 is 1 mm. Moreover, it is assumed that the workpiece material is under plastic behavior, then the damping coefficient for the workpiece is defined as 10 times the damping tool coefficient (c s = 10c). Note that for Ẋt < 0, when the tool mover back from the workpiece, c s vanishes.

F f = H[X t -g] k s (X t -X c ) + H[ Ẋt ]c s Ẋt ) (10)

NUMERICAL SIMULATIONS

The machine tool is represented by the tool holder and it has a prescribed displacement X h defined by the feed velocity v f , the spindle speed w and the run-out distance ρ between the the tool axis and the spindle speed axis.

X h (t) = v f t + ρ sin(wt) (11)
The milling tool is considered as a flexible body and linear elastic behavior is assumed for the end milling tool material and, therefore, the tool parameters m, k and c are calculated by considering solid mechanics principles from the tool geometry and material properties.

The stiffness of the primary system is defined by considering the tool as a circular beam in bending. Under this assumption, there is a relation given by:

k = 3E t πd 4 64L 3 (12)
where tool material Young modulus E t is considered as 200 GPa, d is tool diameter and L is the free tool length.

The workpiece stiffness is evaluated by assuming that the workpiece is a compression beam, therefore:

k s = E c dt 10L ( 13
)
where workpiece material Young modulus E c is considered as 200 GPa and t is the depth of cut.

The feed per tooth for a two flute milling tool is written in equation 13 as z is the number of teeth.

f t = v f wz = v f 2w (14)
The dry fiction is evaluated from the shear strength calculated by the shear area as a function of the depth of cut t and approximated by:

A s = π 0 tf t sinφdφ (15)
Table 1. Tool, workpiece and machine tool properties. The system dynamics is investigated by considering different operational conditions of the milling process. The idea is to represent different aspects related to proper and improper functioning. In order to describe these conditions, different machine tool velocities and workpiece material are used.

Initially, let us consider a condition that represents the ideal functioning during the milling process. The feed velocity and the shear stress limit used for this example was v f = 0.5mm/s and τ s = 100MPa.

Figures 3(a) and 3(b) presents the displacement time history of different parts of the system, namely the tool, the workpiece and the chip. Figure 3 also show the contact force (3c), the tool force (3(d)) and the phase space of the response (3e). It can be noted the moment that the tool holder makes the tool to touch the workpiece between 2 and 3s. After that the contact force varies from zero until 25N, the force that the tool cuts the workpiece. Note that, to achieve this maximum contact force, the force transfered by the tool holder to the tool is approximately 100N.

The uncut chip thickness is not considered as in Martellotti equation [START_REF] Araujo | A model for thread milling cutting forces[END_REF] . The chip thickness is given by the simulation as the displacement steps that the work-piece material let the tool cut. In order to improve the surface quality, the chip steps must be as smaller as possible. Also, the tool cannot touch and leave the work-piece with large displacements, because it causes impacts on the tool.

The increase of system velocities can change the dynamical response of the system and this can be associated with surface quality and operational conditions. In order to show this kind of behavior, different spindle speeds and feed velocities are now considered. Under these conditions, the response is changed. Another important parameter that may cause different dynamical response is related to the material work-piece property, represented by the shear stress limit and its impact on the stick-slip behavior. As when the shear stress limit is achieved, the chip is removed, for a softer material, the velocities can be higher without an impact on chip thickness.

Figure 4(a) and 5(a) presents the time history of the displacements on the system with 0.1mm/s and 2mm/s respectively. It can be noted that for higher feed velocities the chip formation is smoother and the movement described by the tool holder is reflected on the chip movement. Figures 4 and 5 also show the contact force (4(b) and 5(b)), the tool force (4(c) and 5(c)) and the phase space of the response (4(d) and 5(d)).

Figure 6(a) and 7(a) presents the displacement time history on the system with 200M P a and 500M P a, respectively. It can be noted that for higher stress limits the steps are very high and the surface quality is as expected. The contact forces, shown on Figures 6(b This article deals with a dynamical analysis of the milling process. A nonsmooth system is used for the mathematical description of the system dynamics. Numerical simulations show some operational conditions presenting situations related to proper and unproper functioning. The general aspects of the milling process is captured by the proposed model that can be used to identify unproper system response.

  Figure 1. Tool and system modeling

  (b), for stick stage, and Figure 2(c) for slip stage.

Figure 2 .

 2 Figure 2. Contact / Non-contact situations

Figure 3 .

 3 Figure 3. v f = 0.5mm/s, τ s = 100MPa and w = 1800rpm

  Diameter (d) Height (L) Tool Damp. (c) Elastic Mod. (E t ) Run-out (ρ) Spindle Speed (w) Depth of cut (t

  ) and 7(b), are higher for harder materials and the tool forces has an strong transient behavior in Figure7(b) that with softer materials, it is not present. Phase spaces of the response are presented on Figures6(d) and 7(d).

  Figure 4. v f = 0.1 mm/s, τ s = 100MPa and w = 1800rpm

  Figure 5. v f = 2mm/s, τ s = 100MPa and w = 1800rpm

  Figure 6. v f = 0.5mm/s and τ s = 200MPa and w = 1800rpm

  Figure 7. v f = 0.5mm/s, τ s = 500MPa and w = 1800rpm

  Figure 8. v f = 0.5mm/s, τ s = 100MPa and w = 3000rpm

  Figure 9. v f = 0.5mm/s, τ s = 100MPa and w = 6000rpm

20th InternationalCongress of Mechanical Engineering November 15-20, 2009, Gramado, RS, Brazil