Evaluation of 22MnB5 Steel Austenitization Sub-Models for Simulating the Heating Phase of Hot Stamping

The heating phase of hot stamping is the most time-and energy-intensive part of the process, and therefore offers significant potential to improve process efficiency. This requires a robust thermometallurgical model that predicts the blank heating profile and the austenitization progress inside a roller hearth furnace. In this paper, an IM austenitization model for 22MnB5 steel is derived using dilatometry data and evaluated against a JMAK type model using temperature profiles of coupons heated within a laboratory-scale muffle furnace and blanks heated within a roller hearth furnace. The austenitization rates predicted by the two models are then compared to the inferred austenitization processes from those temperature profiles. The result shows that the IM model has potential to capture the trend of austenitization process, especially the rapid pearlite dissolution at the beginning of the transformation, whereas the JMAK type model is too restrictive and fails in this regard.


Introduction
Automakers increasingly rely on ultra-high strength steel (UHSS) parts made through hot stamping for automotive lightweighting. In the direct hot stamping process, blanks are first austenitized, usually in a roller hearth furnace, before they are transferred to a press to be formed and simultaneously quenched into fully martensitic components. Most often the parts are made of 22MnB5 alloy with an Al-Si coating to prevent oxidation and decarburization within the furnace. The coating melts at around 577 ℃ [1], and then transforms into an Al-Si-Fe intermetallic layer that provides long-term corrosion protection. The furnace parameters, consisting of zone temperatures, roller speed, and other ancillary parameters, must be adjusted to ensure full austenite transformation and proper coating transformation while minimizing cycle time, i.e. the time to austenitize one batch of parts. This is often done heuristically based on experience but considering that a typical roller hearth furnace contains more than 15 independently-set heating zones, in addition to the roller speed and batch layout, the outcome is rarely optimal. This is further complicated by the fact that the blank heating rate is largely determined by the radiative properties of the coating, which change dramatically as the coating evolves. Accordingly, there is a need for an accurate thermometallurgical model of the blank transformation within the furnace in order to optimize the heating stage of hot stamping.
Austenitization is a complex multi-stage phase transformation that is affected not only by the temperature of the blank but also its heating rate [2]. Verma et al. [3] recently developed a roller hearth International Deep-Drawing Research Group (IDDRG 2020) IOP Conf. Series: Materials Science and Engineering 967 (2020) 012076 IOP Publishing doi: 10.1088/1757-899X/967/1/012076 2 furnace model that consists of a heat transfer sub-model coupled with an austenitization sub-model. This complete thermometallurgical model can predict both the blank temperature history as it travels through the furnace, and the gradual austenitization of the steel. Verma et al. compared two different austenitization sub-models: a first-order kinetics model proposed by Di Ciano et al. [4] that has two degrees-of-freedom; as well as a more detailed physical model proposed by Li et al. [5], which considers nucleation, growth and impingement and involves 15 degrees-of-freedom. In both cases, the model parameters were identified through nonlinear least-squares regression with dilatometry data. Verma et al. concluded that the model by Li et al. was not only susceptible to overturning but also limited to intercritical annealing as opposed to full annealing, which is the objective in hot stamping. This is also supported by Yan et al.'s work [6], which examined the Li et al. [5] model in more detail and concluded that it performs better at temperatures closer to onset temperature, TAc1, for shorter soaking times, and significantly underpredicts the austenite fraction, fγ, at higher soak temperatures.
Zhu et al. [7] proposed a highly physical model using the 2D Cellular Automation method. This model simulates the nucleation and grain growth of austenite in pearlite, the growth of austenite in ferrite, and the grain coarsening of austenite. However, the calculations are computationally intensive to carry out, with the simulation of a 10-minute austenitization process taking over a full day to compute and is therefore unsuitable for furnace process optimization. Huang et al. [8] developed a non-isothermal austenitization model using the Internal State Variable (ISV) method, which includes an incubation period and a transition period. It is capable of predicting both the starting and finishing temperature of the austenitization process, as well as fγ during the process.
Luo et al. [9] proposed a "model free" approach to predicting the austenitization kinetics for SA508 Gr.3 steel, using the Isoconversional Method (IM). "Model free" means the model does not presume a specific model parameterization. In particular, instead of deriving a set of constant model parameters, activation energy (Ea) and pre-exponential factor (A) for the entire process, these parameters are evaluated at different degrees of conversions, i.e., Ea(fγ), and A(fγ) based on dilatometry data measured for heating rates between 0.02 K/s and 500 K/s. An important advantage of such approach is that it can model complex conversion processes while avoiding oversimplified model assumption.
This paper evaluates the suitability of the IM approach for modelling austenitization of 22MnB5 steel in hot stamping. The model is first derived from a set of dilatometry measurements carried out on inductively heated 22MnB5 specimens. The predictive capability of the model is first assessed by comparing it to the dilatometry data set that was used in the derivation, after which, the IM model is then integrated into the roller hearth furnace model and evaluated using temperature measurements of coupons heated within a laboratory-scale muffle furnace and blanks heated within a roller hearth furnace.

Roller Hearth Furnace Model
The furnace model derived by Verma et al. [3] is applied to predict the heating rate of an unpatched, singe-gauge blank. Normally, ceramic rollers convey the blanks through the furnace. The blanks are heated using natural gas fired radiant tube burners, which may be located above and below the ceramic rollers. The furnace is separated into zones, each of which typically includes multiple burners and a single control thermocouple. The burners in each zone are simultaneously activated or deactivated via a hysteresis control system based on the temperature measured by the control thermocouple relative to the corresponding zone set-point temperature. Each zone is assumed to be large and isothermal relative to the blank, and the blank is assumed to be irradiated equally on both the top and bottom surface. Due to the small thickness, the blank temperature is assumed to be uniform. Since the contact area between rollers and the blank is small, conduction heat transfer from the rollers is excluded from the model. The air within the furnace is assumed to be quiescent [3], therefore the dominant convection mode is natural convection. With these assumptions, the blank heating curve is governed by Equation (1) where ρ and cp denotes the density of the blank and specific heat, respectively [10], ∆hγ is the latent heat of austenitization, dfγ/dt is the instantaneous austenitization rate. Qconv is the total convection obtained from where Tamb is the ambient temperature, Ts is the blank surface temperature, i.e., the blank temperature, and the area A accounts for both the upper and lower blank surfaces. The average convection coefficients for the top (h̅ top) and bottom (h̅ bot) are calculated from natural convection correlations [11], and range approximately from 1 to 10 W/(m 2 K) for the top surface and 1 to 17 W/(m 2 K) for the bottom surface, becoming smaller as the blank temperature approaches that of the surrounding furnace air. The net radiative transfer to the blank, Qrad is the energy received through radiation heat transfer and is given by where σ = 5.67 × 10 -8 W/(m 2 •K 4 ) is the Stefan-Boltzmann constant, and α and ε are the total absorptivity and emissivity. The spectral emissivity is plotted in Figure 1 (a) and the derived total absorptivity and emissivity by Jhajj et al. [12] is shown in Figure 1 (b) . Both the total emissivity and total absorptivity drop at 575 ℃ due to the melting of the Al-Si coating, and then rise as the intermetallic layer evolves and becomes progressively rougher. Changes in surface morphology have great influence on the radiative properties of the material, although the exact relationship remains unclear [1,13,14]. The metallurgical sub-model is used to obtain the austentization rate, dfγ/dt, based on the blank temperature calculated from the heat transfer model. The metallurgical sub-model, in turn, is coupled with the heat transfer model through the latent heat of austenitization. The latent heat of austenitization for 22MnB5 has not been characterized in literature. Verma et al. estimated it to be 30 kJ/kg based on a rule of mixture using ∆hγ values for fully pearlitic and fully ferritic iron [3] and taking the as-received state of 22MnB5 being approximately 80%-ferrite/20%-pearlite. Equation (1) is solved simultaneously with a differential equation that governs the rate of austentization, dfγ/dt. In this work the explicit Euler technique is used to solve these equations.

F1 model
Khawam et al. [15] reviewed many solid-stated kinetics models, including the first-order kinetics model (F1), which is the Avrami (or Johnson-Mehl-Avrami-Kolmogorov, JMAK) equation, with n = 1. Di-Ciano et al. [4] modified the F1 model Equation to account for a non-isothermal heating rate Taking the derivative of Equation (4)  where R is the universal gas constant (8.314J/molK), T(t) is the current temperature of the blank in K and g(fγ) is an intermediate variable The rate of austenitization is then given by The two model parameters are the pre-exponential factor, A, and activation energy, Ea, which Di Ciano derived from dilatometry data to be 3.42  10 18 s -1 and 402 kJ/(molK), respectively, for heating rates between 1-5 K/s.

Isoconversional Method (IM) model
Luo et al. [9] modelled the austenitization process of SA508 Gr.3 steel using an isoconversional method (IM), specifically the Friedman method, which is commonly used in thermal analysis but not for solidstate phase transformation kinetics. The general equation for a solid-state phase transformation is where fγ is the austenitization fraction, and f(fγ) represents the specific austenitic transformation model, such as the F1 model, which is f(fγ) = 1 − fγ according to Khawam et al. [15]. Instead of having a constant activation energy (Ea) and a pre-exponential factor (A) for the entire transformation process, in the IM approach those parameters are evaluated as the reaction progresses, i.e. at different specific fγ values. At a specific fγ, f(fγ) is constant, and Equation (8) can be rearranged into Since A(fγ) and f(fγ) are both constant at this specific fγ value, they can be combined into a modified preexponential factor, A'(fγ), and after taking natural logarithm Equation (9) where β represents the constant heating rate. A total of nine sets of dilatometry data with three repeats per heating rate (10, 100, 200 K/s) are used to derive the model parameters. The dilatometry data is measured with a 1.6 mm thick Al-Si coated 22MnB5 (90 g/m 2 coating weight) steel sample on a DIL805 A/D thermal-mechanical simulator. The sample composition is obtained through optical emission spectrometry (OES) and is listed in Table 1: The austenite fraction fγ is interpolated from the dilatometry curve using a lever rule procedure detailed in Huang et al.'s work [16] and is shown in Figure 2. The derivative of austenite fraction with respect to temperature, dfγ/dT, is calculated numerically from Figure 2 and the result is in Figure 3.    Figure 3 show that the austenitization fraction appears to decrease before rising again at the beginning of the process, which is a measurement artifact caused by induction heating used by the dilatometer. As-received 22MnB5 is a mixture of ferrite and pearlite, which are magnetic phases, whereas austenite is nonmagnetic [17]. Austenitization happens quickly at the beginning stage with pearlite dissolution [16], causing the DIL805 controller to respond rapidly, which result in an overshoot of HF power for a short period of time until the controller stabilizes. In addition to the effect of temperature, the spike in HF power also triggers a magnetostriction effect that contributes to a minor change in length of the sample. Figure 2 and Figure 3 also show that the third dataset of 100 ℃/s is an outlier and is not used to derive the model.
The process for obtaining Ea(f) and A(f) curves is summarized in the flow chart shown in Figure 4. Data pairs of [ln(dfγ/dt), 1/T] at this fγ from each of the eight remaining datasets are plotted and a line is fitted through these points, as shown in Figure 5 for a subset of sampled f values (f = 0.25, 0.51, and 0.75). The actual fitting process follows a much more stringent sampling of f values with different spacings at different ranges of f, considering that the transformation rate is very rapid at the beginning and slows down afterwards, to have a reasonably evenly distributed sampled points. Figure 6 provides an example of how the 100 K/s heating rate dataset No.2 is sampled.   With the derived model parameters Ea(fγ) and A'(fγ), the model prediction can be made by solving differential Equation (9) using the explicit Euler scheme. At each timestep, ti, the austenite fraction fγ,i can be used to linearly interpolate its corresponding effective activation energy, Ea(fγ,i), and preexponential factor, A(fγ,i). Given a small enough change in austenite fraction, the two parameters are assumed to be constant within this timestep, and the transformation rate can be obtained through

IM Sub-model Validation against Dilatometry Data
The temperature history of each dilatometry dataset is incorporated into the IM model to predict the austenitization process, and the model error is calculated by subtracting the data value at this temperature from the model prediction. Figure 8 shows an example of the model validated against the first run of 10℃/s dilatometry measurements; model prediction is in good agreement with the dilatometry data. This is further demonstrated in the error plot in Figure 9. The highest error is under 0.1 and tends to occur at the onset of austenitization. It may be due to the measurement artifact mentioned previously, which causes the model parameter to be missing for the first 5% of austenitization process. In practice, the austenitization fraction is set to fγ,i = 5% when fγ,i ≤ 5% to retrieve the isoconversional values of the two model parameters.

Comparison of F1 and IM sub-models for predicting austenitization in furnaces
The F1 and IM austenitization sub-models are combined with the heat transfer model and are used to predict the heating curves of coupons heated within a laboratory-scale muffle furnace and blanks heated in a roller hearth furnace. The heating curve predictions are compared with thermocouple measurements. The austenitization rate is also inferred from the thermocouple measurements by rearranging Equation (1) to solve for dfγ/dt, which is then compared directly with result obtained from austenitization sub-models.

Muffle Furnace Measurements.
A muffle furnace set to 910 ℃ is used heat Al-Si coated 22MnB5 coupons having dimension of 130  30  2mm, which is mounted in a ceramic holder as described in Ref [3] to minimize conduction heat transfer between the blank and the holder. A K-type thermocouple is welded to the centre of the blank and measured using NI -9213 ® data acquisition unit. The roller hearth furnace model is adapted to this scenario by modeling the muffle furnace as a single heating zone in a roller hearth furnace. The F1 and IM models are activated at the same prescribed TAc1, which is apparent by the drop in heating rate. A comparison of the measured temperature profile and predicted temperature profiles using both the F1 and IM models is shown in Figure 10. Both kinetics sub-models generate temperature profiles that closely match the experimental measurements. Figure 11 shows a more informative comparison between the predicted fγ by both kinetics submodels and the inferred fγ obtained using equation (13). Both the IM and F1 models underpredict the austenite fraction. The austenitization transformation rate is initially very large due to rapid pearlite dissolution, as suggested by both Huang et al. [16] and Li et al. [5]. This is captured by the IM model but not the F1 model. Although inferred, F1 and IM austenite transformation are all set to start at the same onset temperature, there is a delayed start with IM model. This could be related to the missing model parameters due to the measurement artifact. Manually shifting the IM curve forward by 10 K shows a better agreement between it and the inferred result. With more data and a resolution to the measurement artifact, the IM model shows potential for providing an accurate instantaneous austenite phase fraction.

Roller Hearth Furnace
Results. Next, the models are analyzed using temperature measurements made on instrumented blanks heated within a roller hearth furnace. A 625  297  1.6 mm 22MnB5 steel blank with 90 g/m 2 Al-Si coating weight is sent through the roller hearth furnace instrumented with K-type thermocouples and a Datapaq ® unit. The measured temperature profile is compared to the F1 and IM model predictions in Figure 12. Qualitatively the model predictions match the experimental result, although with much less agreement than the muffle furnace results shown in Figure 10. This may be due to the idealization of a constant and uniform zone temperature, while in reality the zones are not isothermal and the average temperature fluctuates about a set-point [3]. Moreover, the radiative properties shown in Figure 1 change with the evolving surface state, which is a complex function of heating rate and is not yet completely understood [13,14].
Inferring the austenite phase fraction using the thermocouple data via equation (13) initially produces a non-physical result shown in Figure 13, because the change in sensible energy inferred from the blank temperature, ΔEsensible, is larger than the heat transfer to the blank predicted by the model, Qnet. This issue likely arises from uncertainties in the radiative properties, which change abruptly over this temperature range as shown in Figure 1 due to coating evolution. In their muffle furnace experiments, Jhajj et al. [12] circumvent this problem by removing the Al-Si coating and applying a protective boron nitride (BN) coating, which remains inert at high temperatures, but this cannot be done for an industrial roller hearth furnace. Instead, the absorptivity value at 700 and 800 ℃ are heuristically increased by 10% and the modified result of net heat transfer is shown Figure 13.     With the modified absorptivity, the inferred fγ is compared to the F1 and the IM prediction in Figure  14, and the rate of transformation comparison is shown in Figure 15. Similar to the muffle furnace result, the F1 model predicts a much lower austenitization rate close to TAc1 compared to the IM model. The IM model, however, still shows a delayed start despite the model starts at the same austenitization onset temperature. Manually shifting the IM curve forward by 5 K shows a better agreement.

Conclusions
Maximizing the profitability of hot stamping operations involves identifying the optimal set of furnace parameters that minimize cycle time while achieving the functional requirements of the heating stage, i.e. austenitization and coating transformation. This, in turn, demands accurate models for both the heat transfer and metallurgical transformations undergone by the blank. This paper derives a new type of austenitization kinetics sub-model for 22MnB5 alloy, using the Isoconversional Method (IM), from dilatometry measurements. The model is incorporated into a thermometallurgical blank model, and the performance is compared with a first-order kinetics model using instrumented Al-Si coated 22MnB5 coupons heated in a muffle furnace and instrumented blanks heated within a roller hearth furnace.
Both the F1 and IM models under-predict the austenite fraction compared to values inferred from thermocouple measurements, although the inferred austenite phase fraction is subject to significant uncertainties due to the uncertain radiative properties of the evolving coating. However, the IM method has recreated the rapid austenitization associated with pearlite dissolution close to TAc1. This result shows that although F1 model is robust and easy to use, it does not have the necessary degrees-of-freedom to capture the physical transformation process. In contrast, the IM method captures the general trend better, although more dilatometry data with a wider range of heating rates and smaller intervals should be used to derive a more robust set of parameters. Lower heating rate data, in particular, would be more representative of conditions within a roller hearth furnace. Other candidate phase transformation n models, such as the ISV model, should also be evaluated as candidates for modeling austenitization within roller hearth furnaces.
Future work should also focus on reducing the uncertainty in radiative properties after the coating melts. The surface morphology and properties evolve as the coating transforms from its as-received Al-Si state to its final intermetallic form, which in turn drastically changes radiative properties. Experiments should be conducted where samples with coatings at different stages of the evolution are preserved and examined to uncover the relationship between coating evolution and changes in radiative properties.