Asymptotic analysis of asymmetric thin sheet rolling

An analytical model for asymmetric rolling is presented, which includes asymmetry in roll friction, roll size and roll speed, for a rigid, perfectly plastic thin sheet deformed with Coulomb friction. This model is solved asymptotically, based on the systematic assumptions that both the roll gap aspect ratio and the friction coefficient are small. While the leading order solution is shown to be consistent with an existing slab model, we are able to derive additional detail by looking to higher orders. We compare our higher order solution and the leading order solution with finite element simulations, and use the results to determine the practical range of validity of our analytical model. Within this region, it gives reasonable quantitative predictions of the force and torque results from finite element simulations and approximates through thickness variation of stress and strain with orders of magnitude shorter computation times. A MATLAB implementation of this solution is included in the supplementary material. & 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Rolling is the process of reducing the thickness of a metal workpiece by passing it between two rotating rolls with separation less than the current workpiece thickness. This process is used to produce sheet metal and other products.
Asymmetry in rolling can arise through inequality in roll radii, roll velocity and interface friction; inhomogeneous or anisotropic workpiece material; or bending end forces. Regardless of the mechanism, asymmetry generally results in curvature of the rolled workpiece. If this is unintentional, perhaps a consequence of tool wear, it can cause compatibility difficulties and machine damage. Intentional asymmetry can be used to reduce the total roll torque and force required to achieve a given reduction and can provide flexibility in machine design.
The mechanical simplicity of driving only a single roll first motivated investigations into asymmetric rolling but it continues to be an active area of research for other reasons. In addition to process efficiency gains, improved workpiece quality and reduced maintenance requirements; curvature can also be desirable if controlled to produce a wider range of products.
Given these attractions, improving online control of asymmetric rolling is an area of interest. For use in control, solution times should not slow the processes when included in the control loop. This prohibits the use of finite element simulations and motivates research into faster, yet still accurate, analytical models.
Early experimentation has been used to quantify roll force and torque for a range of geometries and materials [3,22]. A variety of techniques have also been used to investigate other properties such as contact stress distributions and workpiece curvature [17,1,10].
Analytical work has included a range of approaches; the most popular being modifying one-dimensional 'slab' models of symmetric rolling. These models are constructed without systematic consideration of the physical system so the subsequent assumptions remain questionable. This seems particularly relevant for recent works [21,26,25,11,9] which extend models such as Hwang and Tzou [13] to capture greater asymmetry and predict curvature.
Alternative techniques have also included upper-bound methods [18,12] and slip-line methods [7,5]. While both are able to predict characteristics such as curvature, the roll contact points or the yield region, the solution processes require a priori knowledge or assumptions about the form of the solution. This hinders the development of these models for other geometries and materials, therefore they have seen less attention in recent work.
Finite element simulations have also been applied to provide detailed results for more general configurations. Most studies focus on predictions of roll force, roll torque, and workpiece curvature [27,24,20] with some more recent studies considering microstructure [23]. While impractical for online control, this approach can be used to gain understanding of the processes and facilitate inexpensive exploration of configurations. The results of these publications have also been used to validate some of the previous analytical models. Understanding the operating limits of a model could also be ascertained by running simulations tailored to the model in question.
Another technique has also been used in modelling symmetric rolling: asymptotic analysis. Asymptotic analysis exploits systematic assumptions of scale to find a rigorous yet tractable approximation, as opposed to simplifications through ad hoc assumptions of unknown error and limitation. A series of publications [28,16,8,4,2] develop models of rolling underpinned by an assumption of a small roll gap aspect ratio: the sheet is much thinner than the length of the roll gap. Without breaking this assumption a range of geometries, friction models and material models have been implemented.
This technique has only been applied to asymmetric rolling by Johnson [15] where asymmetries were only considered for the friction coefficients and roll speeds. The friction coefficient was assumed to be an order of magnitude larger than the roll gap aspect ratio which is not representative of many thin sheet processes which are predominantly cold rolling. Experiments [17] and simulations [19,20] also show the sign of curvature can be dependent on geometry indicating that specific account of roll size may be necessary to capture the complete dynamics of the process.
In the interest of developing analytical models with sufficient resolution to potentially make curvature and microstructure predictions, the present work develops an asymptotic model to explicitly include asymmetric roll size with asymmetric roll speed and asymmetric, small magnitude friction. Complete stress and strain fields are achieved with this approach. In Section 2 we present a model of asymmetric rolling assuming a rigid, perfectly plastic workpiece and roll-workpiece interaction driven by slipping Coulomb friction. The choice of material and friction models is for illustration only; by analogy, a solution could be found for any of the friction or material models in the literature [28,16,8,4,2]. This model is non-dimensionalised to find six non-dimensional groups: the aspect ratio, δ, and the friction coefficient μ, which are assumed to be small; the sheet reduction r; and the ratios of roll size, speed and friction, which are considered unrestricted. An asymptotic solution to this model comprises Section 3, and the model is validated against the explicit solver of the commercial finite element package ABAQUS through a range of asymmetries and parameters in Section 5.

Model formulation
We assume a plane strain configuration, which is valid away from the workpiece edges for sufficiently wide workpieces. Hence, Fig. 1 captures the extent of the model.
The rolls are vertically aligned and the workpiece is fed horizontally. The initial workpiece half thickness isĥ 0 and the length of the roll gap isl , giving the aspect ratio as δ =^ĥ l / 0 where δ is assumed to be small. This is most appropriate when considering a thin sheet with large or flattened working rolls. The material model is assumed rigid perfectly plastic; that is, no elasticity and no hardening. We assume that plastic deformation occurs everywhere in the roll gap and adopts vertical boundaries at the entry and exit (marked as the hashed region in Fig. 1). This assumption is typical of existing 'slab' and asymptotic models of rolling.
Although assuming vertical boundaries to the plastic region imposes very specific combinations of bending and shear end conditions for a given asymmetry, it has been shown experimentally that the bending effects from non-extreme end conditions can be neglected [26].
The von Mises yield criteria and associated flow rule, the Levy-Mises equations, are used and slipping Coulomb friction describes the roll-workpiece interaction like Domanti and McElwain [8]. Unlike Domanti and McElwain [8], rather like the first model of Cherukuri et al. [4], the friction coefficient is also assumed to be small, μ⪡1. Asymmetry is introduced into the friction coefficient, μ; roll radius,R; and rotational speeds, Ω, which must all be defined for the top, subscripted t, and bottom, subscripted b, rolls separately. These assumptions are valid for foil rolling, but may also be valid in other regimes that fall within these assumptions.
The velocity on the roll surfaces is restricted by the no-penetration condition. Horizontal and vertical force equilibria on the roll surfaces are combined with the Coulomb friction model to give the shear boundary condition on the top and bottom roll. The model is closed by applying a given force at each end of the roll gap.
Using carets to denote dimensional quantities, we can definep, s ij ,û,v, λ andk as the pressure, ijth deviatoric stress, horizontal velocity, vertical velocity, flow parameter and yield stress respec- is the roll surface, applicable to both top and bottom rolls andF in/out are the end tensions, per unit width, applied to the workpiece, applicable to the upstream and downstream workpiece. We define the upstream velocity of the workpiece asû 0 , although it is not possible to specify it independently of the two roll velocities. Consequently, we shall considerû 0 to be the characteristic velocity for the purpose of non-dimensionalisation then determine its value from the roll velocities.

Non-dimensionalisation
We scale vertical distances with the initial workpiece half thickness,ĥ 0 , and horizontal distances with the length of the roll gap,l . The aspect ratio, δ =^ĥ l / 0 , is assumed to be small. As the friction is also small, , we define a normalised friction coefficient, β μ δ such that δ is the sole small parameter.
Using the scaling choice of Cherukuri et al. [4], the shear stress scales with the friction coefficient and yield stress, δβ =ŝ ks xy xy .   we are now able to determine the non-dimensional governing equations, using, as a consequence of plane strain, −s xx in favour of s yy . Similarly, the velocity boundary conditions are The shear boundary conditions can be simplified by observing that the roll surface speed remains constant while the workpiece surface speed increases throughout the roll gap. This can be deduced by considering a plug-like flow and is verified by the velocity solution in the next section. The locations where the roll surface speed equals the workpiece speed are called neutral points. Coulomb friction forces can hence be more simply considered constant in three regions; Hence, in addition to non-dimensionalising with the bottom roll friction coefficient, the friction coefficients are defined piecewise to eliminate the surface slip from the problem formulation, x x x x : : and 1 : 1 : .
This allows the shear boundary conditions to be expressed as 10b The end force and velocity conditions can also be non-dimensionalised to It is also useful to define the workpiece height throughout the roll gap, , and the total roll friction acting to draw the workpiece through the roll gap, γ γ γ

Solution
We seek a perturbation expansion in δ for each of the variables, u, v, s ij , p and λ of the form Assuming δ is sufficiently small, powers of δ are considered orthogonal so like terms are collected and solved successively, starting from low orders of δ.

Leading order solution
Neglecting all terms of δ ( ) O or higher, the governing equations are reduced to Eq. (17) indicates that the leading order horizontal velocity is vertically homogeneous so the application of conservation of mass to every vertical slice of the workpiece results in Eq. (15) shows that the pressure is homogeneous through thickness, hence applying the stress results to (14) gives xy 0 0 Applying this at = y h t or = y h b and comparing to the boundary shear gives the general forms The Coulomb boundary conditions (20) applied to (28b) produces the differential equation for pressure, The pressure at the entrance and the exit can be determined from the workpiece force end conditions (22), This defines two boundary conditions for the ODE (29). However, the discontinuous nature of γ Δ means that (29) must be solved in three sections, as shown in Fig. 2 . The locations of these neutral points, x nt and x nb , are not known a priori, and must be determined as part of the solution.
x nt and x nb are the locations where, by definition, the surface velocity equals the roll velocity. The pressure, ( ) p 0 , must be continuous at these locations.
Hence, the middle piece of the pressure curve is determined by joining the outer two curves such that the neutral point locations satisfy the correct roll speed ratio; the magnitudes of the values are satisfied by the choice of characteristic velocity,û 0 . This is analogous to choosing the constant of integration marked in panel (c) of Fig. 2.
Once (28b) and (28c) have been solved for ( ) p 0 and K(x), substitution reveals the shear stress to be which completes the leading order solution.

Comparison with Hwang and Tzou [13]
Extracting the surface pressure and shear from this leading order solution is equivalent to employing a 'slab' model. In fact, this particular solution is equivalent to (10) in Hwang and Tzou [13] if we assume that the reduction is small. Hwang and Tzou [13] approximate the horizontal coordinate with an expansion of the tangent function, ω ω ω = ( ) ≃ + x tan /3 3 , to solve much of this in closed form, reducing the neutral point search to numerically inverting an algebraic equation instead of iteratively solving (29). However, this approximation is only valid in the limit of small reductions, whereas our method supports an arbitrary reduction.
Although this agreement validates the assumptions made in Hwang and Tzou [13], the rigour of the present approach is still beneficial as the assumptions are made more explicit and hence the limits are easier to test. For example, neglecting the shear stress contributions in the yield criterion at leading order is a consequence of assuming small friction coefficients so we can, in principle, determine the range of friction coefficients for which accurate predictions can be made.

Correction terms
The same solution process can be applied iteratively to solve for terms with increasing order of δ. The absence of δ ( ) O terms in the governing equations or boundary conditions mean that this order is solved to be identically zero. This suggests why the existing slab models have been generally successful; given their ad hoc assumptions are correct, they achieve accuracy up to terms of δ ( ) O 2 . Further accuracy can still be achieved with the systematic approach presented here by repeating this process for terms of δ ( ) O 2 . The correction terms increase the accuracy of the through-thickness solution. In practice each variable raises an order as a polynomial in y with each correction. Horizontal velocity, pressure, longitudinal deviatoric stresses and the flow parameter become quadratic in y and vertical velocity and shear stress become cubic. Velocity also becomes dependent on the stress distribution, therefore material properties and frictional effects affect the strain field.
For brevity, the derivation of this correction has been relegated to Appendix A and a MATLAB implementation is included in the Supplementary material.

Numerical simulations
Numerical simulations were used to validate the presented model. The commercial finite element package ABAQUS was used with the input modified from the explicit two dimensional rolling model presented in section 1.3.11 of the 'ABAQUS Example Problems Manual' [6]. Symmetry was broken with the addition of a second roll and the initialisation was modified to have the rolls close onto a stationary workpiece instead of feeding the workpiece into the roll gap with a non-zero initial velocity.
The first simulation set was run with a material close to rigidperfect plastic for numerical results as close to the asymptotic modelling assumptions as possible. The yield stress in shear was set to 173 MPa with no hardening effects. As ABAQUS requires some amount of elastic behaviour, the elastic modulus and Poisson's ratio were set as high as stable computation would allow: 200 GPa and 0.45 respectively in this case.
The symmetric base case was a 10 mm strip thinned by 12% with 2.5 m radius rolls; this roll size is realistic when approximating curvature of rolls flattened in thin sheet rolling. The friction coefficients were taken to be 0.1 and roll surface velocities of Given these base values, the top roll was then altered to vary the ratio of the top and bottom friction coefficient, surface speed and roll size. It is worth noting that δ varied with the roll size as the workpiece thickness was held constant.
A second set of simulations was made to observe the performance over a range of parameters: specifically, varying the friction magnitude, aspect ratio and reduction. One dimensionless parameter was varied while the others were held constant. Note that this means two geometric parameters may vary simultaneously; for example, the roll size was reduced as the reduction was increased to ensure the aspect ratio remained constant. Six different sets of initial parameters, as specified in Table 1, were used.
The lower roll surface speed, initial half thickness and yield stress were − 1.2 m s 1 , 0.005 m and 100 MPa respectively. The material used in these simulations was further from perfect plasticity than the previous example. This reflects a more realistic material: Poisson's ratio of 0.35 and an elastic modulus of 100 GPa.

Results and discussion
In this section, we present the comparison between the numerical simulations and asymptotic model for varying asymmetries and several cross sections of non-dimensional parameters.

Numerical comparison over varying asymmetries
Results from the leading order asymptotic solution, our second order corrected asymptotic solution and numerical simulations for the first set of simulations described in Section 4 are presented in Fig. 3.
The trends of the roll force and torque as each ratio is varied are captured well by the asymptotic solution. The median error was 0.85 MN and 0.007 MN m for the force and torque respectively with maximum errors of 2.25 MN and 0.28 MN m occurring for asymmetric speeds where the magnitudes vary the most. Considering characteristic force and torque values of 25 MN and 1.0 MN m, these median values correspond to less than 3.5% error. There is also minimal difference between the leading order and corrected asymptotic solutions.
The discrepancies in this comparison can most likely be attributed to the elastic effects, which are neglected in the asymptotic model but incorporated in the simulations.
The most phenomenologically interesting trend in both cases is the drop in force and transition in direction of torque as the roll Table 1 Default parameter sets for varying non-dimensional parameters. Furthermore, an important consequence of elasticity is that a region of sticking often occurs between the roll and workpiece. In this region the shear stress smoothly changes direction and the workpiece drops below yield. The error bars in Fig. 4 indicate the extent of this sticking region and the mid-point marked as the neutral point. Fig. 4 shows that the plateaus in force and torque occur as the neutral point reaches the end of the roll gap. The asymptotic solution predicted the location of the neutral point in all experiments with similar accuracy to that observed in Fig. 4, although the neutral point varies little while varying friction or roll size.
In the case of large roll speed asymmetry, when a roll is entirely slipping as one of the neutral points reaches an end of the roll gap, the assumption of Coulomb friction means that any further speed asymmetry will not affect the neutral points. The speed of the process, expressed by the velocity scaling, is then entirely controlled by the non-slipping roll. In finite element simulations, elastic effects mean both roll speeds remain important until both neutral points have reached opposite ends of the roll gap. This explains the inaccuracy in Fig. 4.

Numerical comparison over varying non-dimensional parameters
The roll force, roll torque and neutral point for 'Combo 2' of the second set of simulations are presented in Fig. 5. The remaining five parameter sets in Table 1, which exhibit very similar trends, are illustrated by way of absolute error for roll force and torque in Fig. 6.
When δ exceeded 0.3, the finite element simulations failed to reach a steady state solution as the rolls slipped without deforming the workpiece. This seems to be reflected in the asymptotic model as both boundary conditions cannot be satisfied with pressure continuity.
For larger friction coefficients, (typically μ ≥ 0.3 b ), the asymptotic solution breaks down: terms begin to 'jump order'. That is, correction terms become as large as leading order terms, which is a clear sign that the premise of separating orders in the perturbation expansion is incorrect. This is unsurprising considering the violation of the assumption that the friction coefficient and aspect ratio are of similar size.
Within the presented parameter range, however, the asymptotic solutions behave as one would intuitively expect and capture most of the trends exhibited by the simulations; the major discrepancy being the clear deviation in roll torque for increasing friction coefficient. Fig. 6 shows that this large error occurs only in parameter sets that have asymmetric roll speeds. Combined with the widening sticking region around the neutral point in the simulations observed in Section 5.1, this suggests that the cross shear region is influenced significantly by elasticity. This has the effect of smoothing the change in direction of surface shear and hence is likely to reduce the severe changes in roll torque observed by the rigid-plastic asymptotic model.
The increasing error in the roll force is the consequence of the simulation and asymptotic predictions increasing at different rates. This is unsurprising considering the friction coefficient is assumed to be similar in size to the aspect ratio and the imbalance becomes most significant once μ b exceeds δ = 0.1.
Variation due to changes in aspect ratio or reduction is wellcaptured by our model above δ = 0.05 and ≈ r 0.15. The poor agreement for small reduction may result from the workpiece dropping below yield, indicated by the widening sticking zone. The lower reduction rate may be insufficient to produce plastic deformation throughout the roll gap resulting in significant elastic contributions. For small δ, force and torque are generally larger in magnitude so the larger absolute error is not too troubling.

Stress and strain distributions of the asymptotic solution
Referring back to Fig. 3, it is clear that the correction terms make little difference to the force and torque predictions. Nevertheless, there is a gain in qualitative accuracy that stems from the inclusion of higher order terms. This is evident when plotting the  Table 1. Error bars indicate the finite length of sticking between the rolls and workpiece. stress distributions and velocity fields between leading order and our corrected solution: Figs. 7 and 8 respectively. Once again, the 'Combo 2' parameters have been used.
As discussed in Section 3.2, each variable gains additional resolution through thickness by using the corrected asymptotic solution. Although all solutions now exhibit top-bottom asymmetry, this is most pronounced for the horizontal velocity, pressure and longitudinal deviatoric stresses, which are homogeneous in y at leading order. Both velocities gain discontinuities at the neutral points as a consequence of their dependence on the leading order shear stress distribution via (5). These discontinuities are a necessary consequence of Coulomb friction used without elasticity or smoothing at low relative slip speeds. It is interesting to note that these discontinuities are not present in the leading order solution for ( ) u 0 and ( ) v 0 as they are independent of the friction model used.

Computational time comparison
Each asymptotic solution presented here, implemented in MATLAB, was computed in less than 21 seconds CPU time. By contrast, the finite element solutions took between 5 min and 6 h CPU time to complete in ABAQUS.
We note also that the second order correction came at a small additional cost compared to the leading order as the majority of the computation time was spent finding the neutral point. The leading order solution typically required around 20 s whereas less than one additional second was required for the second order correction.

Application to hardening material
Proper treatment of a hardening material would require us to reformulate this model with an alternative yield condition and rework the asymptotic analysis. In the interest of making a quick comparison, an ad hoc approximation is made by selecting the yield stress to correspond to the average effective strain produced by the current model. This is compared to simulations with a material with similar properties to carbon steel. Specifically, we use an elastic modulus of 180 GPa; Poisson's ratio of 0.27; and shear yield stress starting at 186 MPa hardening to 494 MPa at a strain of 1.0.
Over a range of asymmetries, the trends of the roll force and the torque in Fig. 9 appear similar to the ideal case presented in Section 5.1, however, the errors increase. The

Conclusion
A model for asymmetric rolling of rigid perfect plastic sheets with Coulomb friction has been presented and solved asymptotically to a piece-wise ordinary differential equation. This was  Table 1. achieved via the systematic assumptions that the aspect ratio, δ, and the friction coefficient, μ, are small. A MATLAB implementation is included in the Supplementary material.
The leading order asymptotic solution agrees with a 'slab' model [13] for predicting roll force and torque in the limit of small reduction, an assumption not needed by our method. The δ ( ) O 2 correction offers new predictions of the through-thickness variation of each component of stress and strain. This qualitative refinement has a relatively minor effect on the force and torque predictions, however it gains significance when modelling hardening effects (e.g. [16]) or if consideration must be made of the material micro-structure, such as when modelling dynamic recrystallisation.
The asymptotic solution was compared to finite element simulations in the most comprehensive validation of an asymmetric rolling model to date. The asymptotic model captures most trends present in the simulated force, torque and neutral point variation while taking orders of magnitude less time to compute. Specifically, it was confirmed that the model performs well within the region where δ ≤ ≤ 0. ; and asymmetries of roll size, speed and friction between 0.8 and 1.5. Outside these limits, the asymptotic and 'slab' models should be treated with suspicion. In particular, for μ ≥ 0.3 t where the solution was found to 'jump order', indicating that it should be considered invalid.
This regime corresponds to thin sheet rollingfor example, a 4 mm sheet reduced by 25% with a 0.5 m effective roll radius. The assumptions are applicable for materials with minimal hardening and high elastic rigidity compared to the yield stress, such as lead, mild steel, some aluminium alloys and more. Although the tolerance for hardening can be extended, as presented in Section 5.5, including hardening in the model formulation would be recommended to ensure accuracy in these circumstances.
In some cases, it seems that degradation in solution quality stems from the appearance of regions where the workpiece sticks to the roll surfaces. In the simulations, this results in a large subyield cross-shear, which has a strong effect on torque predictions. Although we can capture cross-shear, our rigid-plastic model is incapable of resolving these sub-yield regions correctly.
The numerics, Fig. 5 in particular, also capture an interesting oscillation in the position of the bottom neutral point as the reduction is varied. This may be related to the change in the sign of curvature observed in other studies [3,19,20]. If so, this would indicate that for a model to robustly predict curvature through reduction variations, it would require greater phenomenological detail than ours or the previous 'slab' models.
Future work should focus around more realistic materials. Although work hardening could be approximated by modifying the yield stress in this model based on some mean effective strain estimate, the asymptotic method could be used to provide a rigorous treatment for this or other hardening behaviours, like Smet and Johnson [28], Johnson and Smelser [16], Domanti and McElwain [8], Cherukuri et al. [4]. Further, incorporating elasticity and sub-yield behaviour may capture the trends missed by the present model so would be a desirable addition, although this poses a significant modelling challenge. Incorporating roll deformation could also improve predictions for foil rolling.
Finally, the prediction of curvature has been attempted by several authors [21,26,14] and the same methods could be applied to this asymptotic model and the detail gained here could underpin future curvature predictions to capture the oscillations discussed above. Zealand's kind support through the Cambridge-Rutherford Memorial Scholarship. E.J.B. also gratefully acknowledges the support of a Royal Society University Research Fellowship and a College Lectureship at Gonville & Caius College, University of Cambridge.

Appendix A. Second order corrections
Continuing the expansion in orders of δ from Section 3, we begin by outlining the second order governing equations and boundary conditions. The equations, after removing the δ ( ) O terms that were found to be zero, are