Abstract
Introduced is a new physicsbased threedimensional (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 rollstack deformations due to dynamic disturbances. Formulation of the new 3D model is achieved through a combination of the highly efficient simplifiedmixed finite element method with a Newmarkbeta direct time integration approach to solve the system of differential equations that governs the motion of the rollstack. 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 clustertype 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 4high mill equipped with both workroll bending and workroll crown, a 6high mill with continuously variable crown (CVC) intermediate rolls, and finally a complex 20high cluster mill. Solution times on a single computing processor for the damped 4high and 20high case studies are just 0.37 s and 3.38 s per timestep, respectively.
Introduction
Dimensional quality of coldrolled 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 lightweight 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 rollstack 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 steadystate 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 timehistory responses of the 3D structural dynamics of a rollstack 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 selfexcited 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 quasistatic (steadystate) 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 springmassdamper 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.
Method/model  Investigator  Assumptions  Comments on model characteristics 

Linear lumped parameter: Unimodal  Yarita 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] 


Linear lumped parameter: Multimodal  Yun et al. 1998 [12–14] Hu et al. 2000 [15,16] 


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

 
Brusa et al. 2010 [21] 

 
3D model  Sun et al. 2014 [22] Kapil et al. 2016 [23] Brusa et al. 2009 [24] 


Method/model  Investigator  Assumptions  Comments on model characteristics 

Linear lumped parameter: Unimodal  Yarita 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] 


Linear lumped parameter: Multimodal  Yun et al. 1998 [12–14] Hu et al. 2000 [15,16] 


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

 
Brusa et al. 2010 [21] 

 
3D model  Sun et al. 2014 [22] Kapil et al. 2016 [23] Brusa et al. 2009 [24] 


It is important to note, however, that all linear lumped parameter systems inherently neglect 3D bulkbody deformation behavior, particularly the transverse (widthwise) distributions of bending, nonuniform Hertzian contact flattening, and sheartype rollstack deformations. Indeed, neither unimodal nor multimodel linear lumped systems allow for analysis of the effects of bending and torsional modes on the transient rollstack deformation behavior [26]. Since these “non rollstack” linear lumped model formulations inherently lack 3D bulkbody deformations, vital information is missed regarding the timevarying 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., 20high 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 continuumelementbased 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 4high mill by combining (along the transverse/axial direction) a rolling process model (of the rollbite contact mechanics) with a structural dynamics model. Kapil et al. [23] also developed a 4high 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 rollstack models is the requirement for continuous and constant contact between interacting bodies (i.e., roll–roll and roll–strip contacts). However, since nonuniform 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 continuumelementbased FE formulations, they lack ability to accommodate typical strip profile/flatness control mechanisms, such as workroll 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., 4high). Indeed, even models with some degree of flexibility are still limited to vertically oriented mills such as 2high, 4high, or 6high, i.e., inplane rollstack configurations. Models needed to investigate complex cluster mills require multiple outofplane contact interactions and have seldom been proposed. While Brusa et al. [21] investigated the dynamic behavior of a 20high cluster mill (Zmill), the structural dynamics are represented using a simplified linear lumped parameter system. In an effort toward the design of 20high cluster mills, Brusa and Lemma [24] used a multibody dynamics approach, but assumptions of symmetry and rigid bodies restrict accuracy for any timehistory analyses of the rollstack 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 nonexistent 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 rollbite process parameters under an assumption of plane strain. As a result, the models are not sufficiently formulated to provide realworld time history of the 3D structural dynamics—particularly the transient strip dimensional quality.
Accordingly, this work introduces and demonstrates a generalpurpose, nonlinear 3D dynamic rollstack model formulation that gives an efficient timehistory analysis in cold rolling, including the transient predictions of strip thickness profile, contact force distributions, and rollstack deformations. As is shown, the model is applicable to accommodate various types of flatness control mechanisms and mill configurations, including 20high cluster mills. Note that while chatter prediction with multistand 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 simplifiedmixed finite element method (SMFEM) rollstack 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 ThreeDimensional RollStack Model.
An overview of the main equations to formulate the 3D static rollstack deformation model using the SMFEM 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 crosssection 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. Twinbeamcoupled Winker foundation elements are used to represent nonlinear Hertztype 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 rollbite 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 SMFEM method in both computational validation with largescale continuum FEM and experimental validation against industry data [27,32,33].
The interroll angle of inclination, θ, between bodies 1 and 2 in Eqs. (5)–(8) accounts for both vertical and cluster mill configurations ($\theta =90deg$ for vertical rollstacks as in 2high, 4high, and 6high mills, whereas clustertype 20high mills involve several unique interroll 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 6high 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 timestep 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 timeintegrated solution.
Incorporation of Flatness Control Mechanisms and NonUniform 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 nonuniform roll profiles such as parabolic crowns, tapered profiles, advanced CVC profiles, as well as highfidelity deviations such as roll grinding errors and roll roughness is readily adapted into the general contact interference algorithm as follows:
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 workroll 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 wellknown “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}/2d^{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 rollstack deformation model has been validated against both largescale continuum FE models and industry measurements of strip thickness profile [32,33]. Described next is the adaptation of this global stiffnessbased static SMFEM method into a Newmark implicit time integration approach to obtain a highly efficient, 3D, nonlinear, generalpurpose dynamic rollstack model.
Dynamic Solution in Time.
In Eq. (14), matrix [K] is written concisely here, but it is actually a function of the displacement vector, u, and corresponds to [K_{G}(u)] in Eq. (1) for the nonlinear, static SMFEM model. Similarly, f(t) in Eq. (17) is actually f(u, t) and is obtained from f(u) in Eq. (1).
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, NMbeta is unconditionally stable with secondorder 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 higherorder, nonphysical modes of vibration that need to be suppressed. However, accuracy is then reduced to firstorder [38]. One of the advantages of the formulation in this work is that the static equilibrium achieved by the iterative solution during each timestep (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.
Mill Housing Stiffness Representation.
Note that the static SMFEM rollstack deformation model assumes a quasistatic condition wherein modeling of the mill housing is not required to accurately obtain steadystate 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 inplane rollstack configurations (i.e., 2high, 4high, and 6high 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” forcetype 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” displacementtype 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 SMFEM and NMbeta formulation to address the research gaps in the context of achieving a generalpurpose 3D dynamic rollstack 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 rollstack model using the NMbeta time integration technique are discussed for a 4high mill.
Demonstration of model compatibility with conventional flatness control mechanisms (workroll bending and workroll crown) on a 4high mill.
Model demonstration for 6high mill involving asymmetric, CVCmachined profiles on the intermediate rolls.
Model demonstration for a complex 20high cluster mill.
Figure 6(a) depicts the 4high 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 6high mill for case (3), while a side view of the 20high 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 4high rolling mill configuration is adapted here to demonstrate the 3D structural dynamics solution, i.e., timehistory response. Since the NMbeta 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, steadystate (i.e., longterm) results for the damped case are compared to results of the original (previously validated) static SMFEM formulation in order to validate the 3D dynamic formulation. Mill dimensions and parameters for these 4high 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.
Material  Yield strength (MPa)  Strength coefficient (MPa)  Strain hardening exponent 

SS 301  283.37  2987.49  0.7475 
Material  Yield strength (MPa)  Strength coefficient (MPa)  Strain hardening exponent 

SS 301  283.37  2987.49  0.7475 
Parameter  Value 

Work roll diameter  76.2 mm 
Backup roll diameter  304.8 mm 
Work roll face length  305 mm 
Backup roll face length  305 mm 
Work roll neck diameter  50.8 mm 
Backup roll neck diameter  187.2 mm 
Workroll neck length  152 mm 
Backup roll neck length  152 mm 
Distance between bearings on work roll  457 mm 
Distance between bearings on backup roll  457 mm 
Workroll and backup roll Young's modulus  207 GPa 
Strip width  209.6 mm 
Entry thickness  2.576 mm 
Exit thickness  2.382 mm 
Reduction ratio  7.495% 
Strip modulus  10.14 GPa 
Parameter  Value 

Work roll diameter  76.2 mm 
Backup roll diameter  304.8 mm 
Work roll face length  305 mm 
Backup roll face length  305 mm 
Work roll neck diameter  50.8 mm 
Backup roll neck diameter  187.2 mm 
Workroll neck length  152 mm 
Backup roll neck length  152 mm 
Distance between bearings on work roll  457 mm 
Distance between bearings on backup roll  457 mm 
Workroll and backup roll Young's modulus  207 GPa 
Strip width  209.6 mm 
Entry thickness  2.576 mm 
Exit thickness  2.382 mm 
Reduction ratio  7.495% 
Strip modulus  10.14 GPa 
Spatial and Temporal Convergence.
As mentioned, even though NMBeta is unconditionally stable with average acceleration, inadequate timestep sizes can result in inaccurate solutions, as evident in Fig. 7(a) for a timestep 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 timestep 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 timestep 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 timestep size, the case studies that follow incorporate a 50 − μs timestep. Also, in finitedifference solutions to timedependent, spatially discretized boundary value problems, the required timestep 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 timestep 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 timehistories of the axes of the upper work roll (WR) and upper backup roll (BUR) on the 4high 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 highfrequency 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 timestep.
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 highfrequency 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, timevarying reduction across the strip width as a result of the 3D bulk body deformation and flattening of the rollstack, is shown in Fig. 12(a).
Figure 12(b) shows the timehistory 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 timestep.
Comparison of SteadyState Damped Dynamic Response to Static Model Results.
Comparison of the displacement and contact force results for the steadystate solution of the 3D dynamic response against the original static SMFEM 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 timehistory responses are now provided to demonstrate compatibility of the 3D dynamic rollstack model with conventional flatness control mechanisms, including roll bending and roll crown, on the same 4high 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 steadystate 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 steadystate 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 timestep.
 (3)
6High 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 6high 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 SMFEM model to incorporate CVC profiles was validated. Roll shifting with CVC profiles, while readily applicable, is not included for the demonstration of the timehistory response. As in the previous case, a post steadystate disturbance is introduced (from 175 to 225 ms in this case) via a 100% stepchange 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 steadystate. 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 steadystate 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)
20High Cluster Mill
As mentioned, almost all existing model formulations to investigate rolling mill structural dynamics are restricted to vertical mill configurations, i.e., 2high, 4high, or 6high. Clustertype mills, such as the 20high 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 20high cluster mill is again examined for a disturbance defined by 100% stepchange 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 SMFEM method, which was validated against largescale 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 midplane. The timehistory 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 20high 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 steadystate solution and static SMFEM 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 noniterative version of the SMFEM in Ref. [27], is corrected by NewtonRaphson iteration in the presented formulation and shows zero force (i.e., no contact).
Damped timehistory results for the incoming strip thickness stepchange disturbance on the 20high mill are demonstrated in Figs. 24–27. The average run time (1 CPU, matlab) for the 20high mill is 3.83 s per timestep.
Conclusion
An efficient and numerically stable new 3D nonlinear dynamic rollstack model formulation is presented by coupling the simplifiedmixed finite element method (SMFEM) with a Newmarkbeta direct time integration technique. The combined formulation offers both computational efficiency and millconfiguration flexibility in simulating timehistory 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 SMFEM model during each timestep. 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 4high, 6high, and 20high 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. CMMI1555531).
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.