{'fraction_non_alphanumeric': 0.02666312044948303, 'fraction_numerical': 0.011569533561621064, 'mean_word_length': 4.8566978193146415, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Internal oxidation and formation of Si/Al-enriched oxide bands in the scale of electrical steel grades

Abstract Oxidation during steel hot rolling is responsible for various surface defects. Local enrichments of oxygen-affine alloying elements such as silicon or aluminium can cause such defects by complicating oxide scale removal. In this paper, correlations between the temperature profile during oxidation and enrichment formation in electrical steels are investigated. Diffusion and reaction simulations using the numerical Crank-Nicolson-scheme are evaluated by comparing them to laboratory-scale oxidized samples. Furthermore, analytical approximations (spectral methods) are investigated as an alternative approach, focusing on trade-offs between accuracy and calculation times. A link between heating periods during oxidation and Si/Al-rich bands in the scale was established.


Introduction
During the production of flat steel products, the material is exposed to a series of subsequent rolling steps at high temperatures in air. Especially during hot rolling, the material suffers significantly from the formation of oxides which have to be removed in order to achieve the required excellent surface finish of the final product. Typically, the oxides formed at the steel surface (termed "scale") are removed by hydraulic descaling before rolling. However, certain (electrical) steel grades are known to form strongly adherent scales, resulting in the residual scale being rolled into the metal during the subsequent hot rolling step. One of these critical scale constituents is silicon [1][2][3], the effect of which on descaling motivates the present study. Silicon can increase the adherence between scale and steel substrate due to the formation of fayalite (Fe 2 SiO 4 ), the thermal expansion coefficient of which is closer to the steel substrate than to wustite (Fe 1− x O). Consequently, descalability in the presence of fayalite is deteriorated due to the lower thermal stress between the steel-fayalite boundary during hydraulic descaling [4]. Increased oxidation rates are also commonly observed during the processing of electrical steel grades [5]. Furthermore, Fe 2 SiO 4 forms a eutectic with Fe 1− x O at 1177 • C [6] which leads to the formation of a liquid oxide layer during slab reheating conditions. Even more in the presence of aluminium, a ternary eutectic with hercynite (FeAl 2 O 4 ) may form at a melting temperature as low as 1148 • C [6].
Describing oxidation phenomena, we will distinguish the following: • Internal and external oxidation as conventionally defined, depending on local oxygen partial pressure, temperature and alloy composition • Interior and exterior scale formation on a macroscopic scale, relative to the position of the original metal surface, as identified by marker experiments [7][8][9] Thorough reviews on the oxidation of hot rolled steels, including the effects of alloying elements, atmosphere and gas velocities, can be found in the literature [10,11]. It is reported that scale growth kinetics in dry atmospheres are governed by the ionic transport properties of cations (Fe 2+ , Fe 3+ ) and anions (O 2− ). For the case of pure iron [12,13], the cation transport through the scale is largely dominating, leading to an outward growth of the scale layer, as evidenced by Pt marker experiments [7]. This behaviour is consistent with the stoichiometry of the innermost iron oxide in the scale (wustite, Fe 1− x O, with x ranging between 0.04 and 0.17 [14,15]) showing cation-conducting properties. In the presence of water vapour, the anion transport increases, leading to an inward growth of the scale layer. Consequently, Pt markers are then found inside the (wustite) scale layer [7], implying that H 2 O plays a critical role in the oxidation behaviour of iron-based alloys. A wealth of literature is available on water vapour effects and scale morphology [10][11][12][16][17][18]. It is commonly observed that the addition of water vapour will increase the amount of interior scale formed, although the reasons remain open to debate. Comprehensive descriptions for the oxidation mechanism of dilute Fe-Si alloys can be found in the works of Birks et al. [12] and Sequeira [19]. Under the conditions of high oxygen partial pressure, high temperature and moderate alloying element content, external oxidation of Fe will occur, leading to formation of scale consisting mostly of iron oxides. Si and Al will not migrate into the scale to a significant extent. Instead, they will remain in the steel matrix and diffuse towards the steel/scale boundary whilst selective oxidation reduces the activities of non-oxidised silicon and aluminium before gradually incorporating these oxides into the interior scale. Therefore, local enrichment of these elements will depend on the inward transport of oxygen through the scale layer [20]. Experiments on Fe-1.5Si model alloys, heat treated at 800 • C to 1000 • C under humid atmospheres [8], indicate the formation of an interior Si-rich scale, as opposed to the exterior scale which was mostly composed of plain Fe-oxides. In contrast, oxidation in dry air often leads to the formation of protective layers which will inhibit the inward oxygen transport. Similar observations are reported for other alloy systems within the same temperature range [12]. At higher temperatures, the formation of a Si-rich oxide layer within the (interior) scale is essential to quantify the extent of oxidation of silicon-containing steel slabs prior to hot rolling.
Within this work, the case of a medium alloyed electrical steel grade, which does not form a permanent protective SiO 2 layer, was studied. The formation of the component oxides Fe 1− x O, SiO 2 and Al 2 O 3 as well as Fe 2 SiO 4 and FeAl 2 O 4 was investigated via heat treatments. In the following, experimental results are compared to thermochemical simulations of local phase distribution [21][22][23][24][25]20,26] and differences are critically discussed. The present work aims to provide a combined experimental-theoretical approach for the formation of Si/Al-rich oxide layers during the oxidation of Si/Al-alloyed steels under conditions similar to those in industrial gas-fired slab reheating furnaces.

Experimental procedure
Samples cut from the surface of an industrial grade as-cast steel slab with 2.4 wt.% Si and 1.1 wt.% Al were used for oxidation trials in a horizontal tube furnace. Preparation was performed in a fashion similar to related studies [2,3,27]. For each experiment, up to 10 rectangular samples measuring about 10 mm × 10 mm × 5 mm were ground with SiC paper down to 2000 grit size, resulting in an untarnished metal surface. The samples were then placed on top of fused alumina granulate inside a ceramic furnace boat, with the ground surface facing up. Oxidation was then performed in the furnace inside a mullite tube with an inner diameter of 35 mm using an atmosphere containing 20 % H 2 O (dew point +60 • C), 3 % O 2 and 7 % CO 2 by volume, with the balance being N 2 . This process atmosphere is typical for slab reheating furnaces. In order to prevent water vapour condensation inside the furnace and the gas tubing, the gas lines on the inlet side were kept at 80 • C and the furnace was first purged with 21 min − 1 dry argon (< 2 ppm O 2 , < 3 ppm H 2 O) during initial heating until a temperature of 400 • C was reached, as shown in Fig. 1. After a short holding period at 400 • C with 61 min − 1 Ar flowing through the furnace (typically 15 min), the actual oxidation was performed by flowing the oxidizing process atmosphere into the furnace tube at a rate of 6 l min − 1 . Theoretically, during heating and initial holding at 400 • C, the O 2 and H 2 O traces in the argon gas might cause selective oxidation of Si and Al. However, the effect was deemed to be negligible, since both the amount of oxidizing elements in the gas and the slow diffusion inside the steel at this low temperature were assumed to be severely limiting. During oxidation, all nominal heating and cooling rates were set to 10 K min − 1 . The experiments were performed at isothermal soaking temperatures between 1060 • C and 1160 • C, slightly below the melting point of the FeO-Fe 2 SiO 4 eutectic [6]. Finally, the samples were cooled to room temperature inside the furnace under Ar atmosphere and then removed. Again, the influence of Argon was assumed to be negligible during cooling because of the marginal O 2 and H 2 O contents. At an approximate oxygen partial pressure of 10 − 6 bar in the Ar gas used, only hematite would be reducible by thermal decomposition. While the slow cooling process leads to changes in the scale, e.g. the disproportionation of wustite, the enrichment of Si and Al was assumed to be mostly unaffected. Cooling below 750 • C became increasingly non-linear and the total cooling time was not measured directly, but was observed to be approximately 6.5 h according to experiments. A typical heating profile with a single oxidative heating and isothermal phase can be seen in Fig. 1. In addition, experiments using temperature profiles with an intermediate soaking period during heating as well as multiple heating stages with intermediate cooling were also performed for comparison with simulation results.
Following the oxidation, the whole samples were cold mounted under vacuum in a 10 • taper section, resulting in a magnification factor of about 5.8 in the direction of the angle. Metallographical preparation was then performed using a method based on previous work by Chen and Yuen [27,28]. The scale morphology and microstructure of the samples were investigated using light optical microscopy and SEM (FEI Quanta 200 Mk2), coupled with EDS (EDAX Octane Pro) for elemental analysis.

Mathematical modelling
Theoretical calculations of internal oxidation and local oxide formation were carried out using a two-step based algorithm, containing the diffusion of mobile elements (i.e. oxygen, silicon and aluminium) as well as their thermodynamic stability [24,26]. The diffusion is calculated in small time intervals, with the output serving as starting value to calculate the chemical equilibria for given local concentrations and temperature. Results are displayed as concentration profiles, showing both the nature and the molar fraction of each phase. Implementation of the algorithm was performed in Matlab. The data used for the simulations are shown in Table 1.
The diffusion was calculated via Fick's second law of diffusion, as shown in Eq. (1). It was assumed that the precipitation of oxides can be approximated by the diffusion coefficients of elements in the steel matrix, so that D (T) is constant in space. However, the temperature dependence of the diffusion coefficients was taken into account by assuming an Arrhenius dependency. As the diffusion coefficients change with temperature, the discrete time step size is varied as well in order to maintain accuracy of the diffusion calculation. In particular, to also account for different selections of the spatial step size between simulations, the value of r = Δt Δx 2 D (T) was held constant during each simulation.
∂c ∂t The simulations were performed for a finite steel slab, with the left boundary at the steel surface in contact with the oxidizing atmosphere, and the right boundary inside the steel bulk. Consequently, the concentrations of silicon and aluminium were set to zero at the surface and to the respective concentration for the simulated steel grade in the sample interior [25,26]. For oxygen, the concentration within the domain was set to zero, while at the surface, it was set to the solubility limit of oxygen in pure iron at the respective temperature and oxygen partial pressure at wustite decomposition. The length of the system was chosen so that the concentration of the alloying elements near the right boundary remained virtually constant. The solubility of oxygen in iron can be calculated by using Sieverts' law, where K S denotes Sieverts' constant [32]: The oxygen partial pressure was derived from the equilibrium coefficient for the wustite decomposition, which was calculated using FactSage [33]: The Si-and Al-containing oxides included in the model are Al 2 O 3 , SiO 2 , FeAl 2 O 4 and Fe 2 SiO 4 . For the reaction step, the model relies on a pre-defined reaction sequence based on thermodynamics instead of minimizing the total Gibbs free energy in each calculation step. These reaction equations are sorted in ascending order of Gibbs free energy change per mole of O 2 converted as shown in Table 2 and applied to all nodes, assuming complete reaction. This simplification is permissible for the oxides examined in this paper due to their low solubility products in iron and was deemed necessary, as the model would have been far more computationally intensive otherwise. The Gibbs free energies were calculated using FactSage [33] with the commercial oxide database FToxid. Although the temperature dependence of the Gibbs free reaction energies was included in the model, the reaction order did not change in the temperature region examined.

Numerical approximations using a Crank-Nicolson finite difference method
The partial differential equation for diffusion can be discretised according to the Crank-Nicolson-scheme [34,35]. The resulting linear equations have the form shown in Eq. (4), with the indices j denoting the spatial dimension and n for the time step. It is stable for all values of r, but requires inversion of the resulting tridiagonal coefficient matrix for each calculation step.

Analytical approximations using different functions
Deriving the solution of a diffusion equation can be computationally quite expensive, especially for high spatial accuracy (many meshpoints) as well as for large numbers of equations (i.e. chemical species). Instead of providing a numerical solution for the diffusion step via the Crank-Nicholson-scheme [34], an investigation to approximate the solution via Fourier series, Fast Fourier transforms as well as Chebyshev polynomials has been carried out. Here, the computational effort lies not within the inversion of large matrices but to derive the spatial frequencies of the corresponding sine and cosine waves for each species' concentration c( x → , t). These methods are also commonly referred to as "spectral methods" [36].

Fourier series
The decay of the initial concentration with time can be calculated directly via the analytical solution. For the one dimensional case of a finite slab of length L we obtain The discretized form of Eq. (5) used in the simulations is shown in Eq. (7), allowing the explicit calculation of the concentration profile after a time step Δt by a matrix multiplication. For a given spatial grid and ratio of r, the matrix K is constant, allowing it to be reused and therefore significantly reducing the calculation effort.
The solution c * can thus be approximated by the sum of the individual Fourier terms f * n in Eq. (5) up to the N th order and including the boundary conditions for diffusion at the surface, responsible for the concentration profile c inf after an infinite amount of time. For the chosen system geometry with the surface on the left and the steel bulk as the right boundary, these concentration profiles are simply linear functions between the respective boundary values.

Chebyshev polynomials
Applicable to a wide range of problems, Chebyshev polynomials offer an alternative to Fourier series, especially for non-periodic boundary conditions. They can be defined as the polynomials satisfying Eq. (10) [36]. As proper interpolation requires the spatial mesh to be unevenly spaced, the grid points were chosen according to the Chebyshev-Lobatto points in Eq. (11) [37], where j is the spatial index and J the number of grid points. Similar to the solution using the finite-difference-model, the solution c * (t + Δt) was calculated by applying the boundary conditions and solving the system of linear equations. For the comparisons below, the number of Chebyshev-Lobatto points was chosen so that the minimum distance of the grid is less or equal to that of the corresponding equispaced grid.

Fast Fourier Transform
Eq. (1) can be solved using the fast Fourier transform algorithm to calculate the discrete Fourier transform in space [37][38][39]. An important difference compared to the other methods examined here is that due to periodicity requirements of the FFT, the grid was expanded to be of even symmetry, therefore representing a slab of length L with contact to the oxidizing atmosphere on the left and right boundary. However, since the thickness L of the slab is very large compared to the diffusion lengths of all species, the concentration in the centre region will not change markedly, meaning that this arrangement can still be compared to the geometries chosen for the other methods. For the simulations, Eq. (12) with implicit Euler time stepping was used. As before, for a constant ratio of r, the factors applied to ĉ(t) remain constant as well. The main computational effort here lies in calculating the FFT for diffusion, followed by the inverse FFT for applying the reactions and boundary conditions for each time step.

Method comparison
The comparisons have been carried out for the required computation time and numerical accuracy of the aforementioned models. To keep the comparisons as concise and complete as possible, the root mean square error ϵ of all species' concentrations was formed according to Eq. (14). In order to make the errors of the concentrations of each species comparable, the concentrations were normalized using the maximum concentration of each phase in the reference profile. In the absence of an exact solution, the results of the finite-difference-model using the wellestablished Crank-Nicolson-scheme with a grid spacing Δx = 0.5μm and r = 0.5 were used as the reference c ref . Error summation was performed over the species s and the spatial grid points j, but not over time.
The simulations were performed for the isothermal diffusion and reaction of oxygen and silicon at 1160 • C for 3600 s, with a system length of 500 μm. The comparisons presented hereafter show the errors at the end of the simulation. Due to the different grid spacings, the concentration grids were linearly interpolated to a spacing of 0.25 μm for the comparisons.

Results and discussion
In both the simulations and the experiments, significant local enrichment of silicon and aluminium oxides could be observed. One of the first observations was that the majority of the enrichments was not located directly at the steel-scale-interface, but rather near the original metal surface, suggesting that most of the Si-and Al-oxides were formed during the initial stages of oxidation. The oxide phases formed during the experiments were also in agreement with thermodynamic expectations, with SiO 2 and Al 2 O 3 present as internal oxides inside the steel near the steel-scale-interface and the compound oxides Fe 2 SiO 4 and FeAl 2 O 4 found only embedded into the interior scale.

Comparison of the Crank-Nicolson based model with experimental results
Metallographic sections of samples, oxidised at 1160 • C for varying dwell times, are shown in Fig. 2. In all of the samples, the original steel surface is clearly visible, separating the exterior, brighter scale, containing mostly Fe-oxides, from the interior, darker scale, which consists of mainly wustite, mixed with Si-and Al-containing oxides. The exterior scale exhibited a high degree of porosity, especially for dwell times beyond 15 min. Iron oxides in this region were found to be largely columnar in nature, with a mostly dense oxide layer near the scale surface. A detailed view of the interior scale and interface region is shown in Fig. 3. Inside the interior scale, a darker oxide band can be seen just below the original metal surface. The thickness of this band corresponds to the interior oxide layer formed during heating from 400 • C to 1160 • C without an isothermal holding step, which can be seen in the sample with 0 min dwell time (0 min in Fig. 2). EDS analyses suggest that the darker oxide band is particularly rich in silicon and aluminium oxides.
Performing a simulation with the temperature profile shown in Fig. 1 also results in an increased Si and Al oxide content formed during the heating phase (Fig. 4a). The simulation predicted that altering the heating cycle would result in changes to the Si/Al oxide enriched band.
Introducing an intermediate holding step at 1060 • C for 30 min results in the enriched oxide zone being split by a less Si/Al rich region (Fig. 4b), which corresponds to the intermediate holding step . Fig. 4c shows the results for a sample where an additional heating phase was introduced by first heating to 1160 • C, cooling down to 910 • C without a holding step, heating to 1160 • C again and then holding for 60 min. Compared to the single-ramp heating profile in Fig. 4a, the additional band from the second heating step is clearly visible, while the first band is equivalent in both experiments.
There is a notable discrepancy between the thickness of the interior scale containing the Si-and Al-oxides in the simulations and the experiments. This can be explained by two main contributing factors. First, as mentioned before, experimental samples were prepared using taper section angles, resulting in a magnification factor of about 5.8. Secondly, the simulation did not include wustite, which forms the bulk of the interior scale. The morphology of enrichments, however, is accurately reproduced in the simulation.
For the dwell times examined in the present work, the elemental fluxes during isothermal oxidation remain mostly constant. However, during non-isothermal oxidation, the oxygen solubility and the diffusion coefficients will change significantly. While the activation energies for diffusion of silicon and aluminium are quite similar, they are about twice as high as that of oxygen, which results in their respective diffusion coefficients rising faster with temperature. This causes the transport of silicon and aluminium being temporarily at an advantage, as their diffusion profiles adjust to the changing diffusion distance until a new equilibrium is reached. The oxygen solubility in iron increases with temperature, counteracting the imbalance caused by the changing diffusion coefficients to some degree.

Comparison of models
One of the major differences observed between the finite-difference method and the analytical approximations lay in the effect of local discontinuities. For the analytical approximations, sharp changes in the concentration profile affected all other nodes during diffusion, due to the presence of high wavenumber terms. This typically leads to so called "ringing artefacts", where high frequency oscillations contaminate the concentration profile, which also caused slightly negative concentrations to occur near zero concentrations. This effect was especially noticeable during the initial stages of the simulation, where the discontinuities in concentrations were the largest (e.g. maximum oxygen concentration at the surface node, but zero concentration in all other nodes). While for pure diffusion calculations the smoothening effect of the diffusion can help to mitigate these artefacts, the reaction calculations lead to a continuous generation of such discontinuities. For example, in the case of oxygen, the concentration approximately decreases linearly from the boundary concentration to zero at the reaction front, resulting in a sudden deviation of the profile slope. Finitedifference methods did not suffer from this problem due to the local nature of the approximation, but this comes at the cost of decreased numerical accuracy. Fig. 5 shows a performance comparison of the methods with varying grid distances Δx. The Fourier series method based on the analytical solution produces very similar results compared to the Crank-Nicolson finite-difference-model. For larger grid distances, the calculation times are significantly shorter than for the Crank-Nicolson method, but this situation reverses for very fine grids, with the Fourier series eventually taking about ten times longer for a grid distance of Δx = 0.5 μm. While the Fourier series does not need to perform matrix inversion as the Crank-Nicolson method does, it uses a dense matrix multiplication for these parameters, which grows increasingly more complicated compared to the inversion of the tridiagonal coefficient matrix of the finite-difference-model. The FFT method offers no visible benefits compared to Crank-Nicolson, with the error being roughly twice as large and the calculation times being higher for most simulations. Finally, the Chebyshev polynomials offer the fastest calculation times, but also the largest errors for most values of Δx. It is likely that the reason for both of these phenomena is the uneven grid spacing, resulting in a highly reduced calculation effort, but also decreased accuracy. Comparison with the equispaced grid of the finite-difference-model also causes the irregular behaviour of the error, mainly depending on the coincidence of the grid points containing SiO 2 and Al 2 O 3 at the reaction front. Whereas computation time of all numerical methods can be reduced significantly by decreasing the spatial resolution, it should also be noted that one may end up with wrong results, especially when the grid distance Δx is in a similar range to the range of existence of an (oxide) phase. Fig. 6 compares the methods for different values of r for a grid spacing of Δx = 0.5 μm. Here, the Chebyshev method appears to be the most robust variant for higher values of r and therefore longer time steps Δt. This can be explained by the fact that only at the nodes with the smallest grid spacing the maximum value of r given in the figure is reached, while the other nodes have much lower values, making the accuracy of the diffusion calculation higher for the less densely populated internal nodes. As before, it is also the method with the shortest calculation times due to the smaller number of nodes. However, while the error remains fairly constant for the Chebyshev method, it is still relatively large when compared to the values of the other methods at smaller r. The Fourier series and Crank-Nicolson methods again produce very similar results, but because of the small grid spacing, the calculation times for the Fourier series method are higher by a factor of about ten due to the aforementioned difference in matrix density. The FFT method shows a large growth of the combined error with increasing values of r, likely caused by the weak implementation of the boundary conditions. As for the calculation times, no significant difference in the dependence on the diffusion parameter r could be observed. In particular, the slopes of the curves are roughly the same for all methods.

Conclusions
Oxidation of the 2.4 wt.% Si and 1.1 wt.% Al steel using a water vapour containing atmosphere (H 2 O/O 2 /CO 2 /N 2 (20/3/7, v/v/v)) mimicking the conditions found during slab reheating led to the formation of a finely porous scale with distinct exterior and interior parts separated by the original metal surface. Local enrichments of silicon and aluminium oxides were mostly observed as bands inside the interior  scale directly beneath the original metal surface. These oxide bands could inhibit descaling and thus become a cause of surface defects. Additionally, the internal oxides of silicon and aluminium inside the steel near the steel-scale-surface could also pose a problem during descaling and pickling.
The main insight gained from the simulations was the correlation of the Si-and Al-oxide rich bands with the heating phases of the oxidation process. This finding was substantiated by using different heating profiles and comparing experimental and simulation results, where each heating phase would correspond to an enriched oxide band, while isothermal holding would result in a lower concentration level. For the reheating process, this implies that the enrichment of alloying element oxides is strongest directly after reaching the maximum temperature, making the subsequent descaling more difficult. The simulations showed a good qualitative agreement in regard to the general shape of enrichments formed, despite the fact that the model does not take into account effects such as delamination, formation of wustite or the outwards transport of iron cations. The boundary condition for oxygen at the steel surface for the simulations is a theoretical maximum for internal oxidation. In reality, this concentration will be lower, as oxygen has to pass through the oxide scale first.
Upon comparison of the different methods for simulating the diffusion process, it was found that the Crank-Nicolson scheme delivered the overall best performance under the given parameters. However, for grid distances Δx ≥ 1 μm, the Fourier series method based on the analytical solution of the diffusion problem produced similar results with significantly lower calculation times. It is also likely that different hardware would result in different relative performance of the various methods. For pure diffusion processes, the analytical approximations (spectral methods) are significantly more accurate than finite-difference based models, but they are also far more susceptible to local discontinuities. The reactions will constantly produce such discontinuities during the simulation. The finite-difference-model, based on local approximation, will only be affected locally, while the global analytical approximations will also result in global errors in the solution.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Declaration of Competing Interest
The authors declare no conflict of interest.