Dr Daniel Weisz-Patrault

Thomas Koedinger

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

PhD Siamak Serajzadeh
email: serajzadeh@sharif.edu

PhD Bernard Fedelich
email: bernard.fedelich@bam.de

PhD Dmitrii Legatiuk
email: dmitrii.legatiuk@uni-weimar.de

Takayuki Otsuka
email: otsuka.6gx.takayuki@jp.nssmc.com

Residual stress on the run out table accounting for multiphase transitions and transformation induced plasticity

Keywords: Run out table, Residual stress, Eigenstrain, Multiphysics Run out table, Residual stress, Eigenstrain, Multiphysics

The development of harder and thinner new steel grades necessitates reasonably fast numerical simulations of forming processes in order to optimize industrial conditions through parametric studies. This paper focuses on the evolution of residual stresses of thin strips during cooling on the run out table. Since the problem involves multiphysics and non-linear processes, comprehensive and fully coupled numerical approaches may be difficult to use to design or optimize industrial conditions because of extensive computation times. Therefore, a simplified numerical simulation has been developed and consists in solving first the thermal problem coupled with multiphase transitions and then the mechanical problem accounting for thermal expansion, metallurgical deformation and transformation induced plasticity. Residual stress profiles through the strip thickness are also computed in order to evaluate classic flatness defects such as crossbow and longbow. A postprocessing is also included in order to quantify out of plane displacements that would take place if the strip were cut off the production line. It consists in computing at finite strain the relaxation of residual stresses when the tension applied by the coiler is released. The proposed numerical strategy has been tested on common industrial conditions.

Introduction

Both the constant development of new steel grades and higher quality requirements of flat rolled products (e.g., strip thickness, strip flatness) necessitate a deep understanding of the whole chain of forming processes. Indeed, residual stresses are generated at each stage of the production line and are responsible for flatness defects. Heterogeneous plastic deformations taking place for instance during the rolling process and phase transitions occurring during cooling on the run out table and during the coiling process (see figure 1) are the most significant contributions to the residual stress formation. Therefore, numerical simulations are needed to understand and quantify all mechanisms involved in the different processes. Thus, significant efforts have been made for the simulation of the rolling process [START_REF] Hacquin | Modelisation thermo-mecanique tridimensionnelle du laminage: couplage bande-cylindres [3D thermomechanical modelling of rolling processes: coupling strip and rolls[END_REF][START_REF] Montmitonnet | [END_REF], the run out table [3][4][5] and the coiling process [6][7][START_REF] Weisz-Patrault | International Journal of Solids and Structures[END_REF].

This paper focuses on the evolution of residual stresses during cooling on the run out table (after hot rolling) that is a very non linear and multiphysics process. A cooling pattern is imposed to the strip on the run out table in order to obtain a specific cooling path and phase transitions and thus a targeted microstructure may be obtained [START_REF] Waterschoot | [END_REF]. Several papers [10][11][12][13][14][15] deal with the mathematical modeling of this coupled thermal/metallurgical problem. Indeed, phase transitions are responsible for a source term in the heat conduction problem because of the enthalpy change. Most previously cited papers consider one or two-dimensional unsteady heat conduction problem solved by mixed Finite Differences for the time variable and Finite Elements for space variables. On the other hand phase transitions are computed by using the well-known Avrami's equation proposed in [16]. Classic couplings between heat conduction, phase transition and mechanics are presented in figure 2. Thus, most estimations of residual stresses [17,18,[3][4][5] relies on highly coupled computations even though some couplings are neglected, for instance the effect of strain and stress on thermal and phase transition problems (dotted lines in figure 2). Despite the fact that the cited works present different degrees of details, the most common numerical strategy consists in developing a user material (UMAT) in the Finite Element software Abaqus [START_REF] Abaqus | Abaqus[END_REF] in order to solve simultaneously the heat conduction problem, phase transitions (via Avrami's equation) and the mechanical problem. However, the complete and fully coupled problem is very non-linear and the applicability may be limited because of very long computation times.

For instance [4] can simulate a 6 meters long strip although the run out table is 130 meters long.

The present paper aims at developing a numerical strategy sufficiently fast to enable parametric studies in order to quantify in what extend industrial conditions can be adjusted to optimize flatness defects. The proposed approach relies on the very common assumption in this situation that the mechanical problem has no influence on the thermal problem and phase transitions (dotted lines in figure 2). Indeed, the problem is driven by temperature changes. Phase transitions due to mechanical deformations are negligible as well as self-heating. However, displacements obtained from the mechanical problem can be significant and affect the thermal/metallurgical problem, principally when instabilities occur if mechanical computations are done at finite strain. In this contribution, computations are done under infinitesimal strain assumption and then the effect of geometry changes (on the thermal/metallurgical problem) is neglected. Thus, two numerical calculations can be conducted successively, first the coupled thermal/metallurgical problem and then the mechanical problem. This approach enables us to use different numerical schemes and different geometrical descriptions for each one of the two problems. The global numerical scheme is divided into three different numerical items. First the coupled thermal/metallurgical problem is solved by using TACSI proposed in [START_REF] Iung | [END_REF][21][START_REF] Huin | International Conference on Thermomechanical Processing of Steels-TMP[END_REF]. Thus, temperature and phase proportion fields are extracted at each time step. Thermal expansion, density mismatch between phases and transformation induced plasticity are combined in the form of an eigenstrain increment that is computed by using the free software SCILAB [23].

This eigenstrain being incompatible in general (i.e., it cannot be written as the gradient of a displacement field) an elastic strain should be introduced so that the total strain is compatible. This mechanical problem is solved by using the Finite Element software CAST3M [START_REF]Commissariat a l'Energie Atomique[END_REF].

The paper is constructed as follows. The global numerical strategy is detailed in section 2.

Calculation of the eigenstrain increment is then broached in section 3. In particular transformation induced plasticity is discussed. The mechanical problem of a strip subjected to the eigenstrain increment is then addressed in section 4. Material parameters are extracted from the literature as a function of temperature and representative run out table conditions are selected in section 5.1. Numerical computations are presented and analyzed in section 5.2. A post-processing is also included in section 6 in order to quantify the flatness defect that would take place if the strip were cut off the production line. It consists in computing at finite strain the relaxation of residual stress when the tension applied by the coiler is released. Conclusive remarks are proposed in section 7.

Numerical strategy

The global numerical scheme proposed in this paper is summarized in figure 3. The general problem of the strip cooling down on the run out table is divided into three numerical tasks.

1) The first task consists in solving the coupled heat conduction problem with phase transitions by using TACSI proposed in [START_REF] Iung | [END_REF][21][START_REF] Huin | International Conference on Thermomechanical Processing of Steels-TMP[END_REF]. This program is based on the unsteady heat equation with source terms due to enthalpy changes (phase transitions) given by:

div (λ(T )∇T ) -ρ(T )c p ∂T ∂t = - N p p=2 ∆H p (T ) Ẋp (1)
where T denotes the temperature, λ(T ) the thermal conductivity, ρ(T ) the density of the multiphase mixture, c p the thermal capacity, ∆H p (T ) the enthalpy change during phase transformation, X p the phase proportion of p-th phase and N p the number of considered phases. Time dependent boundary conditions are imposed through various Heat Transfer

Coefficients that describe the known process conditions (water cooling, wet and dry zones).

TACSI relies on solving (1) in 1D along the strip thickness. That is to say that a material particle (described by a segment along the strip thickness) is cooled down by moving to the run out table. Explicit Finite Differences are used and the coupling procedure consists in alternating at each time step the heat conduction solver and the phase transition solver.

Phase transitions are most of the time described by Avrami's equation that relies on statistical treatments and assumptions (e.g., pre-existence of germ nuclei of the product phase with a certain probability of growth) and predicts the overall transformation kinetics for diffusive transformations under isothermal assumption. For non-isothermal conditions, the so-called isokinetic assumption has been extended in [START_REF] Cahn | [END_REF] leading to a modified equation.

This approach is accurate but relies on experimental identifications that should be done for each tested grade. Physical approaches on the other hand, do not require a full experimental fitting procedure for each grade. Nucleation and growth models have been been proposed as reviewed for instance in [START_REF] Pan | Measurement and Modelling of Diffusional Transformation of Austenite in C-Mn Steels[END_REF]. The phase transition model included in TACSI relies on such a physical modeling. The nucleation and growth are modeled by thermodynamical and kinetics approaches. For details the reader is refereed to [START_REF] Iung | [END_REF].

At each time step a coupling with phase transitions is computed and the temperature and phase proportion fields are computed at three positions through the strip thickness (upper and lower surface and mid-plane). Moreover, a steady state is assumed on the run out table, that is to say that the global temperature field does not vary along the strip length even though each material particle undergoes very rapid thermal and structural evolutions.

Since neither the geometrical strip profile (obtained during the rolling process) nor the cooling conditions on the run out table are perfectly constant, no steady state is never strictly reached. However, temperature and phase proportion fields evolve much slower than the strip speed. Therefore, one can assume that by selecting several points along the strip width at one particular position along the strip length, one can construct temperature and phase proportion fields for a strip portion as long as the the run out table. This simplified approach enables us to compute effectively the whole temperature and phase proportion fields just by computing the evolution uni-dimensional material particles (at different position along the strip width) moving to the run out table.

2) The second task consists in computing for each time step the eigenstrain increment (denoted by ∆ε * k where k refers to the time increment) accounting for thermal expansion, phase transformations and transformation induced plasticity. To that end a numerical code has been written by using the free software SCILAB [23].

3) The third task consists in solving at each time step the elastic-plastic problem of the strip undergoing the eigenstrain increment ∆ε * k by using the Finite Element software CAST3M [START_REF]Commissariat a l'Energie Atomique[END_REF]. Indeed, the eigenstrain is not the gradient of a displacement field and an elastic strain should be introduced so that the total strain is compatible. A global plastic strain is also considered in the total strain because the eigenstrain increment can be sufficient to generated macroscopic plasticity. Moreover an initial residual stress profile (obtained at the end of the rolling process) is considered and the traction applied by the coiler is also taken into account. Finite Element computations are done under infinitesimal strain assumption by using quadrangular shell elements in order to obtain short computation times.

It should be noted that the mechanical problem depends on the loading path both because of the fact that the eigenstrain increment depends on the stress state of the previous time step and because of the elastic-plastic behavior. Thus, one cannot compute directly the mechanical problem from the temperature and phase proportion fields obtained on the run out table. Strain and stress history should be considered by applying progressively the temperature and phase proportion fields responsible for the eigenstrain as shown in figure 4. As a result, several FE computations are performed alternatively with computations of the eigenstrain increments.

Eigenstrain and transformation induced plasticity

During cooling on the run out table significant temperature and structural changes take place under mechanical loads (applied tension by the coiler and initial residual stresses). Thus, the strip undergoes thermal expansion, phase transformation and transformation induced plasticity. Indeed, when phase transformations occur under applied mechanical loads that can be much lower than the yield stress, significant plastic strains are observed. This phenomenon called transformation induced plasticity or sometimes super-plasticity received particular attention. The Leblond's model proposed in [START_REF] Leblond | [END_REF][28][29][30] is one of the most fruitful model for industrial implementation. Indeed, because of its solid theoretical basis and the fact that the overall transformation induced plastic strain increment is given by simple explicit formulas, the Leblond's model has been intensively used for various engineering applications and especially for run out table simulations. For instance [4] used a slightly improved version of the Leblond's model proposed in [31]. Within the framework of run out table simulations, a simple formulation for transformation induced plasticity has also been proposed in [32]. However, there is no theoretical derivation from previous mechanical knowledge and the proposed formula is an Arrhenius equation to be fitted with experimental data. This approach has been used for instance in [17].

Most transformation induced plasticity models are dedicated to simple two-phases mixtures. However, multiphase transitions occur in steel during cooling. Residual austenite may be transformed into ferrite, pearlite, bainite and martensite depending on the local temperature kinetics. An attempt to extend the Leblond's model to multiphase transitions has been proposed in [33] but formulas are adapted without theoretical developments as originally made in the Leblond's model. In this paper transformation induced plasticity and classic plasticity related to temperature and stress variations are computed following the recent theoretical work [34] that extends properly the original Leblond's model. Multiphase transitions are taken into account and several assumptions in the original work are relaxed.

Main results useful for this paper are recalled for sake of clarity. It should be noted that expressions are particularized for this study, that is to say during cooling only and for the following phase transitions austenite to ferrite, austenite to pearlite, austenite to bainite and austenite to martensite. For more general equations, the reader is refereed to [34]. At each time step k, the elastic-plastic problem consists in solving the following equation:

                         div (σ k ) = 0 Equilibrium σ k = λ(T k )tr ε e k I + 2µ(T k )ε e k Isotropic behavior ε k = 1 2 ∇u k + ∇u T k Compatibility ε k = ε e k + ε p k + ε k-1 + ∆ε * k Total strain tensor (2)
T k is the temperature field at the k-th increment and λ(T k ) and µ(T k ) are the temperature dependent Lamé's coefficients of the multiphase mixture. In this paper, a simple mixture rule is chosen:

λ(T k ) = N p p=1 X p,k λ p (T k ) and µ(T k ) = N p p=1 X p,k µ p (T k ) (3)
Where λ p and µ p are the Lamé's coefficients of phase p and X p,k denote the phase proportion of phase p and N p is the number of considered phases (p = 1: austenite, p = 2: ferrite, p = 3: pearlite, p = 4: bainite and p = 5: martensite). Moreover, u denotes the displacement vector, σ k is the Cauchy stress tensor, ε k and ε k-1 denote the total strain at the k-th and (k -1)-th increments respectively and ε e and ε p denote the unknown elastic and plastic strain respectively.

It should be noted that the initial residual stress profile can be either set with an initial stress field σ 0 in CAST3M or by imposing that ε 0 corresponds to the elastic counter part of the initial residual stress profile.

In addition, ∆ε * k denotes the eigenstrain increment introduced because of thermal expansion, phase transitions and transformation induced plasticity associated with increment of temperature ∆T k and increment of phase proportion ∆X p,k . The eigenstrain increment ∆ε * k can be decomposed as follows:

∆ε * k = ∆ε thm k + ∆ε tp k + ∆ε cp T,k + ∆ε cp σ,k (4)
Where ∆ε thm k denotes the thermo-metallurgical strain tensor, which corresponds to volume variations due to thermal dilatation and density mismatch between phases, ∆ε tp k is the transformation induced plasticity tensor and ∆ε cp σ,k and ∆ε cp T,k are the classic plasticity tensors induced by equivalent stress and temperature variations respectively. The contribution of ∆ε cp σ,k is neglected.

After [34], the hydrostatic thermo-metallurgical strain increment is shown to be:

∆ε thm = N p p=2 1 3 ρ 1 (T k ) ρ p (T k ) -1 ∆X p,k +         α 1 + N p p=2 X p,k (α p -α 1 )         ∆T k (5)
Where ρ p is the density of phase p and α p is the thermal expansion coefficient of phase p.

Consider X k the total product phase proportion:

X k = 1 -X 1,k = N p p=2 X p,k (6)
and ∆σ Y the following quantity:

∆σ Y = σ Y 1 2 -σ eq 1 2 (7)
where σ Y p is the yield stress, σ eq p the von Mises equivalent stress and ε eq p the von Mises equivalent plastic strain in phase p. Consider the following material parameters:

             ζ = (3λ 1 + 2µ 1 )2µ 1 λ 1 + 2µ 1 ξ p = 2µ 1 (9λ 1 + 14µ 1 ) λ 1 (9µ 1 + 6µ p ) + 2µ 1 (7µ 1 + 8µ p ) (8)
The transformation plastic strain increment reads if σ eq 1 < σ Y 1 :

∆ε tp = N p p=2 s p,k σ Y p σ Y p -σ eq p µ p ξ p ∆X p,k +                                          0 if ∆σ Y ζ ε thm > 1 - 3 ε thm s 1,k σ Y 1 (ε eq 1 ) ln        ∆σ Y ζ ε thm        N p p=2 ∆X p,k >0 ∆X p,k if X k ≤ ∆σ Y ζ ε thm ≤ 1 - 3 ε thm s 1,k σ Y 1 (ε eq 1 ) ln X k N p p=2 ∆X p,k >0 ∆X p,k if X k > ∆σ Y ζ ε thm (9)
Where s p,k is the deviatoric stress tensor of phase p. And if σ eq 1 = σ Y 1 : (10) with the initial value:

∆ε tp = N p p=2 s p,k σ Y p σ Y p -σ eq p µ p ξ p ∆X p,k +                    0 if X k < X min - 3 ε thm s 1,k σ Y 1 ln X k N p p=2 ∆X p,k >0 ∆X p,k if X k ≥ X min
ε tp ini = 2 ε thm 1 -X min X min ln X min (11)
where X min (set to 0.03 in this paper) is the minimal phase proportion that can be produced from pure austenite and where ε thm is the hydrostatic eigenstrain generated in all product phases because of temperature variation and phase transitions:

ε thm = N p p=2 X p,k 1 -X 1 1 3 ρ 1 (T k ) ρ p (T k ) -1 + ρ 1 (T k ) ρ p (T k ) (α p -α 1 )(T k -T ini ) > 0 ( 12
)
It should be mentioned that ε thm given by ( 12) is related to the microstructure and should not be mixed up with the homogenized hydrostatic strain ε thm given by (5). The classic plastic strain rate due to thermal variations reads if σ eq 1 < σ Y 1 :

∆ε cp T =                              0 if ∆σ Y ζ ε thm > 1 -3 αs 1,k σ Y 1 (ε eq 1 ) X k ln        ∆σ Y ζ ε thm        ∆T k if X k ≤ ∆σ Y ζ ε thm ≤ 1 -3 αs 1,k σ Y 1 (ε eq 1 ) X k ln X k ∆T k if X k > ∆σ Y ζ ε thm (13) And if σ eq 1 = σ Y 1 : ∆ε cp T = -3 αs 1,k σ Y 1 X k ln X k ∆T k ( 14
)
where α represents the homogenized difference between thermal expansion coefficient in the product phases and in the austenite:

α = N p p=2 X p,k 1 -X 1 α p -α 1 < 0 ( 15
)
4. Elastic-plastic problem with imposed eigenstrain

Generalized stess and strain

Equation ( 2) is solved at each time step under infinitesimal strain assumption, by using the Finite Element free software CAST3M. Quadrangular linear shell element 'COQ4'is used in order to obtain short computation times. Consider Cartesian coordinates (x 1 , x 2 , x 3 ) corresponding respectively to the direction of strip length, strip width and strip thickness. There are 5 generalized kinematic unknowns (degrees of freedom) for this element: three translations U 1 , U 2 , U 3 and two rotations ϕ 1 , ϕ 2 . Moreover the mechanical state relies on 8 generalized stress components. Thus the in-plane stress resultants N, bending moments M and shear resultants Q are defined as follows:

∀(i, j) ∈ {1, 2} 2 ,                                N i j = h 2 -h 2 σ i j dx 3 M i j = h 2 -h 2 x 3 σ i j dx 3 Q i = h 2 -h 2 σ i3 dx 3 ( 16
)
where h is the strip thickness that depends on x 2 because of the strip geometrical profile. From the computed quantities N i j , M i j , Q i using the FE computation, one should compute an estima-tion of three dimensional Cauchy stresses in order to compute the imposed eigenstrain increment ∆ε * for the next time step (as detailed in section 3). By assuming a linear dependence on the thickness direction x 3 for σ i j and homogenous σ i3 (for (i, j) ∈ {1, 2} 2 ) one obtains:

∀(i, j) ∈ {1, 2} 2 ,              σ i j = 12M i j h 3 x 3 + N i j h σ i3 = Q i h and σ 33 = 0 (17)
The generalized eigenstrain to be imposed in the FE model, is constituted of 8 components (plane and shear eigenstrain and curvature) denoted by ∆E * i j , ∆G * i3 , ∆R * i j . This generalized eigenstrain is given as a function of the three dimensional eigenstrain computed following section 3 as follows:

∀(i, j) ∈ {1, 2} 2 ,                                ∆E * i j = 1 h h 2 -h 2 ∆ε * i j dx 3 ∆G * i3 = 1 h h 2 -h 2 ∆ε * i3 dx 3 ∆R * i j = 12 h 3 h 2 -h 2 x 3 ∆ε * i j dx 3 (18)

Boundary conditions

Boundary conditions are specified as follows. At one end of the strip an average tension is applied by the coiler and displacement U 3 is blocked. Furthermore, displacements U 1 and U 3 and all rotations are blocked at the other end. However, edge effects should be avoided. Indeed, the local distribution of the applied tension is unknown, only the resultant force applied by the coiler is known. Similarly, displacements are not blocked at the other end but a distribution takes place instead. This difficulty is overcome by adding extra pieces of strip (whose length is twice the strip width) at both ends of the strip in order to apply boundary conditions. Thus, mechanical fields unduly affected by simplified boundary conditions does not spread in the real zone of interest. It should be mentioned that in the following results are presented only in this zone of interest.

Material parameters and modeling conditions

A numerical example is proposed in order to present typical outputs of the model. A 140 m long run out table is considered. A steel strip initially at 1150 K and going at 10.25 m/s is cooled down on the run out table with water-jets identically distributed along the strip width.

The cooling path along the strip length is determined by a succession of water, wet and dry zones characterized by Heat Transfer Coefficients. The upper surface is cooled down faster than the lower surface because there are more water jets spraying on the upper surface, as

shown in figure 6. The geometrical profile defines the strip thickness as a quadratic function of the strip width (i.e., x 2 coordinate). In this example, a 3 mm thick strip is considered with a strip crown of 2% (i.e., the strip thickness at the edges is 98 % of the strip thickness at the center). The temperature and phase proportion fields are presented in figure 7a and 7b at three positions along the strip thickness (top, middle and bottom). There are almost no variations along the strip width because the cooling is homogeneous along this direction. However it is clear that the upper surface is cooled down faster than the lower surface. Material parameters are extracted from the literature. The Poisson coefficient is set to 0.3 and the Young modulus is assumed to be identical for all phases and given by [START_REF] Martinez | Jonction 16 MND 5-Inconel 690-316 LN par soudage-diffusion: élaboration et calcul des contraintes résiduelles de procédé[END_REF]:

E(T ) = 2.08 × 10 5 -1.90 × 10 2 T + 1.19T 2 -2.82 × 10 -3 T 3 + 1.66 × 10 -6 T 4 (19)
On the other hands, temperature dependent yield stresses are given for each phase in figure 8a by interpolating numerical values found in [START_REF] Sysweld® | [END_REF].

Results

Numerical results are given at different time steps for the in-plane stress resultant N 11 and the bending moment M 11 in figure 9 and 10. Six other components are calculated but for sake of concision only the most interesting components are presented. The effect of plasticity and transformation induced plasticity is clearly in this example to reduce significantly the residual stress state. This is due to the fact that temperature and phase proportion fields are almost homogenous along the strip width. However, since the cooling path is not symmetric (upper surface cooled down faster than the lower surface), bending moments are generated. It can be noticed that there is a maximum zone in figure 10 that shifts according to time steps. This is only due to the equilibrium because of the quick variation of the in-plane stress resultant N 11 .

The final residual stress profile, that takes place in the strip, is obtained at the last time step at the end of the run out table. Stresses at three positions along the strip thickness are presented in figure 11 as well as the corresponding bending moments. Therefore by comparing initial and final residual stress profiles in figure 8b and 11 , it is clear with the proposed simulation parameters (e.g., cooling path is homogenous along the strip width) that plasticity tends to decrease significantly the stress level. However, bending moments are also introduced and are responsible for flatness defects called longbow and crossbow. It should be noted that the resultant of σ 11 has been subtracted in figure 11a in order to remove the offset due to the applied tension even though at this stage there is no computation of the the residual stress relaxation by properly releasing the tension.

The entire simulation lasts aroud 2 hours and a non-negligible part of the computation time is used for writing files characterizing the mechanical state in order to compute the increment of eigenstrain. Thus, a way of reducing even more computation times is to avoid to write these files by using the flash memory instead.

Relaxation of residual stresses

In this section a post-processing is proposed in order to compute out of plane displacements due the relaxation of the residual stress profile obtained previously. In order to evaluate the influence of the run out table only, the final residual stress profile obtained at the end of the run out table (i.e., x 3 = 140) is applied on the whole 140 m long strip as though there were no subsequent evolution of the residual stress after the run out table (which is obviously not verified). A Finite Element simulation at finite strain is performed in order to quantify how residual stresses are relaxed by out of plane deformation. Thus, U 3 is not blocked any more and the tension is released. In figure 12 the residual stress profile after releasing the applied tension is presented. In comparison with figure 11 there is a clear relaxation of residual stresses through out of plane deformations. The out of plane displacement is presented in figure 13. It corresponds to the displacement that would occur if the strip were cut off just after the run out table. There is a longbow defect at both ends of the strip. Of course, the strip cooling continues during the coiling process and the residual stresses evolve as studied in [6][7][START_REF] Weisz-Patrault | International Journal of Solids and Structures[END_REF].

Conclusion

This contribution deals with a numerical strategy to simulate a non-linear and multiphysics problem: the cooling of a steel strip on the run out table. The proposed approach accounts for heat conduction strongly coupled with metallurgical transformations (multiphase transitions) and an elastic-plastic calculation of residual stress evolution. Thermal expansion, density mismatch between different phases and transformation induced plasticity are included under the form of an imposed increment of eigenstrain. This work enables to obtain reasonable computation times (around 2 hours for a 140 m long strip) because the mechanical problem is solved separately from the thermal/metallurgical problem and authorizing different meshing and time steps.

Numerical results show that homogenous cooling conditions along the strip width tend to decrease significantly the initial residual stress profile, but cooling one surface faster than the other creates an asymmetry that is responsible for bending moments. Therefore, classic flatness defects such as longbow and crossbow are detected.

Acknowledgements

Authors thank Didier Huin (ArcelorMittal) for their expertise and fruitful discussions.

Figure 1 :

 1 Figure 1: Chain of hot forming processes

Figure 2 :

 2 Figure 2: Couplings

Figure 3 :Figure 4 :

 34 Figure 3: Numerical scheme

Figure 5 :

 5 Figure 5: Boundary conditions

Figure 6 :

 6 Figure 6: Prescribed cooling path

Figure 7 :

 7 Figure 7: Coupled thermal/metallurgical problem

Figure 9 : 12 Figure 11 :

 91211 Figure 9: In-plane stress resultant N 11

(a) Relaxed stress σ 11 (b) Relaxed bending moment M 11 Figure 12 :Figure 13 :

 11111213 Figure 12: Relaxed stress profiles

Figure 1 :Figure 2 : Couplings 1 FigureFigure 3 :Figure 4 :Figure 5 :Figure 6 :

 1213456 Figure 1: Chain of hot forming processes

Figure 7 :

 7 Figure 7: Coupled thermal/metallurgical problem

Figure 11 : 11 Figure 12 :Figure 13 :

 11111213 Figure 11: Final stress profiles

Figure 9: In-plane stress resultant N 11

Figure 10: Bending moment M 11


Daniel Weisz-Patrault

Nonlinear and multiphysics evaluation of residual stresses in coils

Keywords: Residual stress, Multiphysics, Contact problem, Coiling process

This article deals with residual stress evaluation within the framework of numerical simulation of the coiling process of steel. Plastic deformations along with multiphase transitions are responsible for large irreversible strain leading to major residual stress issues. A nonlinear mixed analytical/numerical approach is developed to compute residual stresses generated by different contributions of inelastic strain occurring during the coiling process (including both winding and cooling). In particular, transformation-induced plasticity is taken into account. Contact stresses at the interfaces between layers are updated at each time step by a minimization procedure. The analytical part of the proposed strategy is validated by comparison with a finite element computation. Then a coil composed of 196 layers is considered under typical industrial conditions. Numerical results are discussed to present the industrial application of the proposed approach.

Introduction

The current dynamic of steel manufacturing is to regularly develop new stronger grades enabling users to reduce strip thicknesses and thus reduce produced tonnages, which results in energy efficiency by minimizing, for instance, the total mass of vehicles. One of the major issues related to this evolution of steel production is the forming processes (see Figure 1) that lead to serious residual stress problems, which in turn can result in instabilities such as local buckling during the rolling process [START_REF] Abdelkhalek | Coupled approach for flatness prediction in cold rolling of thin strip[END_REF]. This contribution focuses on the coiling process, which consists in winding a strip on itself for storage. Within this framework, residual stress issues can result in local buckling or coil sagging, as demonstrated in Figure 2. Other defects are detailed in [START_REF] Edwards | The mystery of coil winding[END_REF]. When the inner mandrel is removed from the coil, residual stresses may be sufficiently compressive to cause buckling of some layers. Furthermore, the overall inelastic strain due to thermal expansion and phase transitions may be responsible for large displacements leading to coil sagging. In this article I attempt to quantify the overall inelastic strain as well as the associated residual stress by modeling the whole coiling process. The outcomes are necessary for the prediction of some defects such as coil sagging even though additional work is still needed.

Residual stresses generated during this process have been studied, for instance, by Quach et al. [START_REF] Quach | Residual stresses in steel sheets due to coiling and uncoiling: a closed-form analytical solution[END_REF][START_REF] Quach | Residual stresses in press-braked stainless steel sections, ii: Press-braking operations[END_REF], who developed a closed-form analytical solution of a pure bending problem (i.e., disregarding the effect of contacts). The contribution of several successive processes was considered by Moen et al. [START_REF] Moen | Prediction of residual stresses and strains in cold-formed steel members[END_REF] also for a pure bending problem. However, contact pressures are heterogeneous along the coil axis, which results in significant variations of residual stresses. Contact pressures can concentrate near the central axial position and significantly contribute to the formation of plastic deformations and therefore to the final residual stress field. Several coiling models that take into account contacts have been proposed. One can mention purely elastic models developed by Edwards and Boulton [START_REF] Edwards | The mystery of coil winding[END_REF], Hudzia et al. [START_REF] Hudzia | Stress calculation applied to a coil, and optimization of coiling tension[END_REF], and Kedl [START_REF] Kedl | using a two dimensional winding model to predict wound roll stresses that occur due to circumferential steps in core diameter or to cross-web caliper variation[END_REF], where the radial Young's modulus is set as a function of contact pressure to model roughness of contacts. One can also mention a fast simplified coil winding model based on the Airy function proposed by Hinton et al. [START_REF] Hinton | A Study on the Effects of Coil Wedge During Rewinding of Thin Gauge Metals[END_REF][START_REF] Hinton | An airy function to rapidly predict stresses in wound metal strip having asymmetric thickness profile[END_REF]. However, the coiling process consists not only in winding the strip around a mandrel but also in cooling the coil for several hours after removal of the internal mandrel. Thus, contact pressures are completely redistributed and thermal expansion and phase transitions occur. Previous studies dealt only with the isothermal winding phase of the process and ignored the cooling phase. However, it will be demonstrated here that residual stresses are significantly modified during the cooling phase. Recently developed grades undergo significant phase transitions when the coil cools down. Additional inelastic strain occurs because of thermal expansion, phase density mismatch, and transformation-induced plasticity. Thus, this article contributes to the development of a global numerical strategy that solves the nonlinear and coupled problem of a strip winding on itself and then cooling by undergoing phase transitions. Couplings between thermal sciences, metallurgy, and solid mechanics should be considered as well as several types of nonlinearity, such as finite strain computations, elastic-plastic behavior, and contact stress evaluation during the winding phase of the process, and temperature-dependent heat sources because of enthalpy changes and transformation-induced plasticity when the coil cools down.

This article relies on previous contributions. Weisz-Patrault et al. [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF] proposed a nonlinear isothermal elasticplastic model of the winding phase accounting for roughness of contacts and finite strain to consider large rotations. This model was based on previous work by Weisz-Patrault et al. [START_REF] Weisz-Patrault | Imposed curvature of an elastic-plastic strip: application to simulation of coils[END_REF][START_REF] Weisz-Patrault | Non-linear numerical simulation of coiling by elastic finite strain model[END_REF]. Then Weisz-Patrault [START_REF] Weisz-Patrault | Coupled heat conduction and multiphase change problem accounting for thermal contact resistance[END_REF] developed a solution for the cooling phase by coupling heat conduction and multiphase transition problems accounting for thermal contact resistance. The present article aims to establish a numerical strategy to compute residual stresses on the basis of the previously developed models. Each layer of the coil is subjected at each time step to an eigenstrain increment composed of several contributions, namely:

(i) the thermometallurgical hydrostatic strain due to thermal expansion and density mismatch between different phases; (ii) the transformation-induced plastic strain; (iii) classical plastic strain due to temperature variations.

The evaluation of the eigenstrain increment necessitates the computation of temperature and phase proportion rates that are obtained through the method proposed by Weisz-Patrault [START_REF] Weisz-Patrault | Coupled heat conduction and multiphase change problem accounting for thermal contact resistance[END_REF]. Furthermore, the macroscopic stress tensor is also needed for the computation of the eigenstrain increment of the following time step. Thus, the residual stress field should be computed at each time step as an input of the following time step, and one cannot jump directly to the last time step by simply integrating the time-dependent eigenstrain. For the exact same reason, contact stresses are also updated at each time step by use of minimization procedures as detailed in Section 7. The general modeling structure is presented in Figure 3.

Macroscopic plasticity is not considered in this contribution; that is, that the macroscopic equivalent stress does not reach the macroscopic yield stress. Although macroscopic yield stress is not reached, certain phases (e.g., austenite) can locally undergo plastic deformations. This local plasticity is taken into account in the eigenstrain obtained from the homogenization procedure proposed by Weisz-Patrault [START_REF] Weisz-Patrault | Multiphase model for transformation induced plasticity. extended leblond's model[END_REF]. The nonlinear part is therefore limited in this article to the imposed eigenstrain, which should be computed incrementally, and to the contact stress evaluation. Thus, the present work consists in computing all contributions of the overall eigenstrain. Then an analytical solution of the elastic problem subjected to the eigenstrain is sought for each layer of the coil and at each time step. In addition, contact stresses at each interface between layers are also updated at each time step by use of a numerical minimization procedure relying on contact laws accounting for roughness, following the ideas of Sutcliffe [START_REF] Sutcliffe | Surface asperity deformation in metal forming processes[END_REF] or Sheu and Wilson [START_REF] Sheu | Mixed lubrication of strip rolling[END_REF].

Mathematical developments rely, on the one hand, on a series expansion with respect to a certain function sequence and, on the other hand, on Papkovich-Neuber potential theory involving harmonic and biharmonic functions as presented for, instance, by Barber [START_REF] Barber | Elasticity[END_REF]. One can also mention a very elegant potential theory mainly developed by Muskhelishvili [START_REF] Muskhelishvili | Some basic problems of the mathematical theory of elasticity[END_REF] applicable for plane elasticity and relying on complex potentials. An extension in quaternionic algebra dedicated to three-dimensional elasticity was proposed by Weisz-Patrault et al. [START_REF] Weisz-Patrault | Three-dimensional elasticity based on quaternion-valued potentials[END_REF], but is not used in this contribution for simplicity. The main outcome of the proposed approach is the computation of residual stress and inelastic strain in the coil. The work of Weisz-Patrault et al. [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF] and Weisz-Patrault [START_REF] Weisz-Patrault | Coupled heat conduction and multiphase change problem accounting for thermal contact resistance[END_REF] and the present article contribute to the development of an effective tool that enables users to perform parametric studies, especially for recent grades for which residual stress issues arise frequently. To the best of my knowledge, this article presents for the first time a complete approach enabling users to track irreversible strains and residual stresses through the whole coiling process, including the winding phase and the cooling phase. The situation before the coiling process is not considered in this contribution, and can be directly measured through an optical or mechanical measurement system as reviewed by Molleda et al. [START_REF] Molleda | On-line flatness measurement in the steelmaking industry[END_REF]. An inverse method based on conformal mapping techniques was proposed by Weisz-Patrault [START_REF] Weisz-Patrault | Inverse cauchy method with conformal mapping : application to latent flatness defect detection during rolling process[END_REF] to evaluate residual stresses by deflection of the steel strip on a hollow cylinder.

General equations and notation

This section details the global equation system that is solved at each time step (denoted by n). The main notation is summarized in Figure 4. The coil is constituted of N layers. Some assumptions are needed to obtain a convenient analytical solution. Each layer (denoted by (i), where 1 ≤ i ≤ N) is characterized at each time step n by lower and upper surface radii denoted, respectively, by R (i) in f,n (z) and R (i) sup,n (z) (where r and z denote the cylindrical coordinates). They clearly depend on the axial position z. However, the following analytical solution relies on z-independent lower and upper surface radii, denoted by r (i) in f,n and r (i) sup,n , defined as the average values of R (i) in f,n (z) and R (i) sup,n (z) over the coil width:

r (i) in f,n = 1 2L L -L R (i) in f,n (z)dz and r (i) sup,n = 1 2L L -L R (i) sup,n (z)dz (1)
where L denotes the half width. Thus, each layer is described by

Ω (i) n = (r, z) ∈ r (i) in f,n , r (i) sup,n × [-L, L] (2)
In addition, T (i) n (r, z) denotes the temperature field in the ith layer at the current time step n. In this article a multiphase mixture is considered. Phase proportions are denoted by X (i) p,n (with 1 ≤ p ≤ N p , where N p is the number of phases). Thus, Lamé's coefficients (λ (i) p,n (T (i) ), µ (i) p,n (T (i) )) are temperature dependent and therefore inhomogeneous along both the radial direction and the axial direction. The analytical solution relies on constant material properties in each layer (λ (i) p,n , µ (i) p,n ), which are computed as follows:

                     λ (i) p,n = 1 2L r (i) sup,n -r (i) in f,n r (i) sup,n r (i) in f,n L -L λ (i) p,n (T (i) n (r, z))drdz µ (i) p,n = 1 2L r (i) sup,n -r (i) in f,n r (i) sup,n r (i) in f,n L -L µ (i) p,n (T (i) n (r, z))drdz (3)
In addition, Lamé's coefficients of the multiphase mixture (λ (i) n , µ (i) n ) are computed as follows:

λ (i) n = N p p=1 X (i) p,n λ (i) p,n and µ (i) n = N p=1 X (i) p,n µ (i) p,n (4)
At each time step, an additional eigenstrain increment δε * ,(i) n is imposed. Thus, an additional elastic problem should be solved so as to update the previous mechanical state. The global additional problem is written as follows:

∀ 1 ≤ i ≤ N,                      div σ (i) n = 0 (equilibrium) σ (i) n = λ (i) n tr ε e,(i) n I + 2µ (i) n ε e, (i) n
(isotropic behavior)

ε (i) n = 1 2 ∇ u (i) n + ∇ u (i) n T (compatibility) ε e,(i) n = ε (i) n -δε * ,(i) n (elastic strain tensor) (5)
where σ (i) n , ε (i) n , and u (i) n denote the additional stress and strain tensors and the additional displacement vector, respectively. The total stress accumulated from all the previous time steps is denoted by

Σ (i) n-1 = n-1 l=0 σ (i) l ( 6
)
where σ (i) 0 is the initial stress tensor from previous processes as detailed in Figure 3. Furthermore, the eigenstrain increment δε * ,(i) n depends on Σ (i) n-1 as detailed in Section 3. It is decomposed into different contributions:

δε * ,(i) n = δε thm,(i) n + δε tp,(i) n + δε cp,(i) n (7)
where δε thm,(i) n represents the thermometallurgical strain increment, which corresponds to volume variations due to thermal expansion and density mismatch between phases, δε tp,(i) n is the strain increment due to transformationinduced plasticity, and δε cp,(i) n is the classical plastic strain increment induced by temperature variations in the multiphase mixture. Classical plasticity induced by stress variations can also be considered but it is usually negligible.

In addition, the boundary conditions for the additional problem (5) are

∀ 1 ≤ i ≤ N,                        σ (i) n (r (i) sup,n , z).n(r (i) sup,n , z) = T (i) sup,n (z) σ (i) n (r (i) in f,n , z).n(r (i) in f,n , z) = T (i) in f,n (z) σ (i) n (r, ±L).n(r, ±L) = 0 T (i+1) in f,n (z) = -T (i) sup,n (z) T (N) sup,n (z) = T (1) in f,n (z) = 0 ( 8
)
where n is the unit normal vector and T (i) sup,n (z) and T (i) in f,n (z) denote the additional traction (at r = r (i) sup,n and r = r (i) in f,n , respectively), decomposed as follows:

∀ 1 ≤ i ≤ N,      T (i) sup,n (z) = -P (i) sup,n (z)e r + S (i) sup,n (z)e z T (i) in f,n (z) = P (i) in f,n (z)e r + S (i) in f,n (z)e z (9)
where P (i) sup,n (z) = P (i+1) in f,n (z) are the additional contact pressures and S (i) sup,n (z) = -S (i+1) in f,n (z) are the additional contact shear stresses that are evaluated by numerical minimization as detailed in Section 7.

Eigenstrain

This section is dedicated to the computation (for each layer (i) and at the current time step n) of the eigenstrain increment [START_REF] Kedl | using a two dimensional winding model to predict wound roll stresses that occur due to circumferential steps in core diameter or to cross-web caliper variation[END_REF] that contains most of the nonlinear information (along with contact stress evaluation). In all the following, the superscript (i) and the time step index n are omitted for readability. For instance, the displacement vector u (i) n is denoted by u in the rest of this article. Transformation-induced plasticity and classical plasticity are computed following the recent theoretical work of Weisz-Patrault [START_REF] Weisz-Patrault | Multiphase model for transformation induced plasticity. extended leblond's model[END_REF] that extends the well-known model proposed by Leblond et al. [START_REF] Leblond | Mathematical modelling of transformation plasticity in steels i: case of ideal-plastic phases[END_REF]. Multiphase transitions are taken into account, and several assumptions are relaxed. However, expressions are particularized for this study; that is, for cooling only and for the following phase transitions: austenite to ferrite, austenite to pearlite, austenite to bainite, and austenite to martensite. In addition, the following convention is used: p = 1 for austenite, p = 2 for ferrite, p = 3 for pearlite, p = 4 for bainite, and p = 5 for martensite. Equations ( 10)-( 24), extracted from [START_REF] Weisz-Patrault | Multiphase model for transformation induced plasticity. extended leblond's model[END_REF], are recalled here for clarity. The additional stress σ is distributed in different phases as follows:

σ p = E p E σ ( 10
)
where σ p is the additional stress in phase p and

               E = p∈E X p E p + p∈P X p γ p σ Y p E p = E p if p ∈ E γ p σ Y p if p ∈ P (11)
where for each phase p, σ Y p is the yield stress, E p is the Young's modulus, and γ p is the hardening parameter. In addition, E ⊂ 1, • • • , N p is the index set of phases that do not reach the plastic regime and P ⊂ 1, • • • , N p is the index set of phases that reach the plastic regime. The same expression as [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF] holds for the total stress accumulated from all the previous time steps Σ defined by ( 6):

Σ p = E p E Σ (12)
The eigenstrain increment is obtained as follows:

δε * = ε * δt ( 13
)
where δt is the time increment. The hydrostatic thermometallurgical strain rate is shown to be εthm

= N p p=2 1 3 ρ 1 (T ) ρ p (T ) -1 Ẋp +         α 1 + N p p=2 X p (α p -α 1 )         Ṫ ( 14
)
where ρ p and α p are, respectively, the density and the thermal expansion coefficient of phase p. Consider X to be the total product phase proportion:

X = 1 -X 1 = N p p=2 X p (15)
And consider ∆σ Y to be the distance between the yield stress and the von Mises equivalent stress in austenite:

∆σ Y = σ Y 1 2 -Σ eq 1 2 (16)
where Σ eq p is the von Mises equivalent stress in phase p computed from Σ p . Consider the following material parameters:

             ζ = (3λ 1 + 2µ 1 )2µ 1 λ 1 + 2µ 1 ξ p = 2µ 1 (9λ 1 + 14µ 1 ) λ 1 (9µ 1 + 6µ p ) + 2µ 1 (7µ 1 + 8µ p ) (17) If Σ eq 1 < σ Y 1 , the transformation-induced plastic strain rate is εtp = N p p=2 S p σ Y p σ Y p -Σ eq p µ p ξ p Ẋp +                                            0 if ∆σ Y ζ ε thm > 1 - 3 ε thm S 1 σ Y 1 ln        ∆σ Y ζ ε thm        N p p=2 Ẋp >0 Ẋp if X ≤ ∆σ Y ζ ε thm ≤ 1 - 3 ε thm S 1 σ Y 1 ln X N p p=2 Ẋp >0 Ẋp if X > ∆σ Y ζ ε thm (18)
where S p = Σ ptr Σ p 1/3 is the total deviatoric stress tensor of phase p.

If Σ eq 1 = σ Y 1 , εtp = N p p=2 S p σ Y p σ Y p -Σ eq p µ p ξ p Ẋp +                  0 if X < X min - 3 ε thm S 1 σ Y 1 ln X N p p=2 Ẋp >0 Ẋp if X ≥ X min (19)
with the initial value

ε tp ini = 2 ε thm 1 -X min X min ln X min (20)
where X min (set to 0.03 in this article) is the minimal phase proportion that can be produced from pure austenite and where ε thm is the hydrostatic eigenstrain generated in all product phases because of temperature variation and phase transitions:

ε thm = N p p=2 X p 1 -X 1 1 3 ρ 1 (T ) ρ p (T ) -1 + ρ 1 (T ) ρ p (T ) (α p -α 1 )(T -T ini ) > 0 (21)
ε thm given by ( 21) is related to the microstructure and should not be confused with the homogenized hydrostatic strain ε thm given by ( 14). If Σ eq 1 < σ Y 1 , the classical plastic strain rate due to thermal variations is

εcp =                              0 if ∆σ Y ζ ε thm > 1 -3 αS 1 σ Y 1 X ln        ∆σ Y ζ ε thm        Ṫ if X ≤ ∆σ Y ζ ε thm ≤ 1 -3 αS 1 σ Y 1 X ln X Ṫ if X > ∆σ Y ζ ε thm (22) If Σ eq 1 = σ Y 1 , εcp = -3 αS 1 σ Y 1 X ln X Ṫ ( 23
)
where α represents the homogenized difference between the thermal expansion coefficient in the product phases and the thermal expansion coefficient in austenite:

α = N p p=2 X p 1 -X 1 α p -α 1 < 0 (24)

Problem decomposition for each layer

This section deals with the general scheme for establishing an analytical solution. For each layer, boundary conditions are applied to the hollow cylinder, according to [START_REF] Hinton | A Study on the Effects of Coil Wedge During Rewinding of Thin Gauge Metals[END_REF]. The determination of T sup (z) and T in f (z) is done by a minimization procedure detailed in Section 7. The additional problem (5) way be rewritten as a function of displacements instead of the stress tensor. Then for each layer it reduces to the inhomogeneous Navier equation

µ∆ u + (λ + µ)∇div u = div λtr δε * I + 2µδε * (25)
such that the boundary conditions (8) are verified. The Navier equation ( 25) is linear. Thus, nonlinearity aspects are related not to the equation itself but rather to the computation of the eigenstrain increment δε * and the determination of contact stresses T sup (z) and T in f (z). Therefore, a general solution of (25) may be written as the linear combination of a particular solution u and a homogeneous solution u. Corresponding stress and stain tensors are denoted accordingly:

         u = u + u σ = σ + σ ε = ε + ε ( 26
)
The decomposition is summarized in Figure 5. It should be noted that u may be any particular solution of ( 27) without consideration of specific boundary conditions:

µ∆ u + (λ + µ)∇div u = div λtr δε * I + 2µδε * (27)
However, u is a solution of (28) verifying the boundary conditions (29), so the global boundary conditions (8) are verified:

µ∆ u + (λ + µ)∇div u = 0 (28) ∀ 1 ≤ i ≤ N,          σ(r sup , z).n(r sup , z) = -σ(r sup , z).n(r sup , z) + T sup (z) σ(r in f , z).n(r in f , z) = -σ(r in f , z).n(r in f , z) + T in f (z) σ(r, ±L).n(r, ±L) = 0 (29)
T sup (z) and T in f (z) are arbitrary in Sections 5 and 6. Correct values are determined by numerical minimization as detailed in Section 7.

In addition, the third boundary condition at z = ±L in (29) presents more difficulties with respect to the analytical treatment. For simplicity, this edge condition is verified not locally but in a weaker sense; that is, the resultant force of surface traction vanishes at the edges. Moreover, since the analytical development relies on a potential formulation that automatically verifies the general equilibrium, it is sufficient to impose this boundary condition at z = L, and the other condition at z = -L will be satisfied automatically. Therefore, this condition is replaced by 2π 0 r sup r in f σ(r, L).e z rdrdθ = 0 (30)

Particular solution in each layer

This section deals with the analytical determination of a particular solution u of the Navier equation ( 27), without specifying particular boundary conditions. One obtains after basic calculations the right-side term of [START_REF] Martinez | Jonction 16 MND 5-Inconel 690-316 LN par soudage-diffusion: élaboration et calcul des contraintes résiduelles de procédé[END_REF]:

div λtr δε * I + 2µδε * = f * r (r, z)e r + f * z (r, z)e z (31)
where

               f * r (r, z) = (λ + 2µ) ∂δε * rr ∂r + λ ∂δε * θθ ∂r + ∂δε * zz ∂r + 2µ δε * rr r - δε * θθ r + ∂δε * rz ∂z f * z (r, z) = (λ + 2µ) ∂δε * zz ∂z + λ ∂δε * rr ∂z + ∂δε * θθ ∂z + 2µ ∂δε * rz ∂r + δε * rz r (32)
To give an analytical form to the right-side term (32), a function sequence, orthogonal with respect to the following scalar product, is needed:

f, g = r sup r in f r f (r)g(r)dr ( 33
)
where f and g are functions mapping [r in f , r sup ] onto R. Consider the following function sequence:

G m : r ∈ r in f , r sup → J 1 x m r r sup Y 1 x m r in f r sup -J 1 x m r in f r sup Y 1 x m r r sup (34)
And consider F m is a primitive (up to a multiplicative constant) of G m :

F m : r ∈ r in f , r sup → J 0 x m r r sup Y 1 x m r in f r sup -J 1 x m r in f r sup Y 0 x m r r sup (35)
where J η and Y η are Bessel functions of order η ∈ R of the first and second kind, respectively. In addition, x m are successive positive roots of x → J

1 (x) Y 1 x r in f r sup -J 1 x r in f r sup Y 1 (x).
G m vanishes at r = r in f and r = r sup . The following orthogonality property is proven in Appendix A:

G m , G l = 0 if m l M m > 0 if l = m (36)
The proposed approach relies on Fourier series expansions of the right-side terms f * r and f * z , according to the axial direction z:

                   f * r (r, z) = K k=-K f r k (r) exp ikπ L z f * z (r, z) = K k=-K f z k (r) exp ikπ L z (37)
where

f r k (r) = 1 2L L -L f * r (r, z) exp -ikπ L z dz and f z k (r) = 1 2L L -L f * z (r, z) exp -ikπ L z dz (38)
The orthogonality property (36) is used to expand the functions f r k (r) and ∂ f z k (r)/∂r into a series based on the function sequence G m . However, functions G m vanish at r = r in f and r = r sup . Thus, the convergence can be difficult if f r k (r) and ∂ f z k (r)/∂r do not vanish at these points. It is therefore convenient to use the orthogonality property (36) not directly on f r k (r) and ∂ f z k (r)/∂r but on the following auxiliary functions, which vanish at r = r in f and r = r sup :

               g r k (r) = f r k (r) -A r k J 1 r r sup + B r k Y 1 r r sup g z k (r) = f z k (r) -A z k J 0 r r sup + B z k Y 0 r r sup (39) where A r k B r k =         J 1 r in f r sup Y 1 r in f r sup J 1 (1) Y 1 (1)         -1 . f r k (r in f ) f r k (r sup ) (40) and A z k B z k = -r sup         J 1 r in f r sup Y 1 r in f r sup J 1 (1) Y 1 (1)         -1 .        ∂ f z k ∂r (r in f ) ∂ f z k ∂r (r sup )        (41)
The following series expansion is based on an approach very similar to the classic Fourier-Bessel series expansion, detailed in [START_REF] Watson | A treatise on the theory of bessel functions[END_REF]. The Fourier-Bessel series expansion is based on the function sequence H m : x ∈ [0, a] → J η (y η m x), where a > 0 and η ∈ R and where y η m are the successive positive zeros of x → J η (ax). An orthogonality property similar to (36) also holds for the function sequence H m . In addition, all continuous functions on [0, a] can be written in the form +∞ m=1 f m H m (x) (where f m ∈ R). Completeness is proven in [START_REF] Watson | A treatise on the theory of bessel functions[END_REF] for an even more general function space. In this article, such a Fourier-Bessel series expansion is needed on the space of continuous functions on r in f , r sup instead of [0, a]. Thus, the only difference from the classic situation detailed in [START_REF] Watson | A treatise on the theory of bessel functions[END_REF] is the starting point of the definition set, which is not zero. To overcome this difficulty, the Bessel function of the second kind was introduced in the function sequence G m (34). The proof of completeness developed in [START_REF] Watson | A treatise on the theory of bessel functions[END_REF] could likely be adapted to the present situation, which is very similar. However, a simpler approach consists in projecting g r k (r) and ∂g z k (r)/∂r onto the function space of finite dimension, generated by (G m ) 1≤m≤M . Then using the orthogonality property (36) on the auxiliary functions (39), one obtains the projection

                   f * r (r, z) ≈ K k=-K        M m=1 f r m,k G m (r) + A r k J 1 r r sup + B r k Y 1 r r sup        exp ikπ L z f * z (r, z) ≈ K k=-K        M m=1 f z m,k F m (r) + A z k J 0 r r sup + B z k Y 0 r r sup + C z k        exp ikπ L z (42)
where

C z k = f z k (r sup ) - M m=1 f z m,k F m (r sup ) -A z k J 0 (1) -B z k Y 0 (1) (43)
and where the coefficients f r m,k and f z m,k are evaluated numerically as follows:

f r m,k = g r k (r), G m G m , G m and f z m,k = - r sup x m ∂g z k (r) ∂r , G m G m , G m (44)
Furthermore, the inhomogeneous Navier equation ( 27) reads

               (λ + 2µ) ∂ 2 u r ∂r 2 + 1 r ∂ u r ∂r - u r r 2 + µ ∂ 2 u r ∂z 2 + (λ + µ) ∂ 2 u z ∂r∂z = f * r (r, z) (λ + 2µ) ∂ 2 u z ∂z 2 + µ ∂ 2 u z ∂r 2 + 1 r ∂ u z ∂r + (λ + µ) ∂ 2 u r ∂r∂z + 1 r ∂ u r ∂z = f * z (r, z) (45)
A particular solution is sought as follows:

                   u r (r, z) = K k=-K        M m=1 u r m,k G m (r) + U r k J 1 r r sup + V r k Y 1 r r sup        exp ikπ L z u z (r, z) = K k=-K        M m=1 u z m,k F m (r) + U z k J 0 r r sup + V z k Y 0 r r sup + W z k (r)        exp ikπ L z (46)
where

W z k (r) =                  C z k r 2 4µ if k = 0 - C z k (kπ/L) 2 (λ + 2µ) if k 0 (47)
By plugging (46) into (45) and identifying, one obtains

U r k U z k = S k . A r k A z k and V r k V z k = S k . B r k B z k and u r m,k u z m,k = S m,k . f r m,k f z m,k (48)
where

S k =             -(λ + 2µ) 1 r sup 2 -µ kπ L 2 -(λ + µ) 1 r sup ikπ L (λ + µ) 1 r sup ikπ L -µ 1 r sup 2 -(λ + 2µ) kπ L 2             -1 (49)
and

S m,k =             -(λ + 2µ) x m r sup 2 -µ kπ L 2 -(λ + µ) x m r sup ikπ L (λ + µ) x m r sup ikπ L -µ x m r sup 2 -(λ + 2µ) kπ L 2             -1 (50)
A particular displacement field u r , u z , the solution of the inhomogeneous Navier equation, has been established. Therefore, the associated stress field σ can be computed as well, with use of the behavior relation in [START_REF] Moen | Prediction of residual stresses and strains in cold-formed steel members[END_REF]. More precisely, a Fourier series expansion of σ is obtained since the displacement field is known as a Fourier series expansion, so component-wise

σ (..) (r, z) = K k=-K σ (..) k (r) exp ikπ L z (51)
where the symbol (..) should be replaced by rr, rz, zz.

Homogeneous solution

As detailed in Section 4, the problem has been split into a particular solution of the inhomogeneous Navier equation (Section 5) and a homogeneous solution. This section is dedicated to the derivation of the homogeneous solution of the Navier equation ( 28) considering boundary conditions (29). Since the whole analytical solution derived in this article relies on Fourier series expansion along the axial direction, one should write additional contact pressure and shear stress as follows:

P in f (z) = K k=-K P in f k exp i kπ L z S in f (z) = K k=-K S in f k exp i kπ L z P sup (z) = K k=-K P sup k exp i kπ L z S sup (z) = K k=-K S sup k exp i kπ L z (52)
where

P in f k = 1 2L L -L P in f (z) exp -i kπ L z dz S in f k = 1 2L L -L S in f (z) exp -i kπ L z dz P sup k = 1 2L L -L P sup (z) exp -i kπ L z dz S sup k = 1 2L L -L S sup (z) exp -i kπ L z dz (53)
Boundary conditions (29) and ( 9) at r = r sup and r = r in f may be rewritten as

                   σ(r sup , z).n(r sup , z) = K k=-K -P sup k + σ rr k (r sup ) e r + S sup k -σ rz k (r sup ) e z exp i kπ L z σ(r in f , z).n(r in f , z) = K k=-K P in f k + σ rr k (r in f ) e r + S in f k + σ rz k (r in f ) e z exp i kπ L z (54)
In the following, the homogeneous solution u, σ verifying boundary conditions (54) and ( 30) is sought by use of biharmonic and harmonic potentials [START_REF] Barber | Elasticity[END_REF][START_REF] Weisz-Patrault | Three-dimensional elasticity based on quaternion-valued potentials[END_REF][START_REF] Papkovich | Solution générale des équations différentielles fondamentales d'élasticité, exprimée par trois fonctions harmoniques[END_REF]. To simplify the analytical procedure detailed in the following, the homogeneous solution u, σ is decomposed into three contributions:

u = u (a) + u (b) + u (c) σ = σ (a) + σ (b) + σ (c) (55)
The contribution with superscript (a) corresponds to purely biharmonic potentials and the contributions with superscripts (b) and (c) correspond to purely harmonic potentials.

Biharmonic potentials

A homogeneous solution u (a) , σ (a) of the Navier equation ( 28) verifying boundary conditions (54) but not necessarily (30) is sought. A classical formulation using purely biharmonic potentials reads [START_REF] Barber | Elasticity[END_REF]

             2µ u (a) r = ∂ 2 Φ ∂r∂z 2µ u (a) z = -2(1 -ν) ∂ 2 Φ ∂r 2 + 1 r ∂Φ ∂r + (1 -2ν) ∂ 2 Φ ∂z 2 (56)
where ν = λ/(2(λ + µ)) is the Poisson ratio and where

∆∆Φ = 0 (57)
The biharmonic potential Φ is expanded into a Fourier series and modified Bessel functions as done by Weisz-Patrault et al. [START_REF] Weisz-Patrault | Evaluation of contact stress during rolling process, by three dimensional analytical inverse method[END_REF]:

Φ(r, θ, z) = K k = -K k 0 kπ L -3         Φ k kπ L r I 0 kπ L r I 0 kπ L r in f + ϕ k I 0 kπ L r I 0 kπ L r in f + Φ k kπ L r K 0 kπ L r K 0 kπ L r in f + ϕ k K 0 kπ L r K 0 kπ L r in f         exp i kπ L z ( 58
)
where I 0 and K 0 are the zero-order modified Bessel functions of the first and second kind, respectively. Coefficients Φ k , ϕ k , Φ k , and ϕ k are to be determined with the boundary conditions (54). The associated stress field is given by

                                                   σ (a) rr = -ν∆ ∂ ∂z + ∂ 3 ∂r 2 ∂z Φ σ (a) rz = -(1 -ν) ∂ 3 ∂r 3 + 1 r ∂ 2 ∂r 2 - 1 r 2 ∂ ∂r + ν ∂ 3 ∂r∂z 2 Φ σ (a) θθ = -ν∆ ∂ ∂z + 1 r ∂ 2 ∂r∂z Φ σ (a) zz = -(2 -ν) ∂ ∂z ∆ - ∂ 3 ∂z 3 Φ (59)
By inserting (58) into (59), one obtains

                                                                                                               σ (a) rr (r, z) = K k = -K k 0 i          Φ k          (1 -2ν) I 0 kπr L I 0 kπrin f L + kπr L I 0 kπr L I 0 kπrin f L          + ϕ k          I 0 kπr L I 0 kπrin f L - I 0 kπr L kπr L I 0 kπrin f L          + Φ k          (1 -2ν) K 0 kπr L K 0 kπrin f L + kπr L K 0 kπr L K 0 kπrin f L          + ϕ k          K 0 kπr L K 0 kπrin f L - K 0 kπr L kπr L K 0 kπrin f L                   exp i kπ L z σ (a) rz (r, z) = K k = -K k 0          Φ k          - kπr L I 0 kπr L I 0 kπrin f L + 2(ν -1) I 0 kπr L I 0 kπrin f L          -ϕ k I 0 kπr L I 0 kπrin f L + Φ k          - kπr L K 0 kπr L K 0 kπrin f L + 2(ν -1) K 0 kπr L K 0 kπrin f L          -ϕ k K 0 kπr L K 0 kπrin f L          exp i kπ L z σ (a) θθ (r, z) = K k = -K k 0 i          Φ k (1 -2ν) I 0 kπr L I 0 kπrin f L + ϕ k I 0 kπr L kπr L I 0 kπrin f L + Φ k (1 -2ν) K 0 kπr L K 0 kπrin f L + ϕ k K 0 kπr L kπr L K 0 kπrin f L          exp i kπ L z σ (a) zz (r, z) = K k = -K k 0 i          Φ k          2(ν -2) I 0 kπr L I 0 kπrin f L - kπr L I 0 kπr L I 0 kπrin f L          -ϕ k I 0 kπr L I 0 kπrin f L + Φ k          2(ν -2) K 0 kπr L K 0 kπrin f L - kπr L K 0 kπr L K 0 kπrin f L          -ϕ k K 0 kπr L K 0 kπrin f L          exp i kπ L z (60)
Using the boundary conditions (54) at r = r in f and r = r sup , one obtains

               Φ k ϕ k Φ k ϕ k                = a -1 k .                   -P in f k + σ rr k (r in f ) -S in f k + σ rz k (r in f ) -P (i+1) k + σ rr k (r sup ) S (i+1) k -σ rz k (r sup )                   (61)
where a k is detailed coefficient-wise in Appendix B.

Harmonic potentials

Previous biharmonic potentials were expanded into a Fourier series (58), but k = 0 was excluded from the expansion. Thus, the term independent of z was not taken into account in the Fourier series expansion of the stress tensor σ (a) (60). Therefore, boundary conditions (54) cannot be completely verified with the contribution (a) only. In this section, a homogeneous solution u (b) , σ (b) of the Navier equation ( 28) is sought such that σ (b) rr and σ (b) rz do not depend on z. Purely harmonic potentials are used to this end (as detailed in [START_REF] Barber | Elasticity[END_REF]):

           2µ u (b) r = ∂φ ∂r + z ∂γ ∂r 2µ u (b) z = ∂φ ∂z + z ∂γ ∂z -(3 -4ν)γ (62) σ (b) zz = 2 A 0 + 6(1 + ν) 1 + 4ν B 0 z r in f (67)
Using the boundary conditions at r = r in f and r = r sup , one obtains

A 0 C 0 =          1 1 1 r in f r sup 2          -1 . P 0 + σ * ,rr 0 (r in f ) P (i+1) 0 + σ * ,rr 0 (r sup ) (68)
and

B 0 D 0 =       -6(1+ν) 1+4ν 2(1 -ν) -6(1+ν) 1+4ν r sup r in f 2(1 -ν) r in f r sup       -1 .       -S 0 + σ * ,rz 0 (r in f ) S (i+1) 0 -σ * ,rz 0 (r sup )       (69)

Edge conditions

The combination of contributions (a) and (b) satisfies boundary conditions (54) but not boundary conditions (30). In this section, a third (and final) solution u (c) , σ (c) of the homogeneous solution of the Navier equation ( 28) is sought such that (30) is verified. Since the problem is axisymmetric, the boundary condition (30) reduces to

r sup r in f σ (c) zz (r, L)rdr = - r sup r in f σ zz + σ (a) zz + σ (b) zz (r, L)rdr (70)
A formulation using purely harmonic potentials is used as in Section 6.2:

             φ = E 0 ν 1 + ν - r 2 2 + z 2 γ = - zE 0 2(1 + ν) (71)
The only nonvanishing components of the associated displacements and stresses read

                 2µ u (c) r = - E 0 ν 1 + ν r 2µ u (c) z = E 0 1 + ν z σ (c) zz = E 0 (72)
Therefore, by use of (70), the last coefficient is obtained:

E 0 = - 2 r sup 2 -r in f 2 r sup r in f σ zz + σ (a) zz + σ (b) zz (r, L)rdr (73)

Numerical determination of contact stresses

In previous sections, additional contact pressure P sup , P in f and shear S sup , S in f were arbitrary. In this section, the numerical determination of these contact stresses is detailed. It relies on the same principle as in the work of Weisz-Patrault et al. [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF]. That is, contact stresses are evaluated by a minimization procedure using explicit contact laws defined in [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF]. It consists in relating contact stresses to the position discontinuity between two successive layers. This discontinuity of nominal surfaces at each interface is constituted of the interpenetration r C along the radial direction (associated with contact pressures) and the slide z C along the axial direction (associated with shear stresses), where the symbol denotes the jump at the interface considered. Although contact shear is negligible in the numerical study detailed in this contribution, the determination of contact conditions is presented in general. Consider the following contact laws:

P C = f r C and S C = g z C (74)
where f and g are simple functions defining the contact laws. In addition, P C and S C are the contact pressure and the shear stress at the interface considered. The coil consists of many layers, so contact pressure and shear at each interface are gathered in two vectors P and S as well as surface interpenetration r and slides z .

Contact stresses are z-dependent and should be evaluated for each time step, which leads to the solving of many optimization problems. Therefore, the analytical solution proposed for solving the mechanical problem enables us to obtain reasonable computation times.

Validation of the analytical solution

In this section, a numerical validation of the analytical solution developed in Sections 5 and 6 is proposed. Only one layer, subjected to an arbitrary eigenstrain, is modeled for this validation. The combination of several layers and the determination of contact stresses is not included in this section. An axisymmetric tube representing one layer of the coil was simulated by use of the finite element software program Cast3M developed by CEA [START_REF]Cast3m[END_REF]. The geometrical conditions, listed in Table 1, are chosen to be representative with respect to the application. Tractionfree boundary conditions are set at r = r in f and r = r sup , and the only load consists in imposing an eigenstrain ε * as listed in Table 2. These loading conditions necessitate very fine meshes to obtain reliable numerical results, which are presented in Figure 6. The mesh density is defined by the number of elements along the radial direction N r = 100 and along the axial direction N z = 1000. The rather slow convergence of the finite element simulation toward the analytical solution, as a function of the mesh size, is illustrated in Figure 7 (where N r is fixed to 100). The comparison between the analytical solution and the finite element computation (at different axial positions) is presented in Figures 8 and9. Good agreement is observed even though slight discrepancies are noticeable. These discrepancies are due to the lack of accuracy of the finite element computation. It can be observed that the boundary conditions are not verified perfectly: σ rr does not vanish at r = r in f and r sup , and σ rz does not vanish at z = 0 (as it should because of the symmetry). In addition, the traction-free boundary condition at the edges z = ±L is not verified locally in the analytical solution. A weak condition is verified instead; that is, only the resultant forces vanish at the edges. The region affected by the local traction-free boundary condition spreads over a distance roughly equal to three times the thickness. Thus, this edge effect is negligible for applications where layers are very long compared with their thickness.

Numerical results

In this section, numerical results are obtained by use of models developed by Weisz-Patrault et al. [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF] and Weisz-Patrault [START_REF] Weisz-Patrault | Coupled heat conduction and multiphase change problem accounting for thermal contact resistance[END_REF], on the one hand, and in the present contribution, on the other hand. The tested coiling parameters represent classic hot rolling conditions (e.g., for dual-phase steel grades). All parameters used for the simulation are listed in Table 3. Moreover, the Young's modulus as a function of temperature is extracted from [START_REF] Martinez | Jonction 16 MND 5-Inconel 690-316 LN par soudage-diffusion: élaboration et calcul des contraintes résiduelles de procédé[END_REF]:

E(T ) = 2.08 × 10 5 -1.90 × 10 2 T + 1.19T 2 -2.82 × 10 -3 T 3 + 1.66 × 10 -6 T 4 ( 76
)
where T is given in degrees Celsius. The yield stress of each phase as a function of temperature is needed. Grostabussiat-Petit [START_REF] Grostabussiat-Petit | Conséquence mécanique des transformations structurales dans les alliages ferreux[END_REF] performed uniaxial tensile tests at different temperatures and identified the yield stress of austenite and bainite between 400 ˚C and 620 ˚C and between 500 ˚C and 544 ˚C, respectively. Other ferritic phases (ferrite and pearlite) are assumed to have the same yield stress as bainite in this example, even though proper yield stresses for all phases could have been used instead. The martensitic yield stress has been set to 800 MPa. Grostabussiat-Petit [START_REF] Grostabussiat-Petit | Conséquence mécanique des transformations structurales dans les alliages ferreux[END_REF] proposed the following linear interpolation:

σ Y 0 (T ) = -0.147 × T + 183.14 σ Y 3 (T ) = -0.02 × T + 444.8 ( 77
)
where T is given in degrees Celsius. The winding phase of the coiling process is simulated according to [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF], and then the cooling phase is simulated according to [START_REF] Weisz-Patrault | Coupled heat conduction and multiphase change problem accounting for thermal contact resistance[END_REF]. The temperature field obtained is presented at different time steps in Figure 10. Considering the initial temperature, ferrite cannot be produced. Pearlite, bainite, and martensite phase proportion fields are given at different time steps in Figures 11,12, and 13, respectively.

When the coil is wound, the inner mandrel is removed to cool the coil more effectively. Contact pressures before and after mandrel removal are presented in Figure 14. It is clear that before removal of the mandrel, contact tends to concentrate more and more at the center when the number of layers increases. This is classically due to the geometrical profile of the incoming strip, which is assumed to be quadratic with respect to the axial direction. After removal of the mandrel, the contact pressure distribution is completely modified and tends to decrease significantly in amplitude and to concentrate at the center in the whole coil (and not only for the last layers). That is, removal of the internal support of the coil leads to loose contact, especially at the edges, where contact pressures were lower than at the center.

Thermal contact resistance is maximal when contact is lost; therefore, the coil cools more rapidly at the center, where a good contact is ensured, than at the edges, where the distance between two successive layers acts as a significant insulation. However, the first and last layers cool more rapidly at the edges because the air gaps tend to insulate them from the rest of the coil. In other words, the thermal capacity is much lower in these regions than at the center. The proportions of the different phases highly depend on temperature kinetics, and thus the global structures of the phase proportion fields are not surprisingly similar to the temperature variations.

The present contribution enables us to compute displacements and stresses when the coil cools down, considering, on the one hand, the initial state obtained at the end of the winding phase and after removal of the mandrel [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF] and, on the other hand, all nonlinear contributions gathered as a total eigenstrain. When the coil cools down and undergoes phase transitions, contact pressures rapidly decrease, then localize and decrease again, as shown in Figure 15. These results depend on the chosen contact law that is relatively soft in this example. Instabilities are expected when contact pressures vanish all along an interface. In that situation adhesion between the corresponding layers is lost, and one should observe unwinding. It is clear that this issue is critical since contact pressure significantly decreases during cooling of the coil. It can be seen from Figure 15 that the proposed coiling conditions lead to this issue.

Another relevant result is the cumulative inelastic strain exhibited along the circumferential direction in Figure [START_REF] Sheu | Mixed lubrication of strip rolling[END_REF]. Inelastic strain determines classic flatness defects obtained after uncoiling. For instance, one can use the simple purely elastic uncoiling model proposed by Weisz-Patrault et al. [START_REF] Weisz-Patrault | Non-linear simulation of coiling accounting for roughness of contacts and multiplicative elastic-plastic behavior[END_REF] to compute residual stresses on the basis of a flat reference. One of the most interesting ways to present the residual stress profile is in terms of residual curvature. In Figure 17, residual curvatures along the rolling direction denoted by X and along the axial direction denoted by Z are presented for different layers. Residual curvatures are presented before cooling (just after the elastic-plastic winding) and at the end of the cooling process to give insight into the extent to which thermal expansion, phase transitions, and transformation-induced plasticity contribute to the final residual stress. Significant residual stresses are produced during the coiling process and are dramatically increased and redistributed during the cooling phase, as shown in Figure 17.

and the martensitic phase transition is described by X 5 (t) = X ra (t) 1exp (α MS (T (t) -MS ))

where X ra (t) is the retained austenite.

Conclusion

A mixed analytical/numerical strategy has been developed to compute efficiently residual stresses induced by plastic deformation, phase transitions, and temperature variations in a coil. The analytical part of the proposed solution was verified by comparison with a finite element computation. The combination of previous models and the present solution enabled us to simulate the whole coiling process, including both winding and cooling of the coil. Parametric studies can be performed to optimize the process with respect to residual stress issues. In addition, this work also provides a tool for estimating inelastic strain during cooling, which is a significant input for computing coil sagging. Moreover, the example presented in this article shows that contact pressure tends to decrease when the coil cools down, which can lead to the unwinding of the coil if contact is lost. In addition, significant flatness defects are obtained after uncoiling.

k = i        (1 -2ν) K 0 kπrsup L K 0 kπr in f L + kπr sup L K 0 kπrsup L K 0 kπr in f L        a (3,4) k = i        K 0 kπrsup L K 0 kπr in f L - K 0 kπrsup L kπrsup L K 0 kπr in f L        a (

Figure 1 :

 1 Figure 1: Forming processes

Figure 2 :

 2 Figure 2: Instabilities

Figure 3 :

 3 Figure 3: General modeling structure

Figure 4 :

 4 Figure 4: Global system and notation

Figure 5 :

 5 Figure 5: Decomposition in each layer

  Contact stresses are obtained by setting trial values P and S and computing the corresponding interpenetration r and slides z with use of the mechanical model detailed in previous sections. Then the contact law (74) is used for each interface and gives another evaluation of contact stresses denoted by P and S. Thus, contact stresses are determined by the solving of the following optimization problem: P = argmin P≥0 P -P and S = argmin S S -S

  ) = (r -R sup )(r -R in f )/(R sup R in f ).

Figure 6 :Figure 7 :

 67 Figure 6: Finite element computation (condition 1)

Figure 8 :

 8 Figure 8: Validation by comparison with a finite element computation (condition 1). FEM, finite element method.

Figure 9 :

 9 Figure 9: Validation by comparison with a Finite Element computation (condition 2). FEM, finite element method.

Figure 11 :Figure 12 :Figure 13 :Figure 14 :Figure 16 :

 1112131416 Figure 11: Pearlite proportion evolution

Figure 17 :

 17 Figure 17: Residual curvature before and after cooling

  The author gratefully acknowledge Nicolas Legrand (ArcelorMittal Global Research & Development, East Chicago 3001 East Columbus Drive, East Chicago, IN 46312, USA) for fruitfully discussions.

Table 1 :

 1 Modeling parameters

	External radius R sup (mm)	303
	Internal radius R in f (mm)	300
	Half length	L	(mm) 50 (condition 1)
				200 (condition 2)

Table 2 :

 2 Loading conditions

Table 3 :

 3 Coiling conditions

	External mandrel's radius	R ext 381	(mm)
	Internal mandrel's radius	R int 0	(mm)
	Initial temperature	T ini 750	(K)
	External temperature	T ext 300	(K)
	Heat transfer coefficient	H	0.4 × 10 -3 (W mm -2 K -1 )
	Applied nominal tension	σ a	10	(MPa)
	Initial austenite proportion Initial ferrite proportion Initial pearlite proportion Initial bainite proportion Initial martensite proportion	X ini 1 X ini 2 X ini 3 X ini 4 X ini 5	0.4 0.3 0.3 0.0 0.0
	Half length	L	467	(mm)
	Strip thickness at the center	t c	3.3	(mm)
	Strip thickness at the edges	t e	3.261	(mm)
	Number of cycles	N rev 196
	Enthalpy variation for p = 2, 3, 4	∆H p -0.4	(J mm -3 )
	Enthalpy variation for p = 5	∆H 5 -0.3	(J mm -3 )
	Ferritic transition temperature	AE 3 1073	(K)
	Pearlitic transition temperature	AE 1 823	(K)
	Bainitic transition temperature	BS	713	(K)
	Martensitic transition temperature	MS 693	(K)
	Avrami exponents for all phase transitions	n p	1.5
	Avrami coefficients for ferritic transformation	k 2	0.00008
	Avrami coefficients for pearlitic transformation	k 3	0.00015
	Avrami coefficients for bainitic transformation	k 4	0.00015
	Martensitic coefficient	α MS 0.045
	Ferritic, pearlitic, and bainitic phase transitions are described by simple Avrami equations,
	∀p ∈ {2, 3, 4} , X

p (t) = 1exp -k p t n p ,

Figure 15: Contact pressure

where ∆φ = ∆γ = 0 (63)

The associated stress field is given by

Harmonic polynomials according to r and z are sought for φ and γ such that σ (b) rr and σ (b) rz do not depend on z. Thus,

Therefore, by plugging (65) into ( 62) and (64), one obtains The proof of orthogonality relies on several classic indefinite integrals (proven in [START_REF] Watson | A treatise on the theory of bessel functions[END_REF]):

a (1,1) k