45 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 Modeling of Ferrous Metal Diffusion in Liquid Lead Using Molecular Dynamics Simulation Ahmad Anwar Nuris 1,a 1 Berka Semi Strategika, KPP IPB Baranangsiang IV Blok A/35 Tanah Baru, Bogor Utara, Kodya Bogor, Bogor 16154, Indonesia a al.adrislamic@gmail.com Abstract. Modeling of Iron metal diffusion in liquid lead using molecular dynamics simulation has been done. Molecular dynamics simulations are used to predict the value of physical quantities that we want to know based on the designed material model and on the input simulation data. In this research, effect of different geometry of material models was observed to know the diffusion coefficient. The material system was iron (Fe) in liquid lead (Pb). The material models is designed using Packmol software to get the initial configuration of atom's arrangement by inputting the material's characteristics such as mass, density, volume, number of atoms. This work examines the diffusion coefficient of iron in molten lead metal with the geometric shape of the simulation system in the form of iron in molten metal for various simulation models of boxes in a box, balls in a box and balls in balls. To design simulated geometric shapes we use the Packmol program. To calculate the diffusion coefficient we use the molecular dynamics simulation method. To find out which geometry is suitable, we compare the diffusion coefficient of the simulation results with existing references. The diffusion coefficient value of the spherical iron (Fe) system in the spherical liquid lead (Pb) has the best value compared to the other two forms with an accuracy rate of 99.94% because it is influenced by the even distribution of atoms in each part. Keywords: LAMMPS, Molecular Dynamics, Fe-Pb, Packmol, Mean Square Displacement Introduction Computation is an activity to obtain a solution of complex problems following a certain mathematical model [1]. Computer simulation is a tool for studying macroscopic systems by applying microscopic models, especially for prediction of materials properties [2]. One example of research using computer simulations is the prediction of diffusion coefficient of materials that is important data for many applications. Knowing the diffusion coefficient can help us to study the corrosion phenomena as in the field of nuclear reactor design [3,4]. There have been many corrosion experiments in search of superior steel for nuclear applications and determining the appropriate method for corrosion inhibition [5]. The high cost of installation for corrosion experiments and the need for a high safety level is the main constraints today. This is because the metal vapor produced is very toxic, and also, not all experiments can be carried out in an operating reactor. Particularly in Indonesia, this activity seems not possible due to inadequate facilities [6]. Therefore, computation and simulation methods are solutions to overcome these obstacles. One of the computational studies conducted is to use the molecular dynamics method [7-9]. When examining the corrosion of steel materials in the fast reactor, ferrous metal is the largest steel composition element. At the same time, molten lead is a suitable material as a coolant for 46 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 sodium substitute rector. The most relevant property of lead is the massive difference between its melting and boiling points. The boiling and melting points of lead are 601K and 2022K respectively. These properties lead to higher reliability and safety for the reactor installation than the use of sodium [10]. The diffusion coefficient is obtained with a high degree of accuracy when compared with the experimental results. This study using Packmol software for preparing the configuration of coordinate x, y, and z of the atoms of materials. The initial configuration of the arrangement of iron (Fe) and lead (Pb) atoms from Packmol will be running in LAMMPS (Large Scale Atomic / Molecular Massively Parallel Simulator: lammps.sandia.gov), an open-source software that has the advantage of being able to run large-scale computing. With the number of atoms up to millions of grains. LAMMPS software is also widely used for material simulation. This can be seen from the number of journal articles published from the results of LAMMPS-based computational research. Also, this software is always updated [7]. In this study, calculations were carried out to determine the best Mean Square Displacement (MSD) value in the variation of the simulation system model for ferrous metal (Fe) in the liquid lead (Pb) obtained when calculated using the molecular dynamics software LAMMPS. And to find out the value of the diffusion coefficient (D), which is the best in the variation of the system model of iron (Fe) in the liquid lead (Pb) obtained when calculated using the molecular dynamics software LAMMPS. Theoretical Background The interaction model between molecules needed in the simulation is the law of intermolecular force, which is equivalent to the potential energy function between molecules. The selection of the potential energy function must be made before any simulation is carried out. The choice of the interaction model between molecules greatly determines simulation correctness from a physics point of view. Because they are on the atomic scale, interactions must, in principle, be derived quantum, which is where the Heisenberg uncertainty principle applies. However, a classical mechanical approach can be used where the atom or molecule is considered a point mass [11]. For N, the number of atoms in a simulation, the potential energy function is U(R N ), where R N is the set position of the center of mass of the atom or molecule, R N = {R1, R2, R3,…. RN} which can be expressed as: (1) Potential energy is the sum of the interactions between two isolated molecules [12]. There are potential energy models used in molecular dynamics simulations, including the Lennard-Jones potential. This model's use in the simulation provides a reasonably good level of accuracy in describing the interactions between atoms. The main characteristic of Lennard- Jones is that it is very responsive for small r and less attractive for large r [13]. These characteristics can be seen in the image below: 47 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 Figure 1. The atomic interaction model with the Lennard-Jones potential (a) repulsive (b) attractive [14] This potential model equation can be formulated: (2) n and m are positive integers selected n>m, and i, j are the molecule indices, Rij = |Ri-Rj| or the distance between molecules i and j. Whereas σ is the distance parameter, and ε is the parameter that states the interaction's strength. k is the coefficient obtained from the equation: (3) Common choices for m and n are m = 6 and n = 12 [15]. So that we get the equation: (4) The Lennard-Jones potential can be used to determine the characteristics of gases, liquids, clusters, and polycrystalline materials [16]. The Lennard-Jones potential for molecular dynamics simulations has parameters that can be taken from experimental data. Materials and Methods 1. Preparation of Simulation Inputs Based on the flow diagram in the simulation input preparation step, the run's initial simulation is a system consisting of iron (Fe) in the molten lead (Pb). The iron metal being modeled consists of 3527 atoms and 8000 metal atoms of lead (Pb). The temperature used is 1023 K or 750 o C because it is at this temperature that W. M. Robertson experimented calculating the diffusion coefficient of iron in the pure lead. The diffusion coefficient value from the experimental results is 2.8 x 10-9 m 2 /s [17]. The density of the system is calculated using the temperature-dependent density equation, as shown by equation 5 with T in Kelvin [18]. (5) 48 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 2. Determination of the Diffusion Coefficient (D) Determination of the diffusion coefficient (D) of ferrous metal (Fe) can be seen from the data output of simulation from the Mean Square Displacement (MSD) in the format (.txt), which can be processed using Microsoft Office Excel. The MSD value data for iron metal (Fe) is represented in graphical form against simulation time. The diffusion coefficient D value is obtained from the gradient graph of the relationship between the MSD value and the simulation time. After The diffusion coefficient D value of the simulation results is determined, it can be compared with the diffusion coefficient D value of the experimental results so that the simulation results' accuracy level can be seen. Results and Discussion The simulation using the molecular dynamics method basically begins by determining the atoms of the material being studied, namely the iron atoms and the lead atoms. The arrangement of the atomic configuration is made after calculating the volume of the system to be made. Calculation of the iron and liquid lead atoms' volume size is carried out using the services (utilities) available on the Packmol website. The data input provided includes the density of the molecular mass to be made, the number of atoms contained in the molecule, and the molar mass of these atoms. The simulation system made in this study is iron (Fe) in the form of a cube with BCC structure in the lead (Pb) in the form of a cube. This system consists of 3527 iron (Fe) atoms and 8000 lead (Pb) atoms. The visualization of this configuration can be seen in Figure 2 of the following page. Figure 2. (a) cube-shaped iron (Fe) with the structure of BCC in liquid lead (Pb) in the form of a cube visualized with OVITO software, (b) cubic iron (Fe) with the structure of BCC in liquid lead (Pb) in the form of a cube visualized with OVITO with add modification slice 49 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 Figure 3. (a) Spherical iron (Fe) with a BCC structure in the liquid lead (Pb) in the form of a cube visualized with OVITO software, (b) Spherical iron (Fe) with a BCC structure in the liquid lead (Pb) in a cuboid visualized with OVITO with the add modification slice. Besides that, a spherical iron (Fe) simulation system with the structure of BCC in the lead (Pb) is also constructed, consisting of 3527 iron atoms (Fe) and 8000 lead atoms (Pb). Visualization of a spherical iron (Fe) system with a BCC inside structure cuboid of lead (Pb) is shown in Figure 3. While the third simulation system created is a system of 3527 iron atoms (Fe) with a BCC structure in the form of a ball in a spherical liquid lead (Pb) consisting of 8000 atoms. Figure 4. (a) Spherical iron (Fe) with the structure of BCC in the liquid lead (Pb) in the form of a sphere visualized by OVITO software, (b) Spherical iron (Fe) with the structure of BCC in the liquid lead (Pb) in the form of a ball visualized with OVITO with add modification slice Figure 4 is a visualization of a spherical iron (Fe) system with BCC structure in the spherical lead (Pb). These three systems are designed to determine the simulation system's ideal form in the simulation process of the iron (Fe) system in the liquid lead (Pb). The determination of ideal is based on the most diffusion coefficient's value closer to the experimental results. Mean Square Displacement (MSD) System of Iron (Fe) in Liquid Lead (Pb) The diffusion process can affect the crystal structure of a material. Even diffusion can cause crystal defects. In the diffusion process, there is also the movement of atoms to change their position. Each atom that undergoes movement, then the atom has the average square of the atomic movement, commonly referred to as MSD (Mean Square Displacement). The MSD simulation results from LAMMPS against the simulation time are in the form of MSD value data in the format (.txt), which can be processed using Microsoft Office Excel. The MSD relationship curve to the simulation time of the iron (Fe) system in the liquid lead (Pb) is shown 50 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 in Figure 5. In the curve shown in Figure 4.8, it can be seen that the MSD value ranges of three different crystal forms. Figure 5. The curve of the mean square displacement (MSD) relationship with the simulation time of the iron (Fe) system in the liquid lead (Pb) The blue curve shows the simulation results of the iron (Fe) system in the form of a cuboid with BCC structure in the liquid lead (Pb) in the form of a cube, the red curve shows the simulation results of the spherical iron (Fe) system with the structure of BCC in the liquid lead (Pb) in the form of a cube. The green curve shows a spherical iron (Fe) system with BCC structure in the spherical liquid lead (Pb). In the curve shown in Figure 4, it can also be seen that the MSD value of the spherical iron (Fe) system in the liquid lead (Pb) has the highest value compared to the other two forms. The nonlinearity at the start is typical of the MSD curve. This is because the atom has not yet interacted with other atoms. After the nonlinear section at the start of the curve, it is followed by a straight line. The higher the slope of the MSD curve, the more random the atomic arrangement is. Meanwhile, the slope of the MSD curve is directly proportional to the diffusion constant of a material. Diffusion Coefficient (D) of the Iron (Fe) System in Liquid Lead (Pb) The diffusion coefficient can be said to be the most precise physical quantity known from a system. When a system gets disturbed by the influence of the environment, it can cause the movement of the atoms that make up the system randomly so that the coefficients in the system can be known. The movement of the atoms that make up the system can produce atomic trajectories. The atomic trajectories can represent each atom's positions so that each of these atomic positions can cause an interaction force between atoms. According to Newton's law, potential energy is a negative gradient of the interaction force, so from the interaction force between atoms, it can be seen the potential energy between atoms. As explained in the paragraph above, the determination of the diffusion coefficient value of a system made for simulation is also influenced by the potential models used. The potential models used are the Lennard-Jones potential. This potential is recommended for use in a simulation system of iron (Fe) in the liquid lead (Pb) from the previous studies. The Lennard - Jones potential has two parameters, namely the distance parameter and the energy parameter. Cubic Fe in cubic liquid Pb Spherical Fe in cubic liquid Pb Spherical Fe in spherical liquid Pb Simulation time 51 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 These two parameters significantly affect the simulation results used. The distance parameter and the energy parameter have a value that varies from one atom to another. This Lennard - Jones potential will also be an input in the simulation process using the LAMMPS software. The simulation results of the diffusion coefficient value can be shown in Table 1. The most significant diffusion coefficient value is obtained from the spherical iron (Fe) system with BCC structure, in the liquid lead (Pb) in the form of cubes, and the lowest is the spherical iron (Fe) system with BCC structure, in the liquid lead (Pb) in the form of a sphere. This shows that the diffusion coefficient is linear with the MSD calculation results. Table 1. Diffusion Coefficient in the Iron (Fe) System in Liquid Lead (Pb) Crystal Shape Coefficient (m2/s) Accuracy Cubic Fe in cubic liquid Pb 2.91155x 10 -9 96.02% Spherical Fe in cubic liquid Pb 3.00955 x 10 -9 92.52% Spherical Fe in spherical liquid Pb 2.79835 x 10 -9 99,94% Based on Table 1, it is also known that the accuracy value in the spherical iron (Fe) system with BCC structure in the liquid lead (Pb) in the form of a ball is higher than the other two systems. The simulation accuracy level can be determined by comparing the calculation results of the simulation with the experimental results. W. M. Robertson has carried out experiments to calculate the diffusion coefficient of iron in the pure lead. The experimental results' coefficient value is 2.8 x 10-9 m 2 /s [13]. Thus, the spherical shape is an ideal form of a simulation system for simulating liquid iron (Fe) in lead (Pb) systems. This is because the shape of the iron (Fe) simulation system model in the spherical liquid lead (Pb) provides better accuracy values than other forms. Based on table 1, it can also be seen that the forms of simulation system models made on the iron (Fe) system in the liquid lead (Pb) have an effect on the MSD value. This also impacts the diffusion coefficient value obtained because, in theory, the diffusion coefficient value is directly proportional to the MSD. This is inseparable from the distribution of the atoms in each part of the system. The lowest accuracy value is obtained from the spherical iron (Fe) system with a BCC structure in the liquid lead (Pb) in the form of cubes. This is because the lead (Pb) atoms in this system are not evenly distributed in certain parts. The number of lead (Pb) atoms is greater than the side direction in the corner direction. Based on table 2 it is known the number of atomic distributions in each part of the simulation system. Table 2. Comparison of Atomic Distribution in Iron (Fe) System in Liquid Lead (Pb) Crystal Form Number of Atomic Distribution In the corner On the side Cubic Fe in cubic liquid Pb 268 194 Spherical Fe in cubic liquid Pb 308 202 Spherical Fe in spherical liquid Pb 259 238 This can be seen from the visualization results on OVITO, which is shown in the image below. 52 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 Figure 6. Spherical iron (Fe) with a BCC structure in the liquid lead (Pb) in the form of a cubic visualized with OVITO with an add-on modification slice Figure 6 shows the distance between the corner points of the cube and the outer side of the sphere as labeled with (α) has a longer distance than the distance between the cube's side and the outer side of the sphere labeled with (β). Based on this distance factor and the number of atomic distribution causes the interaction force between the atoms to be uneven from each side of the system is created. The accuracy value of the cube-shaped iron (Fe) system with BCC structure in the liquid lead (Pb) is better than the spherical iron (Fe) system with BCC structure in the liquid lead (Pb) in cubic form. Figure 7. Cubic iron (Fe) with BCC structure in the liquid lead (Pb) in cubic shape visualized with OVITO with add modification slice Figure 7 shows that the atomic distribution is better than the atomic distribution. This is because the number of atoms in the cube-shaped iron (Fe) system with the BCC structure in the liquid lead (Pb) is more evenly distributed. If we compare the two, it can be seen that the α / β ratio for a spherical system in a cube is more significant than that for a cube system in a cube. The most significant accuracy value is the ball in the ball model. Because theoretically, the distribution of atoms and energy in the system is almost evenly distributed for each part. 53 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 Figure 8. Spherical iron (Fe) with BCC structure in the liquid lead (Pb) The ball was visualized with OVITO with an add-on slice modification. This is also reinforced by the OVITO visualization results, which shows the value of α = β, as shown in Figure 7. The comparison of α and β distances for the three systems of iron (Fe) in the liquid lead (Pb) can be seen in table 3. Table 3. Comparison of α and β in the iron (Fe) system in the liquid lead (Pb) Crystal Form α (Å) β (Å) Cubic Fe in cubic liquid Pb 30.14 26.10 Spherical Fe in cubic liquid Pb 43.29 12.81 Spherical Fe in spherical liquid Pb 18.69 18.69 Based on table 3, the α / β ratio for the spherical system in a sphere is smaller than the cube system in a cube and the spherical system in a cube. In simple terms, the α / β ratio for the three systems can be denoted in the equation below. (6) Where α is the angular distance of the liquid lead (Pb) from the outer side of the iron (Fe), and β is the angular distance of the liquid lead (Pb) to the outer side of the iron (Fe). Conclusions Based on the research that has been done, it can be concluded that, first, Making the ideal molecular dynamics simulation system for the iron (Fe) system in the liquid lead (Pb) begins by calculating the volume of the system made, namely iron (Fe) in the liquid lead (Pb) by ensuring that the mass density value is maintained. The mass density must also be calculated using the temperature-dependent density equation so that the simulation results can be close to the experimental results. These data become the basis for making the initial configuration. The initial configuration in the form of x, y, z coordinates is made using credible software, one of which is Packmol. Second, the mean square displacement (MSD) value of the spherical iron (Fe) simulation system in the spherical liquid lead (Pb) has the best value compared to the other two forms. Third, the spherical iron (Fe) system's diffusion coefficient value in the spherical liquid lead (Pb) has the best value compared to the other two forms with an accuracy rate of 99.94% because it is influenced by the even distribution of atoms in each part. 54 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 References [1] Warsito, 2005, Komputasi Tomografi dan Aplikasinya dalam Proses Industri, Ohio: Departement of Chemical Engineering The Ohio State University. [2] D. Frenkel and B. Smit, 2002, Understanding Molecular Simulation, Second Edition. London: Academic Press. [3] R. D. Syarifah, Z. Su’ud, K. Basar, D. Irwanto, S. C. Pattipawaej and M Ilham, 2017, Comparison of Uranium and Thorium Nitride Fuel for 500MWth Gas Cooled Fast Reactor (GFR) Longlife without Refueling, International Journal of Energy Research, Special Issue Paper, page 1-7. [4] R. D. Syarifah, Z. Su’ud, K. Basar and D. Irwanto, 2016, The Prospect of Uranium Nitride (UN-PuN) Fuel for 25-100MWe Gas Cooled Fast Reactor Long Life without Refuelling, Journal of Physics: Conference Series 776 012103, page 1-5. [5] Takaya, 2009, the Corrosion behavior of Al-alloying high Cr-ODS steels in the lead-bismuth eutectic, Journal of Nuclear Materials 386–388, page 507–510. [6] A. Maulana, 2006, Aplikasi Paket Program Moldy Untuk Karakterisasi Sifat Bahan Fe, Pb, Bi Dan Pendingin Reaktor Pb-Bi, Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir XVII, page 119-130. [7] A. Arkundato, Z. Su’ud, M. Abdullah, W. Sutrisno, M. Celino, 2013, Inhibition of iron corrosion in high temperature stagnant liquid lead: A molecular dynamics study, Annals of Nuclear Energy 62, page 298-306 [8] A. Arkundato, Z. Su'ud, M Abdullah, W. Sutrisno, 2013, Molecular dynamic simulation on iron corrosion-reduction in high temperature molten lead-bismuth eutectic,Turkish journal of physics (Turk.J.Phys), 37 (1), page 132-144 [9] A. Arkundato, F. Monado and Z. Su’ud, 2017, Effect of temperature on the corrosion inhibition ofiron in liquid lead using oxygen inhibitor: studied by MD simulation, IOP Conf. Series: Journal of Physics: Conf. Series 853 012046 [10] D. Sapundjiev, V. Dyck, and S. W. Bogaerts, 2006, Liquid metal corrosion of T91 and A316L materials in Pb–Bi eutectic at temperatures 400–600 o C, Corrosion Science, Volume 48, page 577–594. [11] J. M. Haile, 1992, Molecular Dynamics Simulation: Elementary Methods, John Wiley & Sons, Inc. [12] A. Witoelar, 2002, Perancangan dan Analisa Simulasi Dinamika Molekul Ensemble Mikrokanonikal dan Kanonikal dengan Potensial Lennard Jones, Tugas Akhir, Bandung: Institut Teknologi Bandung. 55 Computational and Experimental Research in Materials and Renewable Energy (CERiMRE) Volume 2, Issue 1, page 45 – 55 Submitted : February 2, 2019 Accepted : April 5, 2019 Online : Mei 2, 2019 doi : 10.19184/cerimre.v2i1.20561 [13] H. Gould, J. Tobochnik and W. Cristian, 2007, An Introduction To Computer Simulation Methods: Application To Physical Systems, Third Edition, Boston: Addison-Wesley Longman Publishing. [14] A. Arkundato, 2012, Studi penghambatan korosi logam besi dalam lingkungan pendingin reaktor cepat berbahan Pb dan PbBi cair dengan menggunakan metode somulasi dinamika molekul, Disertasi. Bandung: Institut Teknologi Bandung. [15] B. N. Wira and A. Rian, 2012, Simulasi Sifat Fisis Model Molekuler Dinamik Gas Argon dengan Potensial Lennard-Jone, Makassar: Jurusan Fisika FMIPA Universitas Hasanuddin. [16] A. M. Krivtsov and M Wiercigroch, 2002, Molecular Dynamics Simulation of Mechanical Properties for Polycrytal Materials, Mater. Phys. Mech., volume 3, page 45-51. [17] A. Arkundanto, 2013, Study of Liquid Lead Corrosion of Fast Nuclear Reactor and Its Mitigation by Using Molecular Dynamics Method. International Journal ofApplied Physics and Mathematics. [18] A. D. Kirshenbaum, J. A. Cahill and A. V. Grosse, 1961, The Density Of Liquid Lead From The Melting, Journal of Inorganic and Nuclear Chemistry, volume 22, page 33–38.