Muhammad Waqar Nasir
email: muhammad_waqar.nasir@ensam.eu

Hocine Chalal
email: hocine.chalal@ensam.euc

Farid Abed-Meraim
email: farid.abed-meraim@ensam.eu

Prediction of forming limits for porous materials using voidsize dependent model and bifurcation approach

Keywords: Void size effect, Bifurcation approach, Gurson-type model, Forming limit diagram, Ductile damage, Localized necking

The scientific literature has shown the strong effect of void size on material response. Several yield functions have been developed to incorporate the void size effects in ductile porous materials. Based on the interface stresses of the membrane around a spherical void, a Gursontype yield function, which includes void size effects, is coupled with the bifurcation theory for the prediction of plastic strain localization. The constitutive equations as well as the bifurcationbased localization criterion are implemented into the finite element code ABAQUS/Standard within the framework of large plastic deformations. The resulting numerical tool is applied to the prediction of forming limit diagrams (FLDs) for an aluminum material. The effect of void size on the prediction of FLDs is investigated. It is shown that smaller void sizes lead to an increase in the ductility limits of the material. This effect on the FLDs becomes more significant for high initial porosity, due to the increase of void-matrix interface strength within the material.

Introduction

Sheet metal forming, associated with its wide range of products, is undoubtedly a manufacturing process of great importance involving plastic deformation of sheet metals into a desired shape. For safer design perspective of the manufactured parts, it is desirable that sheet metal forming processes should not give rise to localized deformations, which as precursor, ultimately lead to material failure. Thus, prediction of localized deformation is of key importance in selecting the forming process parameters. Forming limit diagram (FLD), originally proposed by Keeler & Backofen [START_REF] Keeler | Plastic instability and fracture in sheets stretched over rigid punches[END_REF] and Goodwin [START_REF] Goodwin | Application of strain analysis to sheet metal forming problems in the press shop[END_REF], is a very useful and practical approach for the characterization of ductility limits in sheet metal forming processes. FLD is the plot of in-plane limiting strains, for a sheet metal subjected to biaxial stretching, at the onset of localized necking. Traditionally, FLDs are determined experimentally, using deep drawing tests with various sheet specimen geometries (covering various strain-path ratios), and measuring the inplane strains at the occurrence of localized necking. However, such experimental determination of FLDs for sheet metals is both costly and time consuming. As an interesting alternative, the development in recent years of constitutive models and material instability criteria for the prediction of plastic strain localization has enabled determining FLDs through numerical simulations.

In addition to the anisotropic plastic behavior of sheet metals, ductile damage, which induces softening, can also be introduced in the constitutive equations, thus providing a more realistic mechanical response for sheet metals. In the past three decades, two damage theories, in this regard, have been developed: Continuum Damage Mechanics (CDM) and micromechanics-based damage theory. The CDM approach treats the damage variable as a scalar or tensor depending upon the softening behavior and anisotropy, as addressed extensively earlier in the literature (see, e.g., [START_REF] Kachanov | On creep rupture time[END_REF][START_REF] Rabotnov | Creep problems in structural members[END_REF][START_REF] Lemaitre | A Course on Damage Mechanics[END_REF][START_REF] Maire | A new formulation of continuum damage mechanics (CDM) for composite materials[END_REF][START_REF] Hambli | Comparison between Lemaitre and Gurson damage models in crack growth simulation during blanking process[END_REF][START_REF] Brünig | Numerical analysis and elastic-plastic deformation behavior of anisotropically damaged solids[END_REF][START_REF] Menzel | A framework for multiplicative elastoplasticity with kinematic hardening coupled to anisotropic damage[END_REF][START_REF] Besson | Non-linear mechanics of materials[END_REF][START_REF] Bouchard | An enhanced Lemaitre model formulation for materials processing damage computation[END_REF][START_REF] Voyiadjis | Advances in damage mechanics: metals and metal matrix composites[END_REF][START_REF] Doghri | Mechanics of deformable solids: linear, nonlinear, analytical and computational aspects[END_REF][START_REF] Lian | A modified Lemaitre damage model phenomenologically accounting for the Lode angle effect on ductile fracture[END_REF]). As to the micromechanics-based damage approach, the damage variable is taken to describe the void growth in porous materials. The first micromechanics-based damage model, which has been subsequently widely used, was proposed by Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]. Later, the Gurson model has been modified and extended in various domains and applications. For instance, Gologanu et al. [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing nonspherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF] extended the Gurson model for oblate and prolate voids, while Madou and Leblond [START_REF] Madou | A Gurson-type criterion for porous ductile solids containing arbitrary ellipsoidal voids-I: Limit-analysis of some representative cell[END_REF][START_REF] Madou | A Gurson-type criterion for porous ductile solids containing arbitrary ellipsoidal voids-II: Determination of yield criterion parameters[END_REF][START_REF] Madou | Numerical studies of porous ductile materials containing arbitrary ellipsoidal voids-I: Yield surfaces of representative cells[END_REF] and Madou et al. [START_REF] Madou | Numerical studies of porous ductile materials containing arbitrary ellipsoidal voids-II: Evolution of the length and orientation of the void axes[END_REF] extended the model for non-spheroidal voids. In the works of Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF], Monchiet and Bonnet [START_REF] Monchiet | A Gurson-type model accounting for void size effects[END_REF] and Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF], variants of the Gurson model have been developed, which take into account the effect of initial void size, while Lacroix et al. [START_REF] Lacroix | Numerical study and theoretical modelling of void growth in porous ductile materials subjected to cyclic loadings[END_REF] and Morin et al. [START_REF] Morin | A Gurson-type layer model for ductile porous solids with isotropic and kinematic hardening[END_REF] proposed a Gurson-type layer model with isotropic and kinematic hardening.

The widely-used Gurson model for describing the mechanical behavior of porous metallic materials involves only one variable for the analysis of void growth. The original Gurson model does not incorporate the void size in the yield surface. However, this model has been applied to materials with micron void sizes as well as submicron void sizes. Due to the mechanism of dislocation motion, plasticity depends upon internal length scale factor, and when the cavity size is equal to or less than the order of internal length scale, application of the original Gurson model is no more justified (see Hutchinson [START_REF] Hutchinson | Plasticity at the micron scale[END_REF]). It is also well known that higher strength is associated with smaller cavity sizes. This observation, in particular, has found applications in many advanced emerging fields, such as Micro Electro-Mechanical Systems (MEMS). In this context, it has been shown in the literature that the yield surface is larger as compared to the original Gurson yield surface for void sizes in the range of micron and submicron (see, e.g., Wen et al. [START_REF] Wen | The modified Gurson model accounting for the void size effect[END_REF]). Experimental investigations have also revealed the void size effects on material behavior (see, e.g., [START_REF] Fleck | Strain gradient plasticity: theory and experiment[END_REF][START_REF] Schlu | Modelling of the damage in ductile steels[END_REF][START_REF] Begley | The mechanics of size-dependent indentation[END_REF][START_REF] Nix | Indentation size effects in crystalline materials: a law for strain gradient plasticity[END_REF][START_REF] Stölken | A microbend test method for measuring the plasticity length scale[END_REF][START_REF] Kawasaki | An investigation of cavity growth in a superplastic aluminum alloy processed by ECAP[END_REF][START_REF] Khraishi | A parametric-experimental study of void growth in superplastic deformation[END_REF][START_REF] Fu | Geometry and grain size effects on the fracture behavior of sheet metal in micro-scale plastic deformation[END_REF][START_REF] Chentouf | Influence of predeformation on microstructure evolution of superplastically formed Al 5083 alloy[END_REF][START_REF] Li | Cavity Behavior of Fine-Grained 5A70 Aluminum Alloy during Superplastic Formation[END_REF]). Concurrently, several numerical and unit cell studies have equally shown the void size effects (see, e.g., Liu et al. [START_REF] Liu | The size effect on void growth in ductile materials[END_REF], Monchiet and Bonnet [START_REF] Monchiet | A Gurson-type model accounting for void size effects[END_REF]). In particular, Molecular-Dynamics analysis performed by Brach et al. [START_REF] Brach | A computational insight into void-size effects on strength properties of nanoporous materials[END_REF], for an aluminum single crystal, shows strengthening of the material with the decrease in void radius. Experiments have also shown that, for nano-porous materials, the decrease in void size causes an increase in material strength (see Biener et al. [START_REF] Biener | Nanoporous Au: A high yield strength material[END_REF][START_REF] Biener | Size effects on the mechanical behavior of nanoporous Au[END_REF], Hakamada and Mabuchi [START_REF] Hakamada | Mechanical strength of nanoporous gold fabricated by dealloying[END_REF]). Atomistic simulations performed by Mi et al. [START_REF] Mi | Atomistic insights into dislocation-based mechanisms of void growth and coalescence[END_REF] showed strong void-size dependency of the stress-strain response and porosity evolution for metallic alloys. They suggested that large void sizes give low values of yield stress as compared to small voids. It has also been found that the porosity evolution for small voids is delayed as compared to large voids before yielding, while it accelerates in plastic regime, causing the rapid decrease in stress carrying capacity for smaller voids as compared to larger voids. They pointed out, through their atomistic simulations, the need for enhancement of the Gurson model to incorporate the void size in the yield function. Their simulations show that, for the same value of initial porosity, smaller voids result in increased material strength. Xu et al. [START_REF] Xu | Size effect affected formability of sheet metals in micro/meso scale plastic deformation: experiment and modeling[END_REF][START_REF] Xu | Geometry and grain size effects on the forming limit of sheet metals in micro-scaled plastic deformation[END_REF] investigated experimentally the size effects on material hardening and forming limits, and proposed a constitutive model accounting for void size effects. Their observations have shown that small grain size, and thus small void size, leads to an increase in strain hardening and forming limits. Later, Xu et al. [START_REF] Xu | An investigation on the formability of sheet metals in the micro/meso scale hydroforming process[END_REF] have implemented their size-dependent constitutive model for the optimization of the hydroforming process of sheet metals.

To account for the interface stresses between the nano-voids and the matrix material, Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] extended the Gurson model to incorporate the void size effect by considering a hypothetical membrane surrounding the voids. A non-dimensional parameter is introduced in the resulting yield function, which depends on the void size and the strength of interface. Dormieux and Kondo [START_REF] Dormieux | Non linear homogenization approach of strength of nanoporous materials with interface effects[END_REF] and Brach et al. [START_REF] Brach | Strength properties of nanoporous materials: a 3layered based non-linear homogenization approach with interface effects[END_REF] further extended the void-size dependent model for nonlinear behavior of porous materials by applying nonlinear homogenization techniques. Their yield surfaces show the increase in material strength for smaller void sizes.

Parallel to the above-discussed constitutive models, which take into account the void size effect, various instability criteria have been developed in the literature for the prediction of diffuse and localized necking. Considère [START_REF] Considère | Memoire sur l'emploi du fer et de l'acier dans les constructions[END_REF] proposed the maximum force criterion as an indicator of diffuse necking. According to this criterion, instability occurs when the maximum load point is yielded during uniaxial tension. The extension of Considère's criterion to in-plane biaxial loading has been proposed by Swift [START_REF] Swift | Plastic instability under plane stress[END_REF]. Using the bifurcation approach, Hill [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF] presented a material instability criterion for localized necking, namely Hill's zero-extension criterion. However, this criterion predicts the ductility limit only for the left-hand side of the FLD. For the complete range of strain-path ratios, some authors suggested combining Hill criterion [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF] with Swift criterion [START_REF] Swift | Plastic instability under plane stress[END_REF]. Adopting the maximum force approach, Hora et al. [START_REF] Hora | A prediction method for ductile sheet metal failure in FEsimulation[END_REF] extended the Swift [START_REF] Swift | Plastic instability under plane stress[END_REF] and Considère [START_REF] Considère | Memoire sur l'emploi du fer et de l'acier dans les constructions[END_REF] criteria for localized necking by introducing the effect of strain path. Marciniak and Kuczyński [START_REF] Marciniak | Limit strains in the processes of stretch-forming sheet metal[END_REF] introduced another approach for the prediction of plastic instability, namely the initial imperfection approach, also referred to as the M-K method. The M-K imperfection approach is widely used in many applications for the prediction of localized necking. In another class of theoretical approaches, based on the condition of loss of uniqueness for the solution of the boundary value problem, the general bifurcation criterion was proposed by Hill [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF]. Similar to the general bifurcation criterion, the limit point bifurcation condition was proposed by Valanis [START_REF] Valanis | Banding and stability in plastic materials[END_REF]. Rudnicki and Rice [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF] and Rice [START_REF] Rice | Localization of plastic deformation[END_REF] extended the bifurcation approach in the form of loss of ellipticity condition. According to this criterion, material instability, in the form of plastic strain localization, occurs when the acoustic tensor becomes singular. Similar to Rice [START_REF] Rice | Localization of plastic deformation[END_REF] bifurcation criterion, Bigoni and Hueckel [START_REF] Bigoni | Uniqueness and localization-I. Associative and non-associative elastoplasticity[END_REF] and Neilsen and Schreyer [START_REF] Neilsen | Bifurcations in elastic-plastic materials[END_REF] proposed the condition of loss of strong ellipticity for the prediction of localized necking. Other types of plastic instability criteria have also been developed in the literature, which are based on the analysis of strain and stress evolutions during the loading by using finite element simulations (see, e.g., Volk and Hora [START_REF] Volk | New algorithm for a robust user-independent evaluation of beginning instability for the experimental FLC determination[END_REF], Situ et al. [START_REF] Situ | A suitable criterion for precise determination of incipient necking in sheet materials[END_REF], Narasimhan and Wagoner [START_REF] Narasimhan | Finite element modeling simulation of in-plane forming limit diagrams of sheets containing finite defects[END_REF]).

In the present work, the void-size dependent constitutive model proposed by Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] is combined with the Rice [START_REF] Rice | Localization of plastic deformation[END_REF] bifurcation criterion to predict forming limit diagrams of sheet metals. Similar approaches have been followed by Mansouri et al. [START_REF] Mansouri | Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis[END_REF] and Chalal and Abed-Meraim [START_REF] Chalal | Hardening effects on strain localization predictions in porous ductile materials using the bifurcation approach[END_REF], in which the Gurson-Tvergaard-Needleman (GTN) model has been coupled with the Rice bifurcation criterion for the prediction of ductility limits. In the latter studies, the effect of void size on the prediction of FLDs was not considered, which is taken into account in the present contribution. The organization of the paper is as follows. First, the constitutive equations related to the Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] model and the equations governing the Rice bifurcation criterion are presented. Then, validation of the resulting model and its numerical implementation is conducted using finite element simulations for aluminum alloy Al5754. Finally, numerical results in terms of FLDs obtained with the proposed approach are presented and discussed along with some concluding remarks.

Constitutive equations

Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] extended the original Gurson model by considering a membrane surrounding a spherical void inside a spherical representative volume element (RVE). The membrane around the void is acted on by a surface tension, which produces surface stresses around the voids. The strength of the membrane follows the von Mises criterion, with int k (N/m) being the cavity interface strength. Performing a limit analysis, the yield function has been derived by Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF], which incorporates the void size. The yield function is expressed in the form of parametric equations:

( ) ( ) ( ) ( ) eq m 1 1 m 2 2 2 2 eq 2 , 0, 6 3 tr 2 sinh ( ) sinh ( ) , 3 / 5 9 1 1 , 5 3 / 5 f f f f Σ - -   Φ Σ ξ Σ ξ =       ξ     Σ = = σ ξ - ξ + Γ       ξ +         Σ = σ + ξ - + ξ + Γ       ξ +   (1)
where m eq 2D D ξ = f is a dimensionless parameter defined as a function of the mean part m D and the equivalent part eq D of the macroscopic strain rate tensor D . In Eq. ( 1), Σ is the Cauchy stress tensor, σ is the yield stress of the fully dense matrix, f is the void volume fraction, eq 3 2

Σ Σ ′ ′ Σ = :
is the macroscopic equivalent stress, with Σ ′ and m Σ the deviatoric and hydrostatic part of the Cauchy stress tensor, respectively, and Γ is a non-dimensional parameter, which depends on the void size a and on the cavity interface strength int k as:

int k a Γ = σ . (2)
For illustration purposes, the above parametric yield function is plotted in Fig. 1 for 0, 0.2 and 0.43 Γ =

, along with the original Gurson model (i.e., without void size effect). It is worth noting that the Gurson model is recovered from the Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] model (D-K model) when 0 Γ = . Note that owing to the symmetry of the parametric yield surface with respect to eq Σ σ axis, only the first quadrant is plotted. From Fig. 1, it can be seen that the yield locus becomes larger for smaller void sizes (i.e., larger values of Γ ). As the void size increases (i.e., smaller values of Γ ), the size dependency of the D-K model becomes negligible, while for 0 Γ = (i.e., very large void sizes), the Gurson model is recovered, as shown in Fig. 1. In addition, two linear stress paths, corresponding to constant stress triaxiality of 0.333 Τ = (i.e., uniaxial tension) and 0.667 Τ = (balanced biaxial tension) are highlighted in Fig. 1. Note that, in sheet metal forming applications, the range of stress triaxiality lies typically between 0.333 and 0.667. In this range of stress triaxiality, it can be observed that the parametric yield locus is sensitive to void size.

Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF] have numerically implemented the D-K model by using the heuristic extension of Gurson's model for isotropic hardening. The yield stress σ , which was considered as a constant value in the D-K model (i.e., without hardening), is actually a function of the cumulated equivalent plastic strain p ε , which allows modeling isotropic hardening. The following isotropic hardening laws are considered in this work:

p p 0 p p 0 0 (ε ) (ε ε ) (Swift's hardening law), ε (ε ) 1 (Power hardening law), / n n K E σ = +      σ = σ +      σ   (3)
where 0 , ε and K n are hardening parameters, while 0 and E σ denote respectively the initial yield stress and Young's modulus. From the equivalence principle of the plastic work rate, the equivalent plastic strain rate p ε ɺ and the macroscopic plastic strain rate tensor p D are related as:

p p (1 ) ε f Σ D -σ = : ɺ . ( 4
)
Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF] have implemented into the D-K model the evolution of the porosity due to growth, based on the matrix incompressibility (see, e.g., Tvergaard, [START_REF] Tvergaard | Effect of yield surface curvature and void nucleation on plastic flow localization[END_REF]):

( ) p g (1 )tr . f f f D = = - ɺ ɺ (5)
The evolution of the void size (i.e., a ɺ ) has been modeled by Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF], by considering that for a spherical void of volume ω , present inside a spherical RVE of volume Ω , the volume constancy of the matrix material implies that ω Ω = ɺ ɺ . As the porosity is defined by ω Ω f = , the evolution equation for the porosity and thus for the void size writes:

( ) 1)

2 2 2 2 d ω Ωω ωΩ ω ω ω Ω 3 , dt Ω Ω ω Ω Ω ω . 3 (
a f f f a a a f f f     -      = = = - = -              = - ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ (6)
Using the definition of the Γ parameter, its time derivative can be expressed as:

( ) ( ) int int 2 d . dt k k a a a a   -   Γ = = σ + σ       σ σ ɺ ɺ ɺ (7)
The time derivative of the yield stress (i.e., σ ɺ ) can be evaluated as:

p p p p ε (ε )ε , ε h ∂σ σ = = ∂ ɺ ɺ ɺ (8)
where p (ε ) h is the hardening slope of the material, which depends on the selected isotropic hardening law.

To determine the elastic-plastic tangent modulus, the plastic multiplier γ ɺ and the yield function Φ can be combined in the form of the following Kuhn-Tucker relation:

0, 0, 0 Φ ≤ γ ≥ Φγ = ɺ ɺ . ( 9
)
The above expression shows that no plastic flow occurs (i.e., 0 γ = ɺ ) when 0 Φ < , while a strict plastic loading (i.e., 0 γ > ɺ ) necessarily implies Φ = 0 ɺ . The latter condition, called consistency condition, is used to derive the elastic-plastic tangent modulus as follows:

( )

V V V 0, f , f, , f Σ Σ V Σ σ Γ Φ Γ σ = : + σ + Γ + = ɺ ɺ ɺ ɺ ɺ (10)
where , V , V and V f Σ V σ Γ are the partial derivatives of the yield function Φ with respect to , , and f Σ σ Γ , respectively. Their expressions are derived as follows:

eq m m eq Φ Φ Φ , Σ V Σ Σ Σ ∂Σ ∂Σ ∂ ∂ ∂ = = + ∂ ∂Σ ∂ ∂Σ ∂ (11)
where

eq 2 m 2 m eq d Φ C , d d Φ C , d Σ ∂ = - ∂Σ ξ Σ ∂ = ∂Σ ξ (12)
with C a real constant, which depends on ξ . By using Eq. ( 12), Σ V can be expressed as:

eq eq 2 2 m m d d C C . d d Σ ∂Σ ∂Σ Σ = - + ξ ∂ ξ ∂ Σ V Σ Σ (13)
By taking the derivative of Eq. ( 1) with respect to parameter ξ , we obtain:

3 2 2 eq 2 2 2 2 3 1 2 2 2 2 2 m 2 2 2 d 9 3 , d 5 5 1 1 d 1 3 3 2 6 . d 3 5 5 1 1 - - -     Σ     ξ ξ Γ ξ     = σ - - ξ +           ξ + ξ + ξ                           Σ σ             = - + Γ -ξ ξ + + ξ +                       ξ + ξ + ξ                 f f f f f f ( 14
)
The remaining partial derivatives in Eq. ( 11) write:

m eq eq 1 , 3 3 , 2 1 Σ Σ Σ ∂Σ = ∂ ∂Σ ′ = ∂ Σ ( 15
)
where 1 is the second-order identity tensor. Finally, Σ V in Eq. ( 13) writes:

N eq 2 2 m eq d d 1 3 C C . d 3 d 2     Σ ′ Σ   = - + =     ξ ξ Σ     Σ Σ Σ V 1 V (16)
By following similar steps, the general expressions of the partial derivatives V , V and V f σ Γ can be obtained as follows:

eq m m eq eq m m eq eq m m eq

Φ Φ Φ V , Φ Φ Φ V , Φ Φ Φ V , f f f f σ Γ ∂Σ ∂Σ ∂ ∂ ∂ = = + ∂σ ∂Σ ∂σ ∂Σ ∂σ ∂Σ ∂Σ ∂ ∂ ∂ = = + ∂ ∂Σ ∂ ∂Σ ∂ ∂Σ ∂Σ ∂ ∂ ∂ = = + ∂Γ ∂Σ ∂Γ ∂Σ ∂Γ (17)
and more specifically:

N N N eq eq 2 2 m m eq eq 2 2 m m eq eq 2 2 m m d d V C C V , d d d d V C C V , d d d d V C C V . d d σ σ Γ Γ   Σ ∂Σ   ∂Σ Σ   = - + =     ξ ∂σ ξ ∂σ       Σ ∂Σ   ∂Σ Σ   = - + =     ξ ∂ ξ ∂       Σ ∂Σ   ∂Σ Σ   = - + =     ξ ∂Γ ξ ∂Γ     f f f f ( 18
)
The macroscopic plastic strain rate tensor p D is defined using the following classical plastic flow rule (normality law):

p . Σ D V = γ ɺ (19)
Substituting Eq. ( 16) in the above equation, the macroscopic plastic strain rate tensor becomes:

N N p 2 N C , Σ Σ D V V = γ = γ ɺ ɺ (20)
where

2 N C . γ = γ ɺ ɺ (21)
The macroscopic Cauchy stress rate tensor is expressed, in the co-rotational material frame, by the following hypoelastic law:

( )

e p ep , Σ C D D C D = : - = : ɺ ( 22
)
where D is the macroscopic strain rate tensor, e C is the fourth-order tensor of elasticity constants, and ep C is the elastic-plastic tangent modulus, which needs to be determined. By using all of the above-described partial derivatives of the yield function, and considering Eqs. ( 20) and ( 22), the consistency condition [START_REF] Besson | Non-linear mechanics of materials[END_REF] can be developed as follows:

N N N N N N e e N V V V 0. f f Σ Σ Σ V C D V C V σ Γ : : -γ : : + σ + Γ + = ɺ ɺ ɺ ɺ (23)
Then, substituting Eqs. (4-8) and ( 20) into the above consistency condition, the expression of the plastic multiplier is derived as follows:

N e N , H Σ V C D γ : : γ = ɺ ( 24
)
where

N N N N N N N N int e int 2 V H V V (1 )( ) V . (1 ) 3 (1 ) f k k h f f a f f a Σ Σ Σ Σ Σ V V C V V 1 Γ γ σ Γ   :           = : : - - -- : -         -σ σ - σ         (25)
Finally, the elastic-plastic tangent modulus is derived as:

( ) ( ) N N e e ep e . H Σ Σ C V V C C C γ : ⊗ : = - (26)

Material objectivity

The constitutive equations described above are implemented into the finite element code ABAQUS/Standard in the framework of large strains. To achieve this, objective derivatives for the tensor variables must be used to ensure the material objectivity (i.e., frame invariance of the constitutive model). A convenient way to maintain the material objectivity is to reformulate the tensor variables in terms of their rotation compensated counterparts. More specifically, if A and B represent a second-and a fourth-order tensor, respectively, their associated rotation compensated expressions can be calculated as:

ij ki l j kl A R R A = ⌢ and ijkl pi qj rk sl pqrs B R R R R B , = ⌢ ( 27
)
where R is a rotation matrix, which is defined using the following rate equation:

T , R R Ω ⋅ = ɺ ( 28
)
where Ω is the skew-symmetric spin tensor. Different objective rates of the Cauchy stress tensor are available in the literature. Among them, the well-known Jaumann stress rate, Green-Naghdi stress rate, and Truesdell stress rate. In the present contribution, Jaumann objective rates are considered for the constitutive equations, which is consistent with the procedure used in the finite element sotware ABAQUS/Standard. The Jaumann objective stress rate J Σ ɺ is expressed as:

J . Σ Σ Ω Σ Σ Ω -⋅ + ⋅ = ɺ ɺ (29)

Bifurcation approach

In the present study, the strain localization criterion proposed by Rudnicki and Rice [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF] and Rice [START_REF] Rice | Localization of plastic deformation[END_REF] is combined with the D-K model to predict the ductility limits of sheet metals. This criterion is based on the bifurcation theory, which predicts the occurrence of strain localization in the form of infinite planar band defined by its normal n (see Fig. 2). The equilibrium conditions along the localization band, which express the continuity of the nominal stress rate vector across the discontinuity surfaces, can be written as:

n N 0 ⋅ = ɺ , ( 30
)
where

G m n ⊗ = ɺ , ( 31
)
where the vector m G n = ⋅ ɺ is the jump amplitude vector, which characterizes the mode of localization (e.g., shear bifurcation mode when m n ⊥ ɺ

).

On the other hand, the nominal stress rate tensor N ɺ is related to the velocity gradient G by the following expression:

N L G = : ɺ , ( 32
)
where L is the corresponding tangent modulus that needs to be determined. By combining Eqs. (30-32), the critical condition, which corresponds to the loss of ellipticity for the boundary value problem, can be expressed as the singularity of the acoustic tensor Q n L n = ⋅ ⋅ :

( ) ( )

det det 0 Q n L n = ⋅ ⋅ = . ( 33
)
As for the sake of material objectivity, Jaumann stress rate has been used in the present work. Within the framework of an updated Lagrangian approach, the expression of the nominal stress rate tensor is given by: ( )

tr N Σ D Σ D Σ Ω Σ = + -⋅ -⋅ ɺ ɺ . ( 34
)
Finally, by combining the above equations, the complete expression of the tangent modulus L is given by: ( )

ep 1 2 3 L C C C C = + --, (35)
( ) 1i jkl ij kl 2ijkl jk il jl ik 3ijkl ik jl il jk C , 1 C , 2 1 C , 2 Σ δ Σ δ Σ δ Σ δ Σ δ = = + = - (36)
where ij δ is the kronecker symbol, which is equal to 1 when i j = , and to 0 otherwise.

It is worth noting that the present work is dedicated to the prediction of localized necking in thin sheet metals. Although the constitutive equations presented in the previous section are formulated in the fully three-dimensional framework, the Rice bifurcation criterion is derived in this work in the framework of plane-stress conditions. To achieve this, the nominal stress rate tensor, defined by Eq. ( 32), is rewritten within the framework of plane-stress theory as follows:

PS N L G , with , , , 1, 2
αβ αβγδ γδ = α β γ δ = ɺ , ( 37
)
where the expression of the plane-stress tangent modulus PS L is deduced from the threedimensional tangent modulus L using the following relationship:

33 PS 33 3333 L L L L L γδ αβγδ αβγδ αβ = - . ( 38
)
Note that, within the plane-stress framework, the emergence of the localization band is searched for in the plane of the sheet, with a variation of the band orientation angle θ (see Fig.

(2)) from 0° to 180°. Localized necking is said to occur when the minimum value, over all possible band orientations, of the determinant of the acoustic tensor (see Eq. ( 33)) becomes zero. The corresponding in-plane strains are considered as the critical strains at localization, which will be used for the FLD plot.

Numerical implementation of the D-K model and validation

Time integration scheme and finite element implementation

The constitutive equations described above have been implemented into the finite element code ABAQUS/Standard in the framework of large strains and three-dimensional approach. In order to reproduce a homogeneous strain state prior to localization, a single finite element with one integration point (specifically, C3D8R solid element in ABAQUS) is used in the simulations, which is subjected to various linear strain paths that are those typically applied to sheet metals under in-plane biaxial stretching. The main motivation behind this choice of loading configuration is to predict plastic instability that is inherent to the 'material' alone, with no interference with structural (geometric) effects. The geometry of the single finite element along with the applied boundary conditions are illustrated in Fig. 3. Directions 1 and 2 represent the major and the minor directions, respectively. The linear strain paths are obtained by varying the strain-path ratio from -0.5, for uniaxial tension, to 1 for balanced biaxial tension. For each loading increment, the stress state is updated using a predictor-corrector approach, following the numerical procedure proposed by Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF]. The resulting numerical algorithm is detailed in Table 1. As to the numerical implementation of the Rice bifurcation criterion, the condition of loss of ellipticity, given by Eq. [START_REF] Stölken | A microbend test method for measuring the plasticity length scale[END_REF], is assessed by computing the determinant of the acoustic tensor for each loading increment and each orientation for the localization band. The numerical detection of the occurrence of localized necking is achieved when the minimum of the acoustic tensor determinant, over all possible orientations for the normal n to the localization band, becomes non-positive. This procedure is repeated for different strain-path ratios ranging from uniaxial tension to balanced biaxial tension. Then, the corresponding in-plane critical strains are plotted in terms of FLD.

Table 1

Outline of the predictor-corrector algorithm and time integration scheme. Σ + and ep C to ABAQUS for solving the equilibrium equations.

Σ T Σ = ) Compute ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 elast elast m eq m n 1 eq n 1 n 1 n 1 n 1 n 1 F Σ Σ Σ Σ + + + + + +     = + - ξ + ξ         Plastic yield condition: • If F 0 < then (elastic

Numerical validation of the D-K model

The above numerical algorithm for the time integration of the constitutive equations of the D-K model has been implemented into the finite element software ABAQUS/Standard via a user-defined material (UMAT) subroutine. This section is devoted to the validation of the numerical implementation of the D-K model. To this end, simulations of the response of an RVE, which is subjected to constant stress triaxialities, are performed with the developed D-K UMAT subroutine and compared with the results provided by Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF]. The powerhardening law, as given by Eq. ( 3), is used for the description of isotropic hardening. The material parameters used in the simulations are summarized in Table 2 (see Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF]). Note that the RVE is subjected to two constant stress triaxialities (i.e., T=1 and T=3). The numerical results in terms of normalized equivalent stress and porosity with respect to equivalent strain, obtained with the present UMAT subroutine, are shown in Figs. 4567along with the reference results provided by Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF].

Table 2

Material parameters used in the simulations. It is recalled here that the original Gurson model is recovered from the D-K model when 0 0 Γ = . In this case, the numerical results obtained with the present UMAT subroutine can be compared with those predicted by the Gurson model, which is available in ABAQUS/Standard. When the void-size independent model is considered (i.e., D-K model with 0 0 Γ = ), all the results predicted using the developed UMAT subroutine coincide with their counterparts obtained with the Gurson model, which is fully consistent. As to the case of void-size dependent model (i.e., D-K model with 0 0 Γ ≠ ), it can be observed a significant difference between the results predicted by the Gurson model and those provided by the D-K model, which includes the void size effect. More specifically, for small initial porosity and moderate stress triaxiality (i.e., softening is observed for small values of 0 Γ , which is due to the rapid evolution of the porosity in this range of small initial non-dimensional parameter 0 Γ (see Fig. 4b). For moderate initial porosity as well as high stress triaxiality (i.e., 0 = 0.01 f and T=3, respectively), the developed UMAT subroutine for the D-K model predicts the same trends as above for both stress and porosity evolutions, which are consistent with the reference predictions provided by Morin et al. [START_REF] Morin | Numerical assessment, implementation and application of an extended Gurson model accounting for void size effects[END_REF].

E (GPa) ν σ 0 (MPa) n 0 f int k (N/m) 0 Γ 200

Prediction of FLDs with the D-K model and bifurcation analysis

In this section, the resulting coupling between the D-K model and the Rice bifurcation criterion is applied for the prediction of FLDs for Al5754 aluminum material. The associated material parameters used in the simulations are taken from Mansouri et al. [START_REF] Mansouri | Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis[END_REF], as reported in Table 3. The additional parameter 0 Γ , which incorporates the void size effect is varied here to analyze its effect on the prediction of localized necking.

Table 3

Material parameters for the Al5754 aluminum material.

Elastic properties Swift's hardening parameters

Initial porosity

E (GPa) ν K ε 0 n 0 f 70 0.33 309.1 0.00173 0.177 0.001
The FLDs predicted with the Rice bifurcation criterion are presented in Fig. 8a for the Al5754 aluminum material for three values of the initial non-dimensional parameter (i.e., 0 Γ 0, 1.5 and 2.5 = ). It is worth recalling that the original Gurson model, which only describes the elastic-plastic behavior with damage due to porosity growth, can be recovered from the D-K model when 0 Γ 0 = . Therefore, for comparison purposes, the FLD predicted with the original Gurson model and the Rice bifurcation criterion is also presented in Fig. 8a. From this figure, it can be seen that the effect of the initial non-dimensional parameter on strain localization is not perceptible when only the void growth mechanism is considered in the D-K model. Moreover, strain localization does not occur in the right-hand side of the FLD (i.e., for positive biaxial strain-path ratios) due to the low values of porosity in this range of strain paths (i.e., right-hand side of FLD, see Fig. 8b). Indeed, in the works of Mansouri et al. [START_REF] Mansouri | Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis[END_REF] and Chalal and Abed-Meraim [START_REF] Chalal | Hardening effects on strain localization predictions in porous ductile materials using the bifurcation approach[END_REF], it has been shown that, using the bifurcation theory coupled with Gurson-type models, the occurrence of strain localization for positive biaxial strain-path ratios requires the critical hardening modulus to be strongly negative.

In order to further investigate the insensitivity of the strain localization predictions to the different values of the initial non-dimensional parameter (see Fig. 8a), the yield surface for Al5754 is illustrated in Fig. 9. In addition, two linear stress paths, corresponding to the extreme left (uniaxial tension) and right (balanced biaxial tension) strain paths of FLD, i.e. constant stress triaxialities of 0.333 Τ = and 0.667 Τ =

, respectively, are highlighted in this Fig. 9. Note that, within this narrow range of stress triaxiality ratio, the yield surface points corresponding to different values of initial non-dimensional parameter are almost indistinguishable, which explains the insensitivity of the predicted FLDs to the initial non-dimensional parameter, as observed in Fig. 8a. To better highlight the effect of the void size on strain localization, nucleation of new voids is considered in what follows, in addition to the void growth mechanism, for both the D-K model and the Gurson model. Accordingly, the evolution of porosity, with an initial value of 0 f , is due to void nucleation and growth, and the porosity rate can then be partitioned as follows:

g n , f f f = + ɺ ɺ ɺ (39)
where n f ɺ represents the contribution to the porosity rate from nucleation. The latter is considered to be strain-controlled following the normal distribution relationship proposed by Chu and Needleman [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF]:

2 p p N N n N N ε ε 1 exp ε , 2 s s 2 f f      -       = -           π       ɺ ɺ (40)
where N N N , s and ε f are the nucleation parameters. The latter are reported in Table 4 for the Al5754 aluminum material. The corresponding FLDs are presented in Fig. 10a, for which the porosity evolution due to nucleation and growth is considered.

Table 4

Nucleation parameters for the Al5754 aluminum material. In contrast to the FLDs obtained with porosity evolution only due to growth (i.e., D-K model and Gurson model, see Fig. 8a), Fig. 10a shows that the effect of void size on strain localization is more perceptible when both growth and nucleation of porosity are considered. More specifically, as observed in Section 3.2, it is expected that smaller values of the initial non-dimensional parameter induce rapid softening (see the stress-strain curves in Figs. 4 and5), thereby promoting early plastic flow localization. This expectation is also confirmed by Fig. 10a, which shows that the predicted FLDs are lowered as initial non-dimensional parameter 0 Γ decreases. This trend is confirmed by Fig. 10b, which also shows the effect of the initial nondimensional parameter on the values of porosity at localization. Note that when 0 Γ 0 = , the FLD predicted by the D-K model coincides with that of the Gurson model, which is also consistent with the theoretical expectation for this particular case (i.e., the Gurson model is recovered from the D-K model when 0 Γ 0 = ).

In addition to growth and nucleation of voids during plastic deformation in ductile materials, coalescence of voids can also be introduced to model the rapid decay of the material stress carrying capacity. Tvergaard and Needleman [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] have modified the original Gurson model to account for the complete kinetics of voids within the material (i.e., nucleation, growth and coalescence), which corresponds to the well-known Gurson-Tvergaard-Needleman (GTN) model. Using the GTN model, the actual void volume fraction f in the expression of the yield surface is replaced by the effective porosity * f , which accounts for the void coalescence mechanism as follows:

* cr GTN cr ( ), f f δ f f = + - (41)
where cr f is the critical porosity, which marks the onset of the coalescence regime, while GTN δ is the accelerating factor that governs the rapid decay of stress carrying capacity at the onset of coalescence. Note that when This approach, which considers the coalescence of voids beyond a critical porosity, is also used for the D-K model in what follows.

The parameters associated with the coalescence regime for the Al5754 aluminum material are reported in Table 5, while the elastic-plastic and nucleation parameters are kept the same as those used in Tables 3 and4. The corresponding FLDs obtained with the developed approach are presented in Fig. 11a. As discussed above, coalescence causes sudden loss of the stress carrying capacity and, therefore, the hardening modulus becomes strongly negative. It has been shown in Mansouri et al. [START_REF] Mansouri | Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis[END_REF] that the prediction of strain localization using the Rice bifurcation criterion requires strongly negative hardening modulus in the right-hand side of the FLD (i.e., for biaxial stretching loading paths), which can be reached by considering the coalescence mechanism. Consequently, it can be observed in Fig. 11a that more localization points are obtained in the right-hand side of the FLD, as compared to the FLDs in Figs. 8a and10a, thanks to the consideration of the void coalescence mechanism. Note that the experimental FLD, which is taken from Brunet et al. [START_REF] Brunet | Analytical and experimental studies of necking in sheet metal forming processes[END_REF], is also shown in Fig. 11a for qualitative comparison. Analogous to the previous FLDs, obtained using only growth and nucleation of voids, the FLDs in Fig. 11a are clearly sensitive to variation of the initial non-dimensional parameter, with lowered FLDs for smaller values of the initial non-dimensional parameter. Note also that for 0 Γ 0 = , the FLD predicted by the D-K model coincides with that of the GTN model, for which the q-parameters for the void interactions are set equal to 1. Furthermore, the yield surfaces for different values of initial non-dimensional parameter are presented in Fig. 12, in which the evolution of porosity is attributed to growth, nucleation and coalescence. In contrast to Fig. 9, with porosity evolution only due to growth, Fig. 12 depicts distinct yield surfaces for different initial values of non-dimensional parameter, and for all possible values of stress triaxiality ratio, which leads in turn to a more perceptible effect on the prediction of FLDs. In Fig. 12, the von Mises yield surface is also illustrated. It is noteworthy that when the initial non-dimensional parameter is increased (i.e., void size is decreased), an increase in the material strength is observed (i.e., expansion of the yield surface). However, the strength of porous material cannot be larger than the strength of the sound material (i.e., undamaged material). This requires that the initial non-dimensional parameter 0 Γ must vary within a physically admissible range, i.e. ).

+ ɺ ɺ ɺ f = f f Gurson model 0 D K model with Γ 2.5 - = 0 D K model with Γ 1.5 - = 0 D K model with Γ 0 - = (a) (b)
It is worth noting that the coalescence criterion used in the present model does not take the effect of void size into account. Other coalescence criteria including the effect of void size can be combined with the present approach. Among them, the constitutive model for nanoporous materials proposed by Gallican and Hure [START_REF] Gallican | Anisotropic coalescence criterion for nanoporous materials[END_REF], in which the effect of void size on coalescence regime has been considered.

The effect of void size on the ductility limits for high initial porosity is also investigated hereafter. Note that, in the previous simulations, the initial porosity for the Al5754 aluminum material has been identified by Brunet et al. [START_REF] Brunet | Analytical and experimental studies of necking in sheet metal forming processes[END_REF] to a value of 0.001. To reveal the effect of void size on strain localization for larger values of initial porosity, Fig. 13 provides the FLDs predicted by the D-K model for different values of initial non-dimensional parameter, when the initial porosity is set to 0.01. The simulation results reported in Fig. 13 are obtained with the D-K model, in which the porosity evolution takes into account all three mechanisms of growth, nucleation and coalescence of voids. As previously discussed for Gurson-type models in general, and for the D-K model in particular, porosity is the most important parameter for the sensitivity of FLDs to the void size. As shown in Fig. 13, for high initial porosity (i.e., 0 0.01 f = ), the effect of the initial non-dimensional parameter on strain localization is more significant than in the case of smaller initial porosity (i.e., 0 0.001 f = , see Fig. 11a). Indeed, increasing the initial nondimensional parameter (i.e., decreasing the void size) will increase the number of voids within the material, in order to keep the same actual porosity, which will increase in turn the number of interfaces between voids and the dense matrix. Therefore, material strain hardening will be increased in this case in the same way as material ductility. This observation has also been found by Wen et al. [START_REF] Wen | The modified Gurson model accounting for the void size effect[END_REF], who concluded that the void size effects on material response is more significant for higher initial porosity.

Conclusions

In the present contribution, the phenomenological model proposed by Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] for ductile porous materials accounting for void size effects has been combined with the Rice bifurcation criterion to predict strain localization in thin sheet metals. The constitutive equations have been implemented together with the Rice bifurcation criterion into the finite element code ABAQUS/Standard within the framework of large strains and fully threedimensional approach. However, the localization bifurcation criterion has been reformulated within the plane-stress framework, as classically done for the prediction of localized necking in thin sheet metals. Linear in-plane biaxial stretching loading paths have been applied for the prediction of forming limit diagrams for the Al5754 aluminum material.

Recall that the Dormieux and Kondo [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF] model, which is an extension of the original Gurson model, has been developed by performing a limit analysis for a spherical void inside a spherical RVE, and by considering a membrane with surface stresses around the void. The resulting yield function depends on the void size through a non-dimensional parameter. Considering at first only the porosity evolution due to void growth, the FLDs predicted with the present approach show no sensitivity to the void size, due to the fact that the porosity at localization remains very small. However, when successively nucleation and coalescence are considered in the evolution of porosity, FLDs exhibit more sensitivity to the void size. More specifically, it is found that the ductility limits are lowered as the initial non-dimensional parameter is decreased (i.e., when the void size is increased). In other words, decreasing the void size has a beneficial effect on formability. These findings are consistent with what has been reported in the literature regarding the void size effects on material response. Nonetheless, there is an upper bound for the initial non-dimensional parameter 0 Γ (i.e., lower bound for the void size), which is physically admissible. This upper bound, which depends on the material parameters, is identified so that the parametric D-K yield surface must lie between the Gurson yield surface and the undamaged elastic-plastic yield surface (i.e., sound material with no void). Moreover, it is found that the FLD predicted using the D-K model with an initial nondimensional parameter 0 Γ 0 = coincides with that obtained with the Gurson model, which is also consistent with the material response for this particular case. Finally, the investigation of the combined effect of void size and initial porosity shows that smaller void sizes lead to an increase in the ductility limits, and this trend becomes more significant for high initial porosity, due to the increase of void-matrix interface strength within the material.

Fig. 1 .

 1 Fig. 1. Parametric yield surface for 0.1 = f and 0, 0.2 and 0.43 Γ = .

Fig. 2 .

 2 Fig. 2. Schematic illustration of localization band in the plane of sheet metal. Maxwell's compatibility equation allows expressing the jump in the velocity gradient G across the localization band as:

3 C

 3 are fourth-order tensors, whose expressions are only function of the Cauchy stress components. Their expressions are:

Fig. 3 .

 3 Fig. 3. Boundary conditions for single finite element (C3D8R) subjected to in-plane biaxial stretching.

For

5 , 0 Γ 1 . 25 , 0 Γ 1 . 25 , UMAT = 0 Γ 2 . 5 , 0 Γ 2 Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 501250125025024567 Fig. 4. Evolution of: (a) normalized equivalent stress, and (b) porosity, for stress triaxiality T=1 and initial porosity 0 = 0.001 f .

Fig. 8 .

 8 Fig. 8. (a) FLDs, and (b) porosity at localization, for Al5754 aluminum.

Fig. 9 .

 9 Fig. 9. Parametric yield surface for Al5754 aluminum and 0 0, 1.5 and 2.5 Γ = .

Fig. 10 .

 10 Fig. 10. (a) FLDs, and (b) porosity at localization, for Al5754 aluminum. Porosity evolution is due to growth and nucleation.

  e., only growth and nucleation of voids are considered), the accelerating factor GTN δ is set to 1, while GTN 1 δ > in the case of coalescence (i.e.,

Fig. 11 .

 11 Fig. 11. (a) FLDs, and (b) effective porosity at localization, for Al5754 aluminum. Porosity evolution is due to growth, nucleation and coalescence.

Fig. 12 .

 12 Fig. 12. Parametric yield surface for Al5754 aluminum and 0 0, 1.5 and 2.5 Γ = . Porosity evolution is due to growth, nucleation and coalescence.

  the parametric D-K yield surface must lie between the Gurson yield surface (i.e., 0 Γ 0 = ) and the undamaged elastic-plastic yield surface. This physical range of initial non-dimensional parameter 0 Γ can be identified by comparing the parametric D-K yield surface for various values of 0 Γ with the particular cases of Gurson yield surface (i.e., 0 Γ 0 = ) and von Mises yield surface (i.e., undamaged isotropic elastic-plastic model), as illustrated in the above Fig.12. This figure reveals that the physical upper bound of the nondimensional parameter 0 Γ is 2.5, for which the ratio between the equivalent stress and the yield stress of the fully dense matrix is lower or equal to 1 (i.e.,

Fig. 13 .

 13 Fig. 13. FLDs for Al5754 aluminum alloy with 0 0.01 f = . Porosity evolution is due to growth, nucleation and coalescence of voids.

Table 5

 5 Coalescence parameters for the Al5754 aluminum material.

	Material	f	cr	GTN δ
	Al5754	0.00284	7

= 0.001 f and T=1, respectively), a softening regime is predicted by the Gurson model, starting from an equivalent strain of about

1, while no softening regime is shown by the D-K model for high values of 0 Γ (i.e., from 0 1.25 Γ = and more in this case, see Fig.4a). However, stress


M Ibrahim

M M Wei

E Deydier

E Manoury

R Poli

P Lecante

K Philippot

Rhodium nanoparticles stabilized by ferrocenyl-phosphine ligands: synthesis and catalytic styrene hydrogenation

Keywords: Rhodium, nanoparticle, ferrocenylphosphine ligand, catalysis, styrene hydrogenation

A series of ferrocenyl-phosphine ligands has been used to stabilize small rhodium nanoparticles that are active in the catalytic hydrogenation of styrene.

Introduction

Metal nanoparticles (NPs) are largely investigated given their unique properties and potential applications in multiple domains such as biology, medicine, microelectronics, optics or catalysis to cite only a few. [START_REF]Physics and Chemistry of Metal Cluster Compounds[END_REF][START_REF]Nanoscale Materials in Chemistry[END_REF][START_REF]Nanoparticles: From Theory to Application[END_REF] In catalysis, the relevance of metal NPs has been known for a long time due to their presence in heterogeneous catalysts. [START_REF] Mondloch | [END_REF]5 Metal NPs are attractive materials because of their high surface/volume ratio. [START_REF]Concepts in Nanocatalysis[END_REF][START_REF] Roucoux | [END_REF] With the recent advances in nanochemistry and particularly in the design of better defined metal NPs, nanocatalysis is now recognized as a full and rich part of catalysis. 8,[START_REF]Nanoparticles in catalysis[END_REF][START_REF]Metal Nanoparticles for Catalysis: Advances and Applications[END_REF] During the past 25 years ligand-stabilized metal nanoparticles have attracted growing attention and numerous investigations are devoted to exploring the effect of the stabilizing agents on the characteristics of NPs such as size, dispersion, morphology and surface properties or also catalytic performance. [START_REF] Tschan | [END_REF] Very recent papers reported on the precise role of ligands and their binding with the metal surface on the NP size control. 12,13 It is important to note that the interaction of the ligands with the metal atoms on the NP surface can be compared to the ligand interaction with the metal centers in homogeneous catalysts, which is of paramount importance for stability and properties like catalytic activity. Thus, ligands also allow tuning the surface properties of metal NPs through electronic or steric effects.

Moreover, ligand-stabilized metal nanoparticles can be applied to catalysis either as stable colloidal suspensions or after deposition on the surface or confinement in the pores of solid supports. 14 The challenge is to find ligands able to stabilize well-defined nanoparticles while controlling accessibility at the metal surface and reactivity. 15 Compared to the investigation of facet dependency, the ligand influence on the catalytic activity has been less intensively studied but recent results illustrate well the interest to do so. 16,17,18,19 One aim lies on the design of more performing catalysts in terms of reactivity and selectivity in order to develop more efficient and eco-compatible chemical production. 20 One efficient method to prepare well-controlled metal nanoparticles is the organometallic approach. This approach has proven powerful to provide controlled nanoparticles with reproducible physical or chemical properties. 21 Metal complexes (organometallic or metal organic complexes) are decomposed under mild conditions to release metal atoms either in the presence of polymers (steric barriers) or ligands (metal surface coordination) as stabilizing agents. The choice of the metal source is important in terms of the resulting metal surface purity. Olefinic metal precursors are privileged because their treatment under hydrogen liberates alkanes as the only reaction byproducts that are inert towards the metal surface. So, metal NPs with clean surfaces capped with the desired ligands are obtained. Numerous ruthenium nanoparticles have been prepared this way, using [Ru( 6 -C8H10)((-C8H12)] as ruthenium source and a large panel of ligands as stabilizers. The ligands are either simple or sophisticated, and mainly contain amines, thiols, carbenes or phosphines that were shown to influence the surface and the catalytic properties of the particles. [START_REF] Catalysis | Organometallic Ruthenium Nanoparticles and[END_REF] Rhodium particles can also be synthesized by the organometallic method, starting from [Rh(acac)(1,5-C8H12)], [START_REF] Phd | [END_REF] [Rh(µCl)(1,5-C8H12)]2 [START_REF] Mavila | [END_REF] , or [Rh( 3 -C3H5)3]. 25 Phosphines ligands were shown to affect the reactivity of Rh NPs in hydrogenation 26,27 and hydroformylation catalysis. 28 Apart from usual ligands to stabilize metal NPs, one can envisage to use metalcontaining ligands. This can be of interest in order to introduce novel properties like redox surface modulation and also more functionality. 29,30 One application could be tandem catalysis for instance. One strategy to achieve this objective is to consider metal-containing ligands with functional groups able to interact with the metal surface such as phosphine groups. 31 To this end, we chose to test ferrocene-based phosphine ligands because these compounds have several attractive properties: they are stable, they may feature planar chirality with potential interest in materials chemistry 32 and especially in asymmetric catalysis, [START_REF] Chen | Phosphorus Ligands in Asymmetric Catalysis[END_REF][START_REF] Gómez Arrayás | [END_REF]35 and they possess a redox-active ferrocene moiety which may be exploited for switchable catalysis. 36,37,38 For instance, 1,2-disubstituted ferrocene phosphines play a major role in several asymmetric catalysis applications. Various examples including bidendate (P,P), 39 (P,N), 40 (P,S), 41 , 42 , [START_REF] Malacea | Phosphorus Ligands in Asymmetric Catalysis[END_REF] (P,O) [START_REF] Mateus | [END_REF]45,46,47 and (P, carbene) 48 compounds have been reported, illustrating the great modularity of this family of ligands. To our best knowledge there is no previous example of the use of such ligands for the stabilization of metal nanoparticles.

In the present work we have focused on a series of racemic ferrocenylphosphine ligands as stabilizing agents for rhodium nanoparticles obtained by a one-pot synthesis from [Rh( 3 -C3H5)3]. As it will be described hereafter, a series of small rhodium nanoparticles in a size range between 1.1 and 1.7 nm have been obtained and these were shown to be active for a model catalytic reaction, namely the hydrogenation of styrene.

Results and discussion

The first objective of this work was to take benefit of previous results in the organometallic synthesis of phosphine-stabilized rhodium nanoparticles, achieved from the [Rh( 3 -C3H5)3] complex as metal source and simple mono-or di-substituted phosphines as stabilizers, 26,28 in order to prepare a new family of Rh NPs by the same approach but using ferrocenyl phosphines. With ferrocene-based-phosphine-capped Rh NPs in hands, a second objective was to test their catalytic activity in a benchmark hydrogenation process, the hydrogenation of styrene, in order to evaluate the influence of the capping ligands at the rhodium surface.

All the ligands employed in this study (Scheme 1) are ferrocenyl diphenyl phosphines bearing different substituted methyl groups at the 2 position. L1 and L2 are both P,N bidendate ligands. 49 In L1 the amine is tertiary with an extra methyl substituent on the nitrogen atom, whereas L2 contains a secondary amine. Ligand L3 50 and L4 35c are respectively P,S and P,O bidendate ligands. L5 is the simplest ligand of the series, the only monodentate one, having a methyl group at the ferrocenyl 2 position.

Scheme 1: Ferrocenyl phosphine ligands used for the synthesis of RhNPs

The synthesis of the NPs was performed in one pot using dichloromethane (CH2Cl2) as the solvent in order to insure the solubility of both the [Rh( 3 -C3H5)3] complex and the ferrocenyl phosphine ligands. Thus, a Fischer-Porter reactor filled with a CH2Cl2 solution containing both [Rh( 3 -C3H5)3] and a given ferrocenyl ligand was pressurized with 3 bar of dihydrogen and then heated at 40°C, under vigorous stirring (Scheme 2). A change of color from a yellow homogeneous solution to a black colloidal suspension was observed over a few hours of stirring, but the reaction was kept running up to 16 h to insure a complete decomposition of the precursor and the formation of a homogeneous population of particles. For the four ligands L1, L2, L3 and L4, a [ligand]/[Rh] ratio of 0.2 molar equivalent was used. These ligands being bidentate, they offer an extra coordination through the second donor group, for an effective [donor]/[metal] ratio of 0.4. For ligand L5, which is monodentate, a molar

[ligand]/[Rh] ratio of 0.4 was applied for comparison purpose. These [ligand]/[Rh] ratios appeared to be sufficient to yield stable colloidal suspensions for several weeks under an argon atmosphere. This is particularly important to note when considering that the L5 ligand has only one -PPh2 coordinating group. Indeed, when PPh3 was used as a stabilizer, a minimum [ligand]/[Rh] ratio of 1 was necessary to reproducibly yield stable colloidal suspensions of small Rh NPs (ca. 1.3 nm). 26 For comparison, the bidentate Ph2P-(CH2)4-PPh2 ligand led to stable Rh NPs at a [ligand]/[Rh] ratio of 0.3 molar equivalent, which amounts to a [donor]/[Rh] ratio of 0.6. Consequently, the bidendate phosphino-ferrocenyl ligands used here behave as efficient stabilizers at lower [ligand]/[Rh] ratios compared to the previously studied similar mono-and diphosphane ligands. The purification of the particles was performed by precipitation after addition of cold pentane to the concentrated colloidal suspensions. The generated black precipitates were isolated by decanting off the solvents, washing and drying under vacuum, leading to dried nanoparticle powders.

Scheme 2:

On-pot synthesis of ferrocenylphosphine-capped Rh NPs TEM images recorded from samples prepared by deposition of a few drops of the crude colloidal suspensions onto a carbon-covered copper grid showed well-dispersed NPs for all the systems studied. They present spherical shapes and quite narrow size distributions, in a size range between 1 and 2 nm (Figure 1), with mean diameters as follows: 1.5 ± 0.6 nm for Rh@L1 NPs, 1.7 ± 0.5 nm for Rh@L2 NPs, 1.3 ± 0.7 nm for Rh@L3 NPs =, 1.1 ± 0.5 nm for Rh@L4 NPs, and 1.2 ± 0.6 nm for Rh@L5 NPs. The L2-and L5-stabilized nanoparticles were investigated in greater details in order to learn more about the effect of a second donor function. Thus, Rh@L2 NPs and Rh@L5 NPs were also characterized by high resolution electron microscopy (HREM). In addition to isolated nanoparticles, the HREM images of the L2-stabilized Rh NPs also revealed the presence of large agglomerates (Figure 2 top), which appeared composed of small and spherical particles gathered together. The isolated particles displayed crystalline regions with interplanar distances of 0.22 nm, in good agreement with the fcc structure (face-centered cubic) previously observed for similarly made Rh nanoparticles and for bulk rhodium. 26 In the case of L5-stabilized Rh NPs no agglomerates were observed besides isolated NPs. These latter were poorly crystallized with only a few of them displaying crystalline plans. Moreover Rh-Rh distances measured did not fit perfectly with the fcc crystalline structure as for the bulk rhodium (Figure 2 bottom). Such a phenomenon was previously observed for PPh3-stabilized RhNPs and explained as the result of a strong coordinating effect of the phosphine ligand. 26 For the Rh@L2 powder, WAXS (Wide-Angle X-ray Scattering) data revealed wellcrystallized particles and confirmed their fcc crystalline structure as for the bulk metal. From the RDF, coherence length reached 3.5-4.0 nm, much more than the average diameter of ca. 1.7 ± 0.5 nm determined from the TEM images (Figure 3 top). However, the envelope of the RDF, related to the size and shape of the objects, is not consistent with single-size spherical NPs but suggests a mixture with a majority of relatively small NPs (ca. 2.0 nm) and a small proportion of much larger ones (up to 4.0 nm). This may be related to a partial coalescence of the small NPs occurring in the agglomerated parts of the sample. For the Rh@L5 sample (Figure 3 bottom), the RDF indicated a much smaller coherence length of ca. 1.3 nm, which is in good agreement with the size measured by TEM (1.2 ± 0.6 nm). These particles, however, appeared less crystalline and displayed in reciprocal space a pattern much different from the fcc network, but in good agreement with the manganese- structure (ICSD 41775, PDF 33-887) already observed for very small Rh NPs 26,51 . These WAXS results suggest that ligand L5 has a stronger influence than ligand L2 on the crystallization state of the particles. Purified samples of L2-and L5-stabilized Rh NPs were also analyzed by Fourier transform infrared (FT-IR) spectroscopy before and after exposure to a CO atmosphere (1 bar; 48h) in the solid state. This investigation provided information on the coordination ability of CO molecules on the rhodium surface. Similar analyses were previously reported for Ru 52 and Rh 53 NPs. The obtained spectra (KBr pellets) are given in figures 4 and 5 for the L2-and L5stabilized Rh NPs, respectively.

In the case of Rh@L5 NPs, the FT-IR spectrum recorded after 48h of reaction with CO (Figure 4b) revealed the presence of complex bands in the region between 2050 cm -1 and 1550 cm -1 in addition to those already present before CO exposure (Figure 4a). The large and intense band at ca. 1990 cm -1 can be attributed to the stretching vibration of CO ligands adopting a terminal coordination mode. The higher frequency shoulder at ca. 2020 cm -1 is attributable to the CO ligands bonded in a geminal or multicarbonyl mode. Finally the band at 1800 cm -1 can be attributed to the standard bridging coordination mode. The high intensity of the band at 1990 cm -1 indicates that terminal CO molecules are abundant. It is generally accepted that this type of coordination occurs on the edges and corners of the particle surface while the bridging coordination mode happens on the extended faces of the particles. 54,55 Given that ligand L5 allows the CO coordination in all possible coordination modes and that the quantity of terminal CO is high, one can conclude that the metal surface remains relatively accessible with this ligand. Conversely, exposure of Rh@L2 to CO yields only one new and very low intensity band at 1924 cm -1 (Figure 5), suggesting a much more limited access for the CO molecules to the metal surface in this case. This effect may result from a tighter L2 coordination, probably induced by the bidentate character of L2 and possibly also by a larger covering of the particle surface by L2 relative to L5. In summary, the FT-IR data suggest that L2 ligand has a stronger steric effect on the RhNPs surface covering than L5 one. This can be attributed to the different coordination modes at the NP surface which can result from the different bonding mode (bidentate vs. monodentate).

Consequently, a greater fraction of Rh free sites are probably available in the case of Rh@L5 NPs than in the case of Rh@L2 NPs.

Rhodium nanoparticles stabilized with ligands L2 and L5 were then evaluated in the catalytic hydrogenation of styrene, in order to assess the ligand effect in this simple model reaction.

Hydrogenation reactions were performed at room temperature in an autoclave under a 30-bar H2 pressure, using isopropanol as a solvent. The activity was determined by GC analysis after 24 h of reaction. The obtained results are summarized in Table 1.

In all cases, the styrene substrate reached full conversion. In the absence of base and with a 0.7% or 0.2% molar percent of L2-stabilized Rh NPs (entries 1 and 2, respectively), ethylbenzene (P1) was the major product but formation of the complete hydrogenation product (ethylcyclohexane, P2) was also observed. The ethylcyclohexane yield increased from 2 to 8 % by increasing the catalyst loading from 0.2 to 0.7 mol. %. The putative product of selective arene hydrogenation, vinylcyclohexane, was never observed. These results suggest a sequential reaction with initial reduction of the vinyl function to ethylbenzene, followed by a slower arene ring reduction to ethylcyclohexane. To confirm this hypothesis, an independent ethyl benzene reduction was carried out under similar conditions (entry 5), yielding 6 % of P2 (similar to the results obtained in entry 1). These reactions were repeated in the presence of excess NaOMe (7-8 equivalents relative to Rh, entries 3, 4 and 6). While the styrene consumption remained complete, the formation of the total reduction product P2 was no longer observed. Thus, the base prevents the arene hydrogenation favoring the selective formation of ethylbenzene (P1). Similar results were observed when using the L5-stabilized Rh NPs. Reactions performed with 0.4 and 0.1 molar % of Rh (entries 7 and 8) gave 10% and 3% of the fully hydrogenated product (P2), and 90% and 97% of ethylbenzene (P1), respectively. The P2 yield increased at higher catalyst loading, as observed with the L2-stabilized Rh NPs. When NaOMe was added (entry 9 and 10) the P2 yield dramatically decreased.

A comparison of entries 5 and 11 shows that Rh@L5 has a greater activity than Rh@L2 for the arene ring hydrogenation, which can be correlated with the FT-IR suggestion of a more accessible surface for the former catalyst. Deeper mechanistic investigations will be necessary to understand how the presence of the base prevents the arene reduction.

Conclusion

Ferrocene-based ligands, namely ferrocenyldiphenylphosphines, have been here reported for the first time as stabilizers for the synthesis of metal nanoparticles. The Rh nanoparticles formed with these ligands have a well-controlled size and are quite small (1-2 nm range). The FT-IR study after exposure to a CO atmosphere has shown that the surface covering effect of the bidentate L2 ligand is greater than that of the monodentate L5 ligand. These NPs were found catalytically active in the hydrogenation of a model substrate, namely styrene, under very mild reaction conditions in spite of the surface coverage by the ligands. In addition, we have discovered that NaOMe greatly inhibits the hydrogenation of the phenyl ring. This inhibition is, to the best of our knowledge, unprecedented. No significant differences in selectivity were observed using mono-and bidentate ligands, but the NPs stabilized by the monodentate ligand L5 are more active than those stabilized by the bidentate ligand L2 in arene hydrogenation of ethylbenzene, in accordance with the lower efficiency of L5 to cover the nanoparticle surface. Deeper investigations will be necessary for a better understanding of the reaction mechanism and product selectivity.

All reagents were used without purification. Reagents and solvents were degassed before use according to a freeze-pump-thaw process. All operations concerning the synthesis of Rh precursor and Rh nanoparticles or the preparation of samples for characterization were carried out using standard Schlenk tubes and Fischer-Porter bottle techniques or in a glove-box (MBraun) under argon atmosphere.

Characterization techniques

Transmission electron microscopy observations low (TEM) and high resolution (HREM) were performed at the "Centre de microcaractérisation Raymond Castaing" in Toulouse (UMS-CNRS 3623). TEM grids were prepared by drop-casting of the crude colloidal suspension in CH2Cl2 onto a holey carbon-coated copper grid under argon atmosphere inside a glove-box. TEM images were obtained using either a JEOL 1100 electron microscope operating at 100 kV with point resolution of Wide-angle X-ray scattering (WAXS) measurements were performed at the CEMES-CNRS in Toulouse. Samples in the solid state were sealed in 1.5 mm diameter Lindemann glass capillaries.

The samples were irradiated with graphite monochromatized molybdenum Kα (0.071069 nm) radiation and the X-ray scattering intensity measurements were performed using a dedicated twoaxis diffractometer. Radial distribution functions (RDF) were obtained after Fourier transformation of the corrected and reduced data.

Infrared Spectroscopy

Fourier Transform Infrared (FT-IR) spectra of L2-and L5-stabilized RhNPs were recorded in the range 4000-400 cm -1 under argon atmosphere using a Perkin-Elmer GX2000 spectrometer installed inside a glovebox. First, the spectra of the purified NPs (KBr pellets prepared in the glove box)

were recorded. Second, powders of puridied nanoparticles were exposed under CO atmosphere (1 bar ; 48h) inside a Fischer porter reactor. After evacuation of excess CO, novel FT-IR spectra were recorded (KBr pellets).

ICP OES and elemental analysis

ICP OES analyses were performed at the "Service d'analyse élémentaire" of LCC-CNRS in Toulouse using Thermo Scientific iCAP 6300 DUO spectrometer with a 3-channel, 12-roller pump, and a 27.12 MHz solid state RF plasma generator. Mineralization of the NPs was performed in a refluxing aqua regia solution.

Synthesis methods

Synthesis of Rh( 3 -C3H5)3] complex

[Rh( 3 -C3H5)3] was prepared according to a slightly improved method from the literature. [START_REF] Beckhaus | Synthetic Methods of Organometallic and Inorganic Chemistry[END_REF] A THF solution (140 mL) of RhCl3.xH2O (1.5 g; 7.2 mmol) was reacted with allylmagnesiumchloride (36 mL of a THF solution at 2 mol.L -1 ; 71.6 mmol) under inert atmosphere and vigorous stirring at -10°C for 1 h and then at 10°C for 16 h. During this time, the solution color changed from red to yellow. Solvent removal under reduced pressure led to a solid, which was dissolved in pentane leading to a yellow solution. After filtration through Celite, pentane was evaporated under reduced pressure and the resulting yellow solid was purified by sublimation at 40°C. Yellow crystals (660 mg, 41% yield) of [Rh( 3 -C3H5)3] were collected and stored under an argon atmosphere at low temperature. [START_REF]Physics and Chemistry of Metal Cluster Compounds[END_REF] H NMR and elemental data were found identical to those from the literature.

Synthesis of phosphino-ferrocenyl ligands

Ligands L1 and L2, 49 L3 50 and L4 44c were prepared according to published methods. Racemic (2-(diphenylthiophosphino)ferrocenyl)methanol (1) was synthesized by optimizing a literature procedure from N,N-dimethylaminomethylferrocene. [START_REF] Hayashi | [END_REF][START_REF] Manoury | the Series: Catalysis by Metal Complexes (CMCO)[END_REF] Synthesis of (2-(diphenylthiophosphino)ferrocenyl)methane 2.

(2-(Diphenylthiophosphino)ferrocenyl)methane 2 was synthesized from alcohol 1 by adaptation of a published procedure: [START_REF] Routaboul | [END_REF] In a Schlenk tube under argon were dissolved 500 mg of alcohol 1 (1.16 mmol) in 6 mL of anhydrous THF. A BH3.THF 1 M solution in THF (5 mL, 5 mmol) was then added by syringe. The reaction mixture was stirred at THF reflux for 70 h. After cooling to r.t., the Schlenk tube was opened to air and water was slowly added (warning: an abundant gas evolution occurred), leading to the formation of an orange precipitate. After THF evaporation on the vacuum line, 10 mL of dichloromethane was added for full dissolution of the orange precipitate. The remaining aqueous solution was then extracted two times with 5 mL of dichoromethane. The combined organic solutions were dried on anhydrous sodium sulfate and purified, after evaporation, by flash chromatography on silicagel using pentane/dichloromethane (50/50, v/v) as eluent to yield, after solvent evaporation, 470 mg (98% yield) of compound 2 as a yellow solid.

Synthesis of ferrocenylphosphine-capped Rh NPs

In a general procedure, 150 mg of [Rh( 3 -C3H5)3] (0.66 mmol) were introduced under an argon atmosphere in a Fischer-Porter After 0.5 h under vacuum at r.t. and cooling to 193 K, a CH2Cl2 solution (150 mL) containing 0.2 molar equivalent of the chosen ligand (L1, L2, L3 or L4) or 0.4 molar equivalent of L5 was added. The reaction medium was then allowed to warm up to r.t. before adding dihydrogen (3 bar) and then introduced in an oil bath preheated at 313 K for 16 h.

After cooling to r.t., the remaining dihydrogen was vented and the volume of the solution was reduced to approximately 10 mL under reduced pressure before its transfer onto a solution of deoxygenated cold pentane (100 mL), generating the product as a black precipitate. After filtration and pentane washings the RhNPs were collected under the form of a fine powder. ICP OES analysis: for L2-stabilized RhNPs : Rh: 29.3 wt% ; for L5-stabilized RhNPs: Rh: 22.7 wt%.

Catalysis studies

In a glove-box, a Schlenk tube was charged with 12 mg of rhodium NPs. Then, 15 mL of dry degassed isopropanol and 500 µL (2.2 mmol) of dodecane (internal reference) were added to generate a stock of colloidal suspension of the catalyst. A portion of the catalyst suspension (1.5 mL) was transferred to a 5-mL glass vial (pre-charged with the desired amount of styrene or ethylbenzene substrate and NaOMe), which was then placed into a stainless steel autoclave. The reaction vessel was pressurized with 30 bar of H2 and stirred with a magnetic bar for 24 h at r.t. .

The reaction was stopped by release of the H2 pressure and quenching with 2 mL of CH2Cl2, followed by filtration through silica gel. The product was finally analyzed by GC (Shimadzu 2014) using a SLB TM 5ms Capillary Column.

Figure 1 :

 1 Figure 1: (a) TEM images and corresponding size distributions of Rh NPs stabilized with L1, (b) L2, (c) L3, (d) L4, (e) L5.

Figure 2 :

 2 Figure 2: HREM images of L2-stabilized Rh NPs (top) and L5-stabilized Rh NPs (bottom) at different magnifications.

Figure 3 :

 3 Figure 3: WAXS measurements (solid line) for L2-stabilized Rh NPs (top) and L5-stabilized Rh NPs (bottom) in reciprocal space (left, with reference diagrams for Rh in the fcc and manganese- structures, dashed line) and real space (right)

Figure 4 :

 4 Figure 4: (a) FT-IR spectrum of L5-stabilized Rh NPs (b) FT-IR spectrum of L5-stabilized Rh NPs after 48h under 1 bar of CO.

Figure 5 :

 5 Figure 5: (a) FT-IR spectrum of purified L2-stabilized Rh NPs (b) FT-IR spectrum of L2stabilized Rh NPs after 48h under 1 bar of CO.

4. 5 Å

 5 or a JEOL JEM 1400 operating at 120 kV with point resolution of 4.5 Å. HTEM observations were carried out with a JEOL JEM 2010 electron microscope working at 200 kV with a point resolution of 2.3 Å and equipped with EDS (Energy Dispersive X-ray Spectroscopy) mode. Statistical size distributions were built by analyzing TEM images with Image J tool software, counting mean diameters of > 200 non-touching particles and assuming a spherical form. Size distributions are quoted as the mean diameter ± the double of the standard deviation (σ) which corresponds to a 95% confidence level. Fourier Fast Transform treatments were carried out with Digital Micrograph Software (Version 1.80.70).

  C{ 1 H} NMR (100 MHz, CDCl3), δ: 14.0 (d, JC-P = 10.0 Hz, CH 3 ), 68.3 (s, subst Cp), 69.7 (s, Cp), 70.4 (d, JC-P = 3.6 Hz, subst Cp), 72.3 (d, JC-P = 3.8 Hz, subst Cp), 75.7 (d, JC-P = 4.3 Hz, quat Cp), 89.4 (d, JC-P = 25.0 Hz, quat Cp), 127.8 (s, Ar), 128.14 (d, JC-P = 7.7Hz, Ar), 128.17 (d, JC-P = 5.9 Hz, Ar), 129.0 (s, Ar), 132.2 (d, JC-P = 18.0 Hz, Ar), 134.9 (d, JC-P = 20.7 Hz, Ar), 137.6 (d, JC-P = 8.8Hz, quat Ar), 139.7 (d, JC-P = 10.9 Hz, quat Ar).

Table 1 :

 1 Styrene and ethyl-benzene hydrogenation using L2-and L5-stabilized Rh NPs under 30 bar of H2 at R.T.

	a: substrate to rhodium molar ratio. b: mol. %.

  Cp), 4.20 (m, 1H, subst Cp), 4.38 (m, 1H, subst Cp), 7.2-7.15 (m, 2H, Ar), 7.3-7.25 (m, 3H, Ar), 7.35-7.45 (m, 3H, Ar), 7.55-7.60 (m, 2H, Ar).13

	P NMR (162 MHz, CDCl3), δ: -21.7. 1 H NMR (400 MHz, CDCl3), δ: 2.08 (s, 3H, CH 3 ), 3.62 (m,
	1H, subst Cp), 4.07 (s, 5H,
	31 P NMR (162 MHz, CDCl3), δ: 42.3. 1 H NMR (400 MHz, CDCl3), δ: 2.17 (s, 3H, CH 3 ), 3.67 (m,
	1H, subst Cp), 4.22 (m, 1H, subst Cp), 4.34 (s, 5H, Cp), 4.43 (m, 1H, subst Cp), 7.35-7.55 (m, 6H,
	Ar), 7.63 (m, 2H, Ar), 7.81 (m, 2H, Ar). 13 C{ 1 H} NMR (100 MHz, CDCl3), δ: 14.3 (CH 3 ), 68.2 (d,
	Synthesis of (2-(diphenylphosphino)ferrocenyl)methane L5.
	(2-(Diphenylphosphino)ferrocenyl)methane L5 was synthesized from compound 2 by deprotection
	of the phosphino group following a published procedure: 60 in a Schlenk tube under argon were
	dissolved 470 mg of compound 2 (1.13 mmol) in 40 mL of toluene. P(NMe2)3 (1 mL, 5.5 mmol, 4.9
	eq) was then added by syringe. The reaction mixture was stirred at toluene reflux overnight and
	purified by flash chromatography on silicagel under argon using pentane/dichloromethane (50/50,
	v/v) as eluent to yield, after solvent evaporation, 421 mg (97% yield) of compound 3 as a yellow
	solid.

JC-P = 10.5 Hz, subst Cp), 70.7 (s, Cp), 74.1 (d, JC-P = 13.2 Hz, subst Cp), 74.34 (d, JC-P = 96.1 Hz, quat Cp), 74.36 (d, JC-P = 9.9 Hz, subst Cp), 88.7 (d, JC-P = 12.4 Hz, quat Cp), 128.0 (d, JC-P = 12.4 Hz, Ar), 128.3 (d, JC-P = 12.4 Hz, Ar), 131.15 (d, JC-P = 2.6 Hz, Ar), 131.18 (d, JC-P = 2.6 Hz, Ar), 131.9 (d, JC-P = 10.9 Hz, Ar), 132.0 (d, JC-P = 10.5 Hz, Ar), 133.7 (d, JC-P = 86.0Hz, quat Ar), 134.8 (d, JC-P = 86.4 Hz, quat Ar). HRMS (ESI+) m/e: 417.0522 (M+H, 100%; calcd. for C23H22FePS: 417.0529).

Acknowledgements

Financial support by the Centre National de la Recherche (CNRS) and the University Paul Sabatier -Toulouse are gratefully acknowledged. M. Ibrahim thanks the University Paul Sabatier -Toulouse for a MESR PhD grant. M.-M. Wei thanks the CNRS and the Midi-Pyrénées Region for a PhD grant.

Conflict of interest

There are no conflicts of interest to declare.

Experimental part 1. General and Materials

Chemicals were purchased as follows: dichoromethane (CH2Cl2) and pentane from SDS; RhCl3.xH2O from Johnson Matthey; [CH2=CHCH2MgCl], P(NMe2)3 and BH3•THF (1 M solution in THF) from Sigma Aldrich ; H2 from Air Liquide. Solvents were purified before use by filtration on an adequate column of a MBraun purification apparatus and handled under argon atmosphere.