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

Kékéli Kpogan, Hamid Zahrouni, Michel Potier-Ferry, Hachmi Ben Dhia

► To cite this version:

Kékéli Kpogan, Hamid Zahrouni, Michel Potier-Ferry, Hachmi Ben Dhia. Buckling of rolled thin sheets under residual stresses by ANM and Arlequin method. International Journal of Material Forming, 2017, 10 (389-404), 10.1007/s12289-016-1288-5. http://01515212

HAL Id: hal-01515212

https://hal.univ-lorraine.fr/hal-01515212

Submitted on 12 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Buckling of rolled thin sheets under residual stresses by ANM and Arlequin method

K. Kpogan<fnref n="1" />,<fnref n="2" /> · H. Zahrouni<fnref n="1" />,<fnref n="2" /> · M. Potier-Ferry<fnref n="1" />,<fnref n="2" /> · H. Ben Dhia<fnref n="3" />

**Abstract** 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 “edges-waves” and “center-waves” defects. The obtained results show the effectiveness of our global approach.

**Keywords** Rolling · Buckling · Residual stresses · Arlequin · Asymptotic numerical method

## 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–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 semi-analytical 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$.

<footnote n="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</footnote>

<footnote n="2">² Laboratory of Excellence on Design of Alloy Metals for low-mAss Structures (DAMAS), Université de Lorraine, Metz, France</footnote>

<footnote n="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</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Simplified schematic of the rolling process and flatness defect</figcaption></figure>

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–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

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 bite 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

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Comparison of longitudinal stress computed with LAM3 code and experimental measurements considered far enough from the bite [1, 4]</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3** New numerical approach for rolling simulation: coupling between a 3D model and a shell model with moving gluing zone</figcaption></figure>

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 applications”, 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 be 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 $\sigma^{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 $\varepsilon$ expressed in terms of $u^{3D}$ as: $\varepsilon = \frac{1}{2}(\nabla u^{3D} + \nabla^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:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_1} {}^t\sigma : \delta\varepsilon \, d\Omega = 0 \\[1em] \displaystyle \sigma = \mathbb{C} : \varepsilon + \lambda_{res}\sigma^{res} \\[1em] \displaystyle \varepsilon = \frac{1}{2}(\nabla u^{3D} + {}^t\nabla u^{3D}) \end{array} \right. \qquad (1) $$

where $\sigma$ denotes the Cauchy stress tensor, $\delta\varepsilon$ is the virtual strain tensor and $\mathbb{C}$ is the elasticity tensor. In the constitutive law, $\sigma^{res}$ denotes residual stresses derived from LAM3 code and $\lambda_{res}$ a scalar parameter ($0 \le \lambda_{res} \le 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:

$$ u^{\text{shell}}(\theta^1, \theta^2, \theta^3) = v(\theta^1, \theta^2) + \theta^3 w(\theta^1, \theta^2) \qquad (2) $$

where $(\theta^1, \theta^2, \theta^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 $\gamma$. 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 $\gamma^c$ and an additional one $\tilde{\gamma}$ so that $\gamma = \gamma^c + \tilde{\gamma}$. 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:

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_2} {}^t S : \delta\gamma^c \, d\Omega = \lambda \int_{\partial\Omega_2} f \delta u^{\text{shell}} \, d\Gamma \\[1em] \displaystyle \int_{\Omega_2} {}^t S : \delta\tilde{\gamma} \, d\Omega = 0 \\[1em] \displaystyle S = \mathbb{C} : (\gamma^c + \tilde{\gamma}) + \lambda_{res} S^{res} \\[1em] \displaystyle \gamma^c = \frac{1}{2}\{\nabla u^{\text{shell}} + {}^t\nabla u^{\text{shell}}\} + \frac{1}{2}{}^t\nabla u^{\text{shell}} \nabla u^{\text{shell}} = \gamma^l(u^{\text{shell}}) + \gamma^{nl}(u^{\text{shell}}, u^{\text{shell}}) \end{array} \right. \qquad (3) $$

## Arlequin coupling

The Arlequin framework consists in partitioning a mechanical system, initially represented by a domain $\Omega$, into two overlapping subdomains $\Omega_1$ and $\Omega_2$ (see Fig. 4) where the intersecting domain is called the gluing zone, $\Omega_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

$$ C(\mu, u) = \int_{\Omega_c(x_0)} \{\kappa\mu u + \varepsilon(\mu) : \mathbb{C} : \varepsilon(u)\} \, d\Omega \qquad (4) $$

where $u$ and $\mu$ are respectively the displacement and the Lagrange multiplier field and $\kappa$ 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 $\Omega_c$ varies during the rolling process. In addition we introduce the weighting functions $\alpha_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:

$$ \left\{ \begin{array}{l} \alpha_1 = 1 \quad \text{in } \Omega_1 \setminus \Omega_2 \\ \alpha_2 = 1 \quad \text{in } \Omega_2 \setminus \Omega_1 \\ \alpha_1 + \alpha_2 = 1 \quad \text{in } \Omega_c \end{array} \right. \qquad (5) $$

Find the global displacement field in the domain $\Omega$ is therefore to weight the displacement $u^{3D}$ and $u^{\text{shell}}$ respectively in the domains $\Omega_1$ and $\Omega_2$. In the following, let us denote $(.)^{3D} \equiv (.)_1$ and $(.)^{\text{shell}} \equiv (.)_2$.
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 $\alpha_i$ in Eqs. 1 and 3,

the Arlequin method allows us to write the following variational form:

$$ egin{cases} 
\int_{\Omega_1} \alpha_1^t \sigma : \delta\varepsilon \, d\Omega + \int_{\Omega_c(x_0)} \{\kappa\mu\delta u_1 + \varepsilon(\mu) : \mathbb{C} : \varepsilon(\delta u_1)\} \, d\Omega = 0 \\
\int_{\Omega_2} \alpha_2^t S : \delta\gamma^c \, d\Omega - \int_{\Omega_c(x_0)} \{\kappa\mu\delta u_2 + \varepsilon(\mu) : \mathbb{C} : \varepsilon(\delta u_2)\} \, d\Omega = \lambda \int_{\partial\Omega_2} f\delta u_2 \, d\Gamma \\
\int_{\Omega_2} \alpha_2^t S : \delta\tilde{\gamma} \, d\Omega = 0 \\
\int_{\Omega_c(x_0)} \{\kappa\delta\mu u_1 + \varepsilon(\delta\mu) : \mathbb{C} : \varepsilon(u_1)\} \, d\Omega - \int_{\Omega_c(x_0)} \{\kappa\delta\mu u_2 + \varepsilon(\delta\mu) : \mathbb{C} : \varepsilon(u_2)\} \, d\Omega = 0 \\
\sigma = \mathbb{C} : \varepsilon + \lambda_{res}\sigma^{res} \\
S = \mathbb{C} : (\gamma^c + \tilde{\gamma}) + \lambda_{res}S^{res}
\end{cases} 
\quad (6) $$

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. 34. 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–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–24]. Let us consider a mixed vector holding all the variables of the nonlinear problem (6):

$$ U = (u_1, u_2, \mu, \tilde{\gamma}, \sigma, S, \lambda, \lambda_{res}) \quad (7) $$

The unknown field $U$ is expanded into power series with respect to a path parameter 'a':

$$ U(a) = \sum_{i=0}^{N} a^i U_i \quad (8) $$

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 $\Omega_1$, the displacement $u_1$ is discretized by using the classical 8-node tri-linear hexahedron elements

and in the shell domain $\Omega_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$:

$$ \begin{bmatrix} K_1 & 0 & C_1 \\ 0 & K_2 & -C_2 \\ {}^t C_1 & {}^t C_2 & 0 \end{bmatrix} \underbrace{\begin{bmatrix} Q_1 \\ Q_2 \\ Q_c \end{bmatrix}}_{\{U\}^k} \underbrace{\begin{bmatrix} 0 \\ F_2 \\ 0 \end{bmatrix}}_{\{K_T\}} = \lambda^k \begin{bmatrix} 0 \\ F_2 \\ 0 \end{bmatrix} + \lambda_{res}^k \begin{bmatrix} F_2^{res} \\ F_2^{nl} \\ 0 \end{bmatrix} \quad (9) $$

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 $\delta$ 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–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 $\lambda_{res}$ and $\lambda$. 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, $\lambda_{res}$ is constant and therefore only $\lambda$ 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
<figure><img src="image_5.png" /><figcaption>**Fig. 4** The Arlequin method in a general mechanical problem. The domain $\Omega$ (a) is split into two superposed domains $\Omega_1$ and $\Omega_2$ that intersect in the gluing zone $\Omega_c$ (b)</figcaption></figure>

solution of the first step to calculate the evolution of the post-buckling of the sheet with respect to the load parameter $\lambda_{res}$. Thus, during this step the load parameter $\lambda$ 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 $\Omega_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 $\Omega_1$ and $\Omega_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.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Coupling procedure: the sheet occupies a volume domain (green). From this domain, we define a surface projection (gray) that represents the shell domain and also the coupling zone (red). A mesh is associated with each of these domains</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Residual stress profile $S_{xx}^{res}$ along the width of the strip, leading to edges-waves flatness defects (profile used by Bush et al. [32])</figcaption></figure>

## 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

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Boundary conditions. Figure (a): Coupled model: edge 1 and 2 in the shell domain are free, edge 4 is under traction $T_0$ and the face 3 in the 3D domain is clamped. Figure (b): Shell model without coupling: the vertical displacement is fixed in the zone located the upstream domain</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Propagation of the edge-waves flatness defects; Figures (a), (b) and (c): coupled model with localization of the roll mill respectively at $x_0 = L/2$, $x_0 = 2L/5$ and $x_0 = L/10$; Figures (d), (e) and (f): shell model without coupling with localization of the roll mill respectively at $x_0 = L/2$, $x_0 = 2L/5$ and $x_0 = L/10$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Parameter $\lambda_{res}$ versus the vertical displacement at point arbitrarily chosen on the strip; Figures (a), (b) and (c) correspond to the positions of the roll mill respectively at $x_o = L/2$, $x_o = 2L/5$ and</figcaption></figure>

$x_o = L/10$; Figure (d): curves for the three positions of the roll mill $x_o = L/2$, $x_o = 2L/5$ and $x_o = L/10$
<figure><img src="image_11.png" /><figcaption>Fig. 10 Propagation of the edge-waves flatness defects for various traction and for $x_0 = L/10$. a case of $T_0 = 10$ MPa, b case of $T_0 = 30$ MPa, c case of $T_0 = 50$ MPa</figcaption></figure>

with the shell code. Let us underline that Abdelkhalek's model recalled above is a pure shell model applied only in the downstream part $x > x_0$, which requires boundary conditions on the line $x = x_0$. In [4], this edge was assumed to be clamped: in the present work this artificial boundary conditions can be avoided.

Let us consider a sheet which geometrical characteristics are given as follows: width $B = 1000$ mm, thickness $h = 0.25$ mm and length $L = 2000$ mm. The elastic material constants are the Youngs modulus $E = 210$ GPa and the

<figure><img src="image_12.png" /><figcaption>Fig. 12 Propagation of local folds flatness defects. Figure (1): position of the roll mill at $x_0 = L/2$. Figure (2): position of the roll mill at $x_0 = 2L/5$. Figure (3): position of the roll mill at $x_0 = L/10$</figcaption></figure>

Poisson's ratio $\nu = 0.3$. In this case, the residual stress field is in the following analytical form:

$$ egin{cases} S^{res}(x, y) = S_{xx}^{res}(y)CDM(x, x_0) \\ S_{xx}^{res}(y) = 60 \left(-5 \left(\frac{2y}{B} - 1\right)^4 + 1\right) \\ CDM(x, x_0) = \frac{1}{2} \left(1 + \erf\left(\frac{x-x_0}{\kappa\sqrt{2}}\right)\right) \\ \erf(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt \end{cases} \quad (10) $$

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 $\kappa$ characterizes this transition. In the calculation, we choose $\kappa = 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 $\Sigma_{min} = -240$ MPa and the maximum tensile stress is $\Sigma_{max} = 60$ MPa (see Fig. 6).

<figure><img src="image_13.png" /><figcaption>Fig. 11 Residual stress parameter $\lambda_{res}$ versus the displacement at a vertical displacement point; The first figure (a) corresponds to all three times of simulation for $T_0 = 10$ MPa and the second figure (b) for $T_0 = 30$ MPa</figcaption></figure>
<table><caption>**Table 1** Characteristics of investigated rolling pass</caption><tr><td>Strip width</td><td>855 mm</td></tr><tr><td>Strip thickness (strip crown)</td><td>0.355 mm (4.81%)</td></tr><tr><td>Exit thickness</td><td>0.252 mm</td></tr><tr><td>Rolling speed</td><td>22 m/s</td></tr><tr><td>Front tension</td><td>100 MPa</td></tr><tr><td>Back tension</td><td>170 MPa</td></tr><tr><td>Type of mill</td><td>4-high</td></tr><tr><td>diameter WR</td><td>555 mm</td></tr><tr><td>diameter BUR</td><td>1300 mm</td></tr><tr><td>Length WR</td><td>1400 mm</td></tr><tr><td>Length BUR</td><td>1295 mm</td></tr><tr><td>WR crown</td><td>0.0322%</td></tr><tr><td>BUR crown</td><td>No crown</td></tr><tr><td>Work roll bending force</td><td>F=480 KN</td></tr><tr><td>Fraction law</td><td>Coulomb: τ =0.03×σ<sub>n</sub></td></tr><tr><td>Young's modulus</td><td>E=210 GPa</td></tr><tr><td>Poisson's ratio</td><td>ν =0.3</td></tr></table>

<figure><img src="image_14.png" /></figure>

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).

## Results and discussions

Figure 8 shows the propagation of the “edge-waves” buckling at three chosen simulation increments, when a tensile stress $T_0 = 50$ MPa is applied and the roll mill is localized at $x_o = L/2$; $x_o = 2L/5$; $x_o = L/10$. The buckling modes and the amplitude of the defects obtained with the

pure shell model and the coupled 3D/shell model look very similar. A more accurate comparison is done by looking at the evolution of the maximal displacement according the loading parameter (cf. 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

<figure><img src="image_15.png" /><figcaption>**Fig. 13** Stress from the LAM3 code; Figure (a): longitudinal component $σ_{xx}$; Figure (b): transverse component $σ_{yy}$</figcaption></figure>
about post-buckling behavior of rolled thin sheets, except in
[8, 34].

Traction plays an important role in rolling process. Dur-
ing the rolling, the sheet is often subjected to a tensile stress
T₀ able to hide all or part of buckling phenomena. From the
work of Fischer et al. [6], an increasing global traction T₀
does not only lead to increase the value of the critical resid-
ual stress, but it also produces shorter buckling waves con-
centrated towards the edge of the strip. This phenomenon
can be simply explained: with a large value of T₀, the buck-
ling 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 longitu-
dinal 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 val-
ues of the tensile stress: T₀ = 10MPa, T₀ = 30MPa
and T₀ = 50MPa (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₀ (see Fig. 10). In par-
allel, 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_{yy}^{res}$. Most of the papers [6, 32]
focus on the longitudinal component $S_{xx}^{res}(y)$ that is relevant
for buckling modes occurring away from the rolls. Here we
choose a profile of residual stress $S_{yy}^{res} = -\frac{4\Sigma}{B^2}y(B-y)$
with $\Sigma = 20$ MPa. We fix the transversal and vertical dis-
placements 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 $\sigma_{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

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Vertical displacement of the sheet buckled under residual stresses; Figure (a): Present model; Figure (b): Abdelkhalek's model</figcaption></figure>
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:

$$ \sigma_0 = (470.5 + 174.4\tilde{\varepsilon}) \times (1 - 0.45e^{(-8.9\tilde{\varepsilon})}) - 175 \quad (11) $$

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 $\sigma_{yy}$ decays and one tends to an uniaxial stress $\sigma_{xx}(y)$ according to Saint-Venant principle. Nevertheless this

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15 Comparison of longitudinal stress profiles computed with the present code, the existing models and experimental measurements (far enough away from the bite); Figure (a): coupled 3D-shell model;</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure (b): Abdelkhalek's model; Figure (c): Counhaye's model; Figure (d): LAM3 code without buckling</figcaption></figure>
transverse stress is compressive and significant (~ -200 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₀ = 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₀) and applying clamped boundary conditions on the line x = x₀. 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₀ is released. Globally,

<figure><img src="image_21.png" /><figcaption>**Fig. 16** Deformed shape of the sheet after releasing tensile stress T₀ for the present model (a) and the model presented by Abdelkhalek (b)</figcaption></figure>
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).

**Acknowledgments** The authors wish to thank the French National Research Agency (ANR) ANR-11-LABX-0008-01, LabEx DAMAS, for its financial support and the partners of the ANR PLAT-FORM, under contract no. 2012-RNMP-019-07, (ArcelorMittal, CEA, CEMEFMINES-PARISTECH, Constellium, Ecole Centrale of Paris, INSA Lyon and Paul Verlaine university of Metz) for the authorization to publish this work.

## 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–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:

$$ 
\begin{pmatrix} u_i \\ \mu \\ \tilde{\gamma} \\ \sigma \\ S \\ \lambda \\ \lambda_{res} \end{pmatrix} = \begin{pmatrix} u_i^0 \\ \mu^0 \\ \tilde{\gamma}^0 \\ \sigma^0 \\ S^0 \\ \lambda^0 \\ \lambda_{res}^0 \end{pmatrix} + a \begin{pmatrix} u_i^1 \\ \mu^1 \\ \tilde{\gamma}^1 \\ \sigma^1 \\ S^1 \\ \lambda^1 \\ \lambda_{res}^1 \end{pmatrix} + a^2 \begin{pmatrix} u_i^2 \\ \mu^2 \\ \tilde{\gamma}^2 \\ \sigma^2 \\ S^2 \\ \lambda^2 \\ \lambda_{res}^2 \end{pmatrix} + \dots + a^N \begin{pmatrix} u_i^N \\ \mu^N \\ \tilde{\gamma}^N \\ \sigma^N \\ S^N \\ \lambda^N \\ \lambda_{res}^N \end{pmatrix} \quad (12)
$$

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
domain. In addition, by considering any vector $\cdot \cdot_i^j$, j
denotes an asymptotic order and i represents the 3D domain

(i = 1) or the shell domain (i = 2). Substituting (12) into
(6), we obtain the following linear problem for order 1:

$$
\left\{
\begin{aligned}
& \int_{\Omega_1} \alpha_1^t \sigma^1 : \varepsilon(\delta u_1) \, d\Omega + \int_{\Omega_c} \{\kappa \mu^1 \delta u_1 + \varepsilon(\mu^1) : \mathbb{C} : \varepsilon(\delta u_1)\} \, d\Omega = 0 \\[1em]
& \int_{\Omega_2} \alpha_2 \{^t S^1 : [\gamma^l(\delta u_2) + 2\gamma^{nl}(u_2^0, \delta u_2)] + ^t S^0 : 2\gamma^{nl}(u_2^1, \delta u_2)\} \, d\Omega \\
& \quad - \int_{\Omega_c} \{\kappa \mu^1 \delta u_2 + \varepsilon(\mu^1) : \mathbb{C} : \varepsilon(\delta u_2)\} \, d\Omega = \lambda^1 \int_{\Omega_2} f \delta u_2 \, d\Gamma \\[1em]
& \int_{\Omega_c(x_0)} \{\kappa \delta \mu u_1^1 + \varepsilon(\delta \mu) : \mathbb{C} : \varepsilon(u_1^1)\} \, d\Omega - \int_{\Omega_c(x_0)} \{\kappa \delta \mu u_2^1 + \varepsilon(\delta \mu) : \mathbb{C} : \varepsilon(u_2^1)\} \, d\Omega = 0 \quad (13) \\[1em]
& \int_{\Omega_2} {}^t S^1 : \delta \tilde{\gamma} \, d\Omega = 0 \\[1em]
& \sigma^1 = \mathbb{C} : \varepsilon^1 + \lambda_{res}^1 \sigma^{res} \\
& S^1 = \mathbb{C} : [\gamma^l(u_2^1) + 2\gamma^{nl}(u_2^k, u_2^0) + \tilde{\gamma}^k] + \lambda_{res}^1 S^{res}
\end{aligned}
\right.
$$

Similarly, we can derive the following linear problem for
order k (1 < k ≤ N)

$$
\left\{
\begin{aligned}
& \int_{\Omega_1} \alpha_1^t \sigma^k : \varepsilon(\delta u_1) \, d\Omega + \int_{\Omega_c} \{\kappa \mu^k \delta u_1 + \varepsilon(\mu^k) : \mathbb{C} : \varepsilon(\delta u_1)\} \, d\Omega = 0 \\
& \int_{\Omega_2} \alpha_2 \{^t S^k : [\gamma^l(\delta u_2) + 2\gamma^{nl}(u_2^0, \delta u_2)] + ^t S^0 : 2\gamma^{nl}(u_2^k, \delta u_2)\} \, d\Omega \\
& \quad - \int_{\Omega_c} \{\kappa \mu^k \delta u_2 + \varepsilon(\mu^k) : \mathbb{C} : \varepsilon(\delta u_2)\} \, d\Omega = \lambda^k \int_{\Omega_2} f \delta u_2 \, d\Gamma \\
& \quad - \int_{\Omega_2} \sum_{j=1}^{k-1} (^t S^j : 2\gamma^{nl}(u_2^{k-j}, \delta u_2)) \\
& \int_{\Omega_c(x_0)} \{\kappa \delta \mu u_1^k + \varepsilon(\delta \mu) : \mathbb{C} : \varepsilon(u_1^k)\} \, d\Omega - \int_{\Omega_c(x_0)} \{\kappa \delta \mu u_2^k + \varepsilon(\delta \mu) : \mathbb{C} : \varepsilon(u_2^k)\} \, d\Omega = 0 \\
& \int_{\Omega_2} {}^t S^k : \delta \tilde{\gamma} \, d\Omega = 0 \\
& \sigma^k = \mathbb{C} : \varepsilon^k + \lambda_{res}^k \sigma^{res} \\
& S^k = \mathbb{C} : [\gamma^l(u_2^0) + 2\gamma^{nl}(u_2^k, u_2^0) + \sum_{j=1}^{k-1} \gamma^{nl}(u_2^{k-j}, u_2^j) + \tilde{\gamma}^k] + \lambda_{res}^k S^{res}
\end{aligned}
\right.
\quad (14)
$$

To solve the system of Eqs. 13 and 14, an additional equa-
tion is needed. In this work, we introduce a condition similar
to the arc-length type continuation condition:

$$
a = \langle \{u_1 u_2\} - \{u_1^0 u_2^0\}, \{u_1^1 u_2^1\} \rangle \quad (15)
$$

where <., > denotes the scalar product for two vectors.

For the continuation procedure, reader can refer to refer-
ences [21, 25, 26].

**References**

1. Abdelkhalek S, Zahrouni H, Potier-Ferry M, Legrand N, Montmi-
tonnet P, Buessler P (2009) Coupled and uncoupled approaches
for thin cold rolled strip buckling prediction. Int J Mater Form
2(1):833–836

2. Marchand H (2000) Modélisation de la planéité en sortie de lam-
image des produits plats (Modelling flatness in flat rolling), Ph.D.
thesis École nationale supérieure des mines de Paris

3. Hacquin A, Montmitonnet P, Guillerault J.-P. (1996) A steady state thermo-elastoviscoplastic finite element model of rolling with coupled thermo-elastic roll deformation. J Mater Process Technol 60(1):109–116

4. Abdelkhalek S (2010) Un exemple de flambage sous contraintes internes: Étude des défauts de planéité en laminage á froid des tôles minces, Ph.D. thesis École nationale supérieure des mines de Paris

5. Fischer FD, Rammerstorfer FG, Friedl N, Wieser W (2000) Buckling phenomena related to rolling and levelling of sheet metal. Int J Mech Sci 42(10):1887–1910

6. Fischer FD, Friedl N, Noe A, Rammerstorfer FG (2005) A study on the buckling behaviour of strips and plates with residual stresses. Steel Res Int 76(4):327–335

7. Rammerstorfer FG, Fischer FD, Friedl N (2001) Buckling of free infinite strips under residual stresses and global tension. J Appl Mech 68(3):399–404
8. Nakhoul R, Montmitonnet P, Potier-Ferry M (2015) Multi-scale method for modeling thin sheet buckling under residual stresses in the context of strip rolling. Int J Solids Struct 66:62–76

9. Kpogan K (2014) Simulation numérique de la planéité des toiles métalliques formées par laminage, Ph.D. thesis, Université de Lorraine

10. Counhaye C (2000) Modélisation et contrôle industriel de la géométrie des aciers laminés à froid (modelling and industrial control of the geometry of cold rolled steels), Ph.D. thesis, Université de Liège

11. Roddeman D, Drukker J, Oomens C, Janssen J (1987) The wrin-kling of thin membranes: Part i—theory. J Appl Mech 54(4):884–887

12. Abdelkhalek S, Montmitonnet P, Potier-Ferry M, Zahrouni H, Legrand N, Buessler P (2010) Strip flatness modelling includ-ing buckling phenomena during thin strip cold rolling. Ironmak Steelmak 37(4):290–297

13. Ben Dhia H (1998) Multiscale mechanical problems: the Arlequin method. C R Acad Sci IIB 326(12):899–904

14. Ben Dhia H, Rateau G (2005) The arlequin method as a flexible engineering design tool. Int J Numer Methods Eng 62(11):1442–1462

15. Ben Dhia H, Rateau G (2001) Mathematical analysis of the mixed Arlequin method. C R Acad Bulg Sci I 332(7):649–654

16. Philippe S (2009) Développement d'une formulation arbitraire-ment lagrangienne eulérienne pour la simulation tridimension-nelle du laminage de produits plats, Ph.d. thesis École nationale supérieure des mines de Paris

17. Büchter N, Ramm E, Roehl D (1994) Three-dimensional exten-sion of non-linear shell formulation based on the enhanced assumed strain concept. Int J Numer Methods Eng 37(15):2551–2568

18. Simo J, Rifai M (1990) A class of mixed assumed strain methods and method of incompatible modes. Int J Numer Methods Eng 37:1595–1636

19. Zahrouni H, Cochelin B, Potier-Ferry M (1999) Computing finite rotations of shells by an asymptotic-numerical method. Comput Methods Appl Mech Eng 175(1):71–85

20. Bischoff M, Ramm E (2000) On the physical significance of higher order kinematic and static variables in a three-dimensional shell formulation. International journal of solids and structures 37:6933–6960

21. Azrar L, Cochelin B, Damil N, Potier-Ferry M (1993) An asymp-totic numerical method to compute postbuckling behavior of

elastic plates and shells. Int J Numer Methods Eng 36(8):1251–1277

22. Tri A, Cochelin B, Potier-Ferry M (1996) Résolution des équations de navier stokes et détection des bifurcations sta-tionnaires par une méthode asymptotique numérique. Revue européenne des éléments finis 5:415–442

23. Boutyour E, Zahrouni H, Potier-Ferry M, Boudi M (2004) Bifur-cation points and bifurcated branches by an asymptotic numer-ical method and padé approximants. Int J Numer Methods Eng 60(12):1987–2012

24. Cadou JM, Potier-Ferry M, Cochelin B (2006) A numerical method for the computation of bifurcation points in fluid mechan-ics. European journal of mechanics B: fluids 25:234–254

25. Cochelin B (1994) A path-following technique via an asymptotic-numerical method. Comput Struct 53(5):1181–1192

26. Cochelin B., Damil N., Potier-Ferry M. Méthode asymptotique numérique. Hermes-Sciences, Paris

27. Tri A, Zahrouni H, Potier-Ferry M (2014) Bifurcation indicator based on meshless and asymptotic numerical methods for nonlinear poisson problems. Numerical methods for partial differential equations 30(3):978–993

28. Rateau G (2003) Méthode Arlequin pour les problémes mécaniques multi-échelles: Applications á des problémes de jonction et de fissuration de structures élancées, Ph.D. thesis, École Centrale de Paris

29. Guidault P.-A., Belytschko T. (2007) On the L² and the H¹ couplings for an overlapping domain decomposition method using lagrange multipliers. Int J Numer Methods Eng 70(3):322–350

30. Roddeman D, Drukker J, Oomens C, Janssen J (1987) The wrin-kling of thin membranes: Part I—theory. J Appl Mech 54(4):884–887

31. Abdelkhalek S, Montmitonnet P, Legrand N, Buessler P (2011) Coupled approach for flatness prediction in cold rolling of thin strip. Int J Mech Sci 53(9):661–675

32. Bush A, Nicholls R, Tunstall J (2001) Stress levels for elastic buckling of rolled strip and plate. Ironmak Steelmak 28(6):481–484

33. Hibbitt K (2001) Sorensen, ABAQUS/standard user's Manual, Vol. 1, Hibbitt, Karlsson & Sorensen

34. Abdelkhalek S, Zahrouni H, Legrand N, Potier-Ferry M (2015) Post-buckling modeling for strips under tension and residual stresses using asymptotic numerical method. Int J Mech Sci 104:126–137

35. Friedl N, Rammerstorfer FG, Fischer FD (2000) Buckling of stretched strips. Comput Struct 78(1):185–190