Modeling of compound layer growth during nitriding of pure iron

This paper presents a physical model which describes the layer growth kinetics and the nitrogen concentration profiles during gaseous or plasma nitriding of pure iron. The model is related to a one dimensional moving boundary value problem where the initial concentration profiles are assumed to be linear and a new boundary condition at the diffusion zone is proposed. The model is solved by using a classical finite difference scheme (FDM-CS) and the Heat Balance Integral Method (HBIM). Due to the proposed boundary condition at the diffusion zone, the numerical solutions can be validated through a simple physical analysis of mass transport theory in the asymptotic time limit. The results obtained are compared, and observed to be in good agreement with available experimental data and other approximate solutions reported in the literature.


INTRODUCTION
A basic framework for modeling the physical processes during nitriding is provided by moving boundary diffusion problems.Management and control of these processes are needed in order to obtain the required and acceptable surface properties such as hardness for tribological purposes of the components (Akhtar et al., 2011;Jacobsen et al., 2015;Park, et al., 2013).From a practical point of view, in order to get a good control of the process during the growth of concomitant layers, it is desirable to have a good estimate of the compact nitrided layers evolution.Several authors have modeled the concomitant growth of compact layers in the Fe-N thermodynamic system (Pye, 2003;Somers andMittemeijer, 1987, 1995;Du, 1993;Du and Agren, 1996;Christiansen and Somers, 2006;Mittemeijer and Somers, 1997;Kooi et al., 1996a,b;Hosseini et al., 2010;Fast and Verrijp, 1955;Ricard et al., 1990;Somers and Christiansen, 2005;Jaoul et al., 2006;Christiansen et al., 2008;Hassani-Gangaraj and Guagliano, 2013;Mufu et al., 2000), employing diffusion mechanisms.The study and solution of such models have been performed mainly by using error function approaches, HBIM, finite difference, and finite element methods.Each of these supply acceptable results for large values of time under *Corresponding author.E-mail: j.a.otero@itesm.mx.
Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License rigorous restrictions, by assuming a parabolic growth for the compound layer (Mufu et al., 2000;Mufu et al., 2003;Torchane et al., 1996;Triwiyanto et al., 2014) and a diffusion process that is already taking place within the compound layer and diffusion zone.In this work, there is not an a priori assumption of the parabolic growth for the compound layer and an incubation period, within which the compound layer and diffusion zone are formed, is taken from the experimental data of Somers and Mittemeijer (1995) and Marciniak (1985).Also, by proposing a diffusion zone of constant width, this work introduces a variant to the physical model first presented in León Cazares et al. (2014).The kinematics of compound layer growth is studied and solved using an HBIM and a finite difference method of second order.The proposal of a diffusion zone of constant width not reported before, brings some interesting insights into the physical nature of the process that are verified with numerical solutions.With the proposed model, it is possible to show that the concentration profile in every layer is exactly linear for large time values and it is also observed to be in good agreement with other models and available experimental data.
From a practical perspective, a mathematical model of the concomitant layers growth is presented, which takes into account the nitrogen concentration on the surface generated as a result of reactions between gas and solid.This is the case with gas nitriding, where mixtures of ammonia or nitrogen are used (Somers and Mittemeijer, 1987).These models take into account the nitrogen potential (or nitrogen concentration on the surface) and the evolution of concentration gradients from the surface generated by the concomitant growth of nitrides.On the other hand, during plasma-assisted processes, some models employ the surface evolution due to sputtering generated by ionized species (Marciniak, 1985;León Cazares et al., 2014).This work is developed on the assumption that growth of concomitant nitrided layers is governed by Fick's laws, which provide a phenomenological picture of nitrogen migration in different iron faces.

METHODOLOGY Mathematical model
Considering a one dimensional specimen where the nitrogen gas has penetrated into the iron so that at time , the , and phases are already formed.The surface concentration is assumed constant and equal to .Also, the effective diffusion coefficients only change with temperature, and they are assumed to be independent of the concentration at each point within the specimen.At each interface, the concentration is determined by the phase diagram of the Fe-N thermodynamic system.Therefore, the concentration profiles in each medium have the following boundary conditions , (1) at the phase, where , has been defined as the sum of the sputtering region , and layer thickness.

,
(2) at the phase, with , where is the layer thickness, and for the diffusion zone, with , where is the diffusion zone thickness.For simplicity the time dependence of each layer has been omitted; however, it should be kept in mind that , , and are functions of time.In general, at each phase, the initial concentration profile of nitrogen is shaped during the incubation period.The shape of these profiles is generally unknown and given as a function of the penetration . (4) Here the function is approximated in order to satisfy the boundary conditions given by Equations ( 1)-(3).
The other equations that model this problem are the diffusion equations in each phase, , where , and are the effective nitrogen diffusion coefficients in pure iron at the nitriding temperature .The velocity of each interface, between the two phases, is given by the Stefan Condition (SC) (8) The above equations are written in a general way, so the motion of the and the interfaces may be described by this condition.The constants that appear in Equation ( 8) have been defined, so that for , , , and .For , , , and .The main difference between the proposed model and the approximate analytical models proposed by other authors (Somers and Mittemeijer, 1995;Hosseini et al., 2010;Mufu et al., 2000Mufu et al., , 2003;;Torchane et al., 1996) is that, this problem is not solved by assuming a semi-infinite bar and it is not assumed beforehand that every layer has a parabolic behavior.Instead, the problem is solved over a finite domain in space, where a diffusion zone of constant width is assumed.The model will be tested against experimental data, and starting from a point in time where each layer is already formed, it is assumed that . ( The problem is completely defined by Equations ( 1)-( 9) and its solution will be approximated using the FDM-CS and the HBIM.

Numerical solutions
Here the finite difference method used to find approximate solutions for the motion of each interface and concentration profiles will be described.The method is generalized in order to consider sputtering and bilayer growth.However, the general description that will be given below can easily be applied to monolayer growth with sputtering, bilayer growth with no sputtering, and monolayer growth without sputtering.In the finite difference method, a continuous region is transformed into a finite number of points (nodes) and an approximate solution is found only at these points which constitute a mesh.For this reason, the differential operators are discretized at the mesh points.

Moving mesh classical finite difference scheme: Second order approximation
The partial time derivative of the concentration is expressed as a first order approximation of the backward difference in time, given by where , is the concentration at each phase.The second partial derivative at node of the concentration, relative to the position variable , is approximated by a central difference of second order . (11) Here and represent the length of the steps in for and , respectively.The discretization of the argument is represented by , and the argument is represented by .Therefore, in this notation, . By substitution of ( 10) and (11) in Equations ( 5)-( 7), the discretized diffusion equation is given by ( 12 .This finite difference method described is a classical moving mesh scheme or (FDM-CS).

Heat balance integral method
The HBIM used in León Cazares et al. (2014) will be described, in order to compare the solutions with the FDM-CS and the analytical method used in Mufu et al. (2000Mufu et al. ( , 2003)).In this approximation, the unknown concentration profiles are assumed to be of the form for the diffusion zone.
In order to satisfy the boundary conditions, the profiles given by Equations ( 16)-( 18) must satisfy the following , Substitution of Equations ( 16)-( 18) in the SC, Equation (8) gives for the motion of the interface, and for the motion of the interface.To obtain an approximate solution, the diffusion Equations ( 5)-( 7) along with the SC at each interface, given by Equations ( 22) and ( 23), must be solved.The central idea behind the HBIM is to assume an adequate profile for the solution which enforces the available boundary conditions.Next, the diffusion equations are transformed into a set of ordinary differential equations in time (ODE's), by averaging the partial differential Equations ( 5)-( 7) over the position variable.The mathematical process is fully described in [26], after which the following set of ODE's is obtained (24) for the layer, where the constant has been defined.Similarly for the layer the averaged diffusion equation leads to  19)-( 21) and Equations ( 22)-( 26) along with Equation (9) a system of six differential equations was obtained in time that must be solved for the and layers.To find a numerical solution to the problem, the FDM-CS and the HBIM start from a point in time after the incubation period with the following set of initial conditions: the thickness of the and layers are taken from the experimental data as and at some time .The concentration profiles at this time are assumed to be linear; therefore, , for .In the results analysis, a comparison between the FDM-CS and the HBIM with experiments reported in Somers and Mittemeijer (1995), Mufu et al. (2000), Marciniak (1985) will be shown.In order to compare with the experimental data of Somers and Mittemeijer (1995) and Marciniak (1985), it is being assumed a diffusion zone of constant thickness.It will be shown that, for the constant diffusion zone model, the solution obtained for the compound layer, converges as the width of the diffusion zone is increased.The FDM-CS and the HBIM presented here are solved numerically for increasing values of until no significant changes on the motion of each interface are observed within the time window considered.

Approximate analytical solution: Semi-infinite one dimensional bar
Here, the approximate analytical method (AAM) used in Mufu et al. (2000Mufu et al. ( , 2003) ) to solve the one-dimensional diffusion problem of nitrogen in pure iron, will be reproduced.Given that this model is used to compare the solutions obtained with the proposed methods, and the experimental results cannot be reproduced with the original equations found in (Mufu et al., 2000(Mufu et al., , 2003)), the analytical model will be described briefly.
Since the original experiments in (Mufu et al., 2000) do not consider sputtering, the AAM assumes beforehand a parabolic behavior of each interface and .With this assumption , ( 27) and where the position of each interface without sputtering was defined earlier as: , and .
The concentration profile at each phase is described by the error functions as where correspond to the epsilon, gamma and alpha phases, respectively.
In the original work (Mufu Yan et al., 2000), the diffusion zone is defined within a region where the minimum concentration is taken as .After substitution of the boundary conditions into Equation (30), the concentration profile at each phase is obtained as and where , and were defined subsequently in the study.By substituting Equations ( 31) and (32) in the SC for the interface, the following is obtained . (34) Similarly, substitution of Equations ( 32) and (33) in the SC for the interface leads to Given that this model is a solution for a semi-infinite region, the concentration at the diffusion zone for large values of , goes asymptotically to zero, therefore , which can be written as Equations ( 27)-( 29) along with Equations ( 34), ( 35) and ( 37) can be solved to obtain the constants for and the concentration profiles at each phase.

RESULTS AND DISCUSSION
Here, the results obtained with the FDM-CS and the HBIM described earlier, their comparison with the AAM subsequently described, and the experimental data of Somers and Mittemeijer (1995), Mufu et al. (2000) and Marciniak (1985) will be shown.Firstly, the convergence of the solution as the size of the diffusion zone is increased will be illustrated.To validate the numerical simulations, the FDM-CS and HBIM solutions will be compared with experimental results for monolayer growth.Next, the numerical simulations are compared with experimental data for bilayer growth due to gaseous nitriding (Somers and Mittemeijer, 1995); also with experimental data due to pulse ion nitriding (Mufu et al., 2000) and the AAM will be revised when bilayer growth in pure iron is discussed.

Numerical experiments for bilayer growth
First, the FDM-CS and the HBIM solutions described earlier are compared.The initial profile is assumed to be linear in the position coordinate, so that equation ( 4 for the diffusion zone.In the above equations the time where the numerical simulations start, is shown explicitly.This concentration profile is used as the starting point for the FDM-CS and the HBIM simulations.The initial values for the functions are obtained from the boundary conditions at each interface.Starting the numerical simulations at , when the , and layers are already formed. For the numerical experiments presented in this section the concentrations at each boundary and diffusion constants at each layer are obtained from Mufu (2000), which are consistent with the phase diagram of the Fe-N thermodynamic system and the effective diffusion constants reported in Somers and Mittemeijer (1995).The concentrations at each boundary are obtained from the following equations ( 41) and ( 42) Hernández et al. 139 where must be in Kelvin.For these experiments, the nitriding temperature considered is and no sputtering is taking place.The effective diffusion constants, concentrations at each interface and the size of each layer at are shown in Tables 1 and 2. The and widths are obtained from the AAM at .These widths are the starting point for the simulations.Figure 1 illustrates the resulting time evolution for the and layers as the width of the diffusion zone is increased from to .The maximum time was taken as ; therefore, the time window used in these simulations is . The obtained solutions from the HBIM and FDM-CS are in almost perfect agreement.Figure 2 also illustrates how these solutions do not show significant changes within the time window considered, for diffusion zone widths of .It is expected that the experimental data can be explained or approximated by using the FDM-CS or the HBIM methods described above for a diffusion zone thickness where the solutions do not change within the time window of the experimental setup.
One of the interesting aspects of the proposed model is that it is possible to find exact analytical expressions for each layer width and concentration profile in the asymptotic time limit.This is possible to achieve due to the boundary conditions at the diffusion zone.Since the concentration at the end of the diffusion zone is assumed to be zero, from basic transport theory the width of each layer for large time values can be obtained as for the layer and for the layer.Here, the and represent the thickness of each layer in this limit.Since the width of the diffusion zone is constant, the flux through this region is constant as well.The layer will grow until the flux through this layer equals the constant flux through the diffusion zone.At this point the net flux through each interface will be zero, which has an interesting physical implication for large time intervals in the concentration profiles as well.
Given that the net flux through any point of the sample is zero, the concentration profile in this limit can be obtained as well.If a particular spot within the sample is chosen, when the incoming flux is equal to the outgoing flux through this point, then for an arbitrary point at some position within the layer  41) and ( 42). , (45) and solving for the concentration profile at the layer is obtained as .( 46) In a similar way, it is straightforward to obtain the concentration profile at the layer where the asymptotic values of the and layers are given by Equations ( 43) and ( 44).A remarkable conclusion from Equations ( 46) and ( 47) may be drawn.In this time limit, the concentration profile at each layer is exactly linear.The result does not depend on the initial profile; this predicted behavior, given by Equations ( 43) for the and layers is shown in Table 1.The asymptotic values predicted by Equations ( 43) and ( 44) are also shown.1.The asymptotic limits predicted by Equations ( 46) and ( 47) are also shown in each case.and (44), is illustrated in Figure 2. Several simulations are shown for different diffusion zone widths, by using the FDM-CS and the HBIM, with the physical parameters given in Tables 1 and 2. The corresponding asymptotic values are also shown in Figure 2, where both methods are observed to be in perfect agreement with the predicted values.
To further validate the numerical methods, in Figure 3, the concentration history at two different positions within the sample, is shown.A point at is chosen so that initially it is immersed in the layer.The concentration at this position will jump, when the layer reaches this position.Similarly, another point at is chosen so initially it is immersed in the diffusion zone.At some time, when the layer reaches this position, the concentration will jump to values that correspond to the phase.
From Figure 3, it can be observed that both methods approach the predicted values given by Equations ( 46) and (47) for large time values.Even if the initial concentration profile is parabolic, the dynamics of the nitrogen transport gradually shapes the profile, until the behavior is exactly linear in this time limit.Even for small time intervals, it is also found that this diffusion mechanism rapidly brings the concentration profiles close to a linear behavior.41) and (42).

Figure 4.
Comparison between the FDM-CS and the HBIM solutions with the experiment reported in Marciniak (1985) for the motion of the layer at different sputtering rates.is the sputtering rate in , with the least error between numerical simulations and the experimental data for a sputtering rate of .

Monolayer growth in pure Iron
The FDM-CS and HBIM are compared with the experiments of Marciniak (1985), where the iron surface is sputtered due to plasma nitriding at .In the experiment, the growth of the layer at a constant sputtering rate is observed.The growth of the sputtering zone is modeled as (48) where is the sputtering rate measured in .The physical parameters used for the experiment can be obtained from Equations ( 41) and ( 42) for the concentrations at each boundary.The effective diffusion constants for the layer and diffusion zone at are and .
Using the model described here, the numerical solution used to compare with the experimental data, is the one that does not show changes within the time window considered in the experiments of Marciniak (1985).The concentrations at each interface are shown in Table 3.
In Figure 4, the FDM-CS and HBIM results for different sputtering rates are shown.As expected, the monolayer growth is not parabolic due to the sputtering of the iron surface.The FDM-CS and HBIM show the best agreement with the experimental data of Marciniak (1985) for a constant rate of .The simulations start from , where the width of the interface is .Experimental data is available up to a time interval of ; therefore, the time window is .(Somers and Mittemeijer, 1995), for .
Figure 5. FDM-CS and HBIM solutions from the constant diffusion zone model proposed in this work and experimental data (Somers and Mittemeijer, 1995) for , and compound layer thickness due to gaseous nitriding of pure iron at .The thickness of each layer, where the simulations start, is obtained from the first data point at .For , the effective diffusion constants at each layer are taken from Mittemeijer and Somers (1997) as , and .

Bilayer growth in pure Iron
The solutions from the FDM-CS and the HBIM are compared with the experimental data of Somers and Mittemeijer (1995) for gaseous nitriding in pure iron at . The concentrations at each interface and effective diffusion coefficients from [12] are shown in Table 4.The constant diffusion zone model is compared with experimental data for bilayer growth.The simulations start at , where the widths of the and layer are and , respectively.The obtained solutions are shifted by and calculated for the next .The effective diffusion coefficients at in each layer are , and for the , and diffusion zone, respectively.Figure 5 shows a comparison between the constant diffusion zone model and experimental data of Somers and Mittemeijer (1995) for the motion of the , and compound layer.The predicted values from the proposed models deviates at most by from the experimental data shown in Figure 5 for the compound layer.
In Figure 6, the concentration profile at is calculated in order to compare it with the experimental data in Somers and Mittemeijer (1995).This profile is the concentration in as a function of the nitrogen penetration along the iron sample.The physical parameters correspond to those used in Figure 5 and Table 4.
At , the FDM-CS and the HBIM use a linear concentration profile as mentioned earlier.The solutions  (Somers and Mittemeijer, 1995) for the concentration profile after 2 h of nitriding.The physical parameters for the FDM-CS and HBIM are obtained from Figure 6 and Table 4.
have been tested as well, by using initial quadratic profiles; however, the result is practically the same as the one obtained by assuming a linear profile.The difference between any solution; obtained with different initial profiles, is not distinguishable within the scale shown in Figure 6.It has been shown that the proposed model is independent of the initial profile in the asymptotic limit; therefore, it is expected that even at smaller time intervals the concentration profile in the compound layer will be almost linear (Somers and Mittemeijer, 1995;Hassani-Gangaraj and Guagliano, 2013).Although this behavior is mentioned in Somers and Mittemeijer (1995), it is not explicitly shown or explained.
Finally, solutions from the AAM reviewed are compared with the experimental data of Mufu (2000) obtained from pulse plasma nitriding at during a time interval of .The model assumes a small jump of in concentration at the end of the diffusion zone.Therefore, in order to compare the solutions from the AAM and the experimental results of Mufu (2000), Equation ( 9) must be modified to consider the concentration jump at the and pure iron interface.Therefore, this interface satisfies a SC given by Equation (9).

. (49)
The concentrations at each interface and effective diffusion coefficients for each layer are given in Tables 1  and 2. In this experiment, the concentration at the end of the diffusion zone is taken as .By using an extra SC, given by Equation ( 49), the growth of the diffusion zone may be obtained.Due to the nature of this experiment, the constant diffusion zone model described earlier, is not used.In Figure 7, an excellent agreement between the FDM-CS, the HBIM proposed in this work and the AAM, is observed.Also, the three obtained solutions show good agreement with the experimental data.

Conclusion
The kinetics of layer growth in pure iron is obtained by a simple diffusion model that assumes a diffusion zone of constant thickness.The model is solved by the FDM-CS and the HBIM.First, the numerical solutions are validated by means of the asymptotic time behavior of the solution for the and layers and the asymptotic behavior of the concentration profile for the model of constant diffusion zone.With various numerical examples, the predicted behavior is shown to be captured by both methods.By means of basic mass transport theory, it is shown that the concentration profile is exactly linear in the asymptotic regime, and therefore, the solution is independent of the  Mufu (2000) with the experimental data for the motion of the compound layer due to pulse plasma nitriding of pure iron at .The thickness of each layer at where the FDM-CS and the HBIM numerical simulations start, is shown in Table 1.initial concentration profile in this limit.Even more, it is observed that the concentration profiles in the and layers are practically linear for small time intervals, which means that even if the initial concentration profile is parabolic, the diffusion mechanism shapes the profile very rapidly to a linear profile, thus approaching the predicted asymptotic limit.Second, the constant diffusion zone model is tested against experimental data, where the greatest percent difference between our model and the experiment in Marciniak (1985) is 4.54% for a sputtering rate of 0.65 µm/h and 6.32% between our model and the experimental data for the compound layer reported in Somers and Mittemeijer (1995).Finally, a diffusion zone model is proposed, where the kinematics of the layer is described by a Stefan condition at the fictitious boundary between the diffusion zone and pure iron.The maximum percent difference between our model and the experimental data of Mufu (2000) is 13.47%.From Figure 7, a good agreement is observed between the mathematical model of Mufu (2000) where parabolic growth is assumed, and the proposal from this work.In this case, both models predict the same values for compound layer growth, up to the fourth decimal place.
are the backward (forward) discretized first order derivatives of the concentration, evaluated at the interface position , andHernández et al. 137

Figure 1 .
Figure 1.FDM-CS and HBIM solutions obtained for the and layers at .Values of the and layers at for the FDM-CS and HBIM computer simulations were obtained from the AAM and are shown in table1.Each frame is a numerical simulation where the thickness of the diffusion zone is kept constant.

Figure 2 .
Figure 2. Time evolution of the and layers obtained using the FDM-CS and the HBIM for different diffusion zone widths.The thickness atfor the and layers is shown in Table1.The asymptotic values predicted by Equations (43) and (44) are also shown.

Figure 3 .
Figure 3. Concentration history for a) and b) obtained with the FDM-CS and the HBIM for a diffusion zone width of .The thickness of the and layers at , from which the FDM-CS and HBIM simulations start are shown in Table1.The asymptotic limits predicted by Equations (46) and (47) are also shown in each case.

Figure 6 .
Figure 6.FDM-CS and HBIM solutions from the constant diffusion zone model and experimental data(Somers and Mittemeijer, 1995) for the concentration profile after 2 h of nitriding.The physical parameters for the FDM-CS and HBIM are obtained from Figure6and Table4.

Figure 7 .
Figure 7.Comparison between the FDM-CS, HBIM from this work and revised AAM from the original work ofMufu (2000) with the experimental data for the motion of the compound layer due to pulse plasma nitriding of pure iron at .The thickness of each layer at where the FDM-CS and the HBIM numerical simulations start, is shown in Table1.

Table 1 .
and layer widths at and effective diffusion constants obtained from the AAM described in the study.

Table 2 .
Nitrogen concentrations at in each interface for , according to Equations (

Table 3 .
Nitrogen concentrations in at each interface for , according to Equations (

Table 4 .
Nitrogen concentrations in at each interface