R Nakhoul

Pierre Montmitonnet
email: pierre.montmitonnet@mines-paristech.fr

N Legrand

Manifested flatness defect prediction in cold rolling of thin strips

Keywords: Rolling, Thin Strips, Finite Element Method, Friction, Flatness Defect, Sheet Buckling, Residual Stress

Flatness defects in thin strip cold rolling are a consequence of roll thermo-elastic deformation, from which heterogeneous strip plastic deformation results. When flatness defects manifest on line, buckling reorganizes the stress field in the pre-and post-bite areas. Comparison with flatness roll measurement requires this effect to be taken into account. A coupled Finite Element Method (FEM) approach is used here to compute stresses and strains in-bite as well as out-of-bite. The detection of buckled (non-flat) areas is demonstrated for a very thin strip cold rolling case ("double reduction"). The model is then applied here to two questions, namely the impact on flatness of the heterogeneous temperature field and the effect of friction on optimal setting of a flatness actuator, Work Roll Bending.

INTRODUCTION

In strip rolling, due to rolls elastic deformation combined with roll grinding crown, thermal crown and wear, the roll gap is not uniform in the width direction, so that strip reduction depends on the width coordinate y. Heterogeneous reduction induces residual stresses, and compressive residual stresses induce buckling (out of plane waviness) if the strip is thin hence easy to bend. Most of the time, the strip tension is high enough to keep the in-plane stress tensile everywhere and avoid buckling on-line. But once the strip tension is relaxed, or the sheet is cut, residual stresses may induce distortions called latent flatness defects, i.e. non developable, out-of-plane displacements forming periodic waves (centre waves, edge waves…). However, in some severe operations such as thin strip cold rolling, strain and stress heterogeneity is exacerbated and on-line strip tension may be insufficient to prevent -or hide -waviness. In such a case, manifested flatness defects may show on-line (Fig. 1).

If defects are only latent, the stress field computed beyond the bite by e.g. a 3D FEM should be correct and consistent with measurements. Buckling upon unloading the strip tension can then be dealt with in an uncoupled way, transferring the stress pattern computed by the rolling model into a buckling model (e.g. shell FEM). The latter then tells if the critical conditions for buckling are met, and in this case computes post-buckling mode, wavelength and amplitude.

On the other hand, it is well known that buckling completely transforms stress fields since the local stiffness of a buckled sheet drops dramatically. Therefore, if defects are manifested online as is most of the time the case for very thin sheets, the solution coming out of the 3D FEM of rolling is completely irrelevant outside of the roll bite. Abdelkhalek et al. [START_REF] Abdelkhalek | Coupled approach for flatness predictions in thin strip cold rolling[END_REF] have developed a fully coupled buckling model in an implicit finite element software devoted to rolling processes. They showed that it makes it possible to correct the out-of-bite stress fields and to come out with values consistent with the σ xx (y) profile measurement by a flatness roll.

The question of the impact of this stress reorganization on the strains and stresses within the bite has also been raised in the same paper. The conclusion is that to the possible exception of temper-rolling, the feedback effect is very small: the out-of-bite stress rearrangement is screened by the immediate bite exit area where complex velocity profile changes take place elastically. In the present paper, it is first recalled how the coupled buckling criterion detects buckled areas. Then applications to practical questions are examined, such as the setting of a flatness actuator (work roll bending), the effect of friction on the stress profile and on the "optimal" work roll bending setting, and the effect of a usual model simplification, isothermal versus thermally coupled.

BRIEF LITERATURE SURVEY

A number of uncoupled models of flatness defects have been presented in the literature. They are based on the standard solution by Timoshenko and Woinovsky-Krieger of Von Karman's equations for a rectangular plate under compressive stress on simply supported edges.

Timoshenko and Woinovsky-Krieger [START_REF] Timoshenko | Theory of plates and shells[END_REF] assumed homogeneous stress and sinusoidal waves in both x and y directions. Bush et al. for instance [START_REF] Bush | Stress levels for elastic buckling of rolled strip and plate[END_REF] extended this solution to the complex fields of strip rolling, allowing only right/left symmetric waves. But as real waves are often anti-symmetric, Fischer and coworkers [START_REF] Fisher | Buckling phenomena related to rolling and levelling of sheet metal[END_REF] chose for the transverse shape either an odd or an even polynomial in y, while keeping the waves sinusoidal in the rolling direction.

Only the longitudinal stress component (σ xx ) was used in these papers. Yet, it is easy to see that not all flatness defects can be handled with this assumption (herringbone defects for instance are excluded). Moreover, the type of defect (centre or edge wave) was determined by these authors a priori from the shape of the stress profile, and ad hoc clamping conditions were applied to get it in the end, to avoid managing too many modes with very similar energies. On the contrary, Yukawa and coworkers [START_REF] Yukawa | Numerical analysis of the shape of rolled strip[END_REF] addressed more complete configurations using shell FEM, again with a non-coupled approach. The bifurcation point was detected as the load parameter making the second variation of the total elastic strain energy (i.e. the stiffness matrix) non definite positive. Post-buckling was computed by introducing a small defect corresponding to the mode, the load step being controlled by a modified Riks method. This is the standard approach of instabilities in structural mechanics. Abdelkhalek et al. [6] also proposed a decoupled approach, importing a FEM-computed stress field in a shell FEM for more generality. In the bifurcation detection stage, they maximized the load step and decreased the computational burden by managing the loading steps by the Asymptotic Numerical Method (ANM) developed initially by Zahrouni et al. [START_REF] Zahrouni | Computing finite rotations of shells by an Asymptotic-Numerical Method[END_REF] and improved by Boutyour et al. [START_REF] Boutyour | Bifurcation points and bifurcated branches by an Asymptotic Numerical Method and Padé approximants[END_REF]. ANM consists in developing the solution (displacement and load parameter) in a power series with respect to a step length parameter, up to a truncation order p. This series is inserted in the non linear equations of the problem, the different orders are identified, giving p linear systems with the same stiffness matrix, the right-hand sides of which are computed from the solutions at lower orders. With this approach, they modelled both on-line shape (under strip tension) and relaxed shape (strip tensions relieved), whereas the previous papers dealt only with the latter.

Counhaye [START_REF] Counhaye | Modélisation et contrôle industriel de la géométrie des aciers laminés à froid[END_REF] was the first who questioned the decoupling of the rolling and the buckling model. Indeed, buckling sets a limit to the allowable compressive stresses: wherever buckling occurs, the loss of rigidity changes the stress locally and global changes occur as the stress field is re-equilibrated. This thorough stress rearrangement on the one hand questions the stress fields obtained from the rolling models, on the other hand might have an impact even on in-bite stress and strain fields, or on roll deformation in strongly coupled cases. Therefore, Counhaye used an approach similar to the one proposed by Roddeman and coworkers [START_REF] Roddeman | The Wrinkling of Thin Membranes: Part I-Theory[END_REF], considering buckling as one more strain rate component, present only in the out-of-bite areas, to be added to the elastic / plastic strain rate splitting. This extra strain rate may be interpreted as the local shortening of a material segment as it becomes wavy due to buckling.

The model presented in [START_REF] Abdelkhalek | Coupled approach for flatness predictions in thin strip cold rolling[END_REF] and used here is an adaptation of Counhaye's model.

COUPLED ROLL -STRIP ROLLING MODEL WITH EMBEDDED BUCKLING

The model has been presented in details elsewhere previously, so that only general lines will be recalled hereafter.

Sheet Rolling model Lam3/Tec3

The rolling model is Lam3/Tec3, a 3D strip / roll stack deformation software developed by Hacquin and coworkers [START_REF] Hacquin | Modélisation thermomécanique tridimensionnelle du laminage -Couplage bande / cylindres (3D thermomechanical model of rolling -Strip-roll coupling[END_REF]. The strip elastic-viscoplastic (EVP) deformation is dealt with by an implicit FEM with a velocity formulation using P1-discretisation on hexahedra. A steady state formulation based on streamlines is implemented. This requires integrating elasticviscoplastic constitutive equations along streamlines determined from the velocity field. In [START_REF] Hacquin | A steady state thermo-elastoviscoplastic finite element model of rolling with coupled thermo-elastic roll deformation[END_REF], they described the heterogeneous time step strategy used for this purpose, called ELDTH (Eulerian-Lagrangian with DT Heterogeneous). In [START_REF] Hacquin | A 3D semi-analytical model of rolling stand deformation with finite element validation[END_REF], they focused on the semianalytical thermo-elastic roll deformation model. Inspired by the work of Berger et al. [START_REF] Berger | Die elastische Verformung der Walzen von Vierwalzengerüsten[END_REF] on the deformation of a cylinder under complex loadings, Hacquin modelled the roll stack thermo-elastic deformation based on Timoshenko beam theory (bending component), Boussinesq solution of a half-space under general loading (flattening part), and Hertz contact mechanics between work rolls and back-up rolls. These analytical approaches were found insufficient near the roll barrel edge, therefore edge effects were corrected by analytical functions calibrated using an extensive FEM roll deformation campaign under various geometries and loading conditions [START_REF] Hacquin | Modélisation thermomécanique tridimensionnelle du laminage -Couplage bande / cylindres (3D thermomechanical model of rolling -Strip-roll coupling[END_REF][START_REF] Hacquin | A 3D semi-analytical model of rolling stand deformation with finite element validation[END_REF].

Equations are discretized by an influence function method (IFM). Writing the global

equilibrium of the roll stack under the combined effects of the roll -strip contact stresses, external (roll load) and internal (e.g. roll bending) forces, a system of equations is found, the unknowns of which are the roll rigid body displacement, contact line displacement field and contact pressure profiles at work roll / back-up roll and work roll / strip contact. This system, non-linear due to unknown contact lines, is solved by Newton-Raphson method.

Both strip and roll temperatures are computed and coupled with mechanics. The strip temperature model uses a straightforward 3D SUPG (Streamline-Upwind Petrov Galerkin) Finite Element scheme. The roll temperature (and roll crown) model is a combination of 2D (r,θ) and (r,z) FEM coordinated by another influence function technique. Details and validation can be found in [START_REF] Hacquin | Modélisation thermomécanique tridimensionnelle du laminage -Couplage bande / cylindres (3D thermomechanical model of rolling -Strip-roll coupling[END_REF].

In [START_REF] Hacquin | Experimental validation of a Rolling Stand Elastic Deformation Model[END_REF], the whole roll stack deformation model was validated under isothermal conditions against experimental measurements performed under diverse loading conditions. Later on in [START_REF] Montmitonnet | Comparison of profile prediction models for strip rolling[END_REF], a successful numerical comparison was carried out with the model of Yanagimoto and Kiuchi [START_REF] Yanagimoto | 3D simulation system for coupled elastic / rigid plastic deformation of rolls and workpieces in strip rolling processes[END_REF], which uses a beam bending model for roll bending and 3D FEM for flattening, and with the model of Kim et al. [START_REF] Kim | An integrated FE process model for the prediction of strip profile in flat rolling[END_REF] where full 3D elastic FEM is used for roll stack deformation and strip deformation altogether; all three models were compared with the experiments reported by [START_REF] Yanagimoto | 3D simulation system for coupled elastic / rigid plastic deformation of rolls and workpieces in strip rolling processes[END_REF]. Therefore, this roll stack deformation model, essential for strip flatness, can be used with confidence for the purpose of the present paper.

Simple buckling model embedded in Lam3/Tec3

Counhaye's model [START_REF] Counhaye | Modélisation et contrôle industriel de la géométrie des aciers laminés à froid[END_REF] has been implemented in Lam3/Tec3 [START_REF] Abdelkhalek | Coupled approach for flatness predictions in thin strip cold rolling[END_REF]. Initially proposed within the membrane theory framework in [START_REF] Roddeman | The Wrinkling of Thin Membranes: Part I-Theory[END_REF], it forbids the appearance of a negative stress: everywhere compression is about to occur, the structure buckles, bringing the stress back to almost zero by providing a stress-free alternative to elastic shortening of a material line. The following critical conditions are introduced: ), the membrane is stiff; if the stress becomes negative, it gets slack and in fact, the corresponding stress is put to 0 (direction 1 n 

0 . . 0 . . 0 . . 2 1 2 2 1 1 = > = n n n n n n       σ σ σ (1)
). The essence of the model consists in determining an extra deformation which elastically brings the stress in the buckled direction back to 0. It may be interpreted as the shortening of a material line due to buckling of the structure. This is more or less analogous to elastic-plastic decomposition, but is activated only out of the roll bite, i.e.

where buckling is liable to occur:

bu el ε ε ε ∆ + ∆ = ∆ (2)
∆ε el is the elastic and ∆ε bu is the "buckling strain" increment. Plane stress is assumed outside of the roll bite. The extra deformation representing buckling is computed in the principal axes, and then transported to the reference frame. Let λ i , i = I, II, be the principal components of this extra strain. It is deduced from σ i , i = I, II as follows:

II I i E c i i , = - = σ σ λ (3)
Moving back to the reference frame (Oxyz), the buckling strain increment is added to the global strain increment. In the following equation, u and v are the two in-plane incremental displacements, α is the angle between principal and reference frames in the plane of the strip, ν is Poisson's ratio and E is Young's modulus:

( )

yy xx zz I II xy I II yy II I xx y v x u y v x u ε ε ν ν ε α α λ λ ε α λ α λ ε α λ α λ ε ∆ + ∆ - - = ∆ - +         ∂ ∂ + ∂ ∂ = ∆ + + ∂ ∂ = ∆ + + ∂ ∂ = ∆ 1 sin . cos 2 1 sin . cos . sin . cos . 2 2 2 2 (4)
This strain increment replaces the standard one fed into the module which solves the constitutive differential equations.

ROLLING OPERATION INVESTIGATED

Definition of rolls, material, strip dimensions

All examples and parametric studies described hereafter refer to the same rolling pass, the last stand of a tinplate sheet mill, with very low thickness. All the characteristics are given in

- - - × + = ε ε σ (5)
σ 0 is the yield stress in tension, ε is the equivalent plastic strain. Coulomb friction is assumed, the friction coefficient is µ = 0.03 unless otherwise stated.

Impact of buckling being accounted for

With a reservation for very small reduction (temper-rolling), it has been found that this rearrangement of out-of-bite stress has significant impact neither on strain and stress fields in the bite, nor on the roll loads and deformation [START_REF] Abdelkhalek | Coupled approach for flatness predictions in thin strip cold rolling[END_REF]. Neglecting buckling is therefore licit if only the behaviour in the bite is to be studied. But the "residual" stress field in the strip, under these conditions, is very far from the flatness roll measurement, as shown in Fig. 4. "Residual" means here post-bite, on-line under tension.

The measurements (dots) show an almost flat stress profile (at the scale of the drawing). FEM results without buckling (dashed line) give an enormous compressive stress at the edge, due to a very large elongation gradient connected with roll flattening and the edge drop defect (see Fig. 3). This very high compressive stress is compensated for by a high tensile stress in the centre, since the resultant must be equal to the sheet tension force (here given by the tension stress of 100 MPa, Table 1). Also to be noted is a stress peak at y ≈ 415 mm (i.e. ca.

10 mm from strip edge). It is much more conspicuous when buckling relaxes partly the stress (as in figure 8 or 12); the same location and amplitude is found whatever the mesh size; it is in fact exactly at the location where the edge-drop starts, i.e. the transition between the flattened and free zones of the roll profile, see figure 2b.

When buckling is accounted for, i.e. equations (1-4) are activated, the computed stress profile (full line) comes close to the experiments. More details can be found in the authors' previous paper [START_REF] Abdelkhalek | Coupled approach for flatness predictions in thin strip cold rolling[END_REF]. Therefore, for out-of-bite stress prediction and comparison with stress-meter measurements, accounting for buckling is essential, at least for thin strips. sides of the roll bite. Hence λ I ≠ 0 corresponds to a x-wave (wavy edge) and λ ΙΙ ≠ 0 to y-wave, the so-called "towel effect". Whereas the former is moving with the strip and constitutes a real flatness defect, the latter is a stationary phenomenon which can be observed sometimes at the exit of the bite. It is due to a flexible membrane being submitted to a tension under a nonstrictly homogeneous tension stress (e.g. due to the clamping effect of the roll bite). It has no practical consequences for the strip most of the time, except if it becomes so large that double thickness rolling occurs; this improbable case is not addressed here.

Both types of waves are detected by the model. It should be emphasized, though, that the present buckling model is embedded in a steady state rolling model, so that it represents a moving defect in a stationary context. It can therefore detect wave formation, but not characterize its geometry.

FIRST APPLICATION: IMPACT OF TEMPERATURE FIELD ON FLATNESS

Most of the time, thermal effects are neglected in steel cold rolling, and flatness defects are studied under isothermal conditions. Here, with all three couplings available: strip -roll, mechanical -thermal, and buckling, many effects may be investigated.

Indeed, at bite exit, the temperature is much larger at the edge than in the centre (Fig. 7), about 438 K (165°C) versus 383 K (110°C). This is due to larger reduction and plastic heating (see Fig. 3). Due to the corresponding differential dilatation, the post-bite stress pattern could be affected. A series of simulations has therefore been run, with roll temperature calculated at steady state (i.e. after a long rolling time). The isothermal case is taken as a reference. In the second case, the strip is allowed to cool after bite exit under the effect of the strip cooling system, with heat transfer coefficient (HTC) H cool = 5 kW.m -2 .K -1 ; strip -roll contact is kept adiabatic (H roll = 0). Temperature increases in the roll bite and slowly decays afterwards (remember 1 m ≈ 0.05 s). In the third case, the roll -strip interface is represented by H roll = 100 kW.m -2 .K -1 . The temperature increase is less in the bite, cooling is similar to the previous case. Finally, in the fourth case, H cool = 5 kW.m -2 .K -1 , H roll = 100 kW.m -2 .K -1 , but rolls are moreover cooled efficiently by water sprays, so that strip temperatures are significantly lower.

Yet, the temperature difference between strip edge and centre is very similar in all 3 cases, 313 to 323 K (40°C to 50°C). Fig. 8 shows that these temperature differences, even in the isothermal case (T = 298 K or 25°C everywhere), have negligible impact on the stress pattern once relaxed by edge buckling. workers have proved in [START_REF] Jiang | Finite element simulation of cold rolling of thin strip[END_REF] that friction impacts strip profile after cold rolling and demonstrated the effect of a parabolic variation of friction in the transverse direction in [START_REF] Jiang | Modelling of the effect of friction on cold strip rolling[END_REF].

It is therefore important to quantify the effects of such friction time variations. Friction is assumed homogeneous here. Work roll bending force is equal to 482 kN

Impact of friction on roll deformation

Impact of work roll bending force at fixed µ

Now, the work roll bending force is varied at constant µ = 0.025 (Fig. 13). The very high work roll bending force (900 kN) again gives a low, yet positive, stress in the centre, which will probably result in a wavy centre, at least after tension is cancelled. The most flat stress profile turns out to be for a work roll bending force equal to 350 kN.

Finally, Fig. 14 summarizes the bending force found "optimal" for each value of the friction coefficient, i.e. giving the most flat stress profile "by eye". This graph gives an idea of how to

preset the work roll bending force as a function of varying friction, in the present rolling operation of Table 1.

Fig. 13 The σ xx (y) stress profile for varying bending force and fixed friction coefficient µ = 0.025 Fig. 14 The relationship between friction and optimal bending force

CONCLUSION

The main outcome of the present paper is to show that a simple buckling model can be coupled to a finite element rolling model and applied to practical questions concerning flatness in cold strip rolling. In-bite stress and strain are not affected by post-buckling stress rearrangement, but the residual stress is, in case of manifested defects expected with very thin strips.

The results presented above suggest that for cold rolling, the heterogeneous temperature field does not change significantly the stress pattern, which is dominated by the rearrangement of the stress field as the elastic strain recovery takes place in the immediate post-bite area. This analysis contributes to shed light on which features of the rolling process may be neglected or not, depending on the objective: in bite stress and strain, roll deformation and strip crown do not need particularly buckling to be coupled; the residual state of the strip is very sensitive to buckling whenever it takes place on-line, but it is not sensitive to thermal transfer.

It has been shown also here that friction has an impact on reduction profile, in particular on the edge thickness, through roll load distribution and roll deformation. This strongly changes residual stress distribution when thin strips are dealt with. For a given work roll bending force, just changing friction may turn the strip shape from flat to wavy edge to wavy centre.

Conversely, for each level of friction, an "optimal" work roll bending force can be proposed based on results of the model.

The buckling model presented in section 3.2 deals with buckling at material point (or finite element) level, which is questionable; but it is strongly coupled at constitutive equations level, which Abdelkhalek et al. have shown to be important for precise residual stress prediction [START_REF] Abdelkhalek | Manifested flatness predictions in thin strip cold rolling using a general rolling FEM model[END_REF]. It cannot predict the geometry of defects quantitatively. These are serious limitations.

Figure captions

Fig. 1

 1 Fig. 1 Schematic view of flatness defects. Note reference frame : x = Rolling Direction (RD),

  the principal directions of the Cauchy stress tensor σ in the buckled structure (hence the third equation). Principal directions are determined first, and the angle they make with (x,y) is called α. When a tension is applied in a direction (here 2 n 

3 4. 2 Fig. 2 Fig. 3

 3223 Fig. 2 Coupled deformation of strip and roll:

Fig. 4

 4 Fig. 4 Comparison of the experimental and numerical σ xx (y) stress profiles, ≈ 1 m after roll

Fig. 5 σ

 5 Fig. 5 σ xx stress maps at the edge, in the vicinity of the roll bite. a -without buckling

Fig. 5 Fig. 6 «

 56 Fig.5shows the σ xx stress map at the edge, slightly before and slightly after the roll bite, i.e.

Fig. 7 Fig. 8

 78 Fig. 7 Effect of thermal boundary conditions on central and edge longitudinal temperature

Figs. 9 -

 9 Figs. 9-12 illustrate a parametric study whereby friction varies between µ = 0.01 (slight

Fig. 9

 9 Fig. 9 Effect of varying friction on the work rolls load profile F(y). Work roll bending force

Fig. 10

 10 Fig. 10 Effect of varying friction on the deformed roll generator shape (bottom). Work roll

Fig. 1 Fig. 2 Fig. 3 Fig. 5 σFig. 6 «

 12356 Fig. 1 Schematic view of flatness defects. Note reference frame : x = Rolling Direction

Fig. 9

 9 Fig. 9 Effect of varying friction on the work rolls load profile F(y). Work roll bending

Table 1 ,

 1 together with the stress-strain curve:

	0	( 470	5 .	. 4 . 175	)	1 (	. 45 . 0	exp(	9 . 8	))	175	(MPa)

Table 1

 1 Characteristics of the rolling pass investigated

	Strip width	855 mm
	Strip thickness	0.355 mm
	(strip crown)	(4.81 %)
	Exit thickness	0.252 mm
	Rolling speed	22 m/s
	Front /	100 MPa /
	back tension	170 MPa
	Type of mill	4-high
	Work rolls diameter	555mm
	Backup rolls diameter	1300mm
	Work rolls length	1400mm
	Backup rolls length	1295mm
	Work rolls crown	0.0322 %
	Backup rolls crown	No crown
	Work roll bending
	force / position

Table caption Table 1

 caption1 Characteristics of the rolling pass investigated

Work is in progress to try and implement more physically founded models which give quantitative information on defect geometry. These efforts will be reported somewhere else.

ACKNOWLEDGEMENT

The authors wish to thank the French National Research Agency (ANR) for its financial support under Contract ANR -2010-RMNP-019-01, and the partners of the Platform Project (ArcelorMittal Maizières research, Constellium CRV, CNRS, CEA, Université de Lorraine, INSA Lyon and ARMINES) for authorization to publish this work. Particular thanks are due to M. Potier-Ferry and H. Zahrouni, from LEM3 laboratory, for most helpful discussions.


Kokou Adokanou
email: akokou.adokanou@mines-paristech.fr

Karim Inal

Pierre Montmitonnet

Frédéric Courtade

Barbara Bonnet

Thermo-mechanical analysis of GaAs devices under temperature-humidity-bias testing

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

For Space applications, devices are generally encapsulated by hermetic packaging to prevent any pollution or moisture penetration during the storage on ground of the equipment that could damage the microelectronic devices. At die level, life tests are performed to validate the use of a chip technology.

In the case of non-hermetic packaging and in addition to standard life tests, damp humidity testing at 85°C and 85% RH for several hours should be performed to access the impact of moisture on the device reliability. This test can be performed with or without DC bias applied to the devices and bias may be direct (a current flows through gate) or reverse (no current).

A usual acceleration law for THB tests was published by Hallberg and Peck [START_REF] Hallberg | Recent humidity accelerations, a base for testing standards[END_REF] to estimate the Mean Time To Failure (MTTF) of the devices. This empirical model takes into account the relative humidity, the bias (applied voltage) and the temperature change by using analytic formulation and experimental results. It has been validated for some silicon technologies but is not adapted for technologies based on GaAs even if others authors have enriched it several times.

Unfortunately this model does not quantify other parameters such as the impact of temperature change on the residual stress in the die and possible microstructure change. Because of CTE mismatch between constituent layers, heating the device during THB (or during use) superimposes thermo-elastic stresses to the residual stresses initially present. The latter are mainly due to the same kind of thermo-elastic stress built up when the device was cooled from fabrication temperature to room temperature. Indeed, the stress at any time is equal to the thermo-elastic stress corresponding to the difference between manufacturing temperature and the temperature at that time. In the present paper, only the thermo-elastic stress due to the heating in THB test is considered.

A vast literature exists on this topic for multilayer structures. While the stress and reliability issues have been widely studied in metallic thin films on substrate and in adhesively bonded assemblies, the thermo-mechanical behavior of devices under aging tests has not been fully studied and understood yet. Most studies focus on the stresses induced during fabrication due to the variation of temperature, the mismatch in the CTE and depending on materials thickness and properties.

This paper presents an evaluation of thermomechanical stresses induced in the active part of nonhermetically packaged Monolithic Microwave Integrated Circuit (MMIC) devices dedicated for Space applications by the temperature raise during THB 85/85 test. First, the thicknesses of the different layers of materials constituting the MMIC must be evaluated, including dielectric layers and gate metallization of the transistors. Results from X-ray Photoelectron Spectroscopy (XPS) analysis and FIB-SEM cross sections (Focused Ion Beam and Scanning Electron Spectroscopy) are presented in section 2. Second, analytical models are summarized in section 3 and the numerical model based on finite element method (FEM) are described in section 4. In section 5, the models developed are used to quantify thermo-mechanical stresses induced by heating for a specific multilayer stack (GaAs/Ti/Al/Au/SiNx/SiOy/SiNx) of the gate finger area. Results are then correlated with cracks studied in [START_REF] Ben Naceur | Evaluation of quasihermetic packaging solutions for active microwave devices and space applications[END_REF] on non-hermetic tested devices with reverse bias.

Layer thickness measurement

The GaAs MMIC used in this study includes two 0,18µm pseudomorphic High Electron Mobility Transistors (pHEMT) and a SiNx/SiOy/SiNx Metal-Insulator-Metal (MIM) capacitor. According to the foundry process description, SiNx and SiOy are respectively 0,15µm and 0,8µm thick.

An investigation has been made to check these thicknesses by using two main methods: XPS with ion etching and FIB-SEM analysis. XPS helps to make a quantitative analysis and extract the atomic ratio of each chemical element at the surface (hydrogen excepted). In the present case, it is combined with ion etching in order to measure the atomic ratio depth profile. The spot size of XPS beam is 200µm. Figure 1 shows the etched spot position on the die; a device-free area has been selected so that the layer sequence is GaAs/SiNx/SiOy/SiNx.

Estimation based on XPS erosion rate

The thicknesses of the superficial layers were obtained by splitting the XPS depth profile (Figure 2) in 3 areas: external SiNx (0-250s), SiOy (250-1500s) and internal SiNx (1500-1800s).

After 1800s of etch time, results show the rise of Ga and As peaks which means that GaAs substrate has been reached. This reveals that the experiment is successful.

With an etch rate around 0,6nm/s (pre-calibrated), the thicknesses of the different layers have been obtained easily with respectively 150nm for the external SiNx (passivation layer), 750nm for the SiOy dielectric and 180nm for the internal SiNx in contact with GaAs. The atomic ratio of silicon (Si) and nitrogen (N) in two SiNx layers is practically the same but it cannot be concluded definitely that the two layers are of the same quality; in particular the hydrogen content could be different, which is very important for stress in the layer and the protection it offers. Ratio N/Si ≈ 0.06 between 250-1500s (SiOy zone), which means either that N is inserted in the SiOy layer or that the etching front is quite rough. Measurements have been repeated at different positions and on other samples, giving the same results.

FIB-SEM measurement

In the non-planar gate area, FIB-SEM is a more adequate technique to measure both metallic and insulating materials. Figures 3a and3b show the cutting plane chosen to measure gate finger metal thicknesses, and a tilted SEM view of the cut. A Platinum (Pt) layer has been deposited to protect the surface during FIB ion milling. The die tilt is 30° in order to observe the bottom of the section. The experimental results are shown in Figure 3c. For the insulating layers, thickness is ~129nm for SiNx (passivation layer), ~698nm for SiOy dielectric and ~216nm for the internal SiNx layer. These values are in reasonable agreement with those deduced from XPS, resp. 150nm, 750nm and 180nm. Note that the areas measured by these two techniques are different and do not show necessarily the same thickness exactly.

The small differences between the predicted thicknesses (foundry datasheet: 800nm for SiOy, 150nm for SiNx) and those measured by XPS and FIB-SEM is explained by the fact that the deposition is chemical (PECVD: Plasma-Enhanced Chemical Vapor Deposition), so that geometrical singularities have a notable impact on local thicknesses. In the rest of this study, the datasheet thicknesses are considered confirmed and are used in the simulations.

Summary of analytic models[4][5] [6]

Residual stresses were first analyzed by Stoney [START_REF] Stoney | The tension of metallic films deposited by electrolysis[END_REF] who considered a bilayer strip consisting of a film and substrate. He derived a simple equation to relate the stress in the film to the curvature of the strip. In the simple case where the film thickness is considered infinitesimal compared to the substrate thickness (𝑡 𝑠 ≫ 𝑡 𝑓 ), Nix [START_REF] Nix | Mechanical Properties of Thin Films[END_REF] has used continuum mechanics laws and Stoney's assumptions to evaluate the thermo-mechanical stress induced in the film far from the edge by die heating from room temperature 𝑇 0 to temperature 𝑇:

𝜎 𝑓 (𝑇) = 𝐸 𝑓 1 -𝜐 𝑓 (𝛼 𝑠 -𝛼 𝑓 )(𝑇 -𝑇 0 ) (1)
Where 𝛼 𝑠 and 𝛼 𝑓 are respectively the CTE of the substrate and the film. 𝐸 𝑓 , 𝜐 𝑓 are the Young's modulus and the Poisson's ratio of the film.

Recently, a more complete, yet simple analytical model for analyzing thermal stresses and deformation in multilayers has been developed by Hsueh [START_REF] Hsueh | Thermal stresses in elastic multilayer systems[END_REF]. The multilayer is subjected to a temperature change ∆𝑇; the CTE of the substrate and the films are 𝛼 𝑠 , and 𝛼 𝑖 , 𝑖 = 1 … 𝑛 respectively. At positions away from the edge of the multilayer, the stresses induced by the CTE mismatch are in-plane (i.e., parallel to the interface) and equi-biaxial for the plane geometry. Both the stress normal to the interface and the interfacial shear stress are zero. The coordinate system is defined such that the layer 1 (𝑖 = 1) and substrate interface is located at 𝑦 = 0, the interface between layers 𝑖 and 𝑖 + 1 is located at 𝑦 = ℎ 𝑖 . Here 𝑡 𝑠 and 𝑡 𝑖 are respectively the thickness of the substrate and the films. The free surfaces of the multilayer are located at 𝑦 = -𝑡 𝑠 and 𝑦 = ℎ 𝑛 . The position in y axis is obtained with ℎ 𝑖 = ∑ 𝑡 𝑗 𝑖 𝑗=1 . The in-plane biaxial stress distributions in the substrate and layers, 𝜎 𝑠 and 𝜎 𝑖 can be expressed as:

𝜎 𝑠 = 𝐸 𝑠 1 -𝜐 𝑠 (𝑐 + 𝑦 -𝑏 𝑟 -𝛼 𝑠 ∆𝑇) (for -𝑡 𝑠 ≤ 𝑦 ≤ 0) (2𝑎) 𝜎 𝑖 = 𝐸 𝑖 1 -𝜐 𝑖 (𝑐 + 𝑦 -𝑏 𝑟 -𝛼 𝑖 ∆𝑇) (for 𝑖 = 1 to 𝑛) (2𝑏)
The three parameters 𝑐, 𝑏, and 𝑟 are respectively the uniform strain component, the position of the bending axis at which the bending strain component is zero and the radius of curvature of the system. They can be obtained by applying the three boundary conditions: (i) the resultant force due to the uniform strain component is zero, (ii) the resultant force due to the bending strain component is zero, (iii) the resultant bending moment due to the stresses described by equation ( 2a) and (2b) is zero. (3𝑐)

𝑐 = (𝐸 𝑠 ′ 𝑡 𝑠 𝛼 𝑠 + ∑ 𝐸 𝑖 ′ 𝑡 𝑖 𝛼 𝑖 𝑛 𝑖=1 )∆𝑇 𝐸 𝑠 ′ 𝑡 𝑠 + ∑ 𝐸 𝑖 ′ 𝑡 𝑖 𝑛 𝑖=1 (3𝑎) 𝑏 = -𝐸 𝑠 ′ 𝑡 𝑠 2 + ∑ 𝐸 𝑖 ′ 𝑡 𝑖 (2ℎ 𝑖-1 + 𝑡 𝑖 ) 𝑛 𝑖=1
Knowing the value of 𝑐, 𝑏 and 𝑟, the general solutions for the stress distributions 𝜎 𝑠 and 𝜎 𝑖 in the multilayer given by equations (2a) et (2b) are complete.

Numerical models

This section focuses on a description of the numerical models. Finite element simulations were performed using Abaqus Standard Software. Two geometrical models have been presented:

 "simple case" considering the whole die as a planar multilayer stack on substrate (Figure 5,6).  the second "complete model" takes into account the gate area and geometrical singularity (Figure 7). A structured mesh composed of linear elements has been built. Because of the in-plane stress simplification, simulation was carried out in 2D and due to geometrical symmetries, only half of the multilayer was modeled. The contacts are assumed to be tie elements. Regarding boundary conditions and symmetries, the node (𝑥 = 0 and 𝑦 = -𝑡 𝑠 ) was locked and the normal displacement of Line A set at zero as shown in Figure 4.

In case gate current is present, a local heat source due to the current flow must be taken into account. However here, since THB tests use the reverse biasing on the transistor, the current flow in the gate is locked and there is no internal heat source. Therefore, the flux is considered homogeneous all around the die (Figure 4). The heating condition was a uniform thermal flux on the external surface, applied from room temperature (~22°C) to test temperature (85°C). In the range 22°C-85°C variations of materials properties are considered negligible and simulations were performed under linear elastic theory. For the input data resumed in Table 1, the values came from bibliography [7] [8].

FEA materials description is more extensive than in previous analytical models: thermal conductivity, heat capacity and density have been added. The contact at the interface is supposed to be perfect. The following thicknesses are used in modeling: 150nm SiNx, 800nm SiOy, 20nm Ti, 140nm Al, 20nm Au and 100µm GaAs substrate.

Simple case

In the simple case, two models have been used for simulations. First, the die is considered as a multilayer with a stack of GaAs/SiNx/SiOy/SiNx because it represents a typical cross-section of the free surface of the GaAs microwave (see Figure 5). Second, calculations are made with a stack of GaAs/Ti/Al/Au/SiNx/SiOy/SiNx representing the surface topography in gate finger areas without geometrical singularities.

Complex case

In reality, geometrical singularities (transistor, capacitor, resistor, interconnections, air bridge…) in the active areas of the chip represent ~30% of the surface topography.

For the complex model, one finger of the 0,18µm pHEMT gate has been defined as a multilayer of Ti/Al/Au stack. Above this, there is a SiNx/SiOy/SiNx stack. Gate length and width details were drawn from the previous FIB-MEB cross-sections. The induced stresses (with ∆𝑇 = 63°𝐶) for Nix and Hsueh models through the thickness of the system (Figure 5,6) are calculated from equations ( 1), (2a) and (2b). In order to examine the accuracy of these results, Abaqus FEA simulations have been performed (Figure 5, 6 and 7) with the same dimensional and materials properties. In the calculation, Abaqus FEA has been used for the two simple configurations represented on Figures 5,6 and compared to Nix and Hsueh methods which are usually used in the literature. Then for the "complex" case shown in Figure 7, only Abaqus FEA has been done in order to complete our analysis. The SiNx/SiOy/SiNx stack directly deposited on the GaAs substrate as shown in Figure 5 represents the configuration about ~70% of the die surface (Figure 1). Each layer has its specificity: the bottom SiNx acts as insulator, the SiOy is the dielectric for MIM capacitors and the upper SiNx layer is the die passivation layer. By considering this cross-section, the thermo-mechanical stresses through the thickness of this stack for Hsueh, Nix and Abaqus FEA models are shown in Figure 8b. The three layers are in tensile stresses (𝜎 𝑥𝑥 component) because for each layer the CTE is lower than for the GaAs. For the same layer, the stress induced at the lower surface is always less tensile than the stress at the upper surface. The results obtained indicate tensile stresses in the internal SiNx, the SiOy dielectric and the SiNx passivation layer respectively 50MPa, 42MPa and 50MPa. In the substrate, induced stress distribution is shown in Figure 8a. According to Hsueh and Abaqus FEA results, the top surface (𝑦 = 0) is subjected to tension and the bottom (𝑦 = -𝑡 𝑠 ) is subjected to compression. As it has been considered in the assumptions for Nix's model (𝑡 𝑠 ≫ 𝑡 𝑖 and 𝜎 𝑠 ≈ 0), the stress induced in the substrate is supposed negligible.

A second stack of GaAs/Ti/Al/Au/SiNx/SiOy/SiNx shown in Figure 6 is used to evaluate the induced stress in the gate metallization of the transistors. In this stack, the inserted Ti/Al/Au layers representing the gate fingers are subjected to compression because they have a higher CTE than GaAs. Mainly among these three layers, Aluminum (Al) is subjected to high compression because its CTE is higher than Ti and Au.

Comparisons of thermomechanical stresses in the gate fingers evaluated by different methods are shown in Figure 8e. It should be noted that the thermo-mechanical stresses in the gate metal (Al), Ti and Au are compressive and respectively around -108MPa, -30MPa and -68MPa. According to Figure 8c, the induced stress in the substrate for this stack is in tension at the top surface and in compression at the bottom.

Figure 9: Thermo-mechanical stresses 𝜎 𝑥𝑥 simulated with Abaqus in the gate finger area (Detail A) for complex case.

For the complex case (Figure 7), simulation results with Abaqus FEA model shown in Figure 9 indicate a tensile stress concentration up to 181MPa in the SiNx passivation layer, and a stress level of 110MPa for the internal SiNx layer. The induced stress in the substrate for this configuration is the same as shown in Figure 8c.

Discussion

Results from the previous section show the thermomechanical stress induced in the GaAs microwave devices by the temperature raise during THB 85°C/85RH testing. Each one of the three methods used has its particularities. Nix's formulation is helpful as a first approximation because of its simplicity (equation 1); Hsueh's equations are a more precise by including the thicknesses and all the details on the films deposited, finally the Abaqus FEA is the only one that can take into account geometrical singularities and include more materials properties.

For the flat layers, the differences are quite small, a few MPa at most except in the substrate where the averaged Nix's stress does not show the gradient. Let us compare for instance Hsueh and Abaqus FEA results for the stack GaAs/Ti/Al/Au/SiNx/SiOy/SiNx (Figure 8c,d) with the outcomes from stack GaAs/SiNx/SiOy/SiNx (Figure 8a,b). With the additional metallic layers (Ti/Al/Au), stress in the GaAs substrate decreases by about 40% compared with (SiNx/SiOy/SiNx). This is explained by the negative sign of ∆𝛼 = (𝛼 𝑠 -𝛼 𝑓 ) for each inserted thin film. However stress in the thick substrate is very small in all cases (<2MPa). A variability of thickness of the thin films (SiNx and SiOy) has been observed in the FIB-SEM cross-sections (Figure 3c). An investigation of its effect on the induced stresses has been made using the model presented in Figure 6. According to Hsueh's model results shown in Figure 10, the stress in the SiNx external passivation layer hardly changes (<3%) when the thicknesses of the other two are decreased even strongly. This low variation is explained by the condition: 𝑡 𝑠 ≫ 𝑡 𝑖 so that the substrate behavior dominates. A comparison of the thermo-mechanical stresses simulated with Abaqus FEA in a simple (Figure 8b) and complex (Figure 9) cases shows similar results in the metallic gate (-108MPa in simple case and -115MPa in the complex case). However the tensile stress in the SiNx passivation layer is strongly affected by geometrical singularities which induce stress concentration areas: from +50 MPa in the simple case, it jump to more than 150MPa near the step, which might have consequences on the integrity of the layer.

Failure analysis of GaAs device under THB 85/85 aging testing has revealed that damages were caused mainly by the degradation of passivation layer (cracks, diffusion…) [START_REF] Hisaka | Degradation mechanisms of GaAs PHEMTs in high humidity conditions[END_REF] [10] which leads to moisture intrusion in the die followed by corrosion. In a previous work on the same technology of GaAs microwave [START_REF] Ben Naceur | Evaluation of quasihermetic packaging solutions for active microwave devices and space applications[END_REF], several cracks have been observed in the SiNx passivation layer (Figure 11) and also in the gate metal (Figure 12). This crack induces a large degradation of drain current monitored during THB test.

Indeed, the temperature change and aggressive environment induce stresses which must be added to the existing residual stresses from fabrication. In our specific case (Figure 7), by evaluating the thermo-mechanical stresses induced by heating from room temperature to 85°C, it should be noticed that the stress concentration areas (Figure 9) are very similar to cracked areas (Figure 11). This result suggests leads us to claim that the failure of the passivation layer is directed linked to architectural failure and thermo-mechanical stresses must be considered in the failure analysis. The author [START_REF] Ben Naceur | Evaluation of quasihermetic packaging solutions for active microwave devices and space applications[END_REF] also explains these failures (Figure 12) by the existing internal stresses added to the induced stresses (thermo-mechanical and hydro-mechanical stresses). Due to the linearity of stress state, any initial stress state, once determined, can simply be added to the results presented.

Conclusions

Coefficients of thermal expansion mismatch between layers induces internal stress in multilayers thin films on substrate. These residual stresses in multilayer systems influence their reliability and can result in cracking.

In the present analysis, the thermo-mechanical stress induced in thin films of GaAs microwave under THB testing is evaluated by heating from room temperature to test temperature. Analytical models [5] [6] have been used to make analysis and the FEA has been used for results validation.

With a model that takes into account the geometrical singularities and the thickness of layers (Figure 7) in the

Figure 1 :

 1 Figure 1: Etched spot position on GaAs MMIC for XPS and ion etching analysis.

Figure 2 :

 2 Figure 2: XPS depth profile for superficial layers on GaAs substrate

Figure 3 :

 3 Figure 3: FIB-SEM analysis on GaAs MMIC.

Figure 4 :

 4 Figure 4: 2D modeling of half die with Abaqus CAE.

Figure 5 :Figure 6 :

 56 Figure 5: Detail A in the simple case.Stack GaAs/SiNx/SiOy/SiNx

Figure 7 :Figure 8 :

 78 Figure 7: Detail A in the complex case. Typical cross-section of GaAs pHEMT gate in contact areas.

  Conference on Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Microsystems (EuroSimE 2015)

Figure 10 :

 10 Figure 10: Thermo-mechanical stresses 𝜎 𝑥𝑥 evaluated with Hsueh's method in the SiNx passivation for different thicknesses of SiNx and SiOy

  thickness (nm) SiN x films thickness (nm) SiN x passivation layer stresses σ xx (MPa) 𝑆𝑡𝑟𝑒𝑠𝑠 𝑐𝑜𝑛𝑐𝑒𝑛𝑡𝑟𝑎𝑡𝑖𝑜𝑛 𝑖𝑛 𝑆𝑖𝑁 𝑥

Figure 11 :

 11 Figure 11: Damages of GaAs die after 3000 hours at 85°C/85%RH without bias [2]: SEM observation of gate finger, drain and source contacts

Table 1 :

 1 Materials properties[START_REF] Ghandhi | Fabrication Principles, Silicon and Gallium Arsenide[END_REF] [START_REF] Massénat | Circuits en couches minces -Couches minces traditionnelles[END_REF]

	Properties	AsGa SiNx	Materials SiOy	Ti	Al	Au
	Young's
	modulus	85.5	150	79.2	114	69	78
	(GPa)
	Poisson's ratio	0.31	0.3	0.22	0.34 0.34 0.42
	CTE (10 -6 /K)	6.86	3	0.55	8.4	23	14.2

16th International Conference on Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Microsystems (EuroSimE 2015)

Acknowledgments

The authors gratefully acknowledge the financial support of CNES (French government space agency) and Thales Alenia Space and the permission to publish.

pHEMT gate fingers (Ti/Al/Au), FE simulations indicate compression of metal gate (Al) reaching -115MPa. This is due to its higher CTE with respect to the surrounding SiNx. Inversely, the SiNx passivation experiences up to 181MPa when its geometric singularities are taken into account. These results have been correlated with the cracks areas observed previously [START_REF] Ben Naceur | Evaluation of quasihermetic packaging solutions for active microwave devices and space applications[END_REF] on non-hermetic tested devices under reverse bias. This proves the usefulness of thermo-mechanical analysis to complete failure analysis under THB aging life test. It should be considered in further analysis of the reliability of active microwave devices.

In view of accurate evaluation, further investigation should determine the main properties of die thin films (SiNx, SiOy…) by using ultrasonic measurements method on dice.


Trong-Son Cao
email: trongson.cao@arcelormittal.com

Pierre-Olivier Bouchard
email: pierre-olivier.bouchard@mines-paristech.fr

Pierre Montmitonnet
email: pierre.montmitonnet@mines-paristech.fr

Trong-Son Cao

Pierre-Olivier Bouchard

T.-S Cao

P.-O Bouchard

Ductile damage prediction in different cold forming processes

Keywords: Ductile damage, metal forming, triaxiality, Lode angle

The purpose of the present paper is to show how and to what extent the introduction of refined, shear sensitive models improves on previous ones, based on triaxiality only, for the phenomenological description of ductile damage in bulk cold metal forming processes. Wire-drawing and wire rolling are taken as examples. A set of mechanical tests has been conducted: round bar, notched bar and plane strain tensions as well as torsion for pure shear deformation. Both constitutive and damage models parameters have been carefully identified, with back-computation of the laboratory tests for validation. Application of the models to the cold forming processes, described here, shows the superiority of the shear-enhanced models for locating maximum damage in flat wire rolling, where a significant amount of shear is present ("blacksmith's cross" deformation pattern). On the contrary, it proves unnecessary for low-shear processes such as wire-drawing. The cavity-growth Gurson-Tvergaard-Needleman model, with a modified formulation for nucleation, seems to be the best basis for damage prediction in patented high carbon steel.

Introduction

In the modern world, products should be produced under low cost and severe environmental constraints, there is thus no longer place for the classical "trial and error" tests. For all industrial cold forming processes, the ability of numerical modeling to predict ductile fracture is crucial. The common issue is to predict ductile damage and fracture in multi-axial and non-proportional loadings involving very large plastic strain, in which the damage may localize away from the maximum critical strain locations. Moreover, the shear effects in several processes are important, damage models thus must be able to capture the associated shear damage mechanism. There is thus real need to develop robust damage models for industrial applications. By robustness, the challenges for a good model are three-fold: (1) physical mechanisms of the onset of damage and fracture must be captured; [START_REF] Cao | Characterization of ductile damage for a high carbon steel using 3d x-ray micro-tomography and mechanical tests -application to the identification of a shear modified gtn model[END_REF] the models have to be suitable for numerical implementation and must be simple enough for real-scale structure simulations; and (3) models parameters should be convenient for identification for massive applications. The purpose of the present paper is to examine the capacity of different types of models in three approaches of ductile damage (namely, fracture criteria, coupled phenomenological and micro-mechanics based models) to predict damage and fracture in two different processes: wire-drawing and wire rolling.

Damage models developments and implementations

Six ductile damage models are compared, which are described in terms of invariants of the stress tensor: the stress triaxiality η = σ m /σ; and the Lode parameter θ = 1

-6θ π = 1 -2 π arccos ( 27 2 det(s)) σ 3
(σ m , σ are the mean and von Mises equivalent stresses; s is the deviatoric stress tensor). Elasto-plastic hardening law was used combined with J 2 plasticity, except the case of Gurson-like models (section 2.3) where pressure-dependent plasticity model was employed. The hardening law used was a modified form of the Voce law, in which the flow stress is defined as:

σ 0 = σ p0 + (σ ps -σ p0 + K 2 p ) (1 -exp(-n p )) (1)
where: σ 0 is the flow stress of material matrix, n and K 2 are material parameters;

σ p0 = σ 0 ( p = 0); σ ps = σ 0 ( p → ∞) if K 2 = 0.
The unit of K 2 , σ p0 and σ ps is M P a. The identified hardening law was validated for different mechanical tests: compression, torsion and tensile tests ( [START_REF] Cao | On the development and identification of phenomenological damage models -Application to industrial wire drawing and rolling processes[END_REF][START_REF] Cao | Characterization of ductile damage for a high carbon steel using 3d x-ray micro-tomography and mechanical tests -application to the identification of a shear modified gtn model[END_REF]).

Fracture criterion

The Bai & Wierzbicki fracture criterion was used ( [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF]), in which the strain to fracture function is defined as:

f (η, θ) = 1 2 D 1 e -D2η + D 5 e -D6η -D 3 e -D4η θ 2 + 1 2 D 1 e -D2η -D 5 e -D6η θ + D 3 e -D4η (2)
where D 1 , D 2 , D 3 , D 4 , D 5 , D 6 are 6 positive material parameters which need to be identified. A linear incremental relationship is assumed between the damage variable D and the equivalent plastic strain p :

D( p ) = p 0 d p f (η, θ) (3)

Coupled phenomenological models

Three models were implemented and used: models in Continuum Damage Mechanics (CDM) framework (Lemaitre [4], Lode-dependent Enhanced Lemaitre-LEL [START_REF] Cao | A Lode-dependent enhanced Lemaitre model for ductile fracture prediction at low stress triaxiality[END_REF]) and a fracture strain-based coupled phenomenological model proposed by Xue [START_REF] Xue | Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading[END_REF]. The governing equations of the LEL model are summarized as:

Ḋ = ˙ p Y S b 1 α1+α2θ 2 , if p > D = D0 exp(-Aη) and η > η cutoff 0 , otherwise (4)
σ = w(D)σ M E = w(D)E M with w(D) =    1 -D , if η ≥ η 1 1 -(1-h)η+hη1-η2 η1-η2 D , if η 1 > η ≥ η 2 1 -hD , if η < η 2 (5)
where α 1 , α 2 , D0 , A, η 1 , η 2 are additional parameters; η cutoff is the cutoff value of the stress triaxiality to be identified, which could be chosen equal to -1/3 from the study of [START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF] for a first approximation (or taken as a function of the Lode parameter as proposed by [START_REF] Lou | Modeling of shear ductile fracture considering a changeable cut-off value for stress triaxiality[END_REF] ). It should be noted that, if α 2 = 0, α 1 = 1 and A = 0, the evolution equation of the original Lemaitre model is retrieved.

Xue ([6]) proposed a phenomenological damage model, which is based on the definition of the equivalent fracture strain f as a function of hydrostatic pressure (p) and Lode angle (θ L ):

f (p, θ L ) = f 0 µ p (p)µ θ (θ L ) , with: µ p = 1 -q ln(1 - p p L ) and µ θ = γ + (1 -γ) 6 |θ L | π k (6)
The scalar damage variable D X was used as an internal variable to represent the material degradation. The weakening function w(D X ) = 1 -D β X is adopted to describe the damage effect on the macroscopic strength:

ḊX = m p f (p, θ L ) m-1 ˙ p f (p, θ L ) and σ = 1 -D β X σ M ( 7
)
where β is the weakening exponent, p is the equivalent plastic strain, σ M is the flow stress of undamaged material.

Micromechanics-based models

The yield function of the GTN model ( [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]) writes:

φ = σ σ 0 2 + 2q 1 f * cosh(- 3q 2 2 p σ 0 ) -1 -q 3 f * 2 = 0 ( 8
)
σ 0 is the flow stress of matrix material; q 1 , q 2 , q 3 = (q 1 ) 2 are material constants; f * is the effective void volume fraction, which accounts for the voids' linkage:

f * = f , if f < f c f c + f * u -fc f f -fc (f -f c ) , if f ≥ f c (9)
where f c represents the critical value of f at which void coalescence begins, f f its value at ductile failure, and

f * u = q1± √ q 2 1 -q3 q3
the corresponding value of f * at failure. The evolution of void volume fraction is described as:

ḟ = ḟnucleation + ḟgrowth ( 10
)
where ḟgrowth is defined as ḟgrowth = (1 -f )trace( ˙ p ) and ḟnucleation is often described by a Gaussian curve introduced by [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF]:

ḟnucleation = f N S N √ 2π exp - 1 2 p -N S N 2 ˙ p = A ( p ) ˙ p (11)
Xue proposed a new damage variable (D) which accounts for void shearing damage:

Ḋ = δ D q 1 ḟ + Ḋshear ( 12
)
where δ D the damage rate coefficient. For the GTN model, which accounts for coalescence (i.e. using f c to define the start of coalescence stage), this coefficient is defined as:

δ D = 1 , if D ≤ D c = q 1 f c f * u -fc f f -fc , if D c < D ≤ 1 (13)
Ḋshear = q * 3 f q4 p ˙ p g θ where g θ = 1 - 6 |θ L | π (14)
In addition to this modification, in order to account for the influence of the stress triaxiality on the nucleation process, Cao et al. [START_REF] Cao | Characterization of ductile damage for a high carbon steel using 3d x-ray micro-tomography and mechanical tests -application to the identification of a shear modified gtn model[END_REF] proposed a formulation for the strain at maximum nucleation N as:

N = N0 exp(-Bη) ( 15
)
where η is the stress triaxiality, which can be taken as the initial stress triaxiality if the loading is nearly proportional; B and N0 are two parameters to be identified. This modification allows the void nucleation process to take place earlier at high stress triaxialities, which is materialized by the experiments of Kao et al. [START_REF] Kao | Tensile fracture and fractographic analysis of 1045 spheroidized steel under hydrostatic pressure[END_REF]. The GTN model with the modification proposed by [START_REF] Xue | Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading[END_REF] and [START_REF] Cao | Characterization of ductile damage for a high carbon steel using 3d x-ray micro-tomography and mechanical tests -application to the identification of a shear modified gtn model[END_REF] is called the modified GTN model in the present study.

Implementation and calibration

All models were implemented in Forge®, an implicit, 3D Finite Element software based on a mixed velocity-pressure formulation using mini-elements ( [START_REF] Cao | A detailed description of the Gurson-Tvergaard-Needleman model within a mixed velocity-pressure finite element formulation[END_REF][START_REF] Cao | Numerical simulation of 3D ductile cracks formation using recent improved Lode-dependent plasticity and damage models combined with remeshing[END_REF]). Identification of the parameters of the models above, plus those describing workhardening, has been performed using a set of experiments covering the plane (η, θ): cylinder compression (-1/3;-1); torsion (0;0); plane-strain, flat grooved specimen (≈ 2/3-1;0); round and notched bar tension (1/3-2/3;1). For Gurson-Tvergaard-Needleman -like models, X-Ray microtomography has been used also to identify nucleation parameters and porosity at fracture ( [START_REF] Cao | Characterization of ductile damage for a high carbon steel using 3d x-ray micro-tomography and mechanical tests -application to the identification of a shear modified gtn model[END_REF]). The calibration of all models were based on inverse analyses using multi-objective automatic optimization program ( [START_REF] Cao | Identification methodology and comparison of phenomenological ductile damage models via hybrid numerical-experimental analysis of fracture experiments conducted on a zirconium alloy[END_REF]).

Examples of applications

In this section, the above-mentioned models are compared in terms of damage prediction in two industrial forming processes: wire drawing and wire flat rolling (high carbon steel). In addition, the Lemaitre and B&W models were also applied to predict fracture in rolling of stainless steel wire. Fig. 1a shows the schematic of the studied wire drawing process. Fig. 1b shows the damage state observed after four passes of drawing, near the wire center (total average deformation ≈ 0.81). As can be observed, deformable inclusions are oriented following the drawing direction. At this stage, inclusion cracking is also observed but negligible.

Wire drawing process

For the wire drawing process, most models predict higher damage in the wire core than in the wire skin ( [START_REF] Cao | A comparative study of three ductile damage approaches for fracture prediction in cold forming processes[END_REF]). In experiment, fracture occured in the wire core after 14 drawing passes. Fig. 2 shows damage distributions on the longitudinal and transverse cross sections, at the second pass, which confirm high localization of damage in the wire core (except the modified GTN model). Moreover, between the Lemaitre and Xue coupled models, a stronger localization in the wire center can be observed with the Lemaitre model, while for the Xue model, damage changes gradually from the wire core to the wire surface. In addition, regarding the LEL model, which accounts for the influence of the Lode parameter, damage accumulation is also predicted on wire surface, although less intense than in the wire core. It means that, for the models based on both the Lode parameter and the stress triaxiality (i.e. B&W, Xue and LEL), the influence of shear, which is maximum on the wire surface, is accounted for with a stronger weight than for the stress triaxality-based model (e.g. Lemaitre). Moreover, in the working zone, since the stress triaxiality in the wire surface is negative, damage accumulation predicted by the stress triaxiality-based Lemaitre model is not noticeable at this location. It should be noted that, in [START_REF] Cao | A comparative study of three ductile damage approaches for fracture prediction in cold forming processes[END_REF], results revealed that only the GTN model gives a satisfactory results regarding the prediction of the moment of fracture (after 13 passes).

Wire flat rolling

Fig. 3a represents the schematic of the studied wire flat rolling process. As the wire passes through the rolls, all material points across the width experience some tendency to expand laterally (transverse direction); this is called "spread". The edges are thus strained in tension, which may lead to edge cracking. The edge center tends to expand laterally more than the upper and lower surfaces, which produces barreled edges similar to those observed in compression of a cylinder. Therefore, theoretically, in this process, the barreled areas are critical zones where cracks tend to initiate first. An experimental study of damage in rolling process was performed in [START_REF] Massé | Mechanical and damage analysis along a flat-rolled wire cold forming schedule[END_REF]. Damage was identified at a microscopic scale through SEM observations, which showed decohesion around non deformable inclusions (matrix-inclusion debonding) and fragmentation of deformable inclusions. The FEM analysis of the strain map was also presented, which showed a localization in form of a cross (the "blacksmith cross"). This strain localization may explain the evolution of microstructure during rolling (see [START_REF] Massé | Mechanical and damage analysis along a flat-rolled wire cold forming schedule[END_REF] for more details). The experimental sketch of damage, superposed on a strain map obtained with FEM analysis, is presented in Fig. 3b. Since this process involves a global vertical compression with a significant transverse flow, the expansion of decohesion thus follows the transverse direction. Moreover, the inclusions located in the branches of the blacksmith cross seem to have a preferential orientation due to rotation associated with shear strain in this area. The authors mentioned a higher void density as well as a higher flattening of cavities in the wire core than in the edge. Therefore, from this experimental observation, the rolling process principally affects the voids in the wire core and on the blacksmith cross. Voids located in the barreled zone were also expanded laterally due to rolling, but the expansion is smaller.

In rolling (Figs. 4 and5), Lemaitre's model gives maximum damage on the edge, with a small secondary peak at the in the working zone, since the stress triaxiality in the wire surface is negative (see Fig. 6a), damage accumulation predicted by the stress triaxiality-based Lemaitre model is not noticeable at this location. Regarding micromechanically-based models, the GTN model (Fig. 7e) predicts damage higher in the wire core than for other models (although at this stage, the void volume 14 center corresponding to the heritage of wire-drawing. On the contrary, its shear-enhanced version ("LEL") depicts the blacksmith's cross. The same difference can be found between Gurson-Tvergaard-Needleman ("Porosity") and shear-enhanced Gurson-Tvergaard-Needleman ("Dx-GTN") (in experiments, the most dangerous zones are the center and the blacksmith cross -see ( [START_REF] Massé | Study and optimization of high carbon steel flat wires[END_REF], or [START_REF] Cao | Modeling ductile damage for complex loading paths[END_REF]). The Bai & Wierzbicki model ("B&W") is close to Gurson-Tvergaard-Needleman, whereas Xue's model resembles the enhanced Lemaitre model

Fracture prediction in wire round rolling

The studied rolling process consists of 5 "passes", each pass involving 2 stands (3+3 rolls or 2+2 rolls). The first pass uses the 3+3 technology (see Fig. 6a), where 2 stands of 3 rolls are set next to each other. For each stand, one roll makes an angle Fig. 7 shows the comparison between cracks observed at the outlet in wire round rolling process on a stainless steel obtained from experiment (Fig. 7a) and from numerical prediction with the Lemaitre-type model (Fig. 7b). For this process, fracture occurs on the wire surface. It should be noted that, in [START_REF] Cao | Modeling ductile damage for complex loading paths[END_REF][START_REF] Cao | Comparison of reduction ability between multi-stage cold drawing and rolling of stainless steel wire -Experimental and numerical investigations of damage[END_REF], the authors showed that, both the Lemaitre-type and B&W models give correct prediction for both fracture location and the instant of fracture for this process.

Closure remark

Six damage models from the recent literature have been implemented, their parameters identified from experiments, and their performance tested on different cold forming schedules. Note that the sets of parameters identified by laboratory mechanical tests have been used "as is", without any modification, to study wire-drawing and wire rolling. The experimental trends in terms of time and location of failure have been checked. In the wire drawing of patented high carbon steel, the Gurson-Tvergaard-Needleman model is the one which best fits measurements. Other models predict a largely premature failure (although critical positions were accurately predicted). Regarding rolling process, experimental void density maps can be reproduced only with shear-emphasizing damage models. It is quite difficult to obtain a good balance between the lateral surfaces (low but non negligible porosity), the center (largest density) and the arms of the blacksmith's cross (intermediate). Again, Gurson-Tvergaard-Needleman is globally the best if used with shear-enhancement (here, the Xue phenomenological formulation). Another application for round-rolling on a stainless steel reveals that both the stress-triaxiality based (Lemaitre) and Lode and triaxiality based models (B&W) provided accurately results in terms of fracture prediction in the zone of high lateral expansion (positive stress-triaxiality zone). These results suggest that, for the process with high triaxiality zone, triaxiality-based models could give correct predictive results. However, regarding processes involving strong shear, models should account for both stress triaxiality and Lode parameter. In addition, as can be observed, voids are elongated in the drawing direction in the wire drawing process and voids are flattened and rotated following the shear band in the rolling process. The GTN model with underlying assumption of spherical cavities shows its limitation in these cases. The use of a purely micromechanical model ( [START_REF] Danas | Numerical modeling of elasto-plastic porous materials with void shape effects at finite deformations[END_REF][START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF]), which accounts for void shape change and void rotation, is of great interest for future studies.

Figure 1 :

 1 Figure 1: (a) Schematic of drawing; and (b) SEM images of microstructure in the drawing direction (horizontal direction) after four passes wire drawing [15]).

Figure 7 :

 7 Figure 7: Damage at steady state of second drawing pass for six damage models: (a) Bai & Wierzbicki, (b) Xue, (c) Lemaitre, (d) LEL, (e) GTN and (f ) modified GTN by Xue. The values are represented on the longitudinal and transverse cross sections.

Figure 2 :Figure 3 :

 23 Figure 2: Damage at steady state of second drawing pass for six damage models: (a) Bai & Wierzbicki, (b) Xue, (c) Lemaitre, (d) LEL, (e) GTN and (f) modified GTN by Xue. The values are represented on the longitudinal and transverse cross sections.

Figure 4 :Figure 5 :

 45 Figure 4: Damage on a cross section perpendicular to the rolling direction with different models: (a) Xue, (b) Lemaitre, (c) B&W, and (d) GTN.

Figure 6 :

 6 Figure 6: The studied rolling mills and mesh used.

Figure 7 :

 7 Figure 7: (a) Experimental observations of cracks on the surface of rolled wire; (b) numerical prediction obtained with the Lemaitre-type model, showing cracks at outlet ([18].

Acknowledgments

The financial support from ArcelorMittal, Cezus-Areva and Ugitech via the METAL project is appreciated.