Imposed curvature of an elastic-plastic strip: application to simulation of coils

– This work is part of the framework of a fast modeling of winding aiming at improving knowledge of residual stress evolution in steel strips and therefore their ﬂatness during the coiling process. An exact analytical solution of an elastic-plastic strip with isotropic hardening at ﬁnite strains under an imposed transformation of curvature is developed. Issues related to ﬂow rules for non-diﬀerentiable yield functions (Tresca) have been broached and a unique solution is obtained. The equivalence for this transformation, between von Mises and Tresca yield functions is demonstrated. This solution contributes to an eﬃcient model by terms of computation times that aims at simulating coiling by taking into account inelastic deformations and enabling parametric studies in order to improve the process


Motivations
The coiling process is very common in the steel making industry.It consists in winding under tension a steel strip (coming from rolling and run out table processes) on a cylindrical mandrel (see Fig. 1).The strip has a residual stress profile, due to heterogeneous plastic deformations during the rolling process and phase changes during cooling on the run out table, that can be sufficiently compressive and lead to out of plane deformations by buckling once the strip is uncoiled and cut.This residual stress profile is therefore called a flatness defect.Many numerical simulations have been developed in order to predict the residual stress state in the strip as a function of rolling process parameters [1,2] (see [3] for a bibliographical review).Inverse methods dedicated to experimental validation of these models have been proposed (see [4][5][6] for the mechanical problem and [7][8][9] for the thermal problem) as well as a recent inverse method dedicated to residual stress evaluation using a hollow cylinder [10].
In contrast, few models enable to simulate the effect of coiling on these flatness defects [11][12][13], and therefore the final quality.A model sufficiently fast to enable parametric studies in order to develop a winding strategy under tension is desirable.Recently, a non-linear elastic model has been published [14].The latter enables to simulate the winding of a strip on itself at finite strains with very short a Corresponding author: weisz@lms.polytechnique.frcomputation times.The strip section is not necessarily rectangular but a geometrical profile can be considered and therefore one can easily model the well-known effects of industrialists of contact length reduction of the strip on itself during winding (the coil has a barrel shape).The main weakness of this model is to rely on a purely elastic behavior, that does not enable to estimate precisely the irreversible plastic deformations causing the evolution of residual stresses during the coiling process.This paper is part of this framework and aims at contributing to the development of the same model relying not on a purely elastic behavior anymore but on an elastic-plastic behavior with isotropic hardening in order to evaluate residual deformations.
The model is based on the idea that for each time-step an infinitesimal strip portion is wound on the rest of the coil by following two distinct steps.
-The first step consists in imposing a simple curvature to the strip (whose mid-plane is initially plane).The trial radius of curvature is denoted by R in the following.The strip mid-plane is transformed into a cylindrical cylinder.Since the transformation is imposed the local equilibrium is not satisfied and body forces are unduly introduced (and can be calculated thanks to the divergence of Cauchy stresses).However a global equilibrium (alike shell theory) through the strip thickness is ensured because the resultant force of body forces compensates the resultant force of surface traction (also unduly introduced).-The second step consists in making contact between the curved (step 1) infinitesimal strip portion and the rest of the coil underneath.The contact pressure depends on the radius of curvature imposed during step 1 because positions of surfaces that should be put in contact are determined by this parameter and control the contact pressure.It should be mentioned that contact pressures depend on the axial direction (strip width) because of the geometrical profile of the strip section.
Finally, stresses due to both successive steps are computed as a function of the radius of curvature (introduced in step 1), thus the resultant force of tensions along the circumferential direction is calculated.The radius of curvature is then optimized so that the latter resultant force matches the force imposed by the user (actually a torque is applied and characterizes the applied force).Although this optimization is performed numerically, computation times are extremely short because each step is solved analytically.The aim of this paper is to propose an analytical solution of step 1 consisting in imposing a transformation of simple curvature (of radius R) with an elastic-plastic behavior with isotropic hardening.An analytical solution of curvature considering an elastic-plastic behavior without hardening under infinitesimal strains assumption and without residual stresses has been proposed within the framework of coiling process [15].This work details an other approach that takes into account large rotations (for very small radii of curvature), isotropic hardening and residual stresses.It should be noted that for the considered applications, the solution at finite strains and the solution under the infinitesimal strains assumption are similar as shown in Figure 2.However, the following developments are not more complicated than under the infinitesimal strains assumption.Moreover, the solution at finite strains is more general and could be used for other applications with very high thickness/radius ratios.Two yield surfaces are studied: von Mises and Tresca.It is significant to consider the latter yield surface because it is reasonable to think that the contact problem of step 2 (which is not broached in this contribution) will be solved analytically more easily with Tresca than with von Mises.However, it is demonstrated that for this specific imposed transformation both yield surfaces are equivalent.

Influence of the loading path
As mentioned above, the strategy is to model the coiling process by dividing the unknown loading path into two distinct steps (curvature and contact).It should be noted that the imposed transformation (8) leads to unwanted body forces whose the integral through the strip thickness vanishes.Ultimately it is relevant to consider (for the curvature step) the integral of the radial stress through the strip thickness instead of the detailed profile.Moreover, the choice of this particular imposed transformation has consequences because the loading path matters due to non-linearity.Other transformations could have been considered and would have led to different results.It should be noted that a priori the most realistic loading path is unknown since in practice curvature and contact happen in the same time.A comparison between two different curvature transformations has been done in order to emphasize the influence of the transformation choice.For instance, plastic deformations and stresses present discrepancies when the transformation (8) is imposed or when heterogeneous tensions along the strip thickness at both ends of the infinitesimal strip portion are applied (in order to have a genuine boundary value problem for the curvature step).In the example presented in Figure 3, performed with the Finite Element software Castem [16], the radius of curvature is around 600 mm for a 2 mm thick strip with a Young modulus E = 210 000 Mpa, a Poisson ratio ν = 0.3, a yield stress σ 0 = 200 MPa and no hardening.Discrepancies are mainly due to the fact that the imposed transformation (8) does not allow the strip thickness to reduce during the curvature, which generates unwanted radial stresses.
The consistency of the proposed approach relies on a comparison between the complete coiling simulation (considering both step 1 and step 2) and a Finite Element modeling.Indeed, the second step, consisting in making contact with the rest of the coil, is expected to slightly reduce discrepancies between different possible loading paths during the curvature step.Since the imposed transformation does not allow variations of the strip thickness, displacements (obtained by using (8) or by applying heterogeneous tensions at both ends of the infinitesimal strip portion) are different.This latter differences are taken into account during the contact, which highly depends on nominal surfaces positions.This idea has been validated by producing, with the Finite Element software Abaqus [17], a complete numerical calculation of the coiling process for a few cycles (considering the very long computation times).Weisz-Patrault et al. [14,18] did a comparison with the presented model (which has been completed in order to take into account the step 2).Good agreement has been observed even though the FEM presents some lack of stability when plasticity is considered, very likely because of the explicit integration scheme.Therefore, the strategy consisting in applying an imposed transformation during the curvature step in order to develop a fast analytical solution seems adapted to the aimed applications.

Imposed transformation
Consider a pre-stressed strip in the reference configuration with Cartesian coordinates (X, Y, Z) where X corresponds to the coiling direction, Y the thickness direction and Z the width direction.Let (e X , e Y , e Z ) denote the basis associated to (X, Y, Z).The component e X ⊗ e X of residual stresses is prevailing and the other components are therefore neglected.The residual stress tensor is denoted by: Following developments are done at finite strains, thus it is convenient to introduce a released configuration without residual stress (which does not correspond to the pre-stressed initial state).An elastic tensor is necessarily responsible for the residual stresses in the reference configuration.This latter tensor is denoted by E (0) with J 0 = det E (0) .The isotropic behavior gives a relationship between this elastic tensor and the residual stress tensor Π (0) : This behavior formulation is obtained by considering the free energy of a compressible isotropic neo-Hookean material and an energy balance.Details are given in [19].Thus, it easily obtained that the elastic part of the transformation gradient is of the form: And then E 0 is the only real root of the following polynomial of degree 3: where: det and: A simple transformation is imposed as shown in Figure 4. Thus, the basis associated to polar coordinates in the actual configuration is defined by: e z = e Z (7) This corresponds to the following imposed transformation: The transformation gradient is calculated: The model is elastic-plastic at finite strains and the following multiplicative decomposition is used (see Fig. 5): where E is the elastic part and P the plastic part of the transformation gradient.Thus: (11) where Consider a simple form for the plastic part of the transformation gradient (which is sufficient for this problem): The plastic transformation is isochoric, therefore: The inverse of the plastic part of the transformation gradient is easily obtained: (15) Therefore, the elastic part of the transformation gradient is determined: The isotropic behavior is written as follows (see [19]): tr (E.and: ) Thus, the Cauchy stress tensor is given component wisely: The only unknowns are the fields P 1 , P 2 and P 3 .In order to determine these unknowns, von Mises or Tresca yield functions are considered.

von Mises yield surface
Let start with the von Mises yield function which is differentiable everywhere and avoids consequently the discussion about flow directions due to vertex plasticity for the Tresca yield surface.It is demonstrated afterward that both solutions obtained with von Mises and Tresca yield surfaces are identical.
Let introduce the plastic stress tensor in the released configuration (bis) defined in Figure 5 which is obtained by locally releasing the elastic tensor E: In general, this stress tensor is not necessarily symmetric, but it is in the present setting.Let note Ψ its deviatoric and symmetric part.That is to say if Thus, the plastic stress tensor in the released configuration (bis) is written component wisely: The von Mises yield surface is defined by: where p cum is the cumulative plastic strain and k(p cum ) is the yield stress assumed to be of the form: where σ 0 is the yield stress before hardening and γ a hardening parameter.In the elastic zone: Besides that, at the elastic/plastic boundary the yield criterion is verified with an equality instead of an inequality.An equation that determines univocally the plastic zones is obtained: It should be noted that for a classic boundary value problem where the transformation is not imposed, plastic zones are unknown and should be determined incrementally.In contrast, if the transformation is imposed plastic zones are determined with the yield function only.
The plastic strain rate is given by: Hence: The cumulative plastic strain rate is given by:  The flow rule is associated that is to say that the direction of the plastic strain rate is normal to the yield surface.Figures 6 and 7 enable to determine the normal vector to the von Mises yield surface (denoted by n C ): where vectors e rr , e θθ and e zz represent respectively directions of the stress components σ rr , σ θθ and σ zz in Figure 6 and e = e rr + e θθ + e zz .Thus: Ṗ2 Hence after integration: Therefore: The sign of Ṗ1 /P 1 is determined in order to ensure continuity of the radial stress at the elastic/plastic interface: Besides that, the following quantity is introduced: Finally, the equality in Equation ( 23) that holds in plastic zones can be written as follows: Therefore P 1 is the only real root of the following polynomial of degree 3: This latter root can be calculated completely analytically and enables to evaluate all the unknowns of the problem P 1 , P 2 and P 3 .

Tresca yield surface
In this section, the problem is solved by considering the Tresca yield surface.This choice is relevant insofar as the second step of the winding model (contact with the rest of the coil) will be much easier with this yield surface.In contrast, for the specific curvature problem presented in this paper, difficulties related to the nondifferentiability at the vertices of the Tresca yield surface (see Fig. 6) should be overcome.This issue is well identified, especially for numerical aspects.Regularization consisting in eroding vertices of non-smooth yield surfaces (Tresca, Mohr-Coulomb etc.) have been proposed [20,21].An other approach consists in considering that the flow direction at the vertex lies in the sub-differential and can be written as a linear combination of normal directions of surfaces that constitute the vertex.This gives rise to efficient numerical schemes [22][23][24] or more recently [25,26].In this paper the latter view is adopted and there is no regularization of the non-smooth Tresca yield surface.The flow direction is sought as a linear combination of normal directions adjacent to the vertex.The following developments show that the only stable flow direction corresponds to the flow direction of the von Mises yield surface.
The Tresca yield surface is written as follows: In the elastic zone: P 1 = P 2 = P 3 = 1, thus Ψ rr = Ψ zz using (22).And the Tresca yield criterion reduces to: Ensuring equality instead of inequality in Equation ( 40) at the elastic/plastic boundary enables to determine univocally the plastic zones: It is exactly the same equation as those obtained with the von Mises yield surface (25).Besides that, at the elastic/plastic boundary plasticity is obtained at one of the two vertices such as σ rr = σ zz .Plastic zones are defined by points where the yield criterion (computed with elastic stresses) exceeds the yield stress in Equation (40), thus in plastic zones the following equation holds: The flow rule is associated that is to say that the plastic strain rate is normal to the yield surface, however as mentioned before, the Tresa yield surface is non-differentiable at the possible points where plastic deformations occur (vertex, see Fig. 8).Therefore, the flow direction can a priori be any vector in the sub-differential.Let n α denote one of these vectors obtained by a rotation of n C with an angle α ∈ − π 6 , π 6 in the plane of deviatoric stresses (normal e = (1, 1, 1)) as shown in Figure 7.The extremal values α = − π 6 and α = π 6 corresponds to classical flow directions of one or the other faces of the Tresca yield surface.It is shown in Figure 7 that: The plastic strain rate being collinear to n α , it is obtained: Hence after integration: This latter expression is compatible with the isochoric condition P 1 P 2 P 3 = 1.Stable flow directions are sought.Clearly enough by considering (22), if α = 0 then Ψ rr = Ψ zz in plastic zones because P 2 = P 3 according to (45).In this case, the vertex of the yield surface is not activated anymore and one of the face adjacent to this vertex is activated instead.Therefore the normal direction of the yield surface becomes immediately α = −π/6 or α = π/6.Thus, the three possible stable flow directions are α ∈ {−π/6, 0, π/6}.It is shown in the following that flow directions α = −π/6 and α = π/6 are also unstable.
Indeed, if α = π/6 (a similar demonstration can be done for α = −π/6), the plastic mechanism that should be activated so that this flow direction is stable is |Ψ θθ − Ψ rr | as shown in Figure 8.By injecting α = π/6 in Equation (45): The plastic mechanism actually activated is not  (see Fig. 8).Thus E0 ≥ 0 (a similar expression would be obtained for the lower vertex of the yield surface i.e., for J 2 E 2 0 ≤ J0 E0 ).Thus, the activated mechanism depends on the computed value P 1 .If P 1 ≥ 1 then the activated mechanism is not |Ψ θθ − Ψ rr | but |Ψ θθ − Ψ zz |, which leads to the instability of the flow direction α = π 6 .In contrast if P 1 ≤ 1 then the activated mechanism is |Ψ θθ − Ψ rr | and the flow direction α = π 6 would be stable.The assumption P 1 ≤ 1 is intuitively impossible because since J 2 E 2 0 ≥ J0 E0 the plastic zone under consideration is stretched.One can propose a very simple proof without hardening.Indeed, in the plastic zone if the activated mechanism was |Ψ θθ − Ψ rr | then the equality in the Tresca yield criterion would give: This is a quadratic equation in P 2 1 whose only positive root is + J0 E0 according to (42).This proves after elementary calculations that P 2  1 ≥ 1 therefore P 1 ≥ 1.Thus, finally the activated mechanism is not |Ψ θθ − Ψ rr | and the flow direction α = π 6 is unstable.If hardening is considered the equation is not quadratic and the proof is less obvious even though the result remains true.
These considerations can be adapted for the lower vertex of the yield surface (see Fig. 8) and for the flow direction α = − π 6 .Therefore directions α = −π/6 and α = π/6 are unstable and the only stable flow direction is α = 0. Indeed, in this case Ψ rr = Ψ zz holds in both elastic and plastic zones.The obtained equations are therefore exactly the same as for the von Mises yield surface.It should be noted that the equality in the Tresca yield criterion in the plastic zones leads to the same Equation (37) as for the von Mises yield criterion because Ψ rr = Ψ zz .Both solutions are therefore identical and given by Equation (38).
A comparison between the analytical solution proposed in this paper and a finite element simulation at finite strains [16] is proposed.Results are perfectly overlapped (see Fig. 9) and this validates the exactness of the solution.For this numerical example the residual stress field and the hardening has been set to zero in order to simplify the finite element simulation.Besides this, E = 210 000 MPa, k 0 = 175 000 MPa, σ 0 = 400 MPa, the strip thickness is 2 mm and the radius of curvature is R = 100 mm.

First results
In this section, the presented solution is tested in order to generate several plastic and elastic zones through the strip thickness.A typical residual stress profile is considered and presented in Figure 10a.It is compressive at both surfaces and the core is under tension.A parametric study is proposed in order to understand the influence of geometry that depends only on the aspect ratio Y/R and material behavior that depends on the yield stress σ 0 and the hardening parameter γ.Material parameters are E = 210 000 MPa and k 0 = 175 000 MPa.The aspect ratio range is [−1/100, 1/100].The advantage to plot the stress and strain profiles versus the aspect ratio Y/R is that for any particular radius of curvature R 0 and thickness th the stress and strain profiles through the strip thickness are inferred by extracting, from the latter plots, the part corresponding to Y/R ∈ [−th/R 0 , th/R 0 ].The tested yield stresses are 200, 400 and 600 MPa and the hardening parameters are 0, 10 and 50.The left Cauchy-Green strain δ = 1/2 F tot .t F tot − 1 depends on the initial residual stress and is plotted in Figure 10b  for higher yield stresses and lower hardening parameters.The elastic zone is also larger for higher yield stresses.The influence of the hardening parameters and the location of the elastic and plastic zones are emphasized in Figure 10f were the yield surface k(p cum ) is presented.
For the global coiling model, during the second step the infinitesimal strip portion is put in contact with the rest of the coil, and is pre-stressed considering these results.An analytical solution is sought and results obtained during step 1 should clearly be interpolated so that pre-stresses at the beginning of step 2 are expressed analytically.Considering results presented in Figure 10 and the non-differentiable points at elastic/plastic boundaries, a good choice seems to interpolate each zone separately with polynomials of low degrees so that the slope changes at the elastic/plastic boundaries are not smoothed, which would be the case if only one polynomial of higher degree were used.

Conclusion
This work extends the approach proposed in Reference [14] to an elastic-plastic behavior with isotropic hardening and finite strains.The former work [14] aims at establishing a fast model of windings in order to evaluate the residual stress evolution during the coiling process.A strip submitted to an imposed transformation of curvature has been considered with von Mises and Tresca yield surfaces.Issues related to the non-differentiability of the Tresca yield surface have been broached.The exactness of the analytical solution has been basically validated using a finite element model.A parametric study has been proposed in order to define a strategy for the following step 2 that takes as inputs the results of the presented solution.The global model based the presented solution and considering the second step has been proposed [18].It should also be mentioned that in this contribution, kinematic hardening has not been broached.However, if the strip is submitted to several cycles of coiling/uncoiling, kinematic hardening may become non-negligible.This is a limitation of the presented approach and subsequent developments are needed for these applications.
. Cauchy stresses are plotted versus the aspect ratio Y/R in Figures 10c and 10d for different yield stresses and hardening parameters.The cumulative plastic strain p cum is presented in Figure 10e.Plastic strains are obviously higher
this leads to an immediate change of the flow direction and α goes from π/6 to −π/6 during the plastic increment.(Reciprocally, if the flow direction is −π/6, the activated mechanism is |Ψ θθ − Ψ rr | instead of |Ψ θθ − Ψ zz |, which leads to change immediately α from −π/6 to π/6).In order to simply demonstrate this, consider a plastic zone where: J 2 E 2 0 ≥ J0 E0 , that is to say that the upper vertex of the yield surface is activated