Abstract

Introduced is a new physics-based three-dimensional (3D) mathematical model capable of efficiently predicting time histories of the nonlinear structural dynamics in cold rolling mills used to manufacture metal strips and sheets. The described model allows for the prediction of transient strip thickness profiles, contact force distributions, and roll-stack deformations due to dynamic disturbances. Formulation of the new 3D model is achieved through a combination of the highly efficient simplified-mixed finite element method with a Newmark-beta direct time integration approach to solve the system of differential equations that governs the motion of the roll-stack. In contrast to prior approaches to predict structural dynamics in cold rolling, the presented method abandons several simplifying assumptions and restrictions, including 1D or 2D linear lumped parameter analyses, vertical symmetry, continuous and constant contact between the rolls and strip, as well as the inability to model cluster-type mill configurations and accommodate typical profile/flatness control mechanisms used in industry. Following spatial and temporal convergence studies of the undamped step response, and validation of the damped step response, the new model is demonstrated for a 4-high mill equipped with both work-roll bending and work-roll crown, a 6-high mill with continuously variable crown (CVC) intermediate rolls, and finally a complex 20-high cluster mill. Solution times on a single computing processor for the damped 4-high and 20-high case studies are just 0.37 s and 3.38 s per time-step, respectively.

Introduction

Dimensional quality of cold-rolled strips or sheets continues to be a primary area of research in the metal rolling industry because downstream applications require increasingly stringent tolerances, particularly in markets driven by light-weight fabrication and fully automated assembly. Rolled strip dimensional quality is broadly addressed in terms of thickness profile and flatness (shape), although the two are interrelated due to mass conservation during plastic deformation. Deviations from desired strip thickness profile and flatness arise due to localized variations in the relative plastic deformation (thickness reduction) across the strip width, which subsequently produce widthwise variations in the longitudinal residual stress within the strip (e.g., Fig. 1(a)). The cause of such deviations is largely attributed to a mismatch between the loaded roll gap profile and the incoming strip profile, while many factors contribute to this mismatch, among the most influential are roll-stack deformations involving bending, shearing, and nonuniform Hertzian contact flattening/indentation of the rolls (Figs. 1(b) and 1(c)), as well as geometric discrepancies on both the rolls and the incoming strip. For reasons of simplicity, contributions to the 3D mathematical modeling necessary to predict thickness profiles (and shape/flatness) in cold rolling have focused almost entirely on steady-state process conditions [14], while such modeling approaches may suffice for setup of the mill parameters and may provide reasonably for prediction accuracy over much of a rolling pass, transient analyses that include time-history responses of the 3D structural dynamics of a roll-stack are needed to predict strip geometric quality throughout an entire rolling pass.

Fig. 1
Depiction of (a) strip flatness (buckling) defects caused by non-uniform residual stresses; (b) 4-high mill; (c) exaggerated roll-stack deformation of 4-high mill in the loaded position, typified by bending, shearing, and Hertzian contact flattening of the rolls in addition to elastic–plastic deformation of the rolled strip; and (d) side view of complex 20-high cluster-type mill
Fig. 1
Depiction of (a) strip flatness (buckling) defects caused by non-uniform residual stresses; (b) 4-high mill; (c) exaggerated roll-stack deformation of 4-high mill in the loaded position, typified by bending, shearing, and Hertzian contact flattening of the rolls in addition to elastic–plastic deformation of the rolled strip; and (d) side view of complex 20-high cluster-type mill
Close modal

Indeed, since the rolling process is dynamic in nature, any disturbance (discrete or periodic) leads to vibratory changes in the roll gap spacing, which may then invoke oscillatory thickness and flatness variations in the exit strip, thus affecting dimensional quality. While disturbances may decay with sufficient damping, they may cause momentary entry or exit tension and speed changes of sufficient magnitude and phase that can once again lead to changes in the roll gap spacing (such vibrations can even become self-excited as chatter). Hence, the prerequisite problem of mathematical analysis in preventing dimensional quality degradation due to dynamic disturbances is the ability to predict and understand the transient 3D thickness profile of the rolled strip, i.e., its time history.

As stated previously, an assumption widely employed in exit strip profile and flatness prediction is that of quasi-static (steady-state) operation; whereas models that have included rolling dynamics are generally simplified as linear, lumped parameter system models in which each roll is assumed to be a single block of mass. Notwithstanding the major disadvantages of this (discussed below), the benefit is that it simplifies the rolling mill structural dynamics into a basic spring-mass-damper system operating in a given dimension/plane.

Table 1 classifies available dynamic models from the literature according to the method of structural dynamics implemented. As indicated in Table 1, linear lumped parameter models can be subdivided into two categories: (1) unimodal (Fig. 2(a)), being a simple 1D linear vibration system, can provide reasonable quantitative predictions on modes of vibration (for the first few frequencies) in the vertical direction [5,25] and may be sufficient to study chatter due to negative damping [25]; and (2) multimodal (Fig. 2(b)), where the coupling of two or more principal modes of vibration is possible since motions of the rolls are permitted in more than one mutually perpendicular direction.

Fig. 2
Depiction of simplified linear lumped parameter systems for dynamic analysis: (a) unimodal schematic of a 4-high mill and (b) 2D multimodal schematic of the cross section of a single roll. Note that neither type can capture 3D dynamic roll-stack deformations, due to lack of transverse (out-of-plane) representation along the strip width direction [25] (Reprinted with permission from Elsevier © 1998).
Fig. 2
Depiction of simplified linear lumped parameter systems for dynamic analysis: (a) unimodal schematic of a 4-high mill and (b) 2D multimodal schematic of the cross section of a single roll. Note that neither type can capture 3D dynamic roll-stack deformations, due to lack of transverse (out-of-plane) representation along the strip width direction [25] (Reprinted with permission from Elsevier © 1998).
Close modal
Table 1

Existing implementations of mill structural dynamics in modeling cold rolling of strip/sheet

Method/modelInvestigatorAssumptionsComments on model characteristics
Linear lumped parameter: UnimodalYarita et al. 1978 [5]
Tamiya et al. 1980 [6]
Tlusty et al. 1982 [7]
Chefneux et al. 1984 [8]
Kapil et al. 2014 [9]
Lin et al. 2003 [10]
Kimura et al. 2003 [11]
  • Neutral point coincides with exit plane

  • Relations of various parameters with rolling force under dynamic conditions remains the same as under steady-state conditions

  • Variation in the entry tension represents negative damping

  • Vertical symmetry

  • Motion of the rolls are in-phase

  • Lacks capability to analyze the effects of roll bending and housing torsional modes

  • A simplified 1D model which can be used to study the effects of various parameters on the roll vibrations

  • Plane-strain slab force model assumption

  • Flow of metal through roll-bite forms the basis of the rolling dynamics

  • Provides a reasonable quantitative analysis of the frequencies of first few modes of vibration

  • Simple and efficient model for investigation of chatter due to negative damping

  • Lacks 3D bulk-body deformation effects of the roll stack

Linear lumped parameter: MultimodalYun et al. 1998 [1214]
Hu et al. 2000 [15,16]
  • Vertical symmetry

  • An improved model to estimate force under dynamic conditions which accounts for the rate of change of roll gap along with variation of the roll gap

  • Relaxed the assumption of neutral point coinciding with exit plane

  • Allows for the roll movement component to be resolved in vertical direction as well as rolling direction, but no 3D representation

Kim et al. 2012 [17]
Lu et al. 2020 [18]
Mehrabi et al. 2015 [19]
Niroomand et al. 2015 [20]
  • Vertical symmetry

  • Movement of the rolls are modeled independently

Brusa et al. 2010 [21]
  • Vertical symmetry

  • Lumped 2D parameter model to analyze the dynamic behavior of 20-high cluster mill

  • No provision for loss of contact

3D modelSun et al. 2014 [22]
Kapil et al. 2016 [23]
Brusa et al. 2009 [24]
  • Vertical symmetry

  • Mathematical formulation in vertical plane only

  • Mill tilting effects are neglected

  • Vertical deflection of backup roll neglected

  • Algorithmic damping

  • Shear deformation of the rolls neglected

  • 4-high mill formulation aimed to study the chatter vibration employing transverse representation of the rolls

  • Formulation cannot handle mode coupling

  • No provision for loss of contact

  • Application of load is simplified for the purpose of chatter analysis, which reduces the accuracy of the strip profile prediction and in turn the actual shape and movement of rolls

  • Formulation does not allow for the phase difference to exist

  • No modeling of roll crowns or inclusion of flatness control mechanisms

Method/modelInvestigatorAssumptionsComments on model characteristics
Linear lumped parameter: UnimodalYarita et al. 1978 [5]
Tamiya et al. 1980 [6]
Tlusty et al. 1982 [7]
Chefneux et al. 1984 [8]
Kapil et al. 2014 [9]
Lin et al. 2003 [10]
Kimura et al. 2003 [11]
  • Neutral point coincides with exit plane

  • Relations of various parameters with rolling force under dynamic conditions remains the same as under steady-state conditions

  • Variation in the entry tension represents negative damping

  • Vertical symmetry

  • Motion of the rolls are in-phase

  • Lacks capability to analyze the effects of roll bending and housing torsional modes

  • A simplified 1D model which can be used to study the effects of various parameters on the roll vibrations

  • Plane-strain slab force model assumption

  • Flow of metal through roll-bite forms the basis of the rolling dynamics

  • Provides a reasonable quantitative analysis of the frequencies of first few modes of vibration

  • Simple and efficient model for investigation of chatter due to negative damping

  • Lacks 3D bulk-body deformation effects of the roll stack

Linear lumped parameter: MultimodalYun et al. 1998 [1214]
Hu et al. 2000 [15,16]
  • Vertical symmetry

  • An improved model to estimate force under dynamic conditions which accounts for the rate of change of roll gap along with variation of the roll gap

  • Relaxed the assumption of neutral point coinciding with exit plane

  • Allows for the roll movement component to be resolved in vertical direction as well as rolling direction, but no 3D representation

Kim et al. 2012 [17]
Lu et al. 2020 [18]
Mehrabi et al. 2015 [19]
Niroomand et al. 2015 [20]
  • Vertical symmetry

  • Movement of the rolls are modeled independently

Brusa et al. 2010 [21]
  • Vertical symmetry

  • Lumped 2D parameter model to analyze the dynamic behavior of 20-high cluster mill

  • No provision for loss of contact

3D modelSun et al. 2014 [22]
Kapil et al. 2016 [23]
Brusa et al. 2009 [24]
  • Vertical symmetry

  • Mathematical formulation in vertical plane only

  • Mill tilting effects are neglected

  • Vertical deflection of backup roll neglected

  • Algorithmic damping

  • Shear deformation of the rolls neglected

  • 4-high mill formulation aimed to study the chatter vibration employing transverse representation of the rolls

  • Formulation cannot handle mode coupling

  • No provision for loss of contact

  • Application of load is simplified for the purpose of chatter analysis, which reduces the accuracy of the strip profile prediction and in turn the actual shape and movement of rolls

  • Formulation does not allow for the phase difference to exist

  • No modeling of roll crowns or inclusion of flatness control mechanisms

It is important to note, however, that all linear lumped parameter systems inherently neglect 3D bulk-body deformation behavior, particularly the transverse (widthwise) distributions of bending, non-uniform Hertzian contact flattening, and shear-type roll-stack deformations. Indeed, neither unimodal nor multimodel linear lumped systems allow for analysis of the effects of bending and torsional modes on the transient roll-stack deformation behavior [26]. Since these “non roll-stack” linear lumped model formulations inherently lack 3D bulk-body deformations, vital information is missed regarding the time-varying strip thickness profile and flatness, as well as the transverse distributions of the transient contact force, contact stiffness, and displacement field. Moreover, all linear lumped parameter systems inherently lack the ability to represent industrial flatness control mechanisms in their analyses, in addition to coupled vertical and horizontal interactions that are characteristic of cluster mill configurations (e.g., 20-high in Fig. 1(d)).

Because of the foregoing deficiencies of linear lumped parameter systems, attempts have been made to develop 3D structural dynamics models, including in the study of chatter vibration. Guo et al. [26] used commercial finite element (FE) software to generate a continuum-element-based 3D model to investigate the effects of separate excitations of the roll chocks, rolls, and strip in chatter. Sun et al. [22] developed a dynamic thickness profile prediction model of a 4-high mill by combining (along the transverse/axial direction) a rolling process model (of the roll-bite contact mechanics) with a structural dynamics model. Kapil et al. [23] also developed a 4-high model that included transverse behavior by model space representation along the roll axes. However, the above two 3D dynamic models employ simplifying assumptions that limit their accuracy; these specific limitations, in the context of the modeling contributions addressed in this research, are discussed next.

One of the major simplifying assumptions in existing roll-stack models is the requirement for continuous and constant contact between interacting bodies (i.e., roll–roll and roll–strip contacts). However, since non-uniform diameter profiles (crowns) are commonly machined onto the rolls, a provision for loss of contact under dynamic conditions (as shown later in the results section) is critical to accurately predict time histories of the contact forces, the motions of rolls, and in turn the exit profile and flatness of the strip. In addition, despite the advancement of proposed 3D dynamic models over linear lumped systems, except for computationally expensive continuum-element-based FE formulations, they lack ability to accommodate typical strip profile/flatness control mechanisms, such as work-roll bending, backup roll bending, and lateral shifting of the work rolls or intermediate rolls.

Another major limitation of existing dynamic models is flexibility and applicability. Most in Table 1 are only applicable to the mill configuration in which the original formulation was developed (e.g., 4-high). Indeed, even models with some degree of flexibility are still limited to vertically oriented mills such as 2-high, 4-high, or 6-high, i.e., in-plane roll-stack configurations. Models needed to investigate complex cluster mills require multiple out-of-plane contact interactions and have seldom been proposed. While Brusa et al. [21] investigated the dynamic behavior of a 20-high cluster mill (Z-mill), the structural dynamics are represented using a simplified linear lumped parameter system. In an effort toward the design of 20-high cluster mills, Brusa and Lemma [24] used a multibody dynamics approach, but assumptions of symmetry and rigid bodies restrict accuracy for any time-history analyses of the roll-stack and exit strip parameter predictions. Moreover, flexibility to accommodate irregular roll profiles such as tapered profiles, parabolic crowns, or more advanced continuously variable crown (CVC) machined roll profiles under dynamic conditions is almost non-existent in the literature.

Apart from the limitations described above, most of the formulations mentioned and listed in Table 1 are primarily aimed at investigating the conditions of dynamic instability based on considering the roll-bite process parameters under an assumption of plane strain. As a result, the models are not sufficiently formulated to provide real-world time history of the 3D structural dynamics—particularly the transient strip dimensional quality.

Accordingly, this work introduces and demonstrates a general-purpose, nonlinear 3D dynamic roll-stack model formulation that gives an efficient time-history analysis in cold rolling, including the transient predictions of strip thickness profile, contact force distributions, and roll-stack deformations. As is shown, the model is applicable to accommodate various types of flatness control mechanisms and mill configurations, including 20-high cluster mills. Note that while chatter prediction with multi-stand analysis is not addressed, the model described can be readily adapted for such a purpose.

Mathematical Model

Described next is a brief derivation of the original, static simplified-mixed finite element method (SM-FEM) roll-stack model and its mathematical coupling with a dynamic time integration technique. The resulting consolidated 3D structural dynamics rolling model is capable of providing highly efficient time histories of the exit strip thickness profile as well as the transverse distributions of contact forces and roll displacements. The case study results that follow the mathematical derivation include model validation and demonstrate its ability to efficiently simulate the 3D dynamic behavior considering important flatness control mechanisms and mill stand configurations.

Static Three-Dimensional Roll-Stack Model.

An overview of the main equations to formulate the 3D static roll-stack deformation model using the SM-FEM method [27] is given briefly next, corresponding to the schematic in Fig. 3. 3D Timoshenko beam finite elements are used to represent bending and transverse shear roll deformations. The domains of these 3D beam elements are represented schematically by lines, although their cross-section properties can represent those of any physical beam; in the formulation here, sections are solid circular discs. In accordance with mesh discretization, the disc widths correspond to the beam element lengths, i.e., distances between consecutive nodes along each beam axis. Twin-beam-coupled Winker foundation elements are used to represent nonlinear Hertz-type elastic contact flattening of the rolls (expressed analytically per [28]). Winkler foundation elements are also used to represent the equivalent load versus plastic deformation relation of the rolled strip, obtainable from either measurement data or a suitable roll-bite model [2931]. A significant advantage of this coupled beam/foundation element formulation is that, in contrast to a continuum FE approach (e.g., with isoparametric hexahedral or tetrahedral elements), the formulation eliminates the need for an enormous number of small elements to accurately capture deformations at the contact interfaces. Excellent prediction agreement has been shown for the static SM-FEM method in both computational validation with large-scale continuum FEM and experimental validation against industry data [27,32,33].

Fig. 3
Schematic of static simplified-mixed finite element method (SM-FEM) formulation: (a) partial front view of roll-stack showing continuous, nonlinear Winkler elastic foundations coupled between 3D shear-deformable (Timoshenko) beam elements (used for roll–strip and roll–roll contact); (b) cross-sectional view of upper work-roll and strip (or sheet); (c) two foundation-coupled Timoshenko beam elements with degrees-of-freedom indicated on left nodes
Fig. 3
Schematic of static simplified-mixed finite element method (SM-FEM) formulation: (a) partial front view of roll-stack showing continuous, nonlinear Winkler elastic foundations coupled between 3D shear-deformable (Timoshenko) beam elements (used for roll–strip and roll–roll contact); (b) cross-sectional view of upper work-roll and strip (or sheet); (c) two foundation-coupled Timoshenko beam elements with degrees-of-freedom indicated on left nodes
Close modal
Equation (1) gives the static SM-FEM nonlinear system in which [KG] is the global stiffness matrix and f is the load vector. Both [KG] and f are dependent on displacement vector u [27].
[KG(u)]u=f(u)
(1)
Nodal displacement vectors, uj, at node j include the 3D translation/rotation degrees-of-freedom (DOF):
uj={ujvjwjθxjθyjθzj}
(2)
Matrix [KG], the global stiffness in Eq. (1), is the superposition of all nonlinear continuous foundation stiffness matrices, [KF], and the combined Timoshenko matrices, [KT], with the latter containing bending and shear stiffness:
[KG]=[KF]+[KT]
(3)
For one twin-coupled beam and foundation element i between two bodies, 1 and 2, which represent, respectively, either the strip and a work roll, or two adjacent rolls (see Fig. 3(a)), the foundation-coupled element contribution to the global stiffness [KG1,2,i] has size 24 × 24 and includes all six translation and rotation DOF at each of the four nodes, as in Eqs. 4(a)4(c)
[KG1,2,i]=[KF1,2,i]+[KT1,2,i]
(4a)
where
[KF1,2,i]=[0likfeq(x)[N11]dx0likfeq(x)[N12]dx0likfeq(x)[N21]dx0likfeq(x)[N22]dx]
(4b)
[KT1,2,i]=[[KT1,i][0][0][KT2,i]]
(4c)
Term li in Eq. 4(b) is the length of the ith twin-coupled element, kfeq(x) is the equivalent Winkler foundation contact stiffness, and shape function matrix subset term [Npq] (for p,q ∈ [1,2]) in the foundation stiffness matrix [KF1,2,i] are defined in Eqs. (5)(8). Terms [KT1,i] and [KT2,i] in Eq. 4(c) are conventional 12 × 12 Timoshenko elements for contacting bodies 1 and 2, i.e., rolls (a zero matrix is assigned to the strip since bending/shear stiffness is negligible)
[N11]=[Nv1]T[Nv1]sin2θ+[Nw1]T[Nw1]cos2θ+[Nv1]T[Nw1]sinθcosθ+[Nw1]T[Nv1]sinθcosθ
(5)
[N12]=[Nv1]T[Nv2]sin2θ+[Nw1]T[Nw2]cos2θ+[Nv1]T[Nw2]sinθcosθ+[Nw1]T[Nv2]sinθcosθ
(6)
[N21]=[Nv2]T[Nv1]sin2θ+[Nw2]T[Nw1]cos2θ+[Nv2]T[Nw1]sinθcosθ+[Nw2]T[Nv1]sinθcosθ
(7)
[N22]=[Nv2]T[Nv2]sin2θ+[Nw2]T[Nw2]cos2θ+[Nv2]T[Nw2]sinθcosθ+[Nw2]T[Nv2]sinθcosθ
(8)

The inter-roll angle of inclination, θ, between bodies 1 and 2 in Eqs. (5)(8) accounts for both vertical and cluster mill configurations (θ=90deg for vertical roll-stacks as in 2-high, 4-high, and 6-high mills, whereas cluster-type 20-high mills involve several unique inter-roll angles, see Fig. 1(d)). Iterative techniques for the solution of the nonlinear static system in Eq. (1) are discussed in Ref. [4]. Note that in addition to eliminating large numbers of elements at roll–roll contact interfaces, which dramatically improves efficiency (e.g., 1610 elements and 6 s run time using matlab on a single processor for a 6-high mill versus 6.75 M elements and ∼5 h on 24 central processing units (CPUs) using a continuum FE approach, giving 8.9% maximum difference [32]), the nonlinear Winkler foundations in the formulation also exploit the concept of a “strip modulus,” which is a linearized equivalent elastic foundation derived from either the tangent or secant2 of the rolling load or plastic strain relationship at the working point [27,34,35] (see Fig. 4).

Fig. 4
Depiction of “strip modulus” as the tangent (or secant) of the unit rolling force versus plastic thickness strain (reduction) at the working point
Fig. 4
Depiction of “strip modulus” as the tangent (or secant) of the unit rolling force versus plastic thickness strain (reduction) at the working point
Close modal

In each iterative calculation with the nonlinear static model (and also at each time-step in the dynamic solution), the strip modulus represents an aggregate constitutive relation for elastic–plastic deformation of the strip; it therefore linearizes material nonlinearity in the system. Another advantage of the strip modulus is in the estimation of energy dissipation from the global system under dynamic conditions. During rolling, damping needs to be reasonably estimated and included in the time-integrated solution.

Incorporation of Flatness Control Mechanisms and Non-Uniform Contact Interference.

Inclusion of roll shifting and crossing mechanisms in the presented model is achieved by adjusting the model geometry and modifying the foundation moduli accordingly, whereas roll bending is accommodated by providing corresponding load values to the forcing vector f. The approach to incorporate specialized non-uniform roll profiles such as parabolic crowns, tapered profiles, advanced CVC profiles, as well as high-fidelity deviations such as roll grinding errors and roll roughness is readily adapted into the general contact interference algorithm as follows:

During elastic contact between two bodies (e.g., two rolls 1 and 2 in a vertical plane), the distance d12 between their axes centers at axial location x is obtained from their relative displacement, v1(x) − v2(x), as well as the respective initial vertical positions, yc
d12=[y1c+v1(x)][y2c+v2(x)]
(9)
The contact interference, δ12, is estimated using the above distance function and the respective local diameter functions, D1(x) and D2(x), as in Eq. (10) [27]
δ12(x)=D1(x)+D2(x)d12(x)
(10)
Roll diameter deviations (low and high fidelity) are readily incorporated into D1(x) and D2(x) for calculation of contact interference. The specific (unit) contact force f(x) is the product of the contact interference and the equivalent foundation stiffness, kfeq(x), where the latter is derived using a series stiffness combination of the elastic foundation stiffness of bodies 1 and 2 (strip/roll or roll/roll). The total contact force, PC, is required to equilibrate the integral of specific contact force
PC=0lf(x)dx
(11)
=0lkfeq(x)δ12(x)dx
(12)
=0l(1kf1(x)+1kf2(x))1δ12(x)dx
(13)
In Eq. (12), note that the equivalent foundation stiffness modulus, kfeq(x), for contact between two rolls can be obtained using the local unit force, f(x), based on either the Matsubara (see Eq. (14)) or Hertz/Föppl [36] plane strain relations, where E1, ν1 and E2, ν2 are the respective elastic moduli and Poisson's ratios of the two rolls. Decomposition of the equivalent foundation stiffness modulus, kfeq(x), into kf1(x) and kf2(x) for the respective rolls (which is needed to determine the contact interface location) can be done by first assuming D1 → ∞ and E1 → ∞, then solving for kf2(x) per Eq. (15), then obtaining kf1(x) from Eq. (16)
B=1ν12E1+1ν22E2
(14a)
kfeq(x)=(Bπ[0.1182+ln(D1+D2)lnBln(f(x))])1
(14b)
B2=1ν22E2
(15a)
kf2(x)=(B2π[0.1182+ln(D2)lnB2ln(f(x))])1
(15b)
kf1(x)=(1kfeq(x)1kf2(x))1
(16)

For the case of contact between the strip and a work roll, assuming the nomenclature 1 = strip, 2 = work roll, then kf1(x) simply becomes the strip modulus (see again Fig. 4). The work-roll foundation modulus, kf2(x), is found using Eq. (15), and the equivalent stiffness modulus, kfeq(x), is subsequently obtained through rearrangement of Eq. (16). To account for the well-known “edge drop” phenomenon based on the transition of foundation stiffness from plane strain to plane stress near the strip edges, the factor (|xx0|2/2d2 + 1) is applied to the strip modulus for |x| ≥ x0, where −w/2 ≤ xw/2 and x0 = w/2 − d correspond to d = 25 mm from either strip edge [27]. This modification reduces the strip modulus parabolically to half its nominal value at either strip edge, beginning 25 mm in from the edges.

The above static 3D roll-stack deformation model has been validated against both large-scale continuum FE models and industry measurements of strip thickness profile [32,33]. Described next is the adaptation of this global stiffness-based static SM-FEM method into a Newmark implicit time integration approach to obtain a highly efficient, 3D, nonlinear, general-purpose dynamic roll-stack model.

Dynamic Solution in Time.

The dynamic solution in time is obtained via a Newmark-beta (NM-beta) method to solve second-order differential equations or systems of equations (note that “Newmark” more correctly refers to a family of methods). NM-beta is a direct integration technique, meaning the second-order equations are solved in their original form at discrete points in time (without the need for transformation to first-order differential equations) and with time derivatives expressed by finite difference approximations. A brief derivation of the method is reiterated here from Ref. [37] for the spatially discretized (by SM-FEM) equation of motion, Eq. (17), where u is the displacement vector, [M], [C], and [K] are mass, damping, and stiffness matrices, respectively, f is the load vector, and t is time
[M]u¨+[C]u˙+[K]u=f(t)
(17)

In Eq. (14), matrix [K] is written concisely here, but it is actually a function of the displacement vector, u, and corresponds to [KG(u)] in Eq. (1) for the nonlinear, static SM-FEM model. Similarly, f(t) in Eq. (17) is actually f(u, t) and is obtained from f(u) in Eq. (1).

In the NM-beta method, explicit expressions for the acceleration and velocity in terms of known values (from the previous time-step) can be obtained as in Eqs. (18) and (19), respectively, where the constants, 0 ≤ γ ≤ 1 and 0 ≤ β ≤ 1/2, are based on the extended mean value theorem
u¨(t+Δt)=u(t+Δt)u(t)βΔt2u˙(t)βΔtu¨(t)(12β1)
(18)
u˙(t+Δt)=γ(u(t+Δt)u(t))βΔt+u˙(t)(1γβ)+Δtu¨(t)(1γ2β)
(19)
The explicit expressions of u˙(t+Δt) and u¨(t+Δt) in Eqs. (18) and (19) are simply substituted into the equation of motion (Eq. (17)) to obtain
[M][1βΔt2(u(t+Δt)u(t))1βΔtu˙(t)(12β1)u¨(t)]+[C][γβΔt(u(t+Δt)u(t))(11β)u˙(t)(1γ2β)Δtu¨(t)]+[K]u(t+Δt)=f(t+Δt)
(20)
which is written compactly as
[K^]u(t+Δt)=F^(t+Δt)
(21)
where [K^] and F^ are the effective stiffness matrix and effective load vector, respectively, defined by Eqs. (22) and (23)
[K^]=[K]+(1βΔt2)[M]+(γβΔt)[C]
(22)
F^(t+Δt)=f(t+Δt)+[M][1βΔt2u(t)+1βΔtu˙(t)+(12β)2βu¨(t)]+[C][γβΔtu(t)(βγ)βu˙(t)(2βγ)Δt2βu¨(t)]
(23)

With known initial conditions, Eq. (21) is solved to get the displacement vector, which is then used to obtain the velocity and acceleration using Eqs. (18) and (19), respectively.

Here, β and γ act as the weights for approximating the acceleration, and they may be adjusted for accuracy and stability. With values of γ = 0.5 and β = 0.25, NM-beta is unconditionally stable with second-order accuracy and is commonly referred to as a “constant average acceleration” method. Using γ > 0.5 induces algorithmic damping in the system. This might be necessary in some cases in nonlinear systems where stiffness is a function of displacement. In such cases, numerical methods generate undesirable higher-order, non-physical modes of vibration that need to be suppressed. However, accuracy is then reduced to first-order [38]. One of the advantages of the formulation in this work is that the static equilibrium achieved by the iterative solution during each time-step (to account for changing contact conditions and nonlinear elastic roll flattening) achieves the piecewise linearity in the time domain; thus, it eliminates the need for algorithmic damping.

Damping.

Note that, being a global stiffness-based system, the model presented here can readily accommodate damping in various DOF that reflects the actual dynamic behavior for specific mills and industrial conditions, for example, the particular strip plastic dissipation as well as frictional behavior between the roll chocks and the mill housing posts, which offers both damping and resistance to vertical movements of the rolls. Moreover, even though a computationally efficient mass and initial stiffness (MIS) proportional Rayleigh damping approach can generate unrealistic “spurious” damping at or near material yield [39], since the system is piecewise linear in the time domain and it uses a strip modulus to effectively make the system elastic or linear at each time-step instead of nonlinear elastic–plastic, the MIS approach can be suitably implemented. The NM-beta method's linearity at each time-step is indeed realized since nonlinearity from changing contact and material behavior is addressed inside the static SM-FEM model that is solved within each time-step calculation, and hence, the dynamic NM-beta solution is “blind” to the system nonlinearity. Equation (24) shows the formulation of the global damping matrix
[C]=c1[Mi]+c2[Ki]
(24)
where [Mi] and [Ki] are initial global mass matrix (using HRZ lumping) and global stiffness matrix, respectively, and c1 and c2 are coefficients of the mass and stiffness proportionalities given by Eq. (25)
{c1c2}=2ζωi+ωj{ωiωj1}
(25)
where ωi and ωj are the respective lower and upper bounds of the frequencies of interest, and ζ is the fraction of critical damping. The value of ζ can be reasonably approximated from empirical data, although ζ equivalent to 12% is applied uniformly to [C] for demonstration purposes with the case studies in the next section.

Mill Housing Stiffness Representation.

Note that the static SM-FEM roll-stack deformation model assumes a quasi-static condition wherein modeling of the mill housing is not required to accurately obtain steady-state results for the strip thickness profile, contact force distributions, and relative displacement profiles of the rolls. Under dynamic conditions, however, mill housing plays an important role and requires implementation in the model. A simple but reasonable means to represent the mill housing for typical in-plane roll-stack configurations (i.e., 2-high, 4-high, and 6-high mills) is via a unilateral spring element located at both ends of the upper backup roll, as shown in Fig. 5 (where the roll “ends” implies bearing locations within the roll chocks).

Fig. 5
Illustration of loading condition and model implementation of mill housing stiffness. Depicted are upper sections of a mill housing in (a) initial (unloaded) position, and (b) loaded position. Also shown in (a) and (b) are (i) the mill posts, (ii) the equivalent discrete stiffness representation via springs, and (iii) alternate equivalent stiffness representation that replicates the hydraulic cylinder movement in (ii), even in the loaded position.
Fig. 5
Illustration of loading condition and model implementation of mill housing stiffness. Depicted are upper sections of a mill housing in (a) initial (unloaded) position, and (b) loaded position. Also shown in (a) and (b) are (i) the mill posts, (ii) the equivalent discrete stiffness representation via springs, and (iii) alternate equivalent stiffness representation that replicates the hydraulic cylinder movement in (ii), even in the loaded position.
Close modal

It is important to note that the loading condition, i.e., the rolling force required for the strip thickness reduction, is applied in such a way as to replicate hydraulic cylinder piston strokes in an industrial scenario (see Fig. 5). But rather than each cylinder's displacement providing equal and opposite force between the housing and backup roll (see Figs. 5(a) and 5(b)(i and ii)), an alternate and equivalent configuration (Figs. 5(a) and 5(b)(iii)) can be readily implemented wherein the displacement boundary condition is applied at the top of the housing spring to replicate the movement of the hydraulic cylinder. Such a boundary condition is more realistic and closer to the actual scenario in which rolling load is applied because it directly accounts for mill housing compliance (i.e., mill “stretching”). Such a boundary condition can also be implemented according to the corresponding equal and opposite force (acting respectively on the housing and backup roll) for a given hydraulic cylinder displacement. In contrast, with a “purely” force-type boundary condition, wherein the rolling load is applied at the ends of the upper backup roll, there is no external resistance via reaction forces from the mill housing; hence, this is analogous to a zero housing stiffness condition in Fig. 5(b)(iii), assuming the stiffness at the junction between the backup roll and housing increases by the presence of a housing stiffness spring. At the opposite extreme, a “purely” displacement-type boundary condition that specifies displacements at both ends of the upper backup roll limits the movement of the rolls and is analogous to an infinite housing stiffness case in Fig. 5b(iii). Accordingly, the housing stiffness used in the case studies presented next aim to emulate the industrial scenarios and assume a (lumped) housing stiffness of 2.7 × 109 N/m. Values for specific mill housings based on empirical data can readily be implemented into the mathematical approach described.

Results and Discussion

The objective of the case studies presented next is to demonstrate the ability of the combined SM-FEM and NM-beta formulation to address the research gaps in the context of achieving a general-purpose 3D dynamic roll-stack deformation model. Note that, while not a model limitation, in the studies presented, the strip modulus is not perturbated dynamically; in other words, the input disturbance/variation to the rolling parameters that leads to dynamic behavior is not considered so significant as to modify the “working point.” This also implies that the sensitivity of the rolling force with respect to strip reduction is maintained. The case studies are presented in four categories:

  1. Model validation in which validity and temporal/spatial convergence solution characteristics for the 3D dynamic roll-stack model using the NM-beta time integration technique are discussed for a 4-high mill.

  2. Demonstration of model compatibility with conventional flatness control mechanisms (work-roll bending and work-roll crown) on a 4-high mill.

  3. Model demonstration for 6-high mill involving asymmetric, CVCmachined profiles on the intermediate rolls.

  4. Model demonstration for a complex 20-high cluster mill.

Figure 6(a) depicts the 4-high mill for cases (1) and (2), wherein work roll (WR) bending (as shown) and WR parabolic crown flatness control mechanisms are applied for (2). Figure 6(b) depicts the 6-high mill for case (3), while a side view of the 20-high is as illustrated earlier in Fig. 1(d). The strip material used in all cases is stainless steel (SS) type 301, the isotropic hardening parameters of which are listed in Table 2.

  1. Model Validation

    A typical 4-high rolling mill configuration is adapted here to demonstrate the 3D structural dynamics solution, i.e., time-history response. Since the NM-beta method is unconditionally stable with average acceleration, to ensure accuracy of the results, both spatial and temporal convergence studies are performed for the undamped dynamic case. Following this, steady-state (i.e., long-term) results for the damped case are compared to results of the original (previously validated) static SM-FEM formulation in order to validate the 3D dynamic formulation. Mill dimensions and parameters for these 4-high mill validation cases are listed in Table 3. The incoming strip is assumed to be of constant rectangular cross section with no prior reduction. No flatness control mechanisms such as roll crowns or roll bending are applied for these initial validation studies.

Fig. 6
(a) Front and side views of 4-high mill, showing WR bending force; (b) front and side views of 6-high mill; (c) illustration of “8 + 8 + 8” spatial mesh designation for beam-coupled foundation elements on upper 4-high mill section; (d) depiction of 6-high mill with intermediate rolls having (exaggerated) machined CVC profiles (opposite lateral shifting of the CVC rolls is used for thickness/flatness profile control of the strip)
Fig. 6
(a) Front and side views of 4-high mill, showing WR bending force; (b) front and side views of 6-high mill; (c) illustration of “8 + 8 + 8” spatial mesh designation for beam-coupled foundation elements on upper 4-high mill section; (d) depiction of 6-high mill with intermediate rolls having (exaggerated) machined CVC profiles (opposite lateral shifting of the CVC rolls is used for thickness/flatness profile control of the strip)
Close modal
Table 2

Material constitutive model (isotropic hardening) parameters of 301 stainless steel strip

MaterialYield strength (MPa)Strength coefficient (MPa)Strain hardening exponent
SS 301283.372987.490.7475
MaterialYield strength (MPa)Strength coefficient (MPa)Strain hardening exponent
SS 301283.372987.490.7475
Table 3

4-High mill dimensions and rolling parameters

ParameterValue
Work roll diameter76.2 mm
Backup roll diameter304.8 mm
Work roll face length305 mm
Backup roll face length305 mm
Work roll neck diameter50.8 mm
Backup roll neck diameter187.2 mm
Work-roll neck length152 mm
Backup roll neck length152 mm
Distance between bearings on work roll457 mm
Distance between bearings on backup roll457 mm
Work-roll and backup roll Young's modulus207 GPa
Strip width209.6 mm
Entry thickness2.576 mm
Exit thickness2.382 mm
Reduction ratio7.495%
Strip modulus10.14 GPa
ParameterValue
Work roll diameter76.2 mm
Backup roll diameter304.8 mm
Work roll face length305 mm
Backup roll face length305 mm
Work roll neck diameter50.8 mm
Backup roll neck diameter187.2 mm
Work-roll neck length152 mm
Backup roll neck length152 mm
Distance between bearings on work roll457 mm
Distance between bearings on backup roll457 mm
Work-roll and backup roll Young's modulus207 GPa
Strip width209.6 mm
Entry thickness2.576 mm
Exit thickness2.382 mm
Reduction ratio7.495%
Strip modulus10.14 GPa

Spatial and Temporal Convergence.

As mentioned, even though NM-Beta is unconditionally stable with average acceleration, inadequate time-step sizes can result in inaccurate solutions, as evident in Fig. 7(a) for a time-step of 1 ms with various spatial mesh refinements (see Fig. 6 for mesh nomenclature description). Accordingly, model tests involving both spatial and temporal convergence are carried out to ensure accuracy. The results are illustrated in Figs. 7(b)7(d), and the satisfaction criterion is considered to be the undamped displacement of the upper backup roll (BUR) axis midpoint converging to within 0.254 µm at 100 ms. As shown in Fig. 7(b), with a decrease in the time-step for a spatial mesh of six elements per roll (or “2 + 2 + 2” spatial discretization per the mesh description in Fig. 6), the normalized BUR midpoint displacement converges, and it is inferred that a time-step of 100 µs (10−4 s) or less is acceptable with this spatial discretization. Since it is well known that computational cost increases dramatically with a decrease in time-step size, the case studies that follow incorporate a 50 − μs time-step. Also, in finite-difference solutions to time-dependent, spatially discretized boundary value problems, the required time-step becomes increasingly smaller with a decrease in spatial element size (i.e., with increasing number of elements). From Figs. 7(c) and 7(d), however, with a time-step of 50 μs, all displacement results with the spatial mesh discretizations shown demonstrate satisfactory convergence. In particular, since spatial convergence is observed in Fig. 7(d) using eight internal and eight external elements (per Fig. 6), to reduce the computational cost a mesh with 24 elements on each roll (“8 + 8 + 8” mesh configuration in Fig. 6) is adapted for the subsequent case studies.

Fig. 7
Undamped convergence plots in 4-high mill: (a) illustration of numerical error in NM-Beta method using an inadequate time-step of 1 ms for various spatial mesh refinements, based on the normalized vertical displacement of upper backup roll (BUR) axis midpoint (see Fig. 6 for mesh interpretation in legend); (b) temporal convergence plot for “2 + 2 + 2” spatial mesh based on the vertical displacement of upper BUR axis for various time-steps at 100 ms total time; (c) spatial convergence based on vertical displacement of upper BUR axis for time-step of 50 µs with various spatial mesh refinements; and (d) spatial convergence plot for mesh specifications of “8 + 8 + 8” and “16 + 16 + 16” (per nomenclature in Fig. 6)
Fig. 7
Undamped convergence plots in 4-high mill: (a) illustration of numerical error in NM-Beta method using an inadequate time-step of 1 ms for various spatial mesh refinements, based on the normalized vertical displacement of upper backup roll (BUR) axis midpoint (see Fig. 6 for mesh interpretation in legend); (b) temporal convergence plot for “2 + 2 + 2” spatial mesh based on the vertical displacement of upper BUR axis for various time-steps at 100 ms total time; (c) spatial convergence based on vertical displacement of upper BUR axis for time-step of 50 µs with various spatial mesh refinements; and (d) spatial convergence plot for mesh specifications of “8 + 8 + 8” and “16 + 16 + 16” (per nomenclature in Fig. 6)
Close modal

Undamped Time History.

Figures 8(a) and 8(b) show the respective normalized displacement time-histories of the axes of the upper work roll (WR) and upper backup roll (BUR) on the 4-high mill for the undamped response based on a stepped displacement boundary condition (of the hydraulic cylinders applied to the upper backup roll) corresponding to Δd = 5.08 × 10−4 m in Fig. 5(b)(iii)). As expected, pseudo harmonic oscillation is observed in Figs. 8(a) and 8(b) although it is noted that points along the roll axes (e.g., roll end points) on the surface plots do not trace smooth sinusoids due to mode coupling. Considering the earlier discussion on modeling the mill housing stiffness, it is important to point out that, in the absence of modeling the housing stiffness (e.g., if fixed displacement boundary conditions were applied to both backup roll endpoints in Fig. 8) constant displacements at these endpoints would result, while the midpoint locus would trace an unrealistically smooth sinusoidal displacement. On the other hand, if constant force boundary conditions corresponding to zero housing stiffness were applied, the resulting BUR endpoint displacements, as well as the midpoint locus of the WR and BUR, would be characterized by relatively smooth sinusoidal curves (since mode coupling due to nonlinearity in the system can still result in multiple peaks). Such effects are much more evident in the movement of the BURs (Fig. 8(b)).

Fig. 8
Undamped time history of the normalized displacements in 4-high mill: (a) upper WR axis displacement and (b) upper BUR axis displacement
Fig. 8
Undamped time history of the normalized displacements in 4-high mill: (a) upper WR axis displacement and (b) upper BUR axis displacement
Close modal

Figure 9 shows the displacements of axis midpoints of both the upper and lower WRs and BURs. As observed in Fig. 9 and experimentally confirmed in Ref. [40], the vertical vibration of the upper and lower halves of the roll stack is typically asymmetrical.

Fig. 9
Undamped transient displacement of axis midpoints of WRs and BURs in a 4-high mill showing the mode coupling and phase difference between upper and lower halves of the mill
Fig. 9
Undamped transient displacement of axis midpoints of WRs and BURs in a 4-high mill showing the mode coupling and phase difference between upper and lower halves of the mill
Close modal

The roll movements are not only asymmetrical, a phase difference also exists, i.e., at some time instances (see line “a”–“a” in Fig. 9) the upper and lower rolls are moving toward each other, while at other time instances the rolls are moving away from one other (line “b”-“b” in Fig. 9). Asymmetry and phase difference can also be noted for BURs in Fig. 9 (see boxes “c” and “c*”). Also, the contribution of high-frequency vibration in both WRs is greater than in the BURs since the WRs experience relatively more freedom in movement as their motions are governed by inertia and contact forces alone. The average run time (matlab using single CPU) for the undamped case is 1.1 s per time-step.

Damped Time History.

Damped transient displacement results (with the same stepped boundary condition) for the upper WR and upper BUR are shown in Fig. 10, wherein mass and initial stiffness proportional damping (classical damping) is incorporated. Note that no high-frequency residual (spurious) vibrations exist. Figure 11 shows the respective normalized dynamic contact force distributions between the strip and the upper WR (Fig. 11(a)) and between the upper WR and the upper BUR (Fig. 11(b)). The dynamic exit thickness profile of the strip, characterized by nonuniform, time-varying reduction across the strip width as a result of the 3D bulk- body deformation and flattening of the roll-stack, is shown in Fig. 12(a).

Fig. 10
Damped time history of the normalized displacements in 4-high mill: (a) upper WR axis displacement and (b) upper BUR axis displacement
Fig. 10
Damped time history of the normalized displacements in 4-high mill: (a) upper WR axis displacement and (b) upper BUR axis displacement
Close modal
Fig. 11
Damped time history of normalized contact force distribution in 4-high mill at (a) strip–upper WR interface force and (b) upper WR–upper BUR interface force
Fig. 11
Damped time history of normalized contact force distribution in 4-high mill at (a) strip–upper WR interface force and (b) upper WR–upper BUR interface force
Close modal
Fig. 12
Damped time history for 4-high mill showing exit strip: (a) semi-thickness (top half) and (b) reduction deviation (used to assess and control flatness). Note: no longitudinal (rolling direction) velocity of the strip is considered; i.e., it is assumed that there is a “new” strip surface at each time-step irrespective of vibration frequency of the work rolls such that repeated reductions are not considered within the time-step even if such WR vibration frequency is sufficiently high. While not a model restriction, the same consideration applies to subsequent case study plots.
Fig. 12
Damped time history for 4-high mill showing exit strip: (a) semi-thickness (top half) and (b) reduction deviation (used to assess and control flatness). Note: no longitudinal (rolling direction) velocity of the strip is considered; i.e., it is assumed that there is a “new” strip surface at each time-step irrespective of vibration frequency of the work rolls such that repeated reductions are not considered within the time-step even if such WR vibration frequency is sufficiently high. While not a model restriction, the same consideration applies to subsequent case study plots.
Close modal
A critical geometric quality metric for rolled strips and sheets is the flatness, which can be either manifested physically by buckles or latently represented as non-uniform widthwise variation in the longitudinal (rolling direction) normal residual stress. Based on the principle of mass conservation, flatness control systems operate by minimizing deviations in the relative reduction (thickness strain) at points across the strip width with respect to the averaged relative reduction. Calculation of reduction deviation at any transverse location x is per Eqs. (26) and (27) where r¯ is the average reduction, and Δr(x) is the deviation from the average reduction
r¯=1nj=1j=n[H(xj)h(xj)]/H(xj)
(26)
Δr(x)=[H(x)h(x)]/H(x)r¯
(27)

Figure 12(b) shows the time-history of the reduction deviation profile (which is used to assess and control flatness) for the exit strip. It is well known that without any active flatness control mechanisms, the strip edges undergo greater thickness reduction than the central region, a phenomenon clearly seen in the flatness plot. In addition, according to the time history, note that the vibration causes different strip flatness defects to be dominant at different time instances. The average run time (single CPU) for the damped case is 0.58 s per time-step.

Comparison of Steady-State Damped Dynamic Response to Static Model Results.

Comparison of the displacement and contact force results for the steady-state solution of the 3D dynamic response against the original static SM-FEM solution (same mesh) is shown in Fig. 13. Very strong agreement is evident (peak difference in normalized displacement is ∼10−7).

  • (2)

    Inclusion of Flatness Control Mechanisms

    Damped time-history responses are now provided to demonstrate compatibility of the 3D dynamic roll-stack model with conventional flatness control mechanisms, including roll bending and roll crown, on the same 4-high mill. Although various roll bending and roll crown conditions can be readily applied to different rolls, WR roll bending (15.5 kN/side) and parabolic diameter crown (7.62 × 10−4 m) on work rolls are applied. As an additional study, once steady-state conditions prevail, a disturbance is introduced from 350 to 400 ms through an instantaneous 100% increase in the incoming strip thickness (e.g., due to lap welded strips).

    Figures 14(a) and 14(b) show the corresponding WR and BUR axes displacement time histories, as well as their relative difference (shown enlarged) to the steady-state response during the stepped thickness change. Figure 15 shows a similar plot for contact force distribution time history at the WR–BUR interface. Time histories of the exit strip thickness profile and reduction deviation are shown in Fig. 16. From Figs. 1416, the benefit of the applied WR bending and WR crown is evidenced by more uniform transverse distributions than in Case (1). During the disturbance (see enlarged relative difference plots in Fig. 14), note that the displacement of the rolls is not linearly proportional to the change in incoming thickness, i.e., with a step increase in the incoming thickness, the increase in roll spacing between the WRs not only displaces the rolls apart but also increases the bending and elastic flattening of the rolls. The combination of both phenomena governs the exit strip thickness and its profile. The presence of the mill housing, i.e., its elastic deformation, allows additional vertical displacement along with bending/flattening of the rolls. Figure 16(a) reveals correspondingly greater thickness variation (greater reduction at strip edges) during this disturbance. The average run time for Case (2) is 0.37 s per time-step.

  • (3)

    6-High Mill With Machined CVC Profile on Intermediate Rolls

    Demonstrated next is the capability of the model to accommodate a machined CVC profile on the intermediate rolls of a 6-high mill. The CVC roll diameter crown at location x defined by Cr(x) in Eq. (28) is applied to the intermediate rolls (see again Fig. 6(d))
    Cr(x)=0.01525x+0.0167318x2+0.0440634x3
    (28)

    Mill configuration and parameters used here are replicated from Zhang et al. [32] wherein the ability of the static SM-FEM model to incorporate CVC profiles was validated. Roll shifting with CVC profiles, while readily applicable, is not included for the demonstration of the time-history response. As in the previous case, a post steady-state disturbance is introduced (from 175 to 225 ms in this case) via a 100% step-change in the incoming strip thickness. Figure 17 shows the displacement time histories of the upper WR and intermediate roll (IR) axes, respectively, as well as their relative difference from steady-state. Figure 18 shows the contact force distribution at the WR–IR of interface and its relative difference for the disturbance period.

    Note that, during the disturbance, vibration of the WR does not occur at very high frequency (see Fig. 17(a)), however, due to the variational phase difference in vibrations of the strip and WR, the movement of the WR relative to the strip occurs at a relatively higher frequency, as clearly visible in the contact force distribution plot (see endpoint locus in Fig. 18). The erratic force spikes in Fig. 18(b) are also noteworthy; the cause of which stems from the CVC profile on the IR. Points on the right end of the IR (point “p” in Fig. 18(a), e.g.) continuously change contact conditions, i.e., once the rolls travel towards each other contact occurs; otherwise, it does not. To better visualize this, when steady-state is reached, point “p” in Fig. 18(a) has a zero contact force (no contact), but during vibration due to the disturbance, contact is made when the WR and strip are squeezed together due to their phase difference. A similar plot for the time history of the IR–BUR interface is shown in Figs. 19 and 20 shows the corresponding exit strip thickness profile.

  • (4)

    20-High Cluster Mill

    As mentioned, almost all existing model formulations to investigate rolling mill structural dynamics are restricted to vertical mill configurations, i.e., 2-high, 4-high, or 6-high. Cluster-type mills, such as the 20-high in Fig. 21, are seldom modeled even for static assumptions because of their inherent geometric and contact interface complexity that includes segmented bearings (rolls 8, 9, 10, and 11 in Fig. 21).

    The structural dynamics time history for the 20-high cluster mill is again examined for a disturbance defined by 100% step-change in the incoming thickness (from 700 to 750 ms). Mill geometry and parameters are taken from the static prediction study by Malik et al. [27] using the SM-FEM method, which was validated against large-scale continuum FE.

    It should be mentioned that in Ref. [27], only the upper half of the mill was considered (assuming symmetry), with rolling load applied via a fixed displacement boundary condition on the bottom surface of the strip mid-plane. The time-history case study here, however, provides the full model with rolling load applied via displacement boundary conditions on top of the segmented bearing rolls (A, B, C, and D in Fig. 21). Therefore, as a model verification, the static case corresponding to Ref. [27] is also presented for comparison. Although readily implementable, the mill housing stiffness is not included; it is widely known that 20-high mills are housed by extremely stiff (“zero crown”) structures, which are assumed to have infinite stiffness in the studies here.

    Figures 22 and 23 show the comparison of contact force distribution between the presented dynamic steady-state solution and static SM-FEM results from Ref. [27]. Note that the negative contact force at the roll edges in Fig. 23 (suggesting loss of contact), as predicted by the non-iterative version of the SM-FEM in Ref. [27], is corrected by Newton-Raphson iteration in the presented formulation and shows zero force (i.e., no contact).

    Damped time-history results for the incoming strip thickness step-change disturbance on the 20-high mill are demonstrated in Figs. 2427. The average run time (1 CPU, matlab) for the 20-high mill is 3.83 s per time-step.

Fig. 13
Comparison of 3D dynamic steady-state solutions for (a) roll axes displacements and (b) contact force distributions, on 4-high mill versus validated static SM-FEM solution using same spatial mesh discretization
Fig. 13
Comparison of 3D dynamic steady-state solutions for (a) roll axes displacements and (b) contact force distributions, on 4-high mill versus validated static SM-FEM solution using same spatial mesh discretization
Close modal
Fig. 14
Time history of normalized displacements on 4-high mill with WR bending and WR crown: (a) upper WR axis displacement and (b) upper BUR axis displacement. To better visualize the disturbance period (from 350 to 400 ms due to 100% step-change in incoming thickness), the lower enlarged images for (a) and (b) show the relative difference between displacement at the current time to that of the dynamic steady-state.
Fig. 14
Time history of normalized displacements on 4-high mill with WR bending and WR crown: (a) upper WR axis displacement and (b) upper BUR axis displacement. To better visualize the disturbance period (from 350 to 400 ms due to 100% step-change in incoming thickness), the lower enlarged images for (a) and (b) show the relative difference between displacement at the current time to that of the dynamic steady-state.
Close modal
Fig. 15
Time history of the normalized contact force distribution at the upper WR–upper BUR interface on 4-high mill (with WR bending and WR crown), as well as the relative difference between the current force and the steady-state dynamic response during the disturbance period (from 350 to 400 ms due to 100% step-change in incoming thickness)
Fig. 15
Time history of the normalized contact force distribution at the upper WR–upper BUR interface on 4-high mill (with WR bending and WR crown), as well as the relative difference between the current force and the steady-state dynamic response during the disturbance period (from 350 to 400 ms due to 100% step-change in incoming thickness)
Close modal
Fig. 16
Time history of exit strip parameters on 4-high mill (with WR bending and WR crown applied) showing (a) semi-thickness (top surface) and (b) reduction deviation. Both plots include disturbance from 350 to 400 ms due to 100% step-change in incoming thickness. Note that oscillation in the strip semi-thickness is evident, as highlighted by the inset image in (a).
Fig. 16
Time history of exit strip parameters on 4-high mill (with WR bending and WR crown applied) showing (a) semi-thickness (top surface) and (b) reduction deviation. Both plots include disturbance from 350 to 400 ms due to 100% step-change in incoming thickness. Note that oscillation in the strip semi-thickness is evident, as highlighted by the inset image in (a).
Close modal
Fig. 17
Time history of normalized displacements on 6-high CVC mill: (a) upper WR axis displacement and (b) upper IR axis displacement. Also shown is the relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 175 to 225 ms due to 100% step-change in incoming thickness).
Fig. 17
Time history of normalized displacements on 6-high CVC mill: (a) upper WR axis displacement and (b) upper IR axis displacement. Also shown is the relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 175 to 225 ms due to 100% step-change in incoming thickness).
Close modal
Fig. 18
Time history of normalized contact force distribution: (a) upper WR–upper IR interface on 6-high CVC mill. Note that point “p” has zero contact force in the steady-state condition, implying no contact between the intermediate roll and work roll at this point. (b) Relative difference of contact force distribution to steady-state force in (a) from 175 to 225 ms. Notice point “p” from (a) is indicated in the magnified time scale. While point “p” was not in contact during steady-state, during the 100% incoming thickness disturbance, the point continuously establishes and loses contact (contact forces periodically change from non-zero to zero contact force. Also, note the high-frequency change in contact force during the disturbance.
Fig. 18
Time history of normalized contact force distribution: (a) upper WR–upper IR interface on 6-high CVC mill. Note that point “p” has zero contact force in the steady-state condition, implying no contact between the intermediate roll and work roll at this point. (b) Relative difference of contact force distribution to steady-state force in (a) from 175 to 225 ms. Notice point “p” from (a) is indicated in the magnified time scale. While point “p” was not in contact during steady-state, during the 100% incoming thickness disturbance, the point continuously establishes and loses contact (contact forces periodically change from non-zero to zero contact force. Also, note the high-frequency change in contact force during the disturbance.
Close modal
Fig. 19
Time history of normalized contact force distribution: (a) upper IR–upper BUR interface on 6-high CVC mill and (b) relative difference of contact force distribution to steady-state dynamic response during the disturbance period (from 175 to 225 ms due to 100% step-change in incoming thickness)
Fig. 19
Time history of normalized contact force distribution: (a) upper IR–upper BUR interface on 6-high CVC mill and (b) relative difference of contact force distribution to steady-state dynamic response during the disturbance period (from 175 to 225 ms due to 100% step-change in incoming thickness)
Close modal
Fig. 20
Time history of exit semi-thickness (upper half) of the strip on 6-high CVC mill. Plot includes disturbance from 175 to 225 ms due to 100% step-change in incoming thickness. Note that slight oscillation can again be seen during the thickness disturbance.
Fig. 20
Time history of exit semi-thickness (upper half) of the strip on 6-high CVC mill. Plot includes disturbance from 175 to 225 ms due to 100% step-change in incoming thickness. Note that slight oscillation can again be seen during the thickness disturbance.
Close modal
Fig. 21
Side view of 20-high cluster mill with roll labeling of the strip and upper roll stack [27] (Reprinted with permission from Elsevier © 2008).
Fig. 21
Side view of 20-high cluster mill with roll labeling of the strip and upper roll stack [27] (Reprinted with permission from Elsevier © 2008).
Close modal
Fig. 22
(a) Contact force distribution reproduced from static SM-FEM solution in Ref. [27] for the roll interfaces shown on a 20-high cluster mill (see Fig. 21) and (b) comparison of contact force distribution between the Ref. [27] and the presented dynamic steady-state solution
Fig. 22
(a) Contact force distribution reproduced from static SM-FEM solution in Ref. [27] for the roll interfaces shown on a 20-high cluster mill (see Fig. 21) and (b) comparison of contact force distribution between the Ref. [27] and the presented dynamic steady-state solution
Close modal
Fig. 23
(a) Contact force distribution reproduced from static SM-FEM solution in Ref. [27] for the roll interfaces shown on a 20-high cluster mill (see Fig. 21); (b) comparison of contact force distribution between Ref. [27] and the presented dynamic steady-state solution
Fig. 23
(a) Contact force distribution reproduced from static SM-FEM solution in Ref. [27] for the roll interfaces shown on a 20-high cluster mill (see Fig. 21); (b) comparison of contact force distribution between Ref. [27] and the presented dynamic steady-state solution
Close modal
Fig. 24
Time history of (a) vertical displacement of upper WR axis (roll 2 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Fig. 24
Time history of (a) vertical displacement of upper WR axis (roll 2 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Close modal
Fig. 25
Time history of (a) vertical displacement of upper first-intermediate roll (FIR) axis (roll 3 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Fig. 25
Time history of (a) vertical displacement of upper first-intermediate roll (FIR) axis (roll 3 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Close modal
Fig. 26
Time history of (a) vertical displacement of upper idler (IDLR) axis (roll 6 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Fig. 26
Time history of (a) vertical displacement of upper idler (IDLR) axis (roll 6 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Close modal
Fig. 27
Time history of (a) horizontal displacement of upper FIR axis (roll 3 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Fig. 27
Time history of (a) horizontal displacement of upper FIR axis (roll 3 in Fig. 21); (b) relative difference between the current displacement and the steady-state dynamic response during the disturbance period (from 700 to 750 ms due to 100% step-change in incoming thickness)
Close modal

Conclusion

An efficient and numerically stable new 3D nonlinear dynamic roll-stack model formulation is presented by coupling the simplified-mixed finite element method (SM-FEM) with a Newmark-beta direct time integration technique. The combined formulation offers both computational efficiency and mill-configuration flexibility in simulating time-history responses of the strip thickness profile, contact force distributions, and roll axes displacements due to dynamic disturbances. Incorporation of the strip modulus concept into the dynamic solution approach enables the realization of piecewise linearity in the time domain since nonlinearity is addressed within the SM-FEM model during each time-step. This coupling technique not only ensures the stability of the numerical solution, but also eliminates “spurious” vibrations typically associated with classical damping. Case studies involving 4-high, 6-high, and 20-high mills reveal the capability of the model to accommodate important mill configurations and profile/flatness control mechanisms while abandoning the limiting simplifications in prior dynamic modeling approaches. Based on the application of a strip thickness dynamic disturbance, results with the new model show that it can provide key insights into the dynamically changing contact conditions, as well as the oscillations that affect critical dimensional quality criteria such as transient deviations in the thickness profile and relative reduction.

Footnote

2

Use of the secant is preferable to the tangent in terms of nonlinear stability, particularly during dynamic analysis.

Acknowledgment

The authors gratefully acknowledge support for this work from the National Science Foundation (Grant No. CMMI-1555531).

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The authors attest that all data for this study are included in the paper.

References

1.
Zhang
,
F.
, and
Malik
,
A.
,
2021
, “
An Efficient Multi-scale Modeling Method That Reveals Coupled Effects Between Surface Roughness and Roll-Stack Deformation in Cold Sheet Rolling
,”
ASME J. Manuf. Sci. Eng.
,
143
(
10
), p.
101005
.
2.
Xie
,
L.
,
He
,
A.
, and
Liu
,
C.
,
2018
, “
A Rapid Calculation Method for Predicting Roll Deformation of Six-High Rolling Mill
,”
J. Iron Steel Res. Int.
,
25
(
9
), pp.
901
909
.
3.
Park
,
H.
, and
Hwang
,
S.
,
2017
, “
3-D Coupled Analysis of Deformation of the Strip and Rolls in Flat Rolling by FEM
,”
Steel Res. Int.
,
88
(
12
), p.
1700227
.
4.
Malik
,
A. S.
, and
Hinton
,
J. L.
,
2012
, “
Displacement of Multiple, Coupled Timoshenko Beams in Discontinuous Nonlinear Elastic Contact, With Application to Rolling Mills
,”
ASME J. Manuf. Sci. Eng.
,
134
(
5
), p.
051009
.
5.
Yarita
,
I.
,
Furukawa
,
K.
,
Seino
,
Y.
,
Takimoto
,
T.
,
Nakazato
,
Y.
, and
Nakagawa
,
K.
,
1978
, “
An Analysis of Chattering in Cold Rolling for Ultrathin Gauge Steel Strip
,”
Trans. Iron Steel Inst. Jpn.
,
18
(
1
), pp.
1
10
.
6.
Tamiya
,
T.
,
Furui
,
K.
, and
Lida
,
H.
,
1980
, “
Analysis of Chattering Phenomenon in Cold Rolling
,”
Proceedings of Science and Technology of Flat Rolled Products, International Conference on Steel Rolling
,
Tokyo, Japan
,
Sept. 29–Oct. 4
, pp.
1191
1202
.
7.
Tlusty
,
J.
,
Chandra
,
G.
,
Critchley
,
S.
, and
Paton
,
D.
,
1982
, “
Chatter in Cold Rolling
,”
CIRP Ann.
,
31
(
1
), pp.
195
199
.
8.
Chefneux
,
L.
,
Fischbach
,
J.-P.
, and
Gouzou
,
J.
,
1984
, “
Study and Industrial Control of Chatter in Cold Rolling
,”
Iron Steel Eng.
,
61
(
11
), pp.
17
26
.
9.
Kapil
,
S.
,
Eberhard
,
P.
, and
Dwivedy
,
S. K.
,
2014
, “
Nonlinear Dynamic Analysis of a Parametrically Excited Cold Rolling Mill
,”
ASME J. Manuf. Sci. Eng.
,
136
(
4
), p.
041012
.
10.
Lin
,
Y.-J.
,
Suh
,
C. S.
,
Langari
,
R.
, and
Noah
,
S. T.
,
2003
, “
On the Characteristics and Mechanism of Rolling Instability and Chatter
,”
ASME J. Manuf. Sci. Eng.
,
125
(
4
), pp.
778
786
.
11.
Kimura
,
Y.
,
Sodani
,
Y.
,
Nishiura
,
N.
,
Ikeuchi
,
N.
, and
Mihara
,
Y.
,
2003
, “
Analysis of Chatter in Tandem Cold Rolling Mills
,”
ISIJ Int.
,
43
(
1
), pp.
77
84
.
12.
Yun
,
I.-S.
,
Wilson
,
W. R. D.
, and
Ehmann
,
K. F.
,
1998
, “
Chatter in the Strip Rolling Process, Part 1: Dynamic Model of Rolling
,”
ASME J. Manuf. Sci. Eng.
,
120
(
2
), pp.
330
336
.
13.
Yun
,
I.-S.
,
Ehmann
,
K. F.
, and
Wilson
,
W. R. D.
,
1998
, “
Chatter in the Strip Rolling Process, Part 2: Dynamic Rolling Experiments
,”
ASME J. Manuf. Sci. Eng.
,
120
(
2
), pp.
337
342
.
14.
Yun
,
I.-S.
,
Ehmann
,
K. F.
, and
Wilson
,
W. R. D.
,
1998
, “
Chatter in the Strip Rolling Process, Part 3: Chatter Model
,”
ASME J. Manuf. Sci. Eng.
,
120
(
2
), pp.
343
348
.
15.
Hu
,
P.-H.
, and
Ehmann
,
K. F.
,
2000
, “
A Dynamic Model of the Rolling Process. Part I: Homogeneous Model
,”
Int. J. Mach. Tools Manuf.
,
40
(
1
), pp.
1
19
.
16.
Hu
,
P.-H.
, and
Ehmann
,
K. F.
,
2000
, “
A Dynamic Model of the Rolling Process. Part II: Inhomogeneous Model
,”
Int. J. Mach. Tools Manuf.
,
40
(
1
), pp.
21
31
.
17.
Kim
,
Y.
,
Kim
,
C.-W.
,
Lee
,
S.-J.
, and
Park
,
H.
,
2012
, “
Experimental and Numerical Investigation of the Vibration Characteristics in a Cold Rolling Mill Using Multibody Dynamics
,”
ISIJ Int.
,
52
(
11
), pp.
2042
2047
.
18.
Lu
,
X.
,
Sun
,
J.
,
Li
,
G.
,
Wang
,
Q.
, and
Zhang
,
D.
,
2019
, “
Dynamic Analysis of Vibration Stability in Tandem Cold Rolling Mill
,”
J. Mater. Process. Technol.
,
272
, pp.
47
57
.
19.
Mehrabi
,
R.
,
Salimi
,
M.
, and
Ziaei-Rad
,
S.
,
2015
, “
Finite Element Analysis on Chattering in Cold Rolling and Comparison With Experimental Results
,”
ASME J. Manuf. Sci. Eng.
,
137
(
6
), p.
061013
.
20.
Niroomand
,
M. R.
,
Forouzan
,
M. R.
, and
Salimi
,
M.
,
2015
, “
Theoretical and Experimental Analysis of Chatter in Tandem Cold Rolling Mills Based on Wave Propagation Theory
,”
ISIJ Int.
,
55
(
3
), pp.
637
646
.
21.
Brusa
,
E.
,
Lemma
,
L.
, and
Benasciutti
,
D.
,
2010
, “
Vibration Analysis of a Sendzimir Cold Rolling Mill and Bearing Fault Detection
,”
Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci.
,
224
(
8
), pp.
1645
1654
.
22.
Sun
,
J.
,
Peng
,
Y.
, and
Liu
,
H.
,
2014
, “
Dynamic Characteristics of Cold Rolling Mill and Strip Based on Flatness and Thickness Control in Rolling Process
,”
J. Cent. South Univ.
,
21
(
2
), pp.
567
576
.
23.
Kapil
,
S.
,
Eberhard
,
P.
, and
Dwivedy
,
S. K.
,
2016
, “
Dynamic Analysis of Cold-Rolling Process Using the Finite-Element Method
,”
ASME J. Manuf. Sci. Eng.
,
138
(
4
), p.
041002
.
24.
Brusa
,
E.
, and
Lemma
,
L.
,
2009
, “
Numerical and Experimental Analysis of the Dynamic Effects in Compact Cluster Mills for Cold Rolling
,”
J. Mater. Process. Technol.
,
209
(
5
), pp.
2436
2445
.
25.
Yun
,
I. S.
,
Wilson
,
W. R. D.
, and
Ehmann
,
K. F.
,
1998
, “
Review of Chatter Studies in Cold Rolling
,”
Int. J. Mach. Tools Manuf.
,
38
(
12
), pp.
1499
1530
.
26.
Guo
,
R.-M.
,
Urso
,
A. C.
, and
Schunk
,
J. H.
,
1993
, “
Analysis of Chatter Vibration Phenomena of Rolling Mills Using Finite Element Methods
,”
Iron Steel Eng.
,
70
(
1
), pp.
29
39
.
27.
Malik
,
A. S.
, and
Grandhi
,
R. V.
,
2008
, “
A Computational Method to Predict Strip Profile in Rolling Mills
,”
J. Mater. Process. Technol.
,
206
(
1–3
), pp.
263
274
.
28.
Matsubara
,
S.
,
Hara
,
K.
, and
Takezoe
,
A.
,
1989
, “
Optimization of Work Roll Taper for Extremely-Thin Strip Rolling
,”
ISIJ Int.
,
29
(
1
), pp.
58
63
.
29.
Roberts
,
W. L.
, “
A Simplified Cold Rolling Model
,”
Iron Steel Eng.
,
42
(
10
), pp.
75
87
.
30.
Fleck
,
N. A.
, and
Johnson
,
K. L.
,
1987
, “
Towards a New Theory of Cold Rolling Thin Foil
,”
Int. J. Mech. Sci.
,
29
(
7
), pp.
507
524
.
31.
Bland
,
D. R.
, and
Ford
,
H.
,
1948
, “
The Calculation of Roll Force and Torque in Cold Strip Rolling With Tensions
,”
Proc. Inst. Mech. Eng.
,
159
(
1
), pp.
144
163
.
32.
Zhang
,
F.
, and
Malik
,
A.
,
2018
, “
A Roll-Stack Contact Mechanics Model to Predict Strip Profile in Rolling Mills With Asymmetric, Continuously Variable Crown Rolls
,”
ASME J. Manuf. Sci. Eng.
,
140
(
1
), p.
011008
.
33.
Zhang
,
F.
,
Malik
,
A. S.
, and
Yu
,
H.
, “
High-Fidelity Roll Profile Contact Modeling by Simplified Mixed Finite Element Method
,”
Proceedings of the ASME 2018 13th International Manufacturing Science and Engineering Conference, Vol. 4: Processes
,
College Station, TX
,
June 18–22
, p.
V004T03A034
.
34.
Guo
,
R.-M.
,
1997
, “
Prediction of Strip Profile in Rolling Process Using Influence Coefficients and Boussinesq’s Equations
,”
ASME J. Manuf. Sci. Eng.
,
119
(
5
), pp.
220
226
.
35.
Guo
,
R.-M.
, and
Malik
,
A. S.
,
2005
, “
Development of a New Crown/Shape Control Model for Cluster Mills
,”
Iron Steel Technol.
,
2
(
8
), pp.
31
40
.
36.
Lee
,
S.-H.
,
Song
,
G.-H.
,
Lee
,
S.-J.
, and
Kim
,
B.-M.
,
2011
, “
Study on the Improved Accuracy of Strip Profile Using Numerical Formula Model in Continuous Cold Rolling With 6-High Mill
,”
J. Mat. Sci. Tech.
,
25
(
8
), pp.
2101
2109
.
37.
Lindfield
,
G.
, and
Penny
,
J.
,
2019
,
Numerical Methods
, 4th ed.,
Academic Press
,
London, UK
, pp.
239
299
.
38.
Cook
,
R. D.
,
2007
,
Concepts and Applications of Finite Element Analysis
,
John Wiley & Sons
,
Hoboken, NJ
.
39.
Lu
,
Y.
, and
Morris
,
G.
,
2017
, “
Assessment of Three Viscous Damping Methods for Nonlinear History Analysis: Rayleigh With Initial Stiffness, Rayleigh With Tangent Stiffness, and Modal
,”
Proceedings of Sixteenth World Conference on Earthquake Engineering
,
Santiago, Chile
,
Jan. 8–13
.
40.
Paton
,
S.
, and
Critchley
,
D. L.
,
1985
, “
Tandem Mill Vibration: Its Cause and Control
,”
Iron Steelmak.
,
12
(
3
), pp.
37
43
.