Buckling of rolled thin sheets under residual stresses by ANM and Arlequin method

We present a numerical technique to model the buckling of a rolled thin sheet. It consists in coupling, within the Arlequin framework, a three dimensional model based on 8-nodes tri-linear hexahedron, used in the sheet part located upstream the roll bite, and a well-suited finite element shell model, in the roll bite downstream sheet part, in order to cope with buckling phenomena. The resulting nonlinear problem is solved by the Asymptotic Numerical Method (ANM) that is efficient to capture buckling instabilities. The originalities of the paper ly, first in an Arlequin procedure with moving meshes, second in an efficient application to a thin sheet rolling process. The suggested algorithm is applied to very thin sheet rolling scenarios involving 1 Laboratoire d’Etude des Microstructures et de Mécanique des Matériaux, LEM3, UMR CNRS 7239, Université de Lorraine, Ile du Saulcy, Cedex 01, 57045 Metz, France 2 Laboratory of Excellence on Design of Alloy Metals for low-mAss Structures (DAMAS), Université de Lorraine, Metz, France 3 Laboratoire de Mécanique des Sols, Structures et Matériaux (MSSMat UMR CNRS 8579), Ecole CentraleSupélec, CNRS, Université Paris-Saclay, Grande Voie des Vignes, Cedex, 92295 Châtenay-Malabry, France “edges-waves” and “center-waves” defects. The obtained results show the effectiveness of our global approach.


Introduction
Rolling of thin sheets generally induces flatness defects due to the thin aspect of the sheet and to thermo-elastic deformation of rolls whose profile in the roll-bite generally does not match perfectly the strip thickness profile. This leads to heterogeneous plastic deformations throughout the strip width and then to out of mid-plane displacements that relax compressive residual stresses (see Fig. 1). The most important flatness defects are "edge-waves" and "center-waves" buckles. These waves are the result of buckling due to self-equilibrating longitudinal residual stresses with a compressive longitudinal membrane stress state in the middle of the strip (center-waves) or in the edge zones (edge-waves), respectively. Theoretically a large strain 3D elastoplastic model must be considered to simulate rolling process. In the literature, there are several codes which can simulate rolling processes. In this work, we focus on the rolling code LAM3 [1][2][3]. This code is involved to compute the behavior of the sheet by three dimensional finite element method (3D FEM), and the roll stack elastic deformation by semianalytical model. LAM3 code gives satisfactory results for cases whose flatness defects do not occur. However, it is unable to model manifested flatness defects, as shown in Fig. 2: the code overestimates the stress field beyond the roll bite. This is mainly due to the buckling phenomenon in thin sheets that is disregarded by LAM3 [2]. Often the width to thickness ratio is of the order of 10 4 . To simulate buckling and post-buckling during rolling, existing models generally treat in a weakly coupled manner the rolling problem and the buckling phenomenon: a profile of the residual stresses is recovered from a first rolling computation step before being injected into a thin sheet buckling analysis procedure, in a second step. Fischer et al. [5][6][7] performed a semi-analytical approach of buckling under residual stresses. The study mainly concerns the state of the plate beyond the critical threshold of residual stresses that can trigger the buckling phenomenon in the sheet. They consider self-equilibrating stress states resulting from a rolling process with non-uniform distribution of the longitudinal residual stress over the width of the strip. One can cite also reduced models in [8,9] where amplitude and wavelength of defects are determined as functions of the residual stresses. In the recent paper [4], Abdelkhalek proposed a full model based on a shell formulation and on an asymptotic numerical method (ANM). In a number of situations, the developed code computes correctly the buckling phenomena under residual stresses, by detecting with significant precision the buckling modes and the amplitudes of defects.
As a matter of fact, the weakly coupled approach is valid only in cases where buckling has a weak feedback on the sheet part located in the roll bite. As reported in [4], for  [1,4] rolling cases with low thickness reduction, for instance, the feedback can no more be ignored. To our best knowledge, up to now strong coupling has been addressed in two approximate ways. The first way has been initiated by Counhaye [10]. A rather "old" idea [11] was used to take into account the influence of wrinkling on membrane behavior: roughly, the membrane constitutive law is modified to induce an upper limit to the compressive stresses due to buckling. Though offering a way of relaxation of compressive stresses, this simple approach suffers drawbacks. First it does not provide the shape of the buckled sheet, an important parameter from an industrial point of view. Second, the convergence of the modified algorithm is not easy to achieve, leading to an inaccurate reduction of compressive stresses due to buckling [1,4]. Abdelkhalek [4,12] proposed another approach based on an iterative coupling between a rolling code and a shell model. It alternates 3D computations to find a first evaluation of the residual stresses generated by rolling and shell computations to relax these residual stresses, downstream the roll bite part of the sheet. It has permitted to find the shape and the distribution of flatness defects as well as the final residual stresses. Nevertheless this procedure is intricate, the considered computation domains are changed at each step, which implies frequent transfers of data. More importantly, an unjustified kinematical boundary condition is imposed on the interface between up and downstream roll bide domains. In this paper, we develop a concurrent coupling between up and the downstream domains with changing sizes, to simulate the buckling phenomenon under residual stresses. A 3D model is used in the first evolving domain, whilst a nonlinear shell model is used in the second one, overlapping the first. The effect of rolling is introduced by the residual stress field computed in beforehand by LAM3 code. The 3D and Shell models are coupled within the Arlequin framework. Other coupling techniques could be used. However, the Arlequin framework, initiated by Ben Dhia [13] and further investigated and developed by Ben Dhia and his collaborators and many other researchers, has proved to be a very flexible tool to handle, not only multi-model (e.g. [14]), but also multi-scale problems (e.g. [15]). Arlequin method relies on a partition of models, each model being valid in the domain on which it is defined. Furthermore, Arlequin method introduces coupling operators and partitions of energy between models. Generally, there is a hierarchy between the two models, one model being a simplified version of the other. As compared with the classical method proposed by Ben Dhia and Rateau [15], we introduce a relevant coupling operator and a simple procedure to build varying meshes (see Fig. 3). Thus, the 3D model and the shell model are applied in sub-domains of variable sizes. By performing this coupling, we consider the same spatial domain in LAM3 code and in the present buckling code. Thus, we avoid boundary conditions problems that are necessary within Abdelkhalek's model described above. The paper is organized as follows. In section "A new numerical model for sheet buckling", we describe the proposed model for sheet buckling under residual stresses. Details are given for 3D and shell models and for Arlequin coupling procedure. We present also in this section resolution and discretization techniques used to compute the global solution of the buckling problem. In section "Numerical applica-tions", numerical results are given to assess the accuracy and the robustness of the coupled 3D/shell model.

A new numerical model for sheet buckling
A new numerical model is presented in order to predict the flatness defects generated by a rolling process. It runs as a complement of a 3D finite element code devoted to rolling. One assumes that the rolling code is able to compute the plastic strain induced by the plastic behavior in the bite, or equivalently the residual stresses. Thus the whole process setting can be taken into account, i.e. rolling forces and friction, rolls bending... The stress calculated in this way will be the data of the new code. The physical domain will be the same in the rolling code and in the new buckling code, which avoids questionable boundary conditions. The rolls will be not represented in this second part and their effect will accounted only via the residual stresses resulting from the rolling code. In the new code, the downstream domain will be discretized by 3D elements and the upstream domain by nonlinear shell elements, both models being coupled by Arlequin method. So the two parts of the domain will be combined in the same finite element simulation.

Preliminary computation (LAM3 code)
LAM3 is a finite element code devoted to the modeling of rolling process. The discretization is performed by classical 8-nodes hexaedra. A first key point is the semi-analytic treatment of the deformation of the rolls, what should require tens of thousands of unknowns by brute force finite elements [16]. Indeed this roll deformation has a crucial role in the appearance of a residual stress field that is heterogeneous in the sheet width.A second original point of LAM3 is the possibility to achieve stationary computations. Because of these two features, the computation time is much smaller than with generic codes for forming processes. This efficiency and many years of industrial use make it very suitable for our preliminary computations. Otherwise LAM3 looks like other commercial codes, including classical elasto-visco-plastic materials or constitutive laws defined by the user, unilateral contact and Coulomb Tresca Anisotropic (CTA) friction law. As explained previously, LAM3 does not permit to predict the buckled shapes of very thin metal sheets. To correct this weak point, we try a non intrusive procedure associating LAM3 to compute the plastic strains generated by the process with a new code to deduce the flatness defects from these residual stresses.

3D model
The 3D domain is represented by a solid subjected to residual stresses distribution σ res , zero in the upstream area and non-zero after the roll bite. In this way, the flatness defects do not manifest downstream. Thus, we consider small strain in this domain and then a linearized kinematic form for the strain tensor ε expressed in terms of u 3D as: ε = 1 2 (∇u 3D + t ∇u 3D ). The residual stresses which come from a full model taking into account the elastoplastic law adapted to rolling (or an analytical form) is introduced in the 3D continuum model. The weak form of equilibrium equation, the constitutive law and the compatibility relation can be written as follows: where σ denotes the Cauchy stress tensor, δε is the virtual strain tensor and C is the elasticity tensor. In the constitutive law, σ res denotes residual stresses derived from LAM3 code and λ res a scalar parameter (0 λ res 1). Note that, in this formulation, we do not apply external forces to the 3D model. In what follows, forces induced by the shell part will be added by the bridging technique.

Shell model
In the literature, many nonlinear shell models are available. In this work, we have chosen a shell model based on a three dimensional formulation proposed by Buchter et al. [17] efficient and easy to implement. It is well adapted to problems involving large displacements and large rotations in a total lagrangian framework. Position of a considered point in the current configuration is given by the displacement of the mid-surface and the variation of the director vector between the reference and the current configurations. The displacement can be expressed as follows: where (θ 1 , θ 2 , θ 3 ) are the local curvilinear coordinates.
Vectors v and w denote respectively the displacement of the mid-surface and updating director vector. In the total lagrangian formulation, it is convenient to express the energy functional in terms of the second tensor of Piola Kirchhoff S and the Green Lagrange strain tensor γ . Furthermore, in the shell model, a linear variation of the strain across the thickness is added via the concept EAS of Simo and Rifai [18]. The strain field is then decomposed in a compatible part γ c and an additional one γ so that γ = γ c + γ . This extra variable, orthogonal to the stress field, allows us to enrich the shell formulation and to use a full 3D constitutive relation without condensation which is not the case in the classical shell formulations. For more details about the proposed shell model, refer to references [17,19,20].
Taking into account residual stresses in the same manner as for the 3D model, the shell formulation can be summarized as follows:

Arlequin coupling
The Arlequin framework consists in partitioning a mechanical system, initialy represented by a domain , into two overlapping subdomains 1 and 2 (see Fig. 4) where the intersecting domain is called the gluing zone, c [14,15]. Thus, Arlequin method allows one to couple two different mechanical states through reliable coupling operators as well as consistent energy distribution between the two subdomains. The central point of the model is the coupling operator C which is selected by analogy with the deformation energy of the shell. In this paper, we limit ourselves to the H 1 coupling, and hence where u and μ are respectively the displacement and the Lagrange multiplier field and κ denotes a parameter related to a length. The variable x 0 is a parameter which varies along the sheet length and allows locating the roll bite. Thus, the gluing zone c varies during the rolling process. In addition we introduce the weighting functions α i in order to share the energy between the two models in the coupling area. These weighting functions are used respectively for strain energy and for work induced by external forces. They satisfy the following relations: Find the global displacement field in the domain is therefore to weight the displacement u 3D and u shell respectively in the domains 1 and 2 . In the following, let us denote We recall that we do not apply external forces to the 3D model. They are only applied to the shell model (domain 2). By introducing the weighting functions α i in Eqs. 1 and 3, the Arlequin method allows us to write the following variational form: To solve system (6), we use asymptotic numerical method which is described in the section that follows.

Resolution algorithm
System (6) is nonlinear. Indeed, if we consider that residual stresses are given for our model, the only nonlinear terms come from Eq. 3 4 . This nonlinear system can be solved by classical iterative technique but in this work we propose to solve it by using the asymptotic numerical method (ANM) which is a useful tool for nonlinear problems involving structural instability [21][22][23][24]. ANM consists in expanding variables of the considered problem into power series truncated at a high order. This allows one to obtain large part of the solution branch by decomposing only one tangent operator leading then to a significant reduction of the computation time. Several applications in fluid and solid mechanics have shown the efficiency and the robustness of ANM [19,[21][22][23][24]. Let us consider a mixed vector holding all the variables of the nonlinear problem (6): The unknown field U is expanded into power series with respect to a path parameter 'a': The parameter N represents the truncation order of the series usually chosen between 10 and 20 for an optimal computation. By introducing the series (8) into the nonlinear problem (6), one obtains a sequence of linear problems admitting the same tangent operator. See Appendix A for more details about this procedure. In this manner a large part of the solution branch is obtained.
After obtaining the linear problems, the displacement field and the Lagrange multiplier are discretized as follows. In the domain 1 , the displacement u 1 is discretized by using the classical 8-node tri-linear hexahedron elements and in the shell domain 2 , the displacement u 2 is discretized by using isoparametric quadrilateral element with eight nodes. Refer to reference [19] for more details of the shell discretization and to reference [14] for Arlequin discretization framework.
Note that the Lagrange multiplier is discretized by using the same interpolation of the shell element. Finally the discretized linear problems can be set in the following form for order k: As the series has an intrinsic convergence radius, a continuation technique is used. It consists in computing a maximum value of the parameter 'a' by considering the difference between the solutions for two consecutive orders must be smaller than a given parameter δ proposed by the user. This technique is very simple to be implemented allowing to obtain a naturally adaptive step length which varies within the local nonlinearity of the response curve [25][26][27]. Note that ANM algorithm does not need any iteration procedure as it is the case for Newton Raphson technique.
The resolution is performed using two steps. Indeed, the system of Eq. 9 has two load parameters λ res and λ. Both parameters do not vary at the same time, i.e. only one load parameter varies during a step and the second will be kept constant. In the first step, λ res is constant and therefore only λ will be developed in the form of power series. During this step, we apply traction T 0 uniformly distributed over an edge of the sheet in the downstream domain. This allows us to put the sheet under traction and to be closer to the actual rolling conditions.
In the second step, a study of post-buckling of the sheet under residual stresses is performed. We start from the Fig. 4 The Arlequin method in a general mechanical problem. The domain (a) is split into two superposed domains 1 and 2 that intersect in the gluing solution of the first step to calculate the evolution of the post-buckling of the sheet with respect to the load parameter λ res . Thus, during this step the load parameter λ is set equal to one. We emphasize that the method of resolution of the two steps is similar to a classical resolution of systems of equations using asymptotic numerical method [19,26].

Mesh management
The main objective of this work is to develop a simplified model of rolling, adapted to actual rolling conditions. As indicated previously, we adopt a coupling procedure which consists in moving the gluing zone c and the residual stresses at every step of simulation. In this procedure, we have to couple fields defined on incompatible meshes. We proceed by a matching technique detailed in [28]. The geometry being very simple, it is not difficult to establish a procedure to define moving meshes for the three parts of the problem: shell part, 3D part and coupling part (see Fig. 5). We start from two underlying meshes: a surface mesh for the shell part and a volume mesh for the 3D part. The Arlequin method permits to couple arbitrary meshes, but for simplicity we limit ourselves to two compatible starting meshes, which means that the projection of 3D nodes are corners of the shell elements. The gluing zone is a rectangle of constant size that can undergo a translation in the Ox direction. Thus the location of the gluing zone depends on a single parameter x 0 . When the position the gluing zone is known, its mesh is defined by a coarsening of the shell mesh: indeed the Lagrange multiplier field has to be discretized from a coarser mesh to avoid a "locking" phenomenon [13,14,29] while a fine mesh is necessary to accurately describe buckling patterns with small wavelength. So we obtain then three zones (zone 1, zone 2 and zone 3) on the two subdomains 1 and 2 . The degrees of freedom of zone 3 in the 3D part are disabled and similarly, the degrees of freedom of zone 1 in the shell model. In this way we do not need boundary conditions at the end of the 3D part and at the beginning of the shell part: these two parts will be coupled only by the Arlequin matching.

Numerical applications
The coupling procedure presented previously will be applied to the numerical simulation of flatness defects generated by a rolling process, i.e. of deformed shapes of the sheet due to residual stress field. We are interested in very thin sheets where flatness defects are of great interest. This procedure is non-iterative: firstly a rolling calculation is performed using LAM3 code; secondly, the residual stress field resulting from this calculation is used as a datum in the present code 3D/shell/Arlequin. The results are compared with Counhaye's [10] and Abdelkhalek's [4] models and with experimental measurements. Counhaye's model is similar as Roddeman's approach to account for the influence of wrinkling on membrane behaviour [30] in which the constitutive law is modified to limit the in-plane compressive stress. The implementation in a forming code is rather easy, this has been done in LAM3 [31] and we shall use this version of LAM3. This permits to relax the excessive compressive stress of Fig. 2, but the convergence can be difficult, which leads to large computation time and to a procedure that is not very robust. In addition, this model cannot compute the deformed shape of the sheet. On the contrary, Abdelkhalek [1,4] uses a shell finite element model so that any deformed shape can be obtained from any residual stress field. The procedure proposed by Abdelkhalek is an iterative coupling between the rolling code LAM3 and a buckling code with shell element. The latter procedure permits also to relax the excessive stresses and moreover to predict any deformed shapes. Nevertheless, when compared to the present procedure, it is very complicated: the studied domain is not the same at the first, second and third iteration and this requires data transfer and questionable boundary conditions.
Two numerical tests will be studied. In the first one, the residual stress field is introduced in analytical form. In the second case, this residual stress field is deduced from residual stresses generated by rolling simulation obtained by LAM3 code, in a configuration where experimental results are available.

Edges-waves flatness defects generated by an analytical profile of residual stress
In this section, the proposed coupling model is applied in a case where the residual stress field is defined a priori. Its shape is designed to generate edges-waves flatness defects. As in actual rolling processes, two types of loading will be applied to the sheet: a self-equilibrating residual stress S res (x, y) and a force that is uniformly distributed on the edge 4, f = T 0 S surf , where S surf is the cross section area. Equivalently, one can say that a uniform tensile stress T 0 is applied on the edge 4 (Fig. 7a). Through this application, we seek to validate the model by comparing it with a pure shell model and to elucidate various techniques proposed in the present code. Note that a pure shell model is simpler than the coupled 3D/shell model discussed in this paper and these two models are relevant in the present algorithm because the 3D field in the bite are not considered directly in the second step. Nevertheless the coupled model will be necessary in the future when the rolling code will be applied concurrently  Poisson's ratio ν = 0.3. In this case, the residual stress field is in the following analytical form: where e is the exponential function. The Cumulative Distribution Function CDM(x, x 0 ), centered at x 0 progressively distributes the residual stresses in the sheet. The length κ characterizes this transition. In the calculation, we choose κ = 1 mm so that the residual stress growths monotonically on a length of about 50 mm. The flatness defect in the sheet appears in the compressive area. In the present case, the maximum compressive stress at the edges is min = −240 MPa and the maximum tensile stress is max = 60 MPa (see Fig. 6).  To avoid rigid-body displacement of the strip, the face 3 of the three-dimensional continuum model is clamped (see Fig. 7). Moreover, to take the influence of the rolls into account, we fix the vertical displacement in the domain 1 (Fig. 7a) and the corresponding zone of the shell model (Fig. 7b).  Fig. 9). This maximum is here located on the edges. So the bifurcation point, the buckling load and the post-buckling behavior are similar within the two approaches which assesses the validity of the coupling technique. Hence the 3D/shell model gives about the same results as a pure shell model. This comparison has been achieved with a specific shell finite element that is known to be equivalent to well established finite elements like the S8R5 element (Shell, 8-node, Reduced Integration, 5 degrees of freedom by node) of the Abaqus code [33]. We underline that there are no papers in the open literature about post-buckling behavior of rolled thin sheets, except in [8,34].

Results and discussions
Traction plays an important role in rolling process. During the rolling, the sheet is often subjected to a tensile stress T 0 able to hide all or part of buckling phenomena. From the work of Fischer et al. [6], an increasing global traction T 0 does not only lead to increase the value of the critical residual stress, but it also produces shorter buckling waves concentrated towards the edge of the strip. This phenomenon can be simply explained: with a large value of T 0 , the buckling will occur when the compressive zone is located on a small region near the edge. Thus a look on plate buckling equations shows a monotonical relation between longitudinal and transverse characteristic lengths. Hence a larger traction implies more edge-localized modes and therefore to shorter longitudinal wavelengths. This behavior has been checked with the coupled numerical model for three values of the tensile stress: T 0 = 10MP a, T 0 = 30MP a and T 0 = 50MP a (see Fig. 10). The number of periods passes from 4 to 5 and 6, which confirms the decrease of the wavelength with an increase of T 0 (see Fig. 10). In parallel, the size of the defect decreases from 4 mm to 3 mm and 2.5 mm, which follows from the increase of the critical value of the bifurcation load λ res (cf. Fig. 11). Hence one recovers the conclusions of Fischer et al. [6,7,35] that were obtained analytically.

Local folds at bite exit
The numerical application proposed in this section allows us to discuss the influence of the transverse component of the residual stress tensor S res yy . Most of the papers [6,32] focus on the longitudinal component S res xx (y) that is relevant for buckling modes occurring away from the rolls. Here we choose a profile of residual stress S res yy = − 4 B 2 y (B − y) with = 20 MPa. We fix the transversal and vertical displacements on the section that locates the roll mill (x o = L/2; x o = 2L/5; x o = L/10), in order to prevent the sheet to widen in width of the bite. The sheet is subjected to a traction T 0 = 20 MPa.
The results obtained are shown in Fig. 12. We observe defects in the form of longitudinal folds located near the bite exit. Far from the bite, the sheet remains almost flat. Hence the transversal component can induce flatness defects as this can be seen in the industrial case discussed in section 3.4. The corresponding modes are longitudinal folds vanishing away from the rolls.

Flatness defects generated by a rolling process and comparison with experimental measurements
We consider a rolling process whose data are reported in Table 1. This corresponds to the last stand of a tinplate sheet mill, with very low thickness. For this case, the experimental stress field σ xx away from the rolls is available and it has been reported in Figs. 2 and 15. It oscillates around the prescribed mean value of 100 Mpa, but the difference between the maximal (∼ 115 MPa) and the minimal stress (∼ 85 MPa) is sufficient to induce buckling during the tension release. The interested reader can refer to [34] for a recent study about tension release. Let us recall that the rolling code LAM3 predicts a huge compressive stress near the edges because the 3D finite elements are not adapted for sheet bending and buckling. The constitutive law relating yield stress and equivalent strain is as follows: The sheet is subjected to a front tension of 170 MPa in the upstream domain and to a back tension 100 MPa in the downstream domain.
A first calculation is performed with the rolling code LAM3 to recover all the components of the stress. To avoid a double counting, the downstream stress is subtracted. Since the meshes of LAM3 and the shell model are different, MLS (Moving Least Squares) method is applied for transferring residual stresses to the new coupled model. Figure 13 shows the stress profiles resulting from the code LAM3. The actual stress field is much more intricate as in the previous analytical case. These predicted stresses are very large in the bite where the transverse stress is larger than the longitudinal one. From the point of view of flatness, the most important information lies in the right part of this figure corresponding to downstream. The transverse component σ yy decays and one tends to an uniaxial stress σ xx (y) according to Saint-Venant principle. Nevertheless this MPa) just after the bite and this can be sufficient to generate transverse buckles. One sees a compressive longitudinal stress that is very large and unrealistic on the sheet sides (∼ −800 MPa).
Next the new nonlinear model described in section "A new numerical model for sheet buckling" is run, first by increasing the tensile load up to its nominal value of T 0 = 100 MPa, second by increasing the residual stress up to the value calculated by LAM3. After the first calculation by LAM3, we recover initially the final stresses from LAM3 and we subtract the elastic stresses to obtain the residual stresses. This procedure is performed after the resolution of the traction problem and before post-buckling calculation. The resulting deformed shapes will be compared with the one by Abdelkhalek [4]. As for the longitudinal stress away of the rolls, it can be also compared with the experimental measurements and with Counhaye model [10]. The computed manifest defects are drawn in Fig. 14. The present model and the one of Abdelkhalek predict longitudinal folds after the bite and they are qualitatively similar, but quantitatively different: for the present model, 0.8 mm of defect amplitude over a 220 mm of length and 430 mm of width, while Abdelkhalek predicts respectively 0.1 mm, 150 mm and 400 mm. The difference seems important, but these deflections are rather small, which means that the bifurcation point is close and, in the neighborhood of a bifurcation, a small deviation of the data can lead to a large difference of the result.
The most important results are reported in Fig. 15, where the distribution of the longitudinal stress away of the bite is plotted. The results of the present model can be compared with those of the two aforementioned models and with online experimental measurements by flatness rolls. The unrealistic starting stress (−800 MPa) is corrected by the three models, but the result of the Counhaye's model remains unsatisfactory: it gives a compressive stress (−60 MPa) that should be tensile according to the experimental evidence and the two other models. Certainly the cause is the bad convergence of the nonlinear algorithm in this procedure [31]. Globally the present model and the one of [4] are in agreement with the experimental profile. One observes strong oscillations of the stress along the width in the range 85 − 120 MPa, whose details are not easily captured because of the resolution of the sensor (∼ 50 mm). The present model manages to follow the oscillations of the stress profile given by experimental measurements better than the method of the thesis [4]. Thus the considered domain is the main difference between the two models, Abdelkhalek working in the subdomain (x > x 0 ) and applying clamped boundary conditions on the line x = x 0 . Hence it seems important to discretize the whole domain at each step of the calculation as with the present approach.
There are slight differences between models and experiments near the edges that can be explained by the resolution of the sensor or by a local inaccuracy in the transfer by MLS.
Finally we present in the Fig. 16 the deformed shape of the sheet when the tensile stress T 0 is released. Globally,

Present model
Abdelkhalek model the present model and the one of [4] lead to a center-wave defect that is consistent with the stress profiles of Fig. 15a, b, with a small quantitative difference. Abdelkhalek predicts a slightly shorter wavelength (315 mm instead of 357 mm) and a slightly larger deflection (1.91 mm instead of 1.69 mm). It would be not reasonable to comment further because the behavior after tension relaxation is very sensitive and does not depend only on the value of the stress, but mainly on the difference between minimal and maximal stresses, which would require very accurate calculations and measurements.

Concluding remarks
We have proposed in this study a numerical model which consists in coupling Arlequin and Asymptotic Numerical Methods to simulate flatness defects observed in rolling process. Our attention was focused on very thin sheets where these phenomena are often observed. Note that in thin sheet rolling, we need three-dimensional model to compute correctly the deformations in the bite and a shell model at the downstream roll to compute buckling modes. Both models are coupled in this paper by using Arlequin method. We emphasize that Arlequin method represents significant analysis possibilities for carrying out the present code. This method allows us to couple the phenomena at the upstream roll mill with the buckling phenomena at the downstream domain. This is a direct coupling using a three dimensional model to represent the upstream domain and a shell model to represent the downstream domain. However, the challenge is to choose appropriately the ingredients of the method, in particular the size of the gluing zone and the weighting functions. In the proposed study, we limit ourselves to the effectiveness and validation of the Arlequin coupling between 3D and shell models in the framework of buckling computation under residual stresses. An originality of this work lies in the coupling area which varies during the process simulation. Since the resulting problem is nonlinear, it is solved using Asymptotic Numerical Method which is an efficient tool for solving problems involving instabilities. We have proposed first a simplified model by introducing residual stresses by analytical formulae. Moreover, we performed an industrial rolling test case taking into account all the components of the residual stresses obtained by the help of the rolling code LAM3. The results were compared with those of the experimental and reference models. We retain that our model gives satisfactory results. It predicts accurately the relaxed stresses after buckling and then the flatness defects. Note that our procedure can be easily used to take into account the presence of several defects localized in thin sheets ( "quarter-waves", "herringbone-waves", etc). Moreover, the model can be introduced directly in the rolling code (LAM3).

Appendix A: Details of asymptotic numerical algorithm
Asymptotic Numerical Method (ANM) is a technique for solving nonlinear partial differential equations based on Taylor series with high truncation order. It has been proved to be an efficient method to deal with nonlinear problems in fluid and solid mechanics [19,[21][22][23][24]26]. This technique consists in transforming a given nonlinear problem into a sequence of linear ones to be solved successively, leading to a numerical representation of the solution in the form of power series truncated at relatively high orders. Once the series are fully determined, an accurate approximation of the solution path is provided inside a determined validity range. Compared to iterative methods, ANM allows significant reduction of computation time since only one decomposition of the stiffness matrix is used to describe a large part of the solution branch without need of any iteration procedure. First, the procedure allows expanding the unknown variables of the problem in the form of power series with respect to a path parameter a and truncated at order N . For the considered problem 6, we set: ⎛ The series thus formed is composed of N sequences of the unknown variables and with the initial state of the problem given for order 0. For simplicity, we assume that there are no forces applied on the boundary of the coupling