MATERIALS TRANSACTIONS
Online ISSN : 1347-5320
Print ISSN : 1345-9678
ISSN-L : 1345-9678
Molecular Dynamics Investigation of the Effect of Lanthanum on the Diffusivity of Niobium in Austenite
Haiyan WangZhaofeng YaoXueyun GaoPeng CuiHuiping Ren
Author information
JOURNAL FREE ACCESS FULL-TEXT HTML

2018 Volume 59 Issue 5 Pages 706-711

Details
Abstract

By employing the force-matching method, an embedded-atom-method potential for Fe–Nb–La dilute solid solution has been constructed by fitting to a set of density-functional theory (DFT) data of reference structures. The developed potential reproduces the properties of the reference structures in good agreement. Using the developed potential, molecular dynamics simulations were performed to investigate the diffusion coefficients of Nb in Fe–Nb and Fe–Nb–La systems, respectively. The results indicate that the addition of La suppresses the diffusion of Nb in fcc Fe, the activation energy of Nb was determined to be 47551 J/mol in the absence of La while 82672 J/mol in the presence of La, and consequently leads to the longer incubation time of NbC precipitation in the La micro-alloyed austenite.

1. Introduction

In recent years, rare earth (RE) microalloying in steel has received an increasing interest.13) Many studies have been carried out to investigate the effect of RE in steels,48) such as the purification, solidification, high-temperature oxidation resistance, corrosion resistance and phase transformation. RE elements are beneficial to the mechanical performance of steels, especially to the toughness by grain refinement and suppression of grain boundary embrittlement,9) which are involved with the segregation of RE to the austenite grain boundaries during thermal deformation process.

It is well known that hot rolled microalloyed steels can be strengthened by a combination of grain refinement and precipitation strengthening, using microalloying elements niobium (Nb), titanium (Ti) and vanadium (V), individually or in combination.1012) Previous experiments prove that the addition of lanthanum (La) suppresses the precipitation of NbC in austenite, and the first-principle calculations indicate that the solubility of niobium and carbon in fcc Fe are increased, and the chemical potential of both are decreased in the presence of La.13) For the case of NbC precipitation in ferrite steels, our diffusion couple experiment reveals that Nb diffuses slightly faster in the presence of La, and consequently leads to the faster precipitation kinetics of NbC.14)

Despite the progress in RE containing steels research so far, it still remains to be investigated, since La generally cannot be detected directly in steels by conventional experimental methods due to its low solubility in steel. Although the advances in computational hardware and algorithms allow the application of first-principle calculations to the related investigation, it is beyond the capability of ab initio method in the case of precipitation and dissolution of NbC processes which involve long-ranged strain fields, long-ranged composition fluctuations and large numbers of atoms. The molecular dynamics (MD) method was demonstrated to be useful for simulating structures with millions of atoms, and studying the mechanical and thermodynamical properties of materials.15) However, to apply MD to the La containing steel, an appropriate interatomic potential is required. Unfortunately, to the best of our knowledge, there is no systematical study of interatomic potential for Fe–Nb–La ternary available in literature.

In this work, we developed an interatomic potential for the Fe–Nb–La system based on the EAM model. The parameters of the potential functions were optimized so that the ab initio calculation results were reproduced with an acceptable accuracy. Then, using the developed potential, we performed a molecular dynamics simulation and investigated the effect of La addition on the diffusivity of Nb in fcc Fe matrix. Finally, we investigated the precipitation kinetics of NbC in austenite region using stress relaxation experiments.

2. Methodology

The embedded atom method (EAM)16) is a widely used atomic level semiempirical model for metals. It implicitly includes many-body interactions by a term which depends on the environment of every atom, and has been used successfully for many studies of metallic materials. The potential energy of a system described with the EAM method can be represented as follows:   

\begin{equation} U = \sum_{i < j}\varphi_{\alpha\beta}(r_{ij}) + \sum_{j}F_{\alpha}(\bar{\rho}_{i})\ \text{with}\ \bar{\rho}_{i} = \sum_{j\neq i}\rho_{\alpha}(r_{ij}) \end{equation} (1)
where φαβ denotes the pair interactions between atom i and j at the distance of rij = |rjri|, α and β represent the element types. The function $F_{\alpha }(\bar{\rho }_{i})$ is the embedding energy of atom i in the effective electronic density $\bar{\rho }_{i}$ which is the sum over the contributions ρα(rij) created by the neighboring atoms situated within the cut-off radius from the atom i.

The force-matching method proposed by Ercolessi and Adams17) is used to develop the Fe–Nb–La EAM potential for a classical MD simulation. This method provides an effective way to construct physically justified interatomic potentials even for complicated multicomponent systems by making use of large amounts of information obtained via ab initio calculations. By adjusting the set of the potential functions φαβ(rij), ρα(rij) and $F_{\alpha }(\bar{\rho }_{i})$, the force-matching method reproduce the energies, forces and stresses per atom computed for a serials references structures which represent various phases existing in the investigated system.

2.1 Fitting procedure

The force-matching algorithm is implemented in the potfit package of Brommer and Gähler,18,19) and has been used to generate effective potentials from ab initio data.20,21) For the present work, all free parameters of the analytic functions were fitted to the ab initio reference database which containing fully relaxed structures at T = 0 K, snapshots from ab initio molecular dynamics simulations at higher temperatures, and configurations containing vacancy point-defect.

The potential functions are represented by analytic model where a set of parameters is given to describe the potential, and the parameters are adjusted to optimally reproduce the quantities calculated from ab initio in a series of reference calculations. The optimization procedure is realized by minimizing the deviations between the reference values (forces, energies and stresses) and the respective calculated values from the EAM potential. The objective function Z to be minimized is defined as   

\begin{equation} Z = Z_{\text{F}} + Z_{\text{C}} \end{equation} (2)
where   
\begin{equation} Z_{\text{F}} = \sum_{i = 1}^{N_{\text{a}}}\sum_{\alpha = x,y,z}W_{i}\frac{(F_{i\alpha}^{\text{EAM}} - F_{i\alpha}^{\text{DFT}})^{2}}{(F_{i\alpha}^{\text{DFT}})^{2} + \varepsilon_{i}} \end{equation} (3)
  
\begin{equation} Z_{\text{C}} = \sum_{i = 1}^{N_{\text{c}}}W_{i}\frac{(A_{i}^{\text{EAM}} - A_{i}^{\text{DFT}})^{2}}{(A_{i}^{\text{DFT}})^{2} + \varepsilon_{i}} \end{equation} (4)
Equation (3) is related to the relative deviation between the EAM forces $F_{i\alpha }^{\text{EAM}}$ and ab initio forces $F_{i\alpha }^{\text{DFT}}$ for Na atoms in the fitting database. Equation (4) is the relative deviation of the EAM energies and stresses from the ab initio values, where Ai denotes the energy of a certain configuration or one of its six stress tensor components, and Nc represents the number of the reference values of the energy and stress. εi is a small value that prevents extremely small denominators, and Wi is the weights of the different terms.

2.2 Reference data

To compute the reference database, we need to create a set of reference structures (so-called configurations). For our case, to develop a potential to describe accurately the interactions between Nb and La in fcc Fe matrix, an extensive set of configurations which represents as many as kinds of interesting structures need to be prepared. The structures used as reference data are listed in Table 1, this set of configurations present the pure Fe element, the binary (Fe–Nb and Fe–La), and the ternary (Fe–Nb–La) system at various temperatures. For the body center cubic and face center cubic structures, the computations were performed within 4 × 4 × 4 and 3 × 3 × 3 supercells, respectively. For the binary and ternary configurations, we substituted 2 to 4 Fe atoms in fcc Fe supercell with solute (Nb and La) atoms. Additionally, we also included several structures containing vacancy point-defect. The target potential for Fe–Nb–La is intended to be used in the simulations of the diffusion of Nb in fcc Fe matrix at elevated temperature, therefore it is important that this potential describes well the high temperature status of Fe–Nb–La with face center cubic structure. Accordingly, we used more high temperature reference configurations with fcc structure to train the fitted potential. Furthermore, we calculated bcc based system at 2000 K to simulate the molten state of that.

Table 1 Configurations used to fit the Fe–Nb–La potential. n denotes the number of configuration of the corresponding structure, Natoms is the atom number in the corresponding configuration, V/V0 represents the ratio of the supercell volume V to the equilibrium volume V0 of fcc structure from the ab initio calculations, T is the absolute temperature of the ab initio MD simulations.

The reference data was calculated using the Vienna ab initio simulation package (VASP) with the projector augmented wave (PAW) method and the generalized gradient approximation of Perdew-Burke-Ernzerhof functional (GGA-PBE).22) The plane-wave cutoff energy was chosen as 350 eV and the Brillouim zones were sampled with 2 × 2 × 2 k-point meshes. We considered the ferromagnetic (FM) and the single layer anti-ferromagnetic (AFM1) states for bcc Fe and fcc Fe supercell, respectively. The reference database contains 92 configurations with a total of 10506 atoms, and consists of 31518 forces, 92 energies and 552 stresses.

2.3 Atomistic simulation

The MD simulations were performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS) package.23) The EAM potential developed in this work was used to describe the interaction between Fe, Nb and La atoms. To investigate the effect of La on the diffusion of Nb in fcc Fe, we constructed Fe–1%Nb and Fe–1%Nb–0.2%La (at%) with a cubic periodic simulation cell of 16 × 16 × 16$a_{0}^{3}$, respectively, a0 is the lattice parameter of fcc Fe, and 0.1 at% vacancy was considered in both alloys. Nb, La and vacancies were randomly distributed over the lattice sites of the simulation cell.

An energy minimization by using the conjugate gradient method was carried out to obtain the equilibrium structure. The atomistic models were relaxed for 50 ps with NPT ensemble which was at the desired temperature under zero-pressure state. Subsequently, the simulations were performed with NVT ensemble under the same temperature and pressure conditions for 2000 ps. A constant integration time step of 1 fs was used for all simulations.

3. Results and Discussion

3.1 Optimized potential

The optimization procedure dealt with 12 functions, i.e. 6 pair potentials φαβ, 3 transfer functions ρα and 3 embedding functions Fα, which define the finally matched potential, and are plotted in Fig. 1. For the investigation of the effect of La on the diffusion of Nb in fcc Fe, the interaction between Nb and La is a crucial factor that ought to be studied.24) Figure 1(b) shows that there is an attractive interaction between Nb and La atoms at the distance exceeding 2.6 Å. From our previous ab initio calculations,13) the interaction of Nb and La is attractive over the range of the second nearest neighbor to the sixth nearest neighbor shell. The obtained Nb–La pair potential is comparable with the ab initio results.

Fig. 1

EAM potential functions of Fe–Nb–La ternary derived through the potfit package. (a) and (b) show the pair potentials φαβ between Fe, Nb and La, (c) and (d) represent the electron densities ρ(r) and embedding functions $F_{\alpha }(\bar{\rho })$, respectively.

Figure 2 shows the distributions of energies, stresses and forces obtained using EAM potential and ab initio density functional theory (DFT) values, respectively. Abscissas and ordinates in Fig. 2 represent the physical quantities obtained by EAM potential and ab initio calculation, respectively. The comparison shows that our potential reproduces the properties of the reference structures in good agreement.

Fig. 2

Scatter plot of energies (a), stresses (b) and forces (c) with the reference data of ab initio density functional theory (DFT) on the vertical axis and the EAM values on the horizontal axis.

By using the obtained potential, we calculated the general characteristics of pure Fe which could be compared with the reported ab initio and experimental data,25,26) and the results are summarized in Table 2. Additionally, we calculated these basic properties using the MEAM potential of Fe–Nb binary developed by I. Sa et al.,27) and the results are listed in the last column of Table 2. It can be seen that the present EAM potential is able to produce better matches with the ab initio calculations and experiments, as well as the results from MEAM model, except the elastic constant C44 that the Fe–Nb MEAM produced a rather lower value.

Table 2 Lattice constants (Å) and elastic constants (GPa) of pure Fe calculated via the developed EAM in comparison with the available data.

3.2 Diffusion of Nb in fcc Fe

The diffusion coefficient of Nb atoms in fcc Fe matrix can be obtained from the MD simulations using the following equation:28)   

\begin{equation} D = \frac{\langle R^{2}(t) \rangle}{6t} \end{equation} (5)
where t denotes the simulation time, ⟨R2(t)⟩ represents the mean square displacement (MSD) of Nb atoms as a function of time and is defined as   
\begin{equation} \langle R^{2}(t) \rangle = \frac{\displaystyle\sum_{j}^{N_{\text{Nb}}}(r_{j}(t) - r_{j}(0))^{2}}{N_{\text{Nb}}} \end{equation} (6)
where NNb is the number of the Nb atoms in the system, rj(t) denotes the location of the j-th Nb atom at time t.

Figure 3 shows the MSD of Nb atom as a function of time in the system of Fe–1%Nb and Fe–1%Nb–0.2%La, respectively. MD simulations were performed at 900°C, 950°C, 1000°C and 1050°C for times up to 2000 ps. As shown in Fig. 3, the MSD of Nb keep increasing over the simulation time, and liner trends can be clearly observed from the beginning of the simulations for these four temperatures, which guarantees a good estimation of the diffusion coefficient of Nb. At the same temperature, Nb has a slower MSD growth rate for the La-containing system, e.g. the Nb MSD of La free system reaches 3.2 Å2 after 2000 ps at 900°C, while that of La-containing system reaches 0.7 Å2.

Fig. 3

MSD of Nb atom vs. time at 900°C, 950°C, 1000°C and 1050°C in the system of Fe–1%Nb (a) and Fe–1%Nb–0.2%La (b), respectively.

The Nb diffusion coefficients of Nb at the certain temperature can be obtained via the slope of the MSD curve versus time. And the temperature dependence of the obtained diffusion coefficient does obey the Arrhenius equation:   

\begin{equation} D = D_{0}\exp \left(-\frac{Q}{RT}\right) \end{equation} (7)
where D is the diffusion coefficient at the temperature T, D0 denotes the pre-factor, Q denotes the activation energy, and R is the universal gas constant.

To obtain the activation energy and the pre-factor of the diffusion coefficient, we calculated the logarithms of the diffusion coefficients versus 1/RT, as shown in Fig. 4. It can be seen that the temperature dependence of the calculated diffusion coefficient is in good agreement with the Arrhenius law, i.e. one can observe a linear dependence of the logarithm of the calculated diffusion coefficient on 1/RT. By taking the double logarithm of eq. (7), we can deduce the following equation:   

\begin{equation} \ln D = \ln D_{0} - \frac{Q}{RT} \end{equation} (8)
This equation shows a linear relationship between ln D and 1/T, and by linearly fitting the Arrhenius plot in Fig. 4, we can obtain the activation energy and pre-factor according to the slope and intercept of the lines, respectively. Therefore, according to the fitting lines shown in Fig. 4, the diffusion equations for Nb in fcc Fe with La free and La presence, respectively, are:   
\begin{equation} D = 3.66\times 10^{-10}\exp \left(-\frac{47551}{RT}\right) \end{equation} (9)
  
\begin{equation} D = 3.81\times 10^{-9}\exp \left(-\frac{82672}{RT}\right) \end{equation} (10)

Fig. 4

Arrihenius plots for Nb diffusion coefficients versus temperature.

The activation energy of Nb diffusion is determined to be 47551 J/mol in the absence of La, while this energy is 82672 J/mol in the presence of La. Therefore, the addition of La suppresses the diffusion of Nb atoms in the fcc Fe matrix.

3.3 Precipitation kinetics of NbC

In the case of precipitation kinetics of NbC, the nucleation rate is given by:29)   

\begin{equation} \frac{dN}{dt} = N_{0}Z\beta'\exp\left(-\frac{\Delta G}{k_{B}T}\right)\exp \left(-\frac{\tau}{t}\right) \end{equation} (11)
where N is the number of precipitates per unit volume, t represents the time, N0 is the number of available sites for heterogeneous nucleation, Z represents Zeldovich factor, β′ is the rate at which atoms are added to the critical nucleus, ΔG is the corresponding activation energy for nucleation, kB denotes Boltzmann constant, T is the temperature, and τ represents the incubation time. And β′ can be obtained by the following:   
\begin{equation} \beta' = \frac{4\pi R_{c}^{2}D_{\text{Nb}}C_{\text{Nb}}}{a_{0}^{4}} \end{equation} (12)
where Rc is the critical radius for nucleation, DNb is the diffusion coefficient of Nb, CNb is the Nb concentration, and a0 is the lattice parameter. The evolution of NbC precipitate radius is expressed by:30)   
\begin{equation} \frac{dR}{dt} = \cfrac{D_{\text{Nb}}C_{\text{Nb}} - C_{\text{Nb}}^{e}\exp \biggl(\cfrac{R_{0}}{R}\biggr)}{RC_{\text{Nb}}^{p} - C_{\text{Nb}}^{e}\exp \biggl(\cfrac{R_{0}}{R}\biggr)} + \frac{dN}{Ndt}(\eta R_{c} - R) \end{equation} (13)
where R is the radius of precipitate, $C_{\text{Nb}}^{e}$ and $C_{\text{Nb}}^{p}$ are the equilibrium concentration of Nb in the system and the Nb concentration in the precipitate, respectively, Rc denotes the capillarity radius, and η represents the size factor of the newly formed nuclei.

It can be seen from eqs. (11)(13) that there is a positive correlation between the Nb diffusion coefficient and the precipitation kinetics of NbC, i.e. the reduction of Nb diffusion coefficient will slow down the precipitation process of NbC.

In the experimental study of the effect of La on the NbC precipitation, we investigated the relation of the La addition with the NbC precipitation kinetics in fcc Fe. Three high-Nb steels were prepared in a 10 kg vacuum induction melted furnace, and the chemical compositions of these three steels are (mass%): (i) Fe–0.043C–0.247Si–1.83Mn–0.085Nb–0.007P–0.005S, (ii) Fe–0.042C–0.243Si–1.86Mn–0.087Nb–0.008P–0.005S–0.0048La, and (iii) Fe–0.043C–0.246Si–1.85Mn–0.086Nb–0.007P–0.005S–0.0052La. The ingots were homogenized at 1250°C for 1 h and then hot-rolled to the thickness of 12 mm. The test specimens of 6 mm diameter and 15 mm length were machined from the hot rolled plates with the long axis parallel to the normal direction of the plate. The stress relaxation method31) was employed on a Gleeble-1500 thermo-mechanical simulator to investigate the precipitation behaviors of NbC. The specimens were reheated at 1250°C for 5 min, subsequently cooled to the deformation temperature (925, 950, 975 and 1000°C) at a rate of 20°C/s, and held for 30 s to ensure the temperature uniformity, followed by a single deformation pass with a strain of 0.3 and a strain rate of 10 s−1. After isothermal holding for 1000 s, the specimens were water quenched to room temperature. The stress relaxation curves were obtained during the isothermal holding. During the stress relaxation, the stress evolution is related to the recover, recrystallization and the strain induced precipitation. Specifically, the recovery and recrystallization process leads to the stress decreasing, while the occurrence of NbC precipitation will slow down the downward trend in stress by suppressing the recovery and recrystallization. Therefore, the point of the stress curve at which the decreasing trends slowing down is corresponding to the start point of precipitation, and the trend of downward comes back at the point which is regarded as the finished point of precipitation.

From the stress relaxation curves, we established the precipitation-time-temperature (PTT) diagrams of the experimental steels, as shown in Fig. 5. It can be seen that the PTT curves all are of ‘C’ shape, and the nose in the start curves correspond to the minimum incubation time of precipitation. For the three experimental steels, the minimum incubation time of No. 1 steel is about 29 s, and this value increases to 41 s for No. 2 steel and 43 s for No. 3 steel, respectively. The incubation time of the La addition steels are longer than the La free one.

Fig. 5

PTT diagrams of the three experimental steels, where Ps and Pf represent the start and finish time of the strain induced precipitation, respectively.

Our previous study13) indicates that there exists a long range attractive interaction between La atom and its neighbor Nb and C atoms, which leads to the increasing of the Nb and C solubility in fcc Fe. And the MD simulations in the present work indicate that the addition of La suppresses the Nb diffusion in fcc Fe, and consequently slows down the nucleation and growth rates of NbC precipitation. These effects of La contribute to the extension of the incubation time for NbC precipitation.

4. Conclusions

In the present work, to investigate the effect of La on the diffusion of Nb in fcc Fe via molecular dynamics method, we developed the EAM potential for Fe–Nb–La system by employing force-matching method to fit to a set of reference structure data from ab initio calculations. The newly fitted potential reproduces the parameters and elastic constants of bcc Fe and fcc Fe in good agreement with the reported data. Using the developed potential, we investigated the mean square displacement of Nb in Fe–1%Nb and Fe–1%Nb–0.2%La systems, respectively. The results indicate that the addition of La suppresses the diffusion of Nb in fcc Fe, and the activation energy of Nb diffusion in fcc Fe is 47551 J/mol without La and 82672 J/mol with La presence. This is one of the reasons that the incubation time of NbC precipitation in the La micro-alloyed γ-Fe longer than that in La free one.

Acknowledgements

The authors are grateful for the financial support of the National Natural Science Foundation of China (No. 51101083), and the Inner Mongolia Autonomous Region Natural Science Foundation of China (No. 2017MS0508).

REFERENCES
 
© 2018 The Japan Institute of Metals and Materials
feedback
Top