Multiscale Modelling of Carburization of Steel Through Fe-110 Surface and Role of Vacancy Defects
Authors: Chakrabarty Aurab , Sahputra Iwan Halim , Restrepo Oscar Antonio , Bouhali Othmane , Mousseau Normand , Becquart Charlotte S. , El Mellouhi Fadwa 
Published: July 2019

Carburization of metal is a catalytic reaction that occurs on metal surfaces exposed to hydrocarbon atmosphere at high temperatures. This reaction is a form of the well-known Fischer-Tropsch synthesis and is immensely important to various industrial aspects, most notably metal dusting corrosion (MDC) [1] and catalytic conversions [2]. On Fe surfaces, carburization occurs at high temperature, initiated by adsorption of gaseous hydrocarbons and is responsible for triggering both MDC, a catastrophic failure of the structural integrity of metals and alloys and a severe threat to the petrochemical industries, and a central reaction in catalytic converters, that are purposefully designed to produce a high yield of carburization to reduce the emission of toxic gasses. One of the most abundant and widely-studied carburizing gas is carbon monoxide (CO) that has been observed to react with Fe surface and dissociate at temperatures 600–900 K. In a hydrocarbon environment, the reaction is given by,CO+H_(2) ( → ⊥ 800K)[C^* + H]_2|OWhere C* is adsorbed on the Fe surface. The reaction is sustainable without the presence of hydrogen [3]. Gaseous hydrogen reacts with O atoms deposited on the surface from the dissociation of CO and removes the O atoms as steam [4].From atomic to continuum scale, the reaction characteristics have been studied widely using computer simulations of various approaches and have produced reliable results and predictions in agreement with relevant experiment results [5–9].From an atomistic modelling perspective, the reaction between CO and Fe surface can be broken into two consecutive processes, adsorption and dissociation. Adsorption describes the process of a CO molecule getting attached to the Fe surface. During dissociation, this molecule decomposes into C and O adatoms directly at the surface, with the two atoms moving into their most stable adatom site. Adsorption of CO is a barrier-less exothermic process. Dissociation, on the other hand, is an endothermic process that requires breaking a C-O triple bond, displaying an 11.2 eV/molecule average dissociation energy in vacuum [10]. The catalytic effect of the metal surface plays therefore a crucial role for breaking this bond at a much lower energy [4]. Earlier density functional studies show that on Fe-110 surface, the CO dissociation barrier is close to 1.5 eV [8]. However, this energy is still high in order to achieve from thermal activation only. A number of experiments suggest that the temperature threshold for CO dissociation can be as low as 380 K [11]. Therefore, it can be assumed that in practice, the reaction pathway is more complex than it has been modelled in previous works [7, 8] and one has to take into account that the metal surface contains impurities and defects that could assist the dissociation process [12]. Furthermore, the interactions originating from the periodic boundary condition of the simulations must be evaluated and dealt with accordingly. Commonly, an unwanted interaction from periodic boundary condition can be minimized by choosing a large simulation box. However, large simulations are computationally expensive and therefore an efficient method should be developed that can describe the CO-Fe surface reaction accurately yet is sustainable regarding computational resources.In this work, these challenges are addressed with a systematic study of the adsorption and dissociation of CO molecule on Fe surface, following further diffusion of C atom to subsurface and bulk using a multiscale technique combining atomistic density functional theory and empirical potential (EP) method. For atomistic scale study, we used density functional theory (DFT) with PBE-GGA pseudopotential technique implemented on the VASP code. Preliminary investigations on different low-index Fe surfaces confirm that the 110 surface is the most densely-packed surface and has the minimum structural reconstruction and minimum surface energy and is chosen as standard for our work. Hence we successfully reproduced CO adsorption energies, energy profiles for dissociation and subsurface diffusion of the C, in good agreement with results computed in earlier works [7, 8, 13, 14]. However, we have noticed that the periodic boundary condition applied to these simulation, which is also used to control the concentration of CO on the surface (coverage), has a significant effect on the energetics of these processes. The CO molecule has a strong dipole moment and it leads to van der Waals interaction between the molecules [15] adsorbed on the surface. Taking van der Waals interaction into consideration, adsorption energy studies on surfaces with CO fractional coverages of 0.25 and 0.0625 monolayers (ML) reveal that adsorption CO on 110 surface is energetically more favourable at dilute coverages. Most importantly we demonstrate here that the dissociation of CO also is energetically favourable in dilute surface as the energy barrier is reduced to half when the CO-coverage is reduced from 0.25 to 0.0625 ML [16].Since one expects surface defects to be present on Fe surfaces and to act as corrosion initiator, here we show computationally using accurate electronic structure calculations that the single vacancy defects on the surface are energetically inexpensive and therefore prone to be naturally abundant. In fact, we find that vacancy defects lower the adsorption energy for CO molecules adsorbed next to a vacancy and allows the C atoms to diffuse to the subsurface layers through the cavity. As a result, the reaction path becomes more complex and in order to examine the role of defect in the carburization process, one must compare the pathway with the combined pathway of dissociation of CO on clean surface, and subsequent surface to subsurface diffusion of C. It is worth noting that in the previously published works on surface to subsurface diffusion of C, the O adatom, as a by-product of the dissociation reaction, is ignored and its effects not taken into account. We take the O adatom adsorbed on the surface into consideration in this work and we carried out electronic structure and charge density analysis. It is demonstrated here that the influence of the O adatom is quite significant: 1. The O atom does not directly bond with the Fe surface but due to its high electron affinity, electronic charge is transferred from the surface to the O atom via the Fe-C bond in order to facilitate the breaking of C-O bond. 2. A surface vacancy defect creates an electron deficit that restricts the C atom from forming strong covalent bonds with Fe atoms. Interestingly, even after the C-O bond is completely dissociated, the O atom influences the surface-to-subsurface diffusion of C. Calculation of surface-to subsurface diffusion paths with and without O adatom on the surface shows that in presence of O, the diffusion barrier is higher and the C penetration depth into Fe is lowered.The dissociation and diffusion pathways discussed above were estimated using the nudged elastic band (NEB) method [18]. In order to find a transition state, the end-point (initial and final) configurations for the reaction are first fully relaxed with high precision. NEB uses an interpolated chain of intermediate configurations between the two end point configurations, connected by springs. The whole chain is then relaxed simultaneously with a fixed spring constant until the total average force minimizes under the tolerance limit of 0.01eV/Å. In this work, 8–16 intermediate configurations or images are considered, depending on the length of the path. It is evident that these calculations with DFT are extremely demanding regarding computational resource and therefore we need to anticipate the best possible dissociation/diffusion path with to minimize computation time. This issue is addressed by employing a combined empirical potentials (EP) – DFT technique where we employ the LAMMPS code, to pre-estimate possible minimum energy paths for dissociation and diffusion using NEB and from empirical Molecular dynamics simulations that are able to easily treat thousands of atoms thus increasing the length scales.Starting with the converged set of configurations from this pre-estimated transition pathways using EP, DFT is used to obtain more precise results. With the use of a potential that has been extensively tested to yield results in agreement with DFT [17], the pre-estimate shortens the computational time significantly. For the calculation of the surface-to-subsurface diffusion path of C, this process has reduced the computation real-time from 33792 to 18464 CPU-hours, which is a 45% drop, still producing identical results within the tolerance of 0.01 eV. However, in order to predict a reliable pre-estimate the potential needs to be tested with DFT. The embedded atom method (EAM) potential used in this work by Becquart et al [18] is rigorously calibrated with several DFT-estimated parameters. Unfortunately, this potential represents only Fe-C systems and for surface calculation with CO, we must consider a potential representing a Fe-C-O system. For this purpose, we have also been characterizing a reactive force-field potential designed by van Duin et al [19], specifically for catalytic reaction on Fe. This potential is anticipated to be used in near future to provide a complete atomistic description of this reaction mechanism. Which efficient multiscale technique are you referring to?In summary, the computational understanding gained in this study on the role of vacancy defects on dissociation of CO molecule, followed by subsurface diffusion of C, is beneficial for predicting the nature of the CO-Fe reaction for a practical scenario such as presence of larger vacancy defects on the surface and whether they act as dusting corrosion initiators. Simulation at different length and time scale are underway towards modelling physical and chemical phenomena governing the carburization of steel. Acknowledgement: The advanced computing facility of Texas A&M University at Qatar is used for all calculations. This work is supported by the Qatar National Research Fund (QNRF) through the National Priorities Research Program (NPRP 6-863-2-355).References[1] S. Lee, M. J. Kim, J. Lee, S-H. Kang, S. J. Lee and Y. Yun, Mater. Testing 57 (5), 409 (2015).[2] J. R. Dunetz and D. Fandrick an H.-J. Federsel, Org. Process Res. Dev. 19, 1325 (2015).[3] C. M Chun, J. D. Mumford and T. A. Ramanarayanan, J. Electrochem. Soc. 149, B348-B355 (2002).[4] D. J. Young, J. Zhang, C. Geers and M. Schütze, Mater. And Corr. 62, 7 (2011).[5] Q. Zhang, B. Han, K. Heier, J. Li, J. Hoffman, M. Lin, A-D. Kovacs and H. Cheng, J. Phys. Chem C 115, 12068 (2011).[6] Q. Zhang, B. Han, K. Heier, J. Li, J. Hoffman, M. Lin, S. Britton, A-D. Kovacs and H. Cheng, J. Phys. Chem C 116, 16522 (2012).[7] S. Booyens, M. Bowker and D. J. Willock, Surf. Sci. 625, 69 (2014).[8] D. Jiang and E. A. Carter, Surf. Sci. 570, 167 (2004).[9] D. Johnsom and E. Carter, J. Phys. Chem 114, 4436 (2010).[10] R. Sanderson, Polar Covalence, Acad. Press Inc. 1983.[11] G. Wedler and H. Ruhmann, Surf. Sci. 121, 464 (1982).[12] M. Elahifard, M. Jigato and J. Niemantsverdreit, Chem. Phys. Chem. 13, 89 (2012).[13] D. Jiang and E. A. Carter, Phys. Rev. B. 71, 045402 (2005).[14] A. Stibor, G. Kresse, A. Eichler and J. Hafner, Surf. Sci. 99, 507-510 (2002).[15] J. J. Victor, B. K. Kimmich, J. C. van der Waal, J. H. Zink, V. Zuzaniuk, J. T. Chapman and L. Chen, Patent US 8461395 B2: Synthesis of acetaldehyde from a carbon source (2013).[16] A. Chakrabarty, O. Bouhali, N. Mousseau, C. S. Becquart and F. El-Mellouhi, submitted to J. Comp. Sci. spl. conf. edn. (2015).[17] G. Henkelman and H. Jönsson, J. Chem. Phys. 113, 9901 (2000).[18] C. S. Becquart, J. M. Raulot, G. Bencteux, C. Domain, M. Perez, S. Garruchet and H. Nguyen, Comp. Mat. Sci. 40, 119 (2007).[19] C. Zou and A. C. T. van Duin, JOM, 64, 1426 (2012).