Simulating the Growth of Dual-Phase Boride Layer on AISI M2 Steel by Two Kinetic Approaches

: Two kinetic approaches (integral method and Dybkov method) have been applied for simulating the boriding kinetics of AISI M2 steel in the range of 1173 to 1323 K, by including the effect of incubation periods. For the integral method, a peculiar solution of the resulting system of differential algebraic equations (DAE) has been employed for assessing the diffusivities of boron in FeB and Fe2B . The boron activation energies in FeB and Fe2B have been deduced from both approaches and compared with the data taken from the literature. Furthermore, to experimentally extend the validity of both approaches, four additional boriding conditions obtained on the boronized samples at 1173, 1223, 1273 and 1323 K for 10 h were then used. The predicted boride layers’ thicknesses were confronted to the experimental values. Consequently, a satisfactory concordance was obtained when comparing the simulated layers’ thicknesses to the experimental values derived from the literature.


Introduction
Boriding treatment is a method of enrichment of the surface of steels or Armco iron by a thermal diffusion of atomic boron from boron-containing media [1] between 800 and 1050 • C for 0.5 to 10 h. This thermochemical treatment results in improvement of tribological and anticorrosion properties of the treated parts. An extremely high surface hardness obtained on ferrous alloys is ascribed to the formation of boride phases [2,3].
As per the Fe-B binary phase diagram [4], two iron borides could be formed at equilibrium condition. The Fe 2 B phase containing 9 wt.% of boron crystallizes in tetragonal system while for the FeB phase, its crystal structure is orthorhombic and having 16.23 wt.% of boron [5]. Practically, different variants of boriding processes can be performed on the basis of physical state of the active components: gas [6], liquid [7], electrochemical [8], paste [9], powder [10], fluidized bed reactor [11], plasma [12] and plasma-paste [13]. However, the most frequently employed process is the powder-pack boriding owing to its ease of implementation and lower economic cost compared to other boriding processes [14].
The AISI M2 steel is employed for manufacturing the cutting tools and punches. But for applications requiring extreme wear conditions, this steel should be borided to enhance its wear resistance.
Kinetically, different kinetic approaches have been established to control the growth kinetics of compact layers composed of Fe 2 B [15][16][17][18][19][20][21][22] or (FeB and Fe 2 B) [23][24][25][26][27][28][29][30][31][32][33][34][35] and formed on ferrous alloys. All these mathematical models help to select the optimum values of boride layers' thicknesses according to the requirements imposed by the utilization of these ferrous materials in the industry. For example, in [23], a comprehensive kinetic model was used for the bilayer growth (FeB and Fe 2 B) on AISI M2 steel by assuming a linear concentration of boron profile with the inclusion of incubation periods. In another paper [24], two kinetic models have been suggested for studying the growth of FeB and Fe 2 B layers on AISI D2 steel between 1223 and 1273 K. The first approach was referred to as the mean diffusion coefficient method and the second one the Dybkov method. In this latter approach, the temperature-dependent parameters for the growth of both layers k FeB and k Fe 2 B were fitted with the experimental data by using the Arrhenius relationships. In other recent study, Keddam and Kulka [25] have modelled the growth of the bilayer (FeB/Fe 2 B) on AISI M2 steel with the integral method and Dybkov model. For the integral method, a peculiar solution of the obtained DAE system was used to assess the values of boron diffusion coefficients in FeB and Fe 2 B and by assuming a parabolic shape for the boron concentration profile in each layer. In other study, Dybkov et al. [28] have suggested a mathematical model consisting in studying the growth of bilayer (FeB/Fe 2 B) on Fe-Cr binary alloys. The developed kinetic approach called the Dybkov model required two fitting parameters in order to be in line with the experimental layers' thicknesses. Delai et al. [29] have used the classical diffusion models to follow the boriding kinetics of 4Cr5MoSiV1 steel grade between 860 and 980 • C by getting a bilayer configuration (FeB/Fe 2 B) and neglecting the existence of incubation period during its formation. The same kinetic approach was also used by Belguendouz et al. [34] to investigate the growth of FeB and Fe 2 B layers on AISI D2 steel by considering the boron diffusivity in the iron phase and by ignoring the influence of boride incubation periods on the kinetics.
Zouzou and Keddam [35] have extended the use of the integral diffusion model for FeB, Fe 2 B and diffusion zone formed on AISI 316 steel with the presence of boride incubation periods.
In this research paper, the integral method [25,26] with a new variable changes and the Dybkov model [25,28] were suggested for simulating the boronizing kinetics of AISI M2 steel in the range of 1173 to 1323 K by using the experimental results [23].
These two kinetic approaches have also been checked empirically by confronting the experimental layers' thicknesses to the simulated values. For the integral method, new relations for boron diffusivities in FeB and Fe 2 B have been proposed based on the new variable changes. For the Dybkov model, the time evolution during of FeB and Fe 2 B layers have been described by the set of ordinary differential equations of first order and using two fitting parameters (k FeB and k Fe 2 B ) with the aim of reproducing the experimental layers' thicknesses. The determined activation energies in AISI M2 steel were contrasted with the literature data. Finally, the mass gain in each iron boride layer was also evaluated in dependence of exposure time that is greater than the boride incubation period.

Integral Diffusion Model
This kinetic model [25,26] has been used to predict the layers' thicknesses of FeB and Fe 2 B formed on the surface of AISI M2 steel in case of a saturated matrix with boron atoms. Figure 1 shows in a schematic way the change in the boron concentrations across the FeB and Fe 2 B layers after a treatment time greater than the incubation period.   The terms C FeB up (=16.40 wt.%B) and C FeB low (=16.23 wt.%B) are the maximum and minimum boron contents in FeB [23,24,31]. C Fe 2 B up (=9 wt.%B) and C Fe 2 B low (=8.83 wt.%B) stand for the maximum and minimum boron contents in Fe 2 B [31]. C ads is the adsorbed concentration of boron [26]. The variable u(t) is the location of the (FeB/Fe 2 B) interface and v(t) that of the second interface (Fe 2 B/substrate). The term C 0 stands for the solubility limit of boron within the matrix (=35 × 10 −4 wt.% B) [3].
The FeB layer thickness u(t) is expressed by Equation (1), see [23,24,30,32], for instance: where k 1 is the kinetic constant at the first interface with a boride incubation time t FeB 0 . The thickness of (FeB + Fe 2 B) layer v(t) is expressed by Equation (2): where k is the kinetic constant at the second interface with an incubation time t 0 . Mathematically, it is possible to re-write the Equation (1) in the following form: The thickness of Fe 2 B layer is the difference (l = v − u) given by Equation (4): The initial conditions imposed by the integral diffusion model are expressed by: The boundary conditions of the same model are summarized in the following equations: C Fe 2 B x[t = t 0 ] = 0} = C Fe 2 B up ; for 8.83wt.%B ≺ C ads ≺ 16.23wt.%B and without FeB phase low for C ads ≺ 8.83wt.%B and without FeB phase (9) The change in the boron concentrations with the two variables x and t are expressed by Equations (14) and (15) based on Goodman's method [36]: The following time-dependent variables a 1 (t), b 1 (t), a 2 (t), b 2 (t), have to verify the boundary conditions. The integral method uses the DAE system (i.e., Equations (16)-(21)) for modelling the growth kinetics of FeB and Fe 2 B layers: 2w 12 and: with: In order to find the relations for boron diffusivities in FeB and Fe 2 B, the variables change provided by Equations (22)-(27) were used. and: The two constants ε 1 and ε 2 stand for the non-dimensional parameters. D FeB and D Fe 2 B are the boron diffusivities in FeB and Fe 2 B. The constants α 1 , β 1 , α 2 and β 2 must be positive, and subjected to the boundary conditions. After substitution of these last equations into the obtained DAE system, and derivation with respect to the time variable, the following set of non-linear equations was deduced: After a numerical resolution of the obtained system of non-linear equations, the parameters α 1 , β 1 , α 2 and β 2 can be determined. The two non-dimensional parameters ε 1 and ε 2 are given by Equations (32) and (33): and: The relations of boron diffusivities in FeB and Fe 2 B are finally expressed by Equations (34) and (35): and: Equations (34) and (35) were obtained by comparing Equations (3) and (4) with Equations (26) and (27).

Dybkov Model
The Dybkov model [25,28] has been used to model the kinetics of boron diffusion through the surface of AISI M2 steel by including the effect of boride incubation periods. At the contact interface (solid/gas), the gaseous phase resulting from the chemical reactions occurred in the powders mixture allows the release of active boron. Afterwards the further growth of boride layers is being controlled by the diffusion stage at the expense of boron element [28]. For the sufficiently compact and continuous layers, Equations (36) and (37) are then used to model the layers' growth in dependence of treatment time at a given boriding temperature: with u(t) and l(t) the respective thicknesses of FeB and Fe 2 B layers. The constant g (=0.60) depends on the molar volume of each phase. The constants p = 1, q = 1, r = 1 and s = 2 are the stoichiometric coefficients of FeB and Fe 2 B [28]. At the further stage, the growth process of FeB and Fe 2 B layers will continue at the considered interface under the two partial chemical reactions [28] over a prolonged time duration: B + Fe 2 B = 2 FeB and FeB + Fe = Fe 2 B The two growth parameters k FeB and k Fe 2 B , which are provided by Equations (38) and (39) [25], can be obtained from the fitting of experimental data [23].
Equations (38) and (39) were derived by substituting Equations (3) and (4) and their derivatives with respect to the time into the ODE system (Equations (36) and (37)). Finally, the obtained ODE system can be numerically solved by using the Runge-Kutta method [37] with the initial conditions u 0 and l 0 at a given incubation period t 0 .

Calculation Results and Discussion
The experimental results [23] were utilized for performing the simulation of diffusion kinetics of boron through the boride layers in case of AISI M2 steel. In this research paper, the solid method of boriding was used to generate a dual-phase boronized layer (FeB + Fe 2 B) on the surfaces of treated AISI M2 steel between 1173 and 1323 K for 4-8 h.
In order to validate their diffusion model related to the continuity equations, the samples were also borided at each temperature for a time duration of 10 h. The AISI M2 steel had the chemical composition of: 0.85% C, 0.30% Mn, 4.5% Cr, 5% Mo, 1.950% V, 6.4% W and Fe balance. The boron source, used during their experiments, was the commercial boriding agent called Durborid. In the present work, two kinetics models have been proposed to predict the kinetics of boronizing of AISI M2 steel when forming a bilayer (FeB + Fe 2 B). For the integral diffusion model, a particular solution of the DAE system was used to assess the boron diffusivities in FeB and Fe 2 B by taking 16.40 wt.% as a maximum boron content in FeB.

Assessment of Boron Diffusivities in FeB and Fe 2 B with the Integral Method
First of all, obtaining the values of boron diffusivities FeB and Fe 2 B are required to simulate the boride layers' thicknesses for the given conditions. The experimental results [23] in terms of thicknesses were exploited to check the regime of growth laws of FeB and (FeB + Fe 2 B). According to [23], the measurements of layers' thicknesses of FeB and (FeB + Fe 2 B) were carried out by using an optical microscope. An arithmetic mean of eighty measurements was taken on different cross-sections of boronized specimens, and the obtained values were compared to the simulated ones. The simulated layers' thicknesses were in accordance with the experimental values, i.e., the simulated values are reliable. Table 1 gives the experimental kinetics constants for the two interfaces with the corresponding boride incubation periods [23]. These values were extracted from the slopes of the straight lines showing the dependencies in time of u 2 (t) and v 2 (t) and fitted with Equations (1) and (2). The boride incubation periods were obtained from a zero layer thickness. It was noticed that the incubation periods decreased with increasing boriding temperatures [23,30] due to the thermal activation of boron diffusion at the material surface. The application of integral method also requires the determination of the new values of experimental kinetics constants for the first interface corresponding to the boride incubation periods of the entire boride layer. In Table 2 are displayed the new kinetics constants at the (Fe 2 B/FeB) interface fitted with Equation (3). In Table 3 are grouped the calculated values of boron diffusivities in FeB and Fe 2 B and those of the non-dimensional parameters (ε 1 and ε 2 ) with the help of the integral method (see Equations (34) and (35)) between 1173 and 1323 K and for C FeB up = 16.40 wt.%. It is noticed that the numerical values of ε 1 and ε 2 are nearly constant. In Figure 2 are plotted the dependencies in temperature of calculated boron diffusion coefficients in FeB and Fe 2 B with the integral diffusion model. By adopting the Arrhenius relationships, Equations (40) and (41) were derived: with R denotes the ideal gas constant (=8.314 J mol −1 K −1 ) and T the process temperature in Kelvin.  In Figure 2 are plotted the dependencies in temperature of calculated boron diffusion coefficients in FeB and Fe2B with the integral diffusion model. By adopting the Arrhenius relationships, Equations (40) and (41) were derived: with R denotes the ideal gas constant (=8.314 J mol −1 K −1 ) and T the process temperature in Kelvin. The activation energies for boron diffusion in FeB and Fe2B are found to be 206.41 and 216.18 kJ mol −1 by means of Equations (40) and (41). These values of energy are found to be higher. It is ascribed to the alloying elements present in AISI M2 steel which limits the diffusion rate of boron atoms leading to a reduced thickness of boronized layer. With the presence of alloying elements, the interface (boride layer/transition zone) flattens in case of high alloy steels [23,30]. The activation energies for boron diffusion in FeB and Fe 2 B are found to be 206.41 and 216.18 kJ mol −1 by means of Equations (40) and (41). These values of energy are found to be higher. It is ascribed to the alloying elements present in AISI M2 steel which limits the diffusion rate of boron atoms leading to a reduced thickness of boronized layer. With the presence of alloying elements, the interface (boride layer/transition zone) flattens in case of high alloy steels [23,30].

Estimation of the Two Temperature-Dependent Parameters k FeB and k Fe2B Using the Dybkov Model
The two fitting parameters k FeB and k Fe 2 B for the Dybkov model were calculated by using Equations (38) and (39) and considering the values of kinetics constants k and k listed in Tables 1 and 2 with respect to the boride incubation periods t 0 . The calculation results are given in Table 4. As a second approach, the Dybkov model was used to estimate the values of boron activation energies in both layers (FeB and Fe 2 B). The two temperature-dependent parameters k FeB and k Fe 2 B were determined by using Equations (38) and (39) in order to match the experimental data [23] with the simulation results. Figure 3 gives the change in temperature of the calculated kinetics constants k FeB and k Fe 2 B based on the experimental data [23]. After a fitting of the obtained values in terms of the squares of kinetics constants with the Arrhenius relations, Equations (42) and (43) were derived: Coatings 2021, 11, x FOR PEER REVIEW 9 of 15 ' k listed in Tables 1 and 2 with respect to the boride incubation periods 0 t . The calculation results are given in Table 4.  were determined by using Equations (38) and (39) in order to match the experimental data [23] with the simulation results. Figure 3 gives the change in temperature of the calculated kinetics constants FeB k and B Fe k 2 based on the experimental data [23]. After a fitting of the obtained values in terms of the squares of kinetics constants with the Arrhenius relations, Equations (42) and (43)

T (K)
where R denotes the ideal gas constant (=8.314 J mol −1 K −1 ) and T the boriding temperature given in Kelvin. Table 5 gives a comparison between the reported values of boron activation energies for some steels with the present values found in this work [9,13,18,[22][23][24]26,32,[38][39][40][41][42][43][44]. At first glance, there are some differences in terms of activation energies. This fact was attributed to many influential factors: (the selection of boronizing method, the selection of the temperature range, the approach used for the computation of boron activation energies, the chemical composition of the substrates and the chemical or electrochemical reactions taking place during the boriding process).
where R denotes the ideal gas constant (=8.314 J mol −1 K −1 ) and T the boriding temperature given in Kelvin. Table 5 gives a comparison between the reported values of boron activation energies for some steels with the present values found in this work [9,13,18,[22][23][24]26,32,[38][39][40][41][42][43][44]. At first glance, there are some differences in terms of activation energies. This fact was attributed to many influential factors: (the selection of boronizing method, the selection of the temperature range, the approach used for the computation of boron activation energies, the chemical composition of the substrates and the chemical or electrochemical reactions taking place during the boriding process). Campos et al. [9] have employed the paste-boriding process to obtain the bilayer (FeB/Fe 2 B) at the surface of AISI M2 steel with two values of boron paste thickness (3 and 4 mm). They reported that the change in boron paste thickness had a significant impact on the growth kinetics of boronized layers. They have used a modified version of Brakman model [31] to assess the boron diffusion coefficients in FeB and Fe 2 B. As a result, the activation energies for boron diffusion were 283 and 239.4 kJ mol −1 , respectively. Gunes et al. [13] have boronized the AISI 8620 steel by using the plasma paste boriding between 973 and 1073 K with two paste mixtures (B 2 O 3 + SiC) and (B 4 C + SiC) and 100% B 2 O 3. This experimental study revealed that the values of activation energies depended on the ratios of either B 2 O 3 or B 4 C in the paste mixtures. Ortiz-Dominguez et al. [18] have utilized the commercial boriding agent (Durborid) to form the Fe 2 B layers on AISI D2 steel with flat morphology between 1123 and 1273 K. A diffusion model has been proposed for evaluating the boron activation energy in Fe 2 B (=201.5 kJ mol −1 ). Ruiz-Trabolsi et al. [22] have pack-boronized the AISI 1018 steel with the commercial powder mixture (Hef-Durferrit) by varying the size of samples (in the form of cylinders with 7 mm height and variable diameter of 1.85 to 12.70 mm). They proposed an equation for estimating the boronized layer thickness that depended on the time duration, the process temperature and the diameter of cylindrical samples. Accordingly, the values of boron activation energies for AISI 1018 steel were in the range 93.8-155.22 kJ mol −1 . Campos et al. [23] have used the commercial powder mixture (Durborid) to get the bilayer (FeB + Fe 2 B) on AISI M2 steel from 1173 to 1323 K. In their work, the mass balance equations were formulated at the two growing interfaces to estimate the boron diffusion coefficients in FeB and Fe 2 B. In addition, empirical equations were also proposed to validate their simulation work. Keddam and Kulka [24] have applied two kinetics approaches (mean diffusion coefficient (MDC) method and Dybkov model) for modelling the growth of FeB and  [32] have proposed a mathematical model that incorporates the effect of incubation periods into the boriding kinetics of AISI M2 steel. In their model, a non-linear distribution of boron element inside each boride layer (FeB or Fe 2 B) was considered. From Table 4, it is noticed that the plasma paste boriding [38] requires a minimum value of activation energy to produce the boronized layers on the surface of steels in comparison with solid boriding method [23,30,43,44].

Experimental Verification of Both Kinetics Approaches
The integral method and Dybkov model have been experimentally validated by utilizing the layers' thicknesses obtained at 1173, 1223, 1273 and 1323 K during 10 h. In the diffusion model based on the integral method, Equations (26) and (27) were used for predicting the layers' thicknesses of FeB and Fe 2 B (u(t) and l(t)) for a maximum boron concentration of 16.40 wt.% in FeB. For the two non-dimensional parameters ε 1 and ε 2 , the mean values were considered as follows: ε 1 = 0.064 and ε 2 = 0.093. For the Dybkov model, a numerical resolution of the obtained ODE system was performed by means of Runge-Kutta method [37] and by taking into account the following initial conditions u 0 = l 0 = 0.10 µm for the boride incubation times listed in Table 2. Tables 6 and 7 compare between the experimental values of FeB and Fe 2 B layers' thicknesses and the simulated values at different temperatures for 10 h. It is seen that the predicted layers' thicknesses coincide quite with the experimental values [23] by using these two kinetics models.  When the maximum value of boron concentration within the substrate is attained, the phenomenon of precipitation of iron borides (FeB and Fe 2 B) occurs at equilibrium condition as per the Fe-B phase diagram [4]. The first phase being formed is Fe 2 B, followed by the FeB phase. In this regard, the mass gain values with respect to each iron boride (FeB or Fe 2 B) [30] can be estimated from Equations (44) and (45), for an exposure time exceeding the corresponding boride incubation period: with ρ Fe 2 B = 7.33 g/cm 3 and ρ Fe = 7.86 g/cm 3 . Figure 4 shows the change in mass gain in FeB and Fe 2 B with the exposure time at increasing temperatures. It is seen that the calculated mass gain followed a parabolic trend with a prolonged treatment time because of the increase in diffusion rate of boron atoms through the double-phase boride layer (FeB + Fe 2 B). When the saturation level is reached in terms of boron solubility within the matrix, the nuclei of crystals start to grow and attain the critical size to become boride needles which come into contact with their neighbours till the coverage of the whole surface of sample. Subsequently, the compact boride layer forms after a time duration exceeding the boride incubation period. So the growth process is directly associated with the recorded mass gain. This mass gain can be regarded as a tool permitting the dimensional control of boronized samples before and after boriding treatment.

Conclusions
In this work, two kinetics approaches have been proposed to analyze the growth of dual-phase boride layer on AISI M2 steel between 1173 and 1323 K for variable time duration. The concluding points can be drawn as follows for this simulation work: (1) A peculiar solution of the obtained DAE system was suggested in order to assess the boron diffusivities in FeB and Fe2B for a maximum boron content in FeB of 16.40 wt.%.

Conclusions
In this work, two kinetics approaches have been proposed to analyze the growth of dual-phase boride layer on AISI M2 steel between 1173 and 1323 K for variable time duration. The concluding points can be drawn as follows for this simulation work: (1) A peculiar solution of the obtained DAE system was suggested in order to assess the boron diffusivities in FeB and Fe 2 B for a maximum boron content in FeB of 16.40 wt.%. (2) The two fitting parameters k FeB and k Fe 2 B for the Dybkov model were estimated as a function of process temperature by utilizing the experimental data of reference [23].