Coupled heat transfer and phase transformations of dual-phase steel in coil cooling

Dual-phase steels are generally used in the car industry due to high tensile strength and good formability, which are obtained by a mixture of bainite and ferrite phases. This microstructure is achieved through slow rate coil cooling. However, the manufacturing of dual-phase steels introduces various challenges such as the instability of the cold rolling process. An important factor affecting this is the non-uniform coil cooling of a hot rolled strip. In coil cooling the cooling rates are not controlled and there are different thermal contacts during coil conveyance causing unequal cooling of the steel coil. Unequal cooling rates lead to non-uniform coil cooling, producing irregular phase transformations on different sides of the coil, which causes periodical variations of the phase fractions in the steel strip. Varying phase fractions cause thickness deviations in the strip during the cold rolling process. A three-dimensional transient heat transfer finite element model was developed and used for modeling the complete coil conveyance chain and coil field cooling of the coil on an industrial scale. A coupled phase transformation model is implemented as a subroutine into the finite element model for calculating the resulting phase fractions. It was found that the different thermal contacts during the coil conveyance produce uneven cooling rates causing length- and widthwise variations in the phase fractions. The heat transfer model is vali- dated by comparing temperature profiles between the simulated and measured coil edges. The phase transformation model is fitted into experimental data and verification is carried out in industrial conditions by comparing the modeled phase fractions and test samples from a cooled and unwound steel coil.


Introduction
Dual-Phase (DP) steels are high-strength automotive steels. They feature a balance of high strength, ductility and formability [1]. Tensile strength ranges from 450 to 1200 MPa and is traditionally controlled by the phase fraction of martensite. In modern DP steels, martensite is replaced by bainite to obtain better formability. In this paper the coil cooling of Ferrite-Bainite Dual-Phase (FBDP) steel is studied.
Several investigations prove how the differential accelerated water cooling on the run-out table can result in residual stresses and flatness defects in a hot strip. Cho et al. [2] developed a thermo-mechanical Finite Element Method (FEM) -model coupled with a phase transformation model to predict edge waves in a hot rolled plate during run-out table cooling. This study showed how significantly differing cooling rates influence strain gradient due to phase transformations. In the current paper, similar phenomena are applied to the modeling of the coil cooling process.
Karlberg [3] described the development of a thermo-mechanically coupled solution to predict coil cooling as well as axial, circumferential and effective stresses in the coil. He claimed excellent agreement between the measured and modeled cooling paths of coil cooling and introduced circumferential stresses within the coil. However, the validation of the temperature model was only carried out indoors in a stationary position. In reality, the coil conveyance line contains several different thermal contacts and possible coil field cooling outdoors. More detailed thermal contact modeling with coil conveyors is required when the phase transformation model is implemented into a coil cooling model. In this study validation is carried out comparing experimental measurements of an ongoing full-scale industrial process to modeled results of the FEM-model.
Wang et al. [4] showed that the phase transformation has a notable effect on the prediction accuracy of residual stress generation on run-out table cooling. In the present study, the phase transformations do not occur during water cooling and thus the fully austenitic strip is assumed to be coiled without any flatness defects. However, the stress state and contact pressure distribution during the coiling process and in post-coiling conveyance are considered in thermal contact conductance.
The residual stresses and flatness defects place challenges on further processing, such as laser cutting and cold rolling, making these processes unsteady. A similar important issue related to further processing is the process of coil cooling, particularly the cooling of hot coiled dual-phase steels. Milenin et al. [5] investigated the longitudinal stresses of a steel strip during hot rolling, laminar cooling and coiling. Witek and Milenin [6] have also studied the effect of coiling temperature and the cooling of coils on residual stresses between upper and lower strip surfaces. In both above-mentioned studies, the phase transformations during the cooling were considered and were found to have a significant impact on the level of residual stresses. However, validation of the temperature model was not presented explicitly over the long-lasting coil cooling and a detailed coil conveyance path corresponding to reality was not considered. According to the results of the current study, the thermal contacts in coil conveyance have a great effect on phase transformation and should always be considered when phase transformations occur in modeling of coil cooling.
Kaputkina et al. [7] experimentally discovered the non-uniformity of mechanical properties in hot rolled low-carbon multi-phase steel produced in conventional process conditions. They concluded that phase transformations under differential coil cooling influence the stability of the mechanical properties of the steel coil. Diverging from these investigations, the model developed in this current study accounts for phase transformations coupled with a complete heat transfer simulation at coil conveyance and coil cooling on the coil field. In this way a reliable cooling paths can be modeled, and close-fitting boundary conditions are produced for the coupled phase transformation model. The focus of this study was on microstructural evolution during thoroughly modeled coil conveyance and heterogeneously formed phase fraction distribution in the steel coil, which allows for quantitative evaluation of the phenomena under different conditions. DP steel strips are generally downcoiled as fully austenitic at ~ 590-700 • C after being hot strip rolled and cooled at a slow rate on a run-out table. Thus, phase transformations take place during coil cooling, which may lead to heterogeneous phase fractions due to nonuniform cooling rates in the coil. This regionally varying microstructure in the coil shows as strength deviations between different phases, making the further cold rolling process unsteady, which finally results in thickness deviations in the cold-rolled strip. As an initial state for the simulation the microstructure was assumed fully austenitic, although some minor fractions of martensite or bainite might have formed on the surface of the water-cooled strip. However, such minor fractions would be similar along the strip length and they would not contribute to the heterogeneous variation of the phase fractions at room temperature.
The steel coil is fixed in a certain position during the conveyance, causing unequal heat flow gradients on the bottom section of the coil compared to the freely cooling (radiation and convection to air) top section. This transient cooling process must be modeled precisely to ensure accurate cooling rates for the implemented phase transformation model. All coil contacts with 13 conveyors are considered including conveyor temperatures and hold times during the industrial process. Temperature measurements of the coil conveyors were carried out experimentally at a hot strip mill during a running process to ensure consistent data.
The heat conduction of superposed coil revolutions is quite different compared to internal heat conduction within the material. Thus, real coil geometries were applied to simulate radial heat conduction between individual revolutions within the coil. Anisotropic contact properties exist on strip interfaces due to impurities like water and oxide as shown by Park et al. in [8] but process deviations such as strip crown, flatness and thickness variations also affect the heat conduction at contact interfaces. Heat transfer on contact interfaces occurs by several mechanisms, such as conduction over actual contact regions, conductions through the air and radiation in contact gaps.
Karlberg [3] used radial contact pressure between adjacent revolutions during the coil conveyance depending on the conveyors' contact points and varying internal stresses. Changes in internal stresses in the coil derive partly from differential cooling and phase transformations. However, the used uncoupled heat transfer analysis does not generate any mechanical contacts and thus this model uses temperature dependent radial post-coiling contact pressure for thermal conductance on strip surface-to-surface contacts. The definition of post-coiling stresses within the steel coil are presented by authors in [9].

Numerical simulations
The main interest of this study was to develop a highly accurate heat transfer model of the coil cooling, including the coupled phase transformation evolution. FEM was chosen as the simulation method to include realistic consideration of the coil geometry and coil eye heat reflection features as well as arbitrarily assignable heat transfer boundary conditions. The simulations utilize the non-linear facilities of a transient uncoupled heat transfer analysis with the standard solver of Abaqus™ software. A phase transformation model was numerically implemented into the Finite Element (FE) calculation using a UMATHT subroutine. Convection from the oxidized surfaces of the steel strip on the outer surfaces of a coil is defined by a convective heat transfer coefficient in a FILM subroutine. Thermal contact conductance conditions were defined in a GAPCON subroutine, available in [10].
Displacements and stresses were not included in the analysis as they cannot be calculated in a pure heat transfer analysis as defined in [10]. Thus, 8-node linear heat transfer brick elements DC3D8 were used in all instances (Fig. 1). This element utilizes full Gaussian quadrature. Due to fully Gaussian integration, only 5 elements were required in the widthwise direction of the symmetrical coil. Coil geometry was created by an Excel-based script to generate 2D DXF geometry of a steel coil and extrusion to 3D was obtained using Catia™. Coil geometry mimics real coiled strips, disregarding the strip crown.

FEM assembly
Uncoupled heat transfer analysis assumes that the thermal and mechanical problems are not dependent on one another and thus displacements and stresses are not calculated. All the conveyors are placed in contact with the coil and only thermal contacts are generated during the analysis, using model change and interaction activation. The model change interaction in Abaqus allows the deactivation and activation of geometries and elements during the analysis but the changes must be implemented in individual steps. This speeds up computing when the analysis contains several steps and certain instances are used only during certain steps. Contact management is also carried out by model changes with reference to [10]. Complete coil conveyance contains 36 coil-conveyor contacts and coil field cooling with a boundary condition coil. An uneven coil field cooling process was achieved by applying an isothermal hot coil next to the simulated coil as a boundary condition. Every contact during coil conveyance was separated into individual steps to control model changes and the initial temperatures of the conveyors. This was required for the conveyor geometries, which go through several contacts with the coil but are placed in different locations in the real process.
Different thermal contacts also produce high temperature gradients on contact regions so time increments must be decreased considerably. The time increment is limited up to one second on all steps to avoid overgrowth of martensite and bainite. Steps begin with a 0.01 s time increment to attain converged thermal contact and then increase to one second with respect to the iterations to achieve the converged time increment. Maximum temperature change per time increment was 1 • C.
All the interactions during 36 contacts as well as the internal heat transfer within the coil, including cavity radiations, are created in the initial step. At the beginning of the first step, all the interactions which are not used during this step are deactivated and then reactivated in the upcoming steps when required. Conveyor geometries are removed by model change and interactions deactivated after every step if they are not reused in upcoming steps. With this optimization the model is computationally less expensive and a large amount of interactions is easier to manage by the user. The conveyance steps with used hold times and conveyor temperatures are shown in Table 1.
The simulation was commenced with step 1, by having the thermal contact between a mandrel of the downcoiler and the strip. The initial temperature state of the strip was isothermal (595 • C). This step does not consider any mechanical coiling process, only thermal interactions occurring during the hot coiling. Mechanical coiling to find out internal contact pressure within the coil was performed in an individual analysis using an explicit solver by authors in [9]. The downcoiler mandrel is water cooled efficiently between coiling processes, but mandrel-strip contact can be assumed to be a pure steel-steel thermal contact without any disrupting factors, such as water. The same steel-steel thermal contact was used in the conveyor-coil contacts for all 36 steps.

Convective and radiant heat transfer
Heat transfer from coil surface to ambient is determined by convection and radiation. An emissivity of 0.8 is defined on an oxidized steel surface in [11] and the actual ambient temperature from experimental measurements is used. Radiation to ambient is calculated with Eq. 1: where θ j and θ i are surface and ambient temperatures, respectively. θ Z is the absolute zero temperature (-273.15 • C), emissivity is multiplied by the Stefan-Boltzmann constant as the radiation constant A as in [10]. The equation for convective heat transfer used in [12], is expressed as follows: where Q conv (J/s) defines heat transferred by convection, h conv (W/ (m 2 K)) is the heat transfer coefficient of convection, A is the area of surface (m 2 ), T p (K) is the steel temperature and T a (K) is the ambient temperature. The Nusselt number (Nu), which is the ratio between the convective and conductive heat transfer, was used to calculate the heat transfer coefficient of convection Nu = (h conv D)/k, where D (m) is the characteristic length and in this case the outer diameter of the coil and k is the thermal conductivity of air (W/mK). Nu depends on the Rayleigh number (Ra), which determines the air flow regime (laminar or turbulent). Dash and Dash [13] defined the formula of Nu for a thick hollow cylinder placed horizontally on the ground. They published a relation between Nu, width to outer diameter ratio (L/D) and the ratio of inner diameter to outer diameter (d/D) of the coil. However, the equation was defined only within the laminar air flow regime. The Rayleigh number (Ra) in our case exceeded 10 9 , which means the air flow turned out to be turbulent due to the large size of the coil as described in [12]. Therefore, Nu is solved by utilizing Ra, the Grashof number (Gr), which depicts the dependence between viscosity and buoyancy within a fluid and the Prandtl number (Pr), which defines the ratio between thermal diffusivity and momentum diffusivity. These are defined as follows: where g (m/s 2 ) is gravity, β (= 1/T a , 1/K), η (m 2 /s) is air kinematic viscosity and μ (kg/ms) is air dynamic viscosity. Since Ra > 10 9 , Nu is defined as Nu = 0.14•Ra 0.33 . Substituting Eq.s 3-5 to Eq. 6, we get: and h conv ' with respect to the steel coil surface temperature used in the FILM subroutine: Nu derived with Ra is applied during coil conveyance when the air Table 1 Conveyor hold times and initial temperatures.
Step inside the rolling mill building is not flowing (no wind). After coil conveyance the coil cools outside on the coil field. Air is flowing due to wind, which enhances the cooling effect. Therefore, Nu is defined using the Reynolds number (Re) instead of Ra. Re depicts the ratio between inertial forces to viscous forces within a fluid which is subjected to relative internal kinetic of different fluid velocities as used in [14]. Nu is presented using Re and Pr as defined by Ozisik in [15].
where A is 12.0, e is Euler's number and V is relative fluid velocity (m/s) (in this case, wind).

Cavity radiation
Reflection symmetry of cavity radiation interaction was exploited to define coil eye reflection and radiation to contiguous coils in the coil field (Fig. 2). The coil center in the coil eye, in particular, radiates profusely between the opposite surfaces and cooling is significantly slower compared to coil edges. Using symmetrical reflection, only half of the width of the steel coil was modeled and finer discretization was applied into the symmetric coil geometry. The three outermost and innermost strip layers on the coil were meshed with extra fine discretization to avoid numerical inaccuracy on contact regions.
Strip thickness in the simulations was 3 mm and a single element was used throughout the thickness. The element length on fine discretized regions was selected to be 5 mm and on coarse regions 20 mm. This ensures good dependability on the coil surface for convection, conduction and radiation interactions. The contact regions with conveyors, in particular, require a very dense element mesh to calculate thermal conduction from the coil surface to the inner revolutions of the coil. The significance of thermal conductance is emphasized with a coupled phase transformation model and fast phase growth of martensite and bainite.
The surfaces which thermally radiate against each other are considered through cavity radiation in Abaqus. The theory of cavity radiation is based on the generally known theories of Siegel and Howell in [16] and Holman in [17]. Cavities are defined in a three-dimensional model by the element faces and they are treated as isothermal and have a constant emissivity. Required cavity radiation elements are automatically generated by Abaqus. The formulation of the cavity radiation flux is briefly shown below and is extensively depicted in [10].
Thermal radiation in cavity radiation is based on the gray body radiation theory and only diffuse reflection is taken into account. In addition, deterioration of the radiation in the cavity medium is not considered. Hereby, the equations for the radiation fluxes per unit area into cavity facets are as follows: where subscripts i and j refer to facets, q c j is the flux into facet j, ε i,j are the emissivities, σ is the Stefan-Boltzmann constant, F ij is the geometrical view factor matrix, θ i.j are the temperatures and δ ij is the Kronecker delta.
Using the temperatures of nodes, the temperature radiation power η i is calculated. Then the average facet temperature radiation power is interpolated from the facet nodal temperatures η N . The subscript i and the superscript N refer to facet and nodal quantities, respectively. P N i is the nodal contribution factor calculated from area integration. Equations for temperature radiation power are the following: where A i is the area and N N i are the interpolation functions for facet i. Now the radiation flux into facet i can be presented as depicted in Eq. 17: And the nodal contributions from the radiation flux on each facet is derived as in Eq. 18: Finally, the total cavity radiation flux at a certain node N is:

Energy balance and constitutive definition
Uncoupled heat transfer analysis is perfectly suited for the coil cooling calculations. It is intended to model temperature-dependent conductivity for a solid body considering latent heat effects, convection and radiation while acknowledging cavity radiation. The basic energy balance needs to be achieved in coil cooling analysis. This is determined by Green and Naghdi in [18] using Eq. 20, written as: where ρ is density, U is the material time rate of the internal energy, V is the volume of solid material, q is the heat flux unit per area of the body, S is surface area and r is the heat generated internally into the body per unit volume. As an uncoupled model, the internal energy is not dictated by the displacements or strains of the body but only by temperature, U = U(θ). Heat conduction is defined by the Fourier law in Eq. 21: where k is the conductivity matrix k = k(θ), f is heat flux and x is position. In this case conductivity is fully isotropic. Thermal conduction coefficients between adjacent strip revolutions are considered according to the theory of thermal contact conductance by Mikić [19] in Eq. 22: where k s is the thermal conductivity of steel by Browne in [20] and is also used by authors in [9], P is the contact pressure defined by authors in [9], E ' is the effective elastic modulus, m s and σ s are mean asperity slope and root mean square value of roughness at the interface, respectively.

Implementation of phase transformation model
In the current study, we apply the differential form of the Kolmogorov-Johnson-Mehl-Avrami equation [21], which was described by Leblond et al. in [22] to calculate the austenite decomposition in to polygonal ferrite and bainite during cooling. The model has been used earlier in the coupled modeling of the transformation behavior during cooling by Pohjonen et al. in [23] and validated for a simulation of transformations which occur during the arbitrary cooling path by Pohjonen et al. [24]. In industrial practice such models have been applied in settings where the cooling is achieved, for example, by using water cooling systems, as described by Donnay et al. in [25] and more recently by Milenin et al. and Ilmola et al. separately in [5] and [26]. In such cases the appropriate method of obtaining the model parameters is by conducting continuous cooling experiments because the aim is to capture the phenomena occurring during short time-scales and rapidly changing temperatures. However, in the case of coil cooling, it is better to fit the model corresponding to isothermal holding experiments because then the model and the parameters reproduce the phenomena which occur during slow cooling.
The phenomena, such as carbon diffusion due to the ferrite and bainite transformations, which are likely to be very different to the case of fast cooling, will then be correctly taken into account. While the transformation model is very similar to the one used in earlier studies for faster cooling by authors as well as other researchers, it is shortly described here in order to provide the obtained model parameters and to show the details of the model, as there are some minor differences in the choice of functions.
The pearlite formation was modeled with the equation applied by Kirkaldy and Venugopalan in [27] for the statistical analysis of transformation kinetics of different steel compositions. A similar model has been applied to describe transformations in welding by Henwood et al. in [28] with the numerical procedure described by Watt et al. in [29]. A similar equation has been used in phase transformation modeling also by Martin and Victor Li et al. separately in [30] and [31]. In the current work, the whole transformation model was written as a subroutine to the Abaqus FE software, which allowed the coupling of heat conduction and phase transformation simulations including the release of the heat from the transformations.
The incubation time for ferrite and bainite formation during cooling is taken into account by calculating the time τ required for the transformation to start (assumed as 1% phase formed) by holding at constant temperature. The beginning of the transformation during cooling can then be estimated by applying the rule of Scheil, as originally described by Scheil in [32] and applied later by many authors, including Pohjonen et al. [33].
where the index i is replaced with f for ferrite or b for bainite referring to the corresponding transformations.
Once the ferrite or bainite transformation has started, the subsequent kinetics [22] are calculated using Eq. 24.
The maximum phase fraction of ferrite, χ max,f is calculated using the equilibrium diagram as first applied by Donnay et al. in [25] and later by Pohjonen et al. in [34]. Below the eutectoid temperature χ max,f = (C max − C 0 )/(C max − C f ), the equilibrium austenite carbon concentration C max is the carbon concentration corresponding to the eutectoid point.
For bainite this is χ max,i = 1, as it was observed in the experiments that all the austenite could transform into bainite (thin austenite films between bainite laths are considered as part of bainitic microstructure in the current work). The rate parameter function k i used in this study is calculated using Eq.s 25 and 26.
for ferrite and for bainite. The carbon concentration of the austenite is calculated as to account for the enrichment of austenite carbon concentration due to the formation of polygonal ferrite and bainite during the slow cooling. The pearlite formation was described in [27] using Eq. 27.

∂ ∂t
where K p is a fitting parameter and the coefficient D was calculated by applying Eq. 28.
where the gas constant R is given in kJ / mol • C. The martensite transformation was modeled using a differential form [30] of the Koistinen-Marburger type equation, Eq. 29.
∂T ∂t (29) where the parameters η and M s were fitted to the experimental data obtained from the continuous cooling experiment. To avoid unphysical formation of new martensite in possible repeated heating and cooling cycle below M s , the martensite formation was only allowed when the temperature dropped below the lowest temperature which the material point had experienced. For this reason, keeping track of the lowest temperature of the material at every position during the simulation was required. The heat capacity c and thermal conductivity k were calculated as a function of the austenite fraction and temperature based on the mixture of ferritic and austenitic phases and the respective temperature dependence of the quantities of those phases, that is, c = (χ f + χ p + χ b + χ m ) c α + χ γ c γ and k = (χ f + χ p + χ b + χ m )k α + χ γ k γ , where the austenite and ferrite heat capacities c γ and c α were calculated according to [20] and [30].
The heat released due to the transformations, Δs, during the timestep Δt was calculated using the values presented by [35] Q f = 5.9 × 10 8 , Q p = 6.0 × 10 8 , Q b = 6.2 × 10 8 and Q m = 6.5 × 10 8 J/m 3 . These values for the latent heat were also used more recently in [36].
where L is the temperature dependent latent heat, calculated as described in [20] and [30].

Experimental methods
Several experimental methods were utilized to verify results of numerical simulations. Gathered experimental data was used to fit the phase transformation model to currently used DP steel. The material was continuously cast on a vertical-curved caster and rough-rolled into a 210 mm thick slab. The nose of the slab was cut after rough rolling and samples of ∅6 × 9 mm were drilled and machined in the direction of the thickness. The composition of the studied steel in wt.% was 0.18C-0.24Si-1.84Mn-0.31Cr-0.15Mo and the steel was aluminum-killed and calcium-treated.
To obtain accurate results, the dilatometry measurements during isothermal holding at different temperatures were performed with the Gleeble thermomechanical simulator (Gleeble 3800, Dynamic systems Inc., New York, USA). A shielding gas environment was used to minimize oxide scale formation. The schematic graph of the thermal schedule is presented in Fig. 3.
In the first stage of the study, the holding temperatures were: 700 • C, 675 • C, 650 • C, 625 • C, 595 • C, holding time 8000 s. 550 • C, 500 • C, 450 • C, 400 • C, holding time 4000 s. For the abovementioned holding times, no pearlite formation was observed for any temperature. However, pearlite was seen in the actual coiled material. It seemed that the holding time in the experiments was too short. For this reason, additional experiments were made, with a 10 -h holding time: 550 • C, 565 • C, 575 • C, 595 • C, 625 • C, holding time 36,000 s. After these holding times, pearlite was observed in the material at temperatures of 550 • C, 565 • C in the sample center and at 575 • C, 595 • C near the edges. These results could be used to fit the pearlite formation model.
To distinguish martensite formation, continuous cooling transformation (CCT) experiments were performed. The heating rate was 10 • C/s to 1100 • C followed by a 2 min holding time. The samples were cooled with a rate of 2 • C/s to 920 • C. Three 0.2 strain deformations were performed with 5-second intervals and finally the samples were cooled to room temperature with nine different rates from 1 • C/s to 70 • C.

Microscopy experiments
Microstructural examinations were carried out using a laser scanning confocal microscope (LSCM, VK-X200, Keyence Ltd, Itasca, USA) and a field emission scanning electron microscope (FESEM, Zeiss Sigma, Zeiss International, Oberkochen, Germany) using 2% nital etching. Because martensite formed during several experiments during the cooling to room temperature (after the holding), it was in some cases difficult to distinguish between bainite and martensite phases. Scanning electron microscopy with a better resolution helped in distinguishing the different types. However, it was difficult to recognize from the microscope images alone which areas corresponded to bainite and which areas corresponded to martensite. For this reason, microhardness measurements were applied. The combined use of different methods and comparing microstructure images from different holding temperatures enabled forming a better understanding of what microstructures were visible in the microscope images. The applied methods and the performed analysis are described in detail in [37].

Hardness measurements
Macrohardnesses were measured as a 5 measurements average at the center of the sample using a hardness tester (Duramin A300, Struers, Copenhagen, Denmark) under a 50 N load (HV5). Microhardnesses tests  were performed using a microhardness tester (MHT-Z-AE, CSM Instruments, Needham Heights, USA). The microhardness measurements were used in distinguishing between the different types of regions observed in the microscopy experiments.

Coil cooling measurements
Experimental coil cooling measurements were carried out for the model validation. Two similar coils were rolled one after another. Both coils went through the same conveyance chain and were cooled side by side on the coil field at a distance from other coils ( Fig. 4a and b). Coil conveyance inside the rolling mill building took about 40 min after which the coils were moved to the coil field for the final cooling. The temperature profile was measured once per hour for eight hours.

Results
The phase transformation model parameters were obtained by fitting the model to the experimental dilatometry data. The experimental procedure of obtaining the dilatometry data is described in the experimental section. For the numerical fitting, the Nelder-Mead algorithm by [38] was applied with the Matlab fminsearch function [39]. The obtained parameters are given in Table 2. The martensite start temperature (M S ) was additionally compared to the dilatometry data assuming linear thermal contraction of the steel at the temperatures above M S , which confirmed that the obtained parameter value was correctly determined. Table 3 below summarizes the results obtained from the microscopy, hardness and image analyses, which were aimed to estimate the fractions of bainite and ferrite which form during the isothermal holding, as well as the fraction of martensite which forms during the subsequent cooling. The fractions of ferrite, bainite, pearlite and martensite were obtained by analyzing the SEM images after the holding and cooling using an advanced image analysis method devised by [37].
For holding at 400 • C, the hardness was higher than at 450 • C because the bainite is finer (lower bainite). For holding at 500-550 • C, less bainite and more martensite is formed during cooling, resulting in higher hardness at room temperature. For holding at 595-650 • C, more polygonal ferrite is formed, which lowers the hardness. When the holding temperatures were higher than 650 • C, less polygonal ferrite and more martensite forms during cooling, resulting again in higher hardness at room temperature.
No pearlite was found in the experiments for the holding times of 4000 or 8000 s, as shown in Table 3. Since pearlite was found in real coil-cooled material, it had to form during longer holding times. For this reason, experiments with 10 -h holding times were conducted. For these long holding time experiments pearlite was observed using microscopy. The comparison of the kinetics calculated using the fitted model to the experimentally determined phase fractions are shown in Fig. 5.
For the purpose of constructing the computational model, this study distinguished between four types of ferritic phases: (polygonal) ferrite F, bainite B, pearlite P and martensite M. Although bainite and ferrite could be further distinguished into different classes (e.g. upper, lower, granular bainite, etc.), here a common class of bainite was used to describe the phase which has a clear acicular morphology and has a distinct C-curve of transformation separate from the C-curve of polygonal ferrite, as shown in Table 3. This broad definition of the bainitic phase includes also the Widmanstätten ferrite, which is shown in Fig. 6. The microstructures for the different holding temperatures 400 • C, 500 • C, 595 • C and 700 • C are shown in Figs. 6a, b, c and d, respectively.
After fitting the ferrite and bainite parameters, the scaling parameter K p for pearlite formation was fitted to the experimentally observed pearlite fraction ( Fig. 7a and b.). The comparison of the phase fractions observed in the 10 -h holding time experiments is shown in Table 4.

Temperature model validation
Measured temperature profiles (vertical line in Fig. 4) are compared to modeled profiles in Fig. 8. The modeled profiles respond well with Table 3 Phase fractions formed during 4000 and 8000 s isothermal holding (ferrite and bainite) and during cooling (martensite) and macro hardness (HV5), as well as the standard deviation of the hardness (HV5 SD   measured profiles. The greatest differences between the results were after 2-4 h of cooling. During this time the model cooled down faster than the experimental profile and the profile shape remained steadier than in the experimental profile. There are a few explanations for this difference. The first is weather. Sunshine, cloudiness and wind cannot be defined exactly as they appears in reality. The second is radial contact pressure inside the coil. In the model a static post-coiling contact pressure distribution within the coil is assumed. In reality the contact pressure varies during cooling due to temperature changes and phase transformations including the effect of latent heat. All these affect the thermal conductance inside the coil and the thermal boundary conditions on the coil surface. However, the temperature profile matches up to the experimental profile after 4 h of cooling and then continues to cool synchronously with the measurements. Also, the profile shape agrees with the measured profile, confirming the validity of the internal thermal conductance.

Cooling paths and phase fractions
The whole of the coil conveyance chain and coil field cooling were simulated. The cooling time was roughly 12 h until microstructure transformations were completed everywhere part within the coil. The coil was treated as isothermal with a coiling temperature of 595 • C. Strip thickness was 3 mm, width 0.583 m (symmetric) and length 808 m. The geometry for the boundary condition coil is equal to the main coil. This coil was used only on the coil field and the initial temperature was copied from the main coil after coil conveyance. The phase transformation model was implemented only in the main coil.
Cooling paths in four different surface points are shown in Fig. 9. Point 1 is at the top section of the coil and it cools freely in the air. Inside the mill building, other hot coils are on the same transportation line and thus cooling is relatively slow before conveyance to the coil field. In the coil eye (Point 2), the contact with the mandrel of the downcoiler causes heavy temperature descent at the strip surface. However, it is noteworthy how the temperature rises in the coil eye when coil conveyance starts. Points 3 and 4 are on conveyor contact surfaces which go through different thermal contacts with the coil conveyors during the conveyance chain. Point 3 is in direct contact with different conveyors and temperature changes are significant between the different conveyors. Point 4 is not directly in thermal contact but heat is transferred through the radiation and conduction from adjacent conveyor contact interfaces. Temperatures are collected at the end of every conveyance step.
Coil cooling is noticeably different on the bottom contact surface compared to the contactless top surface of the coil. Figs. 10a, b, c and d depict cooling paths for eight different positions (not the points in Fig. 9) on opposite sides of the coil. The letters B, F, M and P in Figs. 10a, b, c and d refer to phases bainite, ferrite, martensite and pearlite, respectively. After the phase grade and underscore, T refers to freely cooling top section and B to the bottom section. The last number defines the position for the data curve within the coil. This number portrays the amount of revolutions starting from the outermost layer. Differences in the cooling paths in opposite positions are clear even underneath the coil outer surface.
Temperatures of positions near the conveyor contacts descend effectively compared to the freely cooling top region of the coil which cools slowly. In internal revolutions the temperature differences are not that clear and generated latent heat due to phase transformations can be noticed. Phase transformations depend on cooling rates and do not occur simultaneously in the coil and thus latent heat "waves" can be seen in Fig. 10d.
The transformed phase fractions for bainite, pearlite, martensite and ferrite regarding the cooling paths are presented in Figs. 10a, b, c and d. The outermost revolutions cool down rapidly, simultaneously forming bainite. In regions with slower cooling rates, mainly ferrite transforms up to 71-% due to carbon partitioning and diffusion. The remaining austenite transforms into pearlite and bainite. The high cooling rate required for martensite growth is obtained only on conveyors and mandrel contact interfaces. Therefore, martensite is formed only at certain contact areas and reaches only a few revolutions below contact surfaces.
Figs. 10a and b prove how the different phases generate heterogeneously within the DP steel coil during the cooling. Phase fractions in the cooled and unwound coil with finished microstructure are depicted in Figs. 11a, b, c and d. The phase fraction distribution in a lateral direction is shown in the charts which are equally divided into six slices between symmetry plane and coil edge. Phase fraction data is gathered from the node points of the FE-model. The averaging threshold value was set to 75 percent. In which case nodes common to two or more elements receiving multiple contributions will average contributions if the relative difference between them is greater than 75-%.
The outermost revolutions are mainly bainitic apart from the contact interfaces which are partially martensitic. After 20 m (four revolutions) of outer strip layers, the phase fraction of ferrite increases linearly, coinciding with the diminished bainite phase fraction. An important observation is that the phase fractions of ferrite and bainite alternate periodically from 500 m to 811 m. This is a consequence of unequal cooling rates on the opposite sides of the coil. In the body of the coil, the phase fractions are steadier on a single strip revolution. Minor periods can be separated between bainite and pearlite phase fractions while the ferrite fraction remains constant.
Near the strip head that is the coil eye, phase fractions begin to vary again but the periods of variation are slighter in comparison to the strip tail. The same effect in phase fractions can be seen in the lateral direction of the coil. The cooling rates are naturally higher at the coil edges than at the symmetry plane (centerline of the coil) in pure air-cooling conditions. The different phase fractions are specified in individual figures; Figs. 11. a) bainite, b) ferrite, c) martensite and d) pearlite.

Comparison to industrial data
Regionally inhomogeneous microstructure in rolling or transverse direction results as flatness defects and thickness deviations due to strength differences in different phases. In this case the strip is almost fully bainitic at the outermost revolutions until 20 m from strip tail. The

Table 4
Observed phase fractions formed during 36,000 s isothermal holding (ferrite, bainite, pearlite) and during subsequent cooling (martensite, bainite). On the other hand, phase fractions vary considerably on outer revolutions where the cooling effect of conveyors is present and unequal cooling rates produce inhomogeneous microstructures on the bottom and top sections of the coil as well as in the coil eye. Phase fraction variations during a single revolution in the strip cause thickness deviations in the high-speed cold rolling process.
Thickness deviations (tdev) measured from an industrial cold rolling process of DP steel compared to phase fractions and hardness of bainite, ferrite and pearlite from a modeled and unwound strip are presented in Figs. 13-16. Strip lengths are normalized due to the unequal final lengths of a hot strip and a cold-rolled strip. A comparison to averaged hardness distribution (Fig. 13) in the strip was calculated from the modeled phase fractions on the symmetry plane. The measured HVhardnesses of individual phases were 184, 261, 322 and 422 for ferrite, pearlite, bainite and martensite, respectively. The periodical variations of ferrite and bainite at both ends of the strip correspond to thickness deviations in the cold-rolled DP steel strip. The hard ends of the strip are emphasized due to a major fraction of bainite at both ends. This also correlates with the thickness deviations of the strip, which makes the thickness deviations clearly explicable with phase fraction variations in Fig. 14. Fig. 14 is divided into four sections (Fig. 15 a, b, c and d) to compare phase fraction variations to thickness deviations. In the coil eye (section 1), phase fractions of bainite and ferrite vary significantly in a relatively short distance which can be seen as large thickness deviations. Thickness deviations are notably reduced when phase fractions become uniform.
In section 2 all the phase fractions remain at a certain level and there are no significant variations. The smallest thickness deviations occur in section 2. Increasing thickness deviations and periodical phase fractions reappear in section 3. In this section, a fraction of ferrite remains and bainite increases while simultaneously decreasing the fraction of pearlite. The phase fraction oscillations of bainite and pearlite become stronger in section 3, concurrently increasing thickness deviations.
At the outermost revolutions in section 4, the fractions of bainite and ferrite vary heavily while the fraction of bainite increases towards the strip end and the fraction of pearlite fades completely. The microstructure of the outermost revolutions is almost fully bainitic (Fig. 16 a) with a minor fraction of ferrite (Fig. 16 b).
The rolling force of the first (S1) and fourth (S4) cold rolling stands of a four-stand cold rolling mill are accountable for the averaged phase fraction of bainite as well as for averaged hardness (Figs. 17a and b). The rolling force deviation level of S1 is almost directly comparable with the bainite phase fraction, whereas the roll force of the fourth stand does not increase from the normalized length of 0.5-1.0 like at the first stand.

Discussion and conclusions
In this study a coupled heat transfer and phase transformations model for complete coil conveyance and coil field cooling is developed and validated with industrial measurements. The model includes all the steps of a conveyance chain for the modeled process configuration. The phase transformation model is fitted in experimental laboratory thermomechanical experiments and the effect of latent heat is included. Real coil geometry is used in order to obtain an accurate heat transfer model between the strip revolutions within the coil. Static post-coiling contact pressure distribution is implemented in the calculation of thermal contact conductance. Boundary conditions and model configuration are easily convertible for arbitrary simulations. Validation simulation is accountable for industrial data and the measurements prove the model dependability.
The required ferrite-bainite microstructure for a dual-phase steel is achieved by hot coiling and slow rate coil cooling. In this cooling process, the cooling rates are not managed in any way. The coil is fixed in a certain position and conveyed through the different conveyors to the coil field. This produces unequal cooling conditions on the opposite sides of the coil. In the above-mentioned results, these differential cooling rates lead to significant phase fraction variations in the coil and especially within a single revolution of the steel strip. Thermal contacts with coil conveyors accelerate the transformation of bainite on the bottom section of the coil. Conveyors are changed frequently depending on the production rate, making cooling at these contact interfaces more effective. Thus, the bainitic phase fraction is higher on the bottom section of the coil. Correspondingly the ferrite fraction is minor on the bottom section and major on the top section of the coil. Bainite phase fraction diminishes when moving inward the coil form outer shell. Simultaneously ferrite fraction increases. Though, oscillating phase fractions variations are clear on both ends of the strip. Emphasized fraction of bainite on the bottom section of the coil produces this oscillating phase fraction distribution. The outer and inner shells of the coil are principally bainitic and the body is ferritic with minor fractions of pearlite.
Because the mechanical properties between ferrite and bainite are different, the cold rolling parameters need to be adjusted accordingly for achieving a similar outcome. Periodical phase fractions of ferrite and bainite in an unwound steel strip cause unstable cold rolling process. This means that varying rolling forces on the oscillating phase fraction regions which leads to thickness deviations in the strip.
The main results of this work are: 1 According to the model results, the thickness deviations and the unsteady cold rolling process of a dual-phase steel are caused by heterogeneous phase fractions.
2 Phase fraction inequality is caused by the fixed position of the coil during the coil conveyance which leads to unequal cooling rates on a bottom contact region compared to the freely cooling top region of the coil. This causes different phase fractions on the opposite sides of the coil. These forming phases are mainly bainite, pearlite and ferrite and their fractions vary considerably during a single revolution in the coil, especially on the inner and outer shells of the coil.  3 The strength and hardness properties as well as load required for plastic deformation between bainite and ferrite phases are remarkably different and thus likely cause thickness deviations in the highspeed cold rolling process.
These variations are very difficult to take into account in the cold rolling process parameters in order to avoid thickness deviations, even if they could be modeled accurately. Instead, the coil cooling process should be modified to ensure uniform cooling throughout the coil. A possible solution to reduce the phase fraction differences could be rotating the coil during the coil conveyance so that the same section of the coil does not go through all the conveyor contacts. This creates challenges with the strip tail, which is usually appointed below the coil. Another suggestion for this problem would be an alternative cooling path to minimize the several thermal contacts in coil conveyance. Coiling temperature and chemical composition could be changed when any mechanical changes of the process are not required. Possible solutions for this issue will be discussed and simulated in further studies.

Data availability statement
The raw/processed data to reproduce these results cannot be shared because they include process data that are potential business secrets.

Author contributions
The corresponding author developed the FE model used in this study. The theoretical part of the manuscript as well as the experimental coil cooling measurements and modeling results were written by the  corresponding author. The phase transformation model was developed and introduced in this paper by Dr. Aarne Pohjonen. Sami Koskenniska assisted with experimental microstructure investigations to define phase transformation model parameters. All authors contributed to the final draft of this manuscript and gave approval for publication.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to affect the work reported in this paper.