## 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 [1–4], 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.

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.

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 [29–31]. 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].

*u*_{j}, at node

*j*include the 3D translation/rotation degrees-of-freedom (DOF):

*K*_{G}], the global stiffness in Eq. (1), is the superposition of all nonlinear continuous foundation stiffness matrices, [

*K*_{F}], and the combined Timoshenko matrices, [

*K*_{T}], with the latter containing bending and shear stiffness:

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

*a*)

*b*)

*c*)

*l*

_{i}in Eq. 4(

*b*) is the length of the

*i*th 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)

The inter-roll angle of inclination, *θ*, between bodies 1 and 2 in Eqs. (5)–(8) accounts for both vertical and cluster mill configurations ($\theta =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 secant^{2} of the rolling load or plastic strain relationship at the working point [27,34,35] (see Fig. 4).

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:

*d*

_{12}between their axes centers at axial location

*x*is obtained from their relative displacement,

*v*

_{1}(

*x*) −

*v*

_{2}(

*x*), as well as the respective initial vertical positions,

*y*

_{c}

*D*

_{1}(

*x*) and

*D*

_{2}(

*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,

*P*

_{C}, is required to equilibrate the integral of specific contact force

*f*(

*x*), based on either the Matsubara (see Eq. (14)) or Hertz/Föppl [36] plane strain relations, where

*E*

_{1}, ν

_{1}and

*E*

_{2}, ν

_{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

*D*

_{1}→ ∞ and

*E*

_{1}→ ∞, then solving for $kf2(x)$ per Eq. (15), then obtaining $kf1(x)$ from Eq. (16)

*a*)

*b*)

*a*)

*b*)

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 (|*x* − *x*_{0}|^{2}/2*d*^{2} + 1) is applied to the strip modulus for |*x*| ≥ *x*_{0}, where −*w*/2 ≤ *x* ≤ *w*/2 and *x*_{0} = *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.

**is the displacement vector, [**

*u***], [**

*M***], and [**

*C***] are mass, damping, and stiffness matrices, respectively,**

*K***is the load vector, and**

*f**t*is time

In Eq. (14), matrix [** K**] is written concisely here, but it is actually a function of the displacement vector,

**, and corresponds to [**

*u*

*K*_{G}(

**)] in Eq. (1) for the nonlinear, static SM-FEM model. Similarly,**

*u***(**

*f**t*) in Eq. (17) is actually

**(**

*f***,**

*u**t*) and is obtained from

**(**

*f***) in Eq. (1).**

*u**γ*≤ 1 and 0 ≤

*β*≤ 1/2, are based on the extended mean value theorem

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.

*M*_{i}] and [

*K*_{i}] are initial global mass matrix (using HRZ lumping) and global stiffness matrix, respectively, and

*c*

_{1}and

*c*

_{2}are coefficients of the mass and stiffness proportionalities given by Eq. (25)

*ω*

_{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 [

**] for demonstration purposes with the case studies in the next section.**

*C*#### 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).

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 × 10^{9} 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:

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.

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

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

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.

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.

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

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

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.

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

*x*is per Eqs. (26) and (27) where $r\xaf$ is the average reduction, and Δ

*r*(

*x*) is the deviation from the average reduction

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. 14–16, 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*C*_{r}(*x*) in Eq. (28) is applied to the intermediate rolls (see again Fig. 6(d))(28)$Cr(x)=\u22120.01525x+0.0167318x2+0.0440634x3$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. 24–27. The average run time (1 CPU, matlab) for the 20-high mill is 3.83 s per time-step.

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

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

*Sendzimir*Cold Rolling Mill and Bearing Fault Detection