Finite Difference Modeling of the Interstand Evolutions of Profile and Residual Stress during Hot Strip Rolling

Elastic recovery and viscoplastic stress relaxation occur in the interstand of hot rolling, impacting the evolutions of strip profile and residual stress, which are major concerns for obtaining high-quality flat products. A better understanding of the evolutionary mechanisms would help develop shape control strategies. Therefore, a quasi-3D steady-state elasto-viscoplastic rolling model is developed based on the finite difference method. Predictions of spread, profile, and residual stress are validated through comparisons with a two-stand finite element model. The new model is also complemented with a roll stack model and with a viscoplastic constitutive model calibrated by hot compression tests to simulate a seven-stand hot rolling industrial experiment with low carbon steel. Comparisons between the predicted and measured profiles show a satisfactory accuracy. The simulation costs approximately a minute of CPU time, enabling the new model to run massive parametric campaigns for process optimization. It is found that during the interstand elastic recovery, the transverse compressive stress releases and the strip velocity tends to be uniform, revealing residual stress after a significant change of stress pattern. The stress relaxation mainly occurs at the edge near the roll bite and therefore increases the edge drop of the profile; it also decreases the center crown by changing the distribution of the rolling pressure in the roll bite.


Introduction
In hot strip rolling, the thickness profile and the residual stress (which leads to flatness defects) across the strip width are major quality concerns that adhere to increasingly tight tolerances. Excessive residual stress may cause strip breakages and limit the rolling speed and productivity; profile and flatness defects make subsequent processing difficult and degrade product performance. Profile and residual stress are usually attributed to the transverse variations in reduction and elongation in the roll bite. Moreover, two phenomena occurring in the interstand of a hot tandem mill can reshape profile and residual stress: one is the elastic recovery after the compression in the roll bite, and the other is the viscoplastic stress relaxation under high temperature. Therefore, to predict the entire evolution of the profile and residual stress during rolling, a reliable modeling strategy considering these interstand phenomena is needed. It would help to understand the evolutionary mechanisms and to develop corrective measures for the defects. However, most modeling efforts have been devoted to the deformation in the bite, as pointed out by Montmitonnet et al. [1] in a review.
Usually, a strip model provides the rolling pressure distribution and is coupled with a roll stack model to calculate the thickness reduction in the vertical direction. It also predicts the elongation and the residual stress in the rolling direction. However, metal flow in the transverse direction (spread) is active in hot rolling and has a substantial influence on the transverse variations in reduction and elongation because the metal flows in all three directions are highly coupled. Therefore, plane strain models that ignore the spread cannot give a reliable prediction of profile and residual stress in hot rolling; a 3D model is necessary.
The finite element method (FEM) is widely used in a complex 3D rolling analysis because it has the advantage of being physically rich. In addition, some studies have noticed the effects of the interstand (pre-bite or post-bite) deformation. Kim et al. [2] developed an FEM model to predict strip profile and reported that insufficient lengths of pre-bite or post-bite zones in the model leads to an erroneous prediction. Moazeni and Salimi [3] studied the effect of nonuniform reduction on residual stresses using a commercial finite element package (ABAQUS/explicit). Similarly, they suggested that the simulation regime should extend far enough from the roll bite to obtain a steady solution for residual stress. Hacquin et al. [4] built a steady-state FEM model and found that residual stress significantly reduces after changing the material behavior from elastic-plastic to elasto-viscoplastic (i.e., introducing stress relaxation). Based on an FEM hot rolling simulation of thick plate, Zaepffel et al. [5] showed that the viscoplastic deformation could significantly modify the profile in the post-bite region. Legrand et al. [6] predicted strip width variation with an FEM model and highlighted the increase of width due to post-bite elastic recovery. The above observations suggest that the interstand deformation is not negligible and there is a strong interaction between the interstand and the roll bite.
The computational costs of FEM models are generally high. Compared to the standard incremental Lagrangian formulation, the steady-state Eulerian and arbitrary Lagrangian-Eulerian techniques in rolling modeling can save CPU time, as reviewed by Fourment et al. [7]. However, as the interstand regions in multi-stand hot rolling are to be modeled, the total computation is too heavy for the FEM even with current computational power. Therefore, an efficient modeling strategy is needed. A fast model can not only analyze tandem rolling but also run massive parametric campaigns for industrial process optimization.
Different methods have been used to develop fast models for the roll bite deformation, which are the basis of the interstand analysis. They usually neglect elasticity or simplify variations in the vertical direction (often referred to as 'quasi-3D'). Johnson [8] reduced the governing equations for thin strip rolling by an asymptotic analysis and derived an analytical solution of spread. Dixon and Yuen [9] developed a real-time model based on the asymptotic techniques to predict the tension and transverse strain at the exit of the roll bite. Liu et al. [10] obtained the stress fields using a hybrid model of the upper bound method and finite difference method (FDM), assuming the transverse flow had a polynomial pattern along the rolling direction. Ngo [11] built an online width variation model based on the upper bound method with the assumed polynomial velocity field. Tozawa et al. [12] solved the stress field and metal flow problems based on the FDM, assuming that stresses are constant through the thickness; however, the convergence of the solution is sensitive to the strip width and the initial guess of the spread. Later, they simplified the model by introducing a linear relationship between the transverse and longitudinal strains for application to a wider strip [13]. Yao et al. [14] proposed an FDM model with a global iterative solution to predict the distributions of spread, tension, and rolling pressure; the model is robust for typical wide strip hot rolling without the transverse strain simplification.
Some efforts have been made to develop a fast interstand model. The conventional ribbon model proposed by Shohet and Townsend [15] is widely used as a complement to the roll bite model; the post-bite strip is treated as a system of elastic narrow ribbons interconnected at the ends to estimate the residual stress caused by nonuniform elongation. However, the stress in a ribbon segment is assumed to be constant. Lee et al. [16] modeled the plastic yielding of the ribbon segment under excessive stress to predict the corresponding changes of profile and residual stress in the interstand.
Cresdee et al. [17] introduced a viscoplastic creep model without yield stress into the ribbon segment and divided the interstand into several zones along the rolling direction where the increments of creep strain and the resulting stress relaxation were calculated. Milenin et al. [18] simulated the stress relaxation in a similar way to study the residual stress of hot rolled strips during laminar cooling. Nevertheless, only longitudinal stress may occur in these ribbon models, neglecting transverse stress and shear stress. Domanti et al. [19] derived a purely elastic plane stress formulation for a moving strip in the interstand and found an analytical solution for the stress decay with assumed boundary conditions. If coupled with a roll bite model, this model has the potential to predict the residual stress evolution under elastic recovery. Kim et al. [20] considered the pre-bite region in a reduced model with parameters calibrated by the FEM and suggested that a small amount of pre-bite plastic deformation may affect the distributions of tension and rolling pressure. In the online shape control system of hot rolling, Zhao et al. [21] calculated the interstand change of shape and crown with an empirical coefficient for stress relaxation; the value of the coefficient relies on experimental data or simulation results.
Generally, in fast models, the roll bite module neglects the elasticity and therefore cannot provide reliable elastic and plastic strain for the residual stress calculation. The interstand module is not sophisticated enough to describe both elastic recovery and stress relaxation. Moreover, the interaction between the two modules is insufficiently modeled, neglecting the influence of the postbite deformation on the roll bite.
The lack of suitable models limits investigations on the interstand phenomena. Therefore, the aim of the present research is: to develop a fast model that can predict the interstand evolutions of profile and residual stress and that can handle the heavy simulations of hot tandem rolling with different configurations; to analyze the evolutionary mechanisms under the effects of interstand elastic recovery and stress relaxation, and to examine how these effects change from the upstream to the downstream stands.
In Section 2, to reduce the calculation complexity while preserving the essential physics of the problem, the governing equations of elastic-viscoplastic rolling deformation are determined based on the quasi-3D approximation; for efficiency and robustness, a global iterative solution with linearization is carried out based on the FDM. In Section 3, simplified tandem rolling conditions (two stands) are simulated by both the new strip model and an FEM model for validation, and the effects of elastic recovery are studied. In Section 4, the new strip model is complemented with both a roll stack model and a calibrated material model to simulate an industrial experiment of hot rolling (seven stands), and the predicted and measured profiles are compared; the effects of stress relaxation on the evolutions of profile and residual stress are analyzed. Discussion of the results is in Section 5.

Basic Assumptions
We assume that the rolling deformation has reached the steady state and is symmetric with respect to planes Oxz and Oxy. As shown in Figure 1, x, y, z are the coordinates along the rolling direction (RD), the transverse direction (TD), and the vertical direction (VD), respectively; O is the origin at the center of the cross-section at the roll bite exit. We assume that the strip velocity becomes uniform in the middle of the interstand (i.e., no new deformation occurs). This assumption is based on the observation that elastic recovery and stress relaxation in hot rolling cause rapid stress decay out of the roll bite so that the residual stress in the middle of the interstand is low and stable; correspondingly, viscoplastic deformation under low stress is slow enough to be ignored. As a result, the tandem rolling can be decoupled. A simulation is executed in each stand successively by passing the strip velocity, profile, and residual stress in the middle of the interstand to the downstream stand, as shown in Figure 1a. Further discussion on this assumption is in Section 5.
The calculation domain of a single stand is divided into five zones, as shown in Figure 1b: the upstream interstand zone (A), the entrance loading zone (B), the central roll bite zone (C), the exit unloading zone (D), and the downstream interstand zone (E). Zones B and D are characterized by sharp gradients of pressure: the transformations between the low-pressure interstand and the highpressure roll bite take place in short distances. A very fine mesh is required in these two zones to trace the sharp gradients.
In the roll bite (Zones B, C, and D), we introduce the quasi-3D approximation of strip rolling based on the asymptotic analysis by Johnson [8]. When the ratio of the half entry thickness to the roll bite length is small (which is the case in a typical hot rolling finishing mill), the strip velocity components vx, vy, the strain rates x ε , y ε , z ε , xy γ , and the stresses σx, σy, σz, τxy can all be considered constant along VD, and the stresses τxz and τyz can be considered to have a linear distribution along VD. This approximation is extended out of the roll bite, making the interstand deformation (Zones A and E) a 2D plane stress problem (σz vanishes when the strip is free from the contact pressure at the strip/roll interface). Furthermore, we assume that the strip remains flat (i.e., no bucking occurs) in the interstand to study the latent residual stress.

Friction and Force Equilibrium Equations
The velocity components of metal flow are shown in Figure 1c with a streamline trajectory. The relative slide velocity of the strip with respect to the roll surface, vs, can also be expressed by its components [10] h h y ′ = ∂ ∂ ) and vr is the line speed of the work roll. Then, based on the Coulomb model, the components of the frictional stress are given by where ks is the shear yield stress, p the rolling pressure (contact pressure at the strip/roll interface), and ( ) g s v the regularization function to eliminate the discontinuity of the frictional stress at the neutral point where the regularization factor η is set to 0.01 (more details in [14]). Based on the analysis of an element in the roll bite with infinitesimal length and width, as shown in Figure 2, equilibrium equations can be obtained as

Constitutive Relationship
The strain rates are given by ε ε ε ε ε ε γ γ γ For the elastic constitutive behavior, the objective stress rate is simplified to the material time derivative of the stress. This is because the rotation rate of the material is close to zero in most of the calculation domain except for the roll bite. However, the roll bite is dominated by plastic deformation. The material time derivative of the stress is given by where approximation is also made by ignoring the item with vy because the value of vy is generally two or three orders of magnitude smaller than vx in typical hot rolling. Then, based on Hooke's law, we have where E is Young's modulus, G the shear modulus, Kv the bulk modulus, ν Poisson's ratio, v ε the volumetric strain rate, σm the hydrostatic pressure, σm = (σx + σy + σz)/3. For the plastic behavior, we use the Sellars and Tegart viscoplastic constitutive relationship [22] ( ) ( ) where A, α, n, and Q are material constants, R the gas constant, σ and ε  the equivalent stress (or flow stress) and the equivalent strain rate respectively, where ij σ ′ is the stress deviator. In addition, the Levy-Mises flow rule is introduced as where λ is the plastic multiplier.

Solution
The above nonlinear governing equations are solved by an iterative method: the calculation domain is meshed; the governing equations are linearized based on an initial guess of the solution; then the linearized equations are discretized by the finite difference method (FDM) to build a global system of equations which is solved iteratively by updating the initial guess and the coefficients of the linearized equations, and the exact solution of the original nonlinear equations is approached.

Meshing and Initialization
As shown in Figure 3a, the calculation domain is meshed along the rolling and transverse direction, and the mesh is refined near the strip edge and the entrance and exit of the roll bite where large gradients of stress exist. The illustration exaggerates the spread, but the elements are virtually rectangular because the spread is two or three orders of magnitude smaller than the strip width in the finishing hot rolling. As shown in Figure 3b, nodes (solid points) are set at the midpoints of the transverse edges of the elements.
The initial guess of the solution is based on the conventional plane strain theory for the roll bite deformation (no spread or vy = 0), including roll speed, lengths of Zones B and D (more details in Appendix A), velocity, strain, and stress fields.

Linearization of Governing Equations
Based on the initial guess of the solution, the nonlinear terms of the governing equations can be linearized by substituting some variables with their initial value. This operation is based on the observation that: during the iteration, some variables such as vx and λ change less than the other variables in the nonlinear terms, so they can be fixed to serve as coefficients, remaining the essential relationship in the equation. The processed variables are updated in the iteration and are denoted by a 'hat' in the following equations. The transformed governing equations are based on velocity and stress.

Finite Difference Discretization and Iteration
The linearized partial differential equations can be discretized in each element into the finite difference forms with discretized variables at the nodes. In the following, u represents any variable in the partial differential equations, i and j the serial number of the element along RD and TD, Δx and Δy are the length and width of the element, c the common ratio for edge refinement (Figure 3b). Then, the discrete forms are given by , The discretized boundary conditions are: ( ) at the upstream boundary, , 0, , where vup and σx,up(y) are the uniform velocity and residual stress distribution at the upstream boundary, respectively; at the downstream boundary, 0, 0; at the center, , 0, , , 0, ; at the edge, 0, 0.
After the boundary conditions are substituted into the discretized governing equations, a sparse system of linear equations can be assembled, expressed as , where u is the vector of unknown variables at all nodes and I the iteration index ( 0 u represents the initial guess of the solution), and A is the square matrix of coefficient, which depends on the initial guess or the updated solution, b the coefficient vector. This linear system can be solved by LU factorization with the SuperLU library (v4.3, University of California, Oakland, CA, USA) [23]. The calculation flow chart is shown in Figure 4. After solving the system of equations, we obtain the calculated values of the variables, based on which the initial guess of the solution and the corresponding coefficients of the equations can be updated with relaxation factors. Then, the updated system of equations will be solved again. This process will be iterated until the residual errors of the variables meet the convergence criterion (10 −4 m/s for the velocities and 0.1 MPa for the stresses). The update of the solution is given by where a is the relaxation factor with a range of 0.2-0.5 to ensure iterative stability.
where kvT is the influence coefficient of roll speed to downstream tension, which can be obtained numerically from Equations (A1) and (A2) in Appendix. The length of Zone D (lD) is fine-tuned in the iteration to satisfy the boundary condition: the rolling pressure at the bite exit pout is zero where klp is the influence coefficient of lD to the rolling pressure at the exit, which can be estimated from the gradients of the rolling pressure.
Based on the streamline method [24], the coordinate y of each node will be adjusted based on the transverse displacement of metal flow, maintaining the node on the streamline. The transverse displacement V is obtained by integration of vy along the streamline

Finite Element Validation and Effects of Elastic Recovery
As the details of the interstand evolutions of profile and residual stress are difficult to obtain in experiments, a numerical validation with the finite element method is performed. The wellestablished finite element method is widely used in the rolling simulation and can be considered as a reference for the new FDM model. A simplified two-stand tandem rolling simulation is designed to evaluate the predictions of the new model with respect to the interstand evolutions of spread, profile, and residual stress. The material behavior is assumed to be elastic perfectly plastic in this section to analyze the effects of elastic recovery alone (excluding stress relaxation).

Finite Element Model
A two-stand hot rolling model is built in the ABAQUS/explicit finite element package (incremental formulation, v6.13, Dassault Systèmes, Vélizy-Villacoublay, France). Due to the symmetry, only a quarter of the system is modeled, as shown in Figure 5. The initial strip thickness is 10.0 mm, and the reductions in both stands are 20%; therefore, the interstand thickness is 8 mm, and the final thickness is 6.4 mm. The initial strip width is 1 m, and the distance between the two stands is 4 m. The strip is bitten into the stands successively, the initial length of the strip is 12 m, and a length of 10 m is rolled to simulate the steady state. The unrolled strip has no crown and residual stress. The rolls are assumed to be rigid and are structured by rigid shell elements (type R3D4); roll radii are 400 mm. The roll contour of the first stand is a parabolic curve with a crown (C1) varying from −10 μm (concave) to 10 μm (convex) to introduce an uneven reduction: the crown of −10 μm corresponds to a higher reduction at the strip edge, the crown of 10 μm to less reduction at the strip edge. The roll contour of the second stand is flat.
In the FEM model, the strip consists of around 520,000 hexahedron elements (type C3D8R) to ensure sufficient refinement at the roll bite and the edge. The material is isotropic and elastic perfectly plastic with Young's modulus of 80.0 GPa, Poisson's ratio of 0.3, and flow stress of 220 MPa. The Coulomb friction coefficient μ is 0.3 at the strip/roll interface. No tension is applied at either end of the strip, and the roll speed of the second stand is fine-tuned so that the average interstand tension is approximately zero.
The same rolling parameters and boundary conditions are used in the simulation by the new FDM model. Figure 6 shows the variation in the half strip width (spread). The pre-bite shrink and post-bite expansion are due to transverse elastic compression and recovery, respectively. The variation weakens as the strip moves away from the roll bite: beyond the distance of 1.5 m, there is an interstand steady region without an evident variation.     As the roll crown of the first stand decreases, in addition to the direct impact on the profile (Figure 7), the spread increases after the first stand ( Figure 6) and the compressive residual stress at the edge becomes larger (indicating an edge wave trend, Figure 8). This is because more reduction at the edge increases the metal flows in the transverse and rolling directions. Moreover, after the rolling of the second stand, the differences between the three cases-with respect to spread, profile, and residual stress-are significantly reduced. This is due to the counteraction of the reverse reduction pattern in the second stand: Case 1 (concave roll) has more reduction at the strip edge in the first stand, but less in the second stand; Case 3 (convex roll) has a contrary situation.

Results
In the above results, the spread and residual stress predicted by the FEM model are the average value of the nodes through the thickness. The new model produces a vertical average due to its quasi-3D nature. The longitudinal evolutions and the transverse distributions of the results predicted by both models agree in general. However, the downstream profile of the second stand predicted by the FEM model cannot reflect the subtle differences induced by the various roll crowns of the first stand ( Figure 7). This is because the profile is studied in the micrometer scale and is sensitive to small fluctuations of the results caused by elements' entering and exiting the roll bite in the transient incremental FEM model.
To understand the evolutionary mechanism of residual stress, details of stress and velocity in the first stand are extracted from the results of the new model. From the longitudinal and transverse stress fields shown in Figure 9a,b, we can see that most areas in the roll bite are under compressive stresses, except the low-pressure region near the lateral free surface (strip edge). Correspondingly, longitudinal tensile stress is observed at the edge near the roll bite. After the strip leaves the roll bite, the transverse compressive stress in the center gradually releases within a range of around 300 mm (causing the post-bite spread in Figure 6). As a result of the Poisson effect of the transverse elastic recovery, metal in the center tends to shrink in the longitudinal direction, corresponding to the decreasing longitudinal velocity at the strip center, as shown in Figure 9c. After the transverse elastic recovery, the uneven strip velocity along the width gradually becomes uniform to satisfy the boundary condition in the interstand steady region: the central velocity increases with the metal elongation, so tensile stress will rise; conversely, compressive stress will increase at the edge. Moreover, the interstand evolution of the thickness profile is attributed to the combined Poisson effects of the longitudinal and transverse stress changes.
Substituting the boundary conditions σy,dn(y) = 0 and vx,dn(y) = vdn, the downstream residual stress can be expressed by the state at the roll bite exit As shown in Figure 9, the evolutionary patterns of stress in the upstream and the downstream interstand regions are symmetric to some extent with respect to the roll bite (excluding the downstream residual stress introduced by rolling). This is because both regions are constrained by similar boundary conditions: the uniform-velocity steady region at one end and the high-pressure roll bite at the other end.
The new model is implemented by C++ programming on a PC (Intel Core i5, 2.60 GHz). There are 20 elements along the transverse direction and 60 along the rolling direction (20 for Zone A or E; 2 for Zone B; 15 for Zone C; 3 for Zone D) with sufficient accuracy. It takes approximately 0.1 s to solve the linear equation system in an iteration; 20-30 iterations are needed to obtain the final solution for one stand. Therefore, the two-stand simulation cost around 5 s. The incremental FEM model requires several days on a workstation (Intel Xeon E5, 2.20 GHz).

Experimental Validation and Effects of Stress Relaxation
An industrial experiment was carried out in a 1580 mm seven-stand finishing hot rolling mill to obtain the strip profile after each stand [14]: the rolling process was stopped in a short time by the emergency button, and the strip profiles at the interstand steady regions were measured by an ultrasonic thickness gauge ( Figure 10). The resolution of the gauge was 0.01 mm and three measurements were performed at each point to produce an average. The new strip model is complemented with a roll stack model to simulate the experiment and compare the predicted and measured profiles for further validation. The viscoplastic material model used is calibrated by hot compression tests. The effects of stress relaxation on profile and residual stress are analyzed by comparisons with purely elastic behavior in the interstand.

Material Model Calibration and Coupling with Roll Stack Model
The rolled material is low carbon steel HP295 (a Chinese standard for making pressure vessels) whose composition is shown in Table 1. The rolling parameters are listed in Table 2. In simulation, the residual stress before the finishing mill is assumed to be zero. The strip temperature is calculated by the online model in the Level 2 control system. The friction coefficients are calibrated by matching the predicted and measured rolling forces, which are within the range of 0.15-0.45 encountered in the literature [25]; the relatively low magnitude in F2 could be attributed to the lubrication effect of the oxide scale formed in the interstand with high temperature and long passing time. To obtain the viscoplastic constitutive relationship of the material, hot compression tests were carried out using a Gleeble 3500 thermomechanical simulator (Dynamic Systems Inc., NY, USA). Cylindrical specimens of 10 mm diameter and 15 mm long were tested (Figure 11a). Figure 11b shows the experimental procedure: the specimens were deformed at different temperatures of 900-1050 °C and strain rates of 0.01-10 s −1 . The reduction was set to 30% (true strain −0.36) to suit the typical finishing hot rolling conditions. Based on the flow stresses obtained under different conditions (Table  3), Equation (14) is fitted, giving A = 2.521 × 10 11 , α = 1.153 × 10 −2 , n = 4.087, and Q = 3.065 × 10 5 . These parameters are obtained by minimizing the average absolute relative error of the flow stress. The correlation between the measured and predicted flow stresses is shown in Figure 12; the maximum error does not exceed 10%.
where the temperature T is in K, Young's modulus E in GPa. These functions are linear regressions of the experimental data from Fukuhara et al. [26], in which the elastic modulus of low carbon steel under different temperatures (300-1500 K) was obtained by the ultrasonic pulse sing-around method.
The new model is coupled with a previously developed roll stack deformation model [27]. In the roll stack model, roll bending and shear based on classic beam theory are calculated by the finite difference method, while roll flattening is calculated by the influence function method. Figure 13 illustrates the coupled calculation: first, with an initial guess of the strip/roll interface geometry, the strip model calculates the rolling pressure distribution for the roll stack model; then, the roll model calculates the new interface geometry for the strip model. These steps will iterate until the results converge. Moreover, the interface geometry and rolling pressure are updated with relaxation factors to achieve stability.  Figure 14 compares the measured and predicted profile at the interstand steady region after each stand. From F1 to F7, the transverse difference of thickness gradually decreases, but the edge drop becomes increasingly apparent. The prediction has captured these trends. To analyze the effects of stress relaxation on profile and residual stress, the stress relaxation is suppressed in the model by setting the plastic multiplier λ in the interstand elements to zero. It is shown that, compared to the original prediction, the thickness difference decreases without stress relaxation in F1 and F2, and the edge drop becomes gentler in the downstream stands. Figure 15a shows the root mean square error (RMSE) of the prediction at the positions of measurements to evaluate the overall agreement across the full width. Asymmetric wedges exist on the measured profiles, while the model is assumed to be symmetric. To reduce its effect on the model evaluation, the measured values used in the evaluation are the averages of pairs of measuring points symmetric to the center. It is shown that when stress relaxation is ignored, the prediction error increases in general: the error in F1 increases by about 15 μm, corresponding to the large deviation at the edge; an exception occurs in F3, which will be discussed below.   Figure 15b shows the relative error of the crown C40 (the difference between the center thickness and the average thickness at 40 mm from both edges, a typical profile index used in industry, helping describe the error introduced by ignoring stress relaxation). A pattern can be found: in F1 and F2, C40 decreases when stress relaxation is ignored, causing negative errors; however, the effect gradually reverses, and C40 becomes larger in F5, F6, and F7, causing positive errors exceeding 20%. The transition from the negative error to the positive occurs in F3, explaining the small error in F3. Figure 16a-c compares the distributions of exit tensions and downstream residual stresses with/without considerations of stress relaxation (the data of F1, F4, and F7 are selected as representatives). The residual stress vanishes in the interstand steady region after F1 (only a uniform tension that is imposed as boundary condition exists), while a compressive residual stress of about −30 MPa remains at the edge (edge wave trend) after F7. The discrepancy indicates that the high temperature and long passing time in the upstream stands may cause nearly complete stress relaxation. When stress relaxation is ignored, a tensile stress of about 50 MPa appears at the edge after F7, contrary to the original prediction. This is because the large tensile stress at the edge near the roll bite is relaxed by plastic elongation, leading to a trend of compressive residual stress. The plastic elongation corresponding to stress relaxation will also trigger metal flows in the transverse and vertical directions, leading to shrinkage in width and thickness. This is the direct cause for the larger edge drop when considering stress relaxation (most evident in F1 and F2, as shown in Figure 14). It is worth noting that the interstand stress relaxation also affects the roll bite deformation, which indirectly shapes the profile. Figure 16d-f compare the transverse distributions of rolling pressure (obtained by integration along the contact arc). The correlation between the distributions of rolling pressure and exit tension is evident; the lower the tension is, the higher the rolling pressure. Therefore, if the exit tension at the edge is relaxed, the overall distribution of rolling pressure will shift towards the edge. As a result, the roll stack deflection will be reduced, and the center crown (excluding the edge drop) at the roll bite exit will decrease, as shown in Figure 16g-i. In addition, from F1 to F7, the exit profile inherits increasingly more from the entrance profile because the transverse flow weakens with the decreasing thickness; this heredity effect also contributes to the smaller center crown in downstream stands when considering stress relaxation.

Results
The coupling model converges in around 40-70 iterations (downstream stands require more due to stronger coupling between the strip model and the roll model). The fast roll model takes approximately 0.01 s to calculate a new strip/roll interface geometry, which is a small cost compared to the strip model. As a result, it takes about 5-10 s to obtain a solution for a stand. The total calculation time of the seven-stand tandem rolling is 56 s.

Discussion
We have developed an elasto-viscoplastic finite-difference rolling model that can predict the entire evolutions (including the interstand regions) of profile and residual stress. The calculation time of a seven-stand hot tandem rolling is approximately 1 min, making multi-pass and multi-case simulation feasible. To achieve a balance between accuracy and efficiency, the key modeling strategies are as follows. 1. Quasi-3D approximation. Simplification along the vertical direction allows reducing the governing equations and meshing only in the rolling and transverse directions, greatly alleviating the computational complexity. On the one hand, this approximation has a theoretical foundation from asymptotic analysis [8]. On the other hand, similar approximations have been adopted by many fast models for hot and cold strip rolling since the development of the classic slab method [28]. The results in the present research also validate the quasi-3D approximation for predicting profile and the transverse distribution of residual stress in the hot finishing mill. However, the quasi-3D approximation will fail in a hot roughing mill or a cold temper mill where thickness reductions are small and the deformation is highly nonuniform throughout the thickness, as suggested by Montmitonnet [29]. 2. Steady-state formulation. Compared to the standard incremental formulation in a commercial FE package, the steady-state formulation allows mesh refinement only near the roll bite, so the number of elements is significantly reduced. This is important when the long interstand regions are considered. 3. Decoupling of multi-stand rolling. The assumption of the interstand steady region with uniform velocity implies that the uneven elastic stress field caused by the roll bite in one stand cannot reach beyond the steady region due to decay (elastic recovery), and therefore, it will not affect the upstream stand. As a result, the system of multi-stand rolling, which is too large to be solved as a whole, can be decoupled. The classic Saint Venant's principle suggests the decay length is of the order of the strip width, and Domanti et al. [19] agreed based on their analytical solution.
Based on the results in Figure 6 and Figure 9, the decay length is about 1.5 times the width. Therefore, this is a proper assumption for the typical hot tandem rolling mill studied in Section 4, in which the length between two adjacent stands is 5.5 m, long enough for stress decay. 4. Global iterative solution with linearization. In the conventional finite difference model proposed by Tozawa et al. [12], the nonlinear governing equations were solved successively from the edge elements to the center elements; this procedure might produce false values and cause failure. In the previous roll bite model [14], the governing equations were linearized to form a global system of equations that considered most boundary conditions simultaneously and could be solved robustly and efficiently in iterations. Based on the same solution technique, the present model can be seen as an extension of the previous model, from the rigid plastic behavior to the elasto-viscoplastic, from the roll bite region to the interstand. A similar solution is found in a meshless model for bar rolling [30], where the material properties are linearized during each iteration. Allwood [31] also proposed an efficient global solution of nonlinear equations in a 2.5D rolling model; however, no mathematical details about the algorithm were given for further understanding.
It is found that elastic recovery is the transition from the high-pressure roll bite to the interstand steady region, during which the transverse compressive stress releases and longitudinal velocity tends to be uniform. During this transition, the residual stress is gradually revealed from the exit tension after a significant pattern change (140 MPa to −30 MPa at the edge and −40 MPa to 10 MPa at the center for the flat roll case in Figure 8). Correspondingly, the residual stress can be expressed as a function of the longitudinal stress, the transverse stress, and the longitudinal velocity at the roll bite exit (Equation (37)) if only elastic deformation is considered in the interstand (in cold rolling for example). Equation (37) may serve as a simplified interstand module for the roll bite module to estimate the residual stress. Compared to the ribbon model proposed by Shohet and Townsend [15], Equation (37) provides richer elastic physics.
It is worth noting that, even if the reduction in the roll bite is uniform along the width, there is residual stress of −30 MPa at the edge in the interstand steady region (the flat roll case in Figure 8). Similarly, Moazeni and Salimi [3] found that there is an edge wave trend with the uniform reduction in their FEM simulation. However, the conventional plane strain theory will suggest no residual stress for the same condition, as the effects of transverse metal flow and elastic recovery are neglected.
It is found that the stress relaxation increases the edge drop of the profile directly in the interstand region in hot tandem rolling but decreases the center crown indirectly by changing the distributions of tension and rolling pressure in the roll bite. The indirect effect indicates a strong coupling relationship between the interstand and the roll bite; fast models often neglect this relationship and calculate these two regions separately. To the best knowledge of the authors, few studies have been conducted to simulate hot tandem rolling and analyze the effect of stress relaxation on the profile.
If stress relaxation is ignored, the profile error increases (from 10% to 25% in crown C40 in F7, as shown in Figure 15b), and the pattern of residual stress changes (a false tensile residual stress appears at the edge in the downstream stands instead of the original compressive stress, as shown in Figure  16a-c). These errors indicate that it is necessary to consider stress relaxation in the online shape control model in industry.
From the purely elastic analyses ( Figure 9) and those considering viscoplasticity (Figure 16), we can see that elastic recovery sets the basic pattern for the interstand stress field, while the stress relaxation works conditionally: its effect is more evident under high stress (existing at the edge region near the roll bite) and under high temperatures and long interstand passing times (corresponding to the upstream stands). Both phenomena are necessary to obtain reliable evolutions of profile and residual stress in hot strip rolling.
In addition to elastic recovery and stress relaxation, elastic buckling (manifest flatness defects) may also occur in the interstand under excessive compressive residual stress and insufficiently imposed tension. Abdelkhalek et al. [32] introduced a buckling criterion into an FEM model for cold rolling. They found that post-bite buckling has little influence on the roll bite deformation because it occurs beyond the region (whose length is about half the strip width) adjacent to the roll bite. On the contrary, elastic recovery and stress relaxation are intense near the roll bite (causing stress decay), so buckling should have limited influence on them. As a result, manifest flatness defects could be evaluated based on the predicted interstand residual stress in a decoupled way.
A combined validation based on both the FEM model and the measured profiles in the industrial experiment is carried out in this research. Comparisons with more experimental or industrial data (especially for the residual stress) are required for further assessment of the model. The Sellars and Tegart constitutive model is used in this research because it can cover both the low-stress behavior in the interstand (particularly the stress decay region near the roll bite) and the high-stress behavior in the roll bite. However, the new model is open to various constitutive models as long as they provide a relationship between the equivalent stress and the equivalent strain rate. Sophisticated physical-based models, which consider microstructure evolution including recrystallization and recovery, can be coupled with the new model to study the interactions between stresses, geometry, and microstructure.
The iterative algorithm with relaxation factor is used for stability when solving the coupled system of strip and roll stack; however, it will cost more iterations and CPU time. One possible solution is to linearize and discretize the governing equations of roll deformation similar to the strip model and solve the strip-roll system as a whole. It should further increase efficiency, making online applications possible.

Conclusions
An elasto-viscoplastic finite-difference rolling model is developed to analyze the evolutionary mechanisms of profile and residual stress under the effects of the interstand elastic recovery and stress relaxation. A combined validation, based on both the FEM model and the industrial experiment measurements, shows satisfactory accuracies, with a maximum error of 11% in the predicted crown C40. A seven-stand tandem rolling simulation including the interstand regions costs approximately 1 min of CPU time, allowing the new model to run massive parametric campaigns for industrial process optimization and giving it the potential for online applications.
Model reduction techniques-including quasi-3D approximation, steady-state formulation, and decoupling of tandem rolling-are found to be effective. The proposed global iterative solution with linearization is the key to efficiency and robustness.
The results obtained from the analysis can be summarized as follows: • The interstand elastic recovery extends from the roll bite exit for a length of about 1.5 times the strip width, and it sets the basic pattern for the interstand stress field. • During elastic recovery, the transverse compressive stress releases, and strip velocity tends to be uniform. As a result, the residual stress is gradually revealed from the exit tension with a significant change of stress pattern, and the strip profile is reshaped.

•
Stress relaxation is more evident in the high-tension edge region near the roll bite and in the high-temperature low-speed upstream stands.

•
Stress relaxation increases the edge drop of the strip profile directly in the interstand but decreases the center crown indirectly by changing the distribution of the rolling pressure in the roll bite. When stress relaxation is ignored, the error of the predicted crown C40 will increase from 10% to 25% at the exit of the tandem mill, and a false tensile residual stress will appear at the edges in the downstream stands. A better understanding of these interstand phenomena and the evolutionary mechanisms of profile and residual stress would help develop simplified online models or process strategies for shape control.
The present model can be further developed to adapt to the conditions in multi-pass or multistand cold rolling. In addition, since the model calculates the transverse metal flow, if the effects of thermal expansion are considered, it could be used to predict width variation during rolling.