2017 Volume 58 Issue 1 Pages 26-32
We have evaluated the diffusion coefficient of phosphorus in α-iron transported in the form of a mixed dumbbell using a kinetic Monte Carlo simulation based on first-principles calculations. The evaluated diffusion coefficient has been compared with both that for the migration mode via octahedral interstitial sites and that for the migration mode of the vacancy mechanism estimated previously using a first-principles-based kinetic Monte Carlo simulation. As a result, we have confirmed that the phosphorus diffusion by the two interstitial migration modes is much faster than that of the vacancy migration mode. In addition, by applying the rate-theory model incorporating the evaluated coefficients to the simulation of irradiation-induced grain-boundary phosphorus segregation, it was made clear that modifications are required to the previous model for GB P segregation and for evaluating GB P coverage.
In order to secure the integrity of nuclear reactor pressure vessel (RPV) steels, the hardening and embrittlement that are brought about by neutron irradiation and thermal aging need to be assessed by experiments and numerical simulations. In particular, grain boundary (GB) embrittlement is a crucial problem because it may lead to degradation of the steel toughness. Phosphorous (P) is known to be an element that causes such GB embrittlement in RPV steels. So far, the experimental results from neutron-irradiated RPV steels have revealed that P atoms tend to be segregated on the GBs1), and first-principles calculations confirm that P atoms segregated on α-iron GBs reduces its cohesive energy2). However, the influences of dose, dose rate, and irradiation temperature on GB P segregation are not clear because of the lack of an experimental database. Hence, estimations by computer simulations are desired. Thus, several papers report theoretical estimations of irradiation-induced GB P segregation by rate-theory models based on the kinetics of point defects3–7). However, the previous estimations are still unsatisfactory. In Ref. 3), the model for P diffusion in FCC nickel4) is applied to the estimation of P segregation in RPV steel. In Ref. 5), the P migration mode by vacancy (V) is ignored. Besides, all or part of the transport coefficients that are evaluated by the analytical method10,11) including several approximations are employed in the model of Refs. 6, 7), as described below.
So far, we have applied a rate-theory model based on a first-principles calculation to neutron-irradiation-induced GB P segregation considering the effect of carbon atoms6) and the V drag effect7) on the GB P segregation. The model uses the transport coefficient of P atoms for the V mechanism, which is evaluated by kinetic Monte Carlo (kMC) simulation8) based on the model that is determined by the first-principles calculation9). Meanwhile, as for the dumbbell mechanism, its model is also obtained from the first-principles calculation5), but the transport coefficients for the model are obtained from the analytical method10,11) using several approximations. Because our previous model7) shows that the transport coefficient of the V mechanism, which is evaluated by the kMC simulation, indicates the V drag effect that the analytical method10,11) cannot treat, it was suggested that we have to obtain more appropriate coefficients by a kMC simulation based on the first-principles calculation also in the case of a dumbbell mechanism. Furthermore, our previous rate-theory model does not include the effect of the interstitial P atoms that are generated by self-interstitial atoms (SIAs) kicking out substitutional P atoms and that migrate via octahedral interstitial sites5). We refer to such an interstitial P atom as an octahedral P atom hereinafter. In the neutron-irradiation environment in which SIAs are abundant, it is necessary to consider octahedral P atoms in the rate-theory model for irradiation-induced P segregation.
In this study, we developed a kMC simulation code based on the dumbbell migration model reported in Ref. 5), and evaluated the diffusion coefficient of both an SIA and a mixed dumbbell composed of an iron (Fe) atom and a P atom. The evaluated diffusion coefficient of the mixed dumbbell was compared with both that of the octahedral P atom, which was calculated analytically using the migration activation energy5), and that of the vacancy mechanism, which was evaluated previously7). The evaluated diffusion coefficient of SIA was also compared with the analytically evaluated value. In addition, we modified our rate-theory model by incorporating the octahedral P atoms and applied it to the simulation of irradiation-induced GB P segregation. As a result, we found that the previously adopted model of the GB is not adequate in the case of including the octahedral P atoms.
We considered two types of the interstitial P atom: an octahedral P atom and a <110> mixed dumbbell. We employed their migration models that are determined from first-principles calculations5). Figures 1(a) and 3 show each model. Figure 2 shows the migration model of an SIA, which is similar to that of a mixed dumbbell because SIAs almost always migrate in the form of a dumbbell. The transition rate in each figure is obtained from the barrier energy $E_i$, which is evaluated by first-principles calculation, through $W_i = \nu_{I} {\rm exp} \left( - \frac{E_i}{k_{B}T} \right)$, where $\nu_{I}$ is the jump frequency, $k_{B}$ is the Boltzmann constant, and $T$ is the absolute temperature. As for the formation and the annihilation of an octahedral P atom, there is no barrier energy, so that in Figs. 1(c) and 1(d), the configuration of the left-hand side always becomes the configuration of the right-hand side.
(a) Migration model of an octahedral P atom. (b) Transition between an octahedral P atom and a mixed dumbbell. (c) Formation of an octahedral P atom. (d) Annihilation of an octahedral P atom with a V. In the figure, the cube is the unit cell of bcc iron. The black dot, gray dot, black diamond, gray diamond, and white square represent a substitutional P atom, an Fe atom at a normal lattice site, an octahedral P atom, an interstitial Fe atom, and a vacancy, respectively. The mixed dumbbell is represented by MID. $W_{i}$ represents the transition rate, and $E_{i}$ is the barrier energy of the corresponding transition.
SIA migration model (a) rigid translation (RT) and translation rotation (TR); (b) jump between second-nearest–neighbor sites; (c) 60° rotation without translation. The meanings of the symbols are the same as those in Fig. 1.
According to Fig. 1(b), an octahedral P atom and a mixed dumbbell are interchangeable. The migration process of an octahedral P atom is shown in Fig. 1(a), and a mixed dumbbell migrates by the dumbbell mechanism shown in Fig. 3. Therefore, we estimate P transport by the interstitial mechanism from evaluating the diffusion coefficient for both migration modes. We evaluate the diffusion coefficient of octahedral P atoms from the analytical form
| \[ D_{int} = \nu_{I} \frac{a_{0}^{2}}{6} {\rm exp} \left( - \frac{E_{oct \to oct}}{k_{B}T} \right),\] | (1) |
Mixed dumbbell migration model (a) RT and TR; (b) jump between second-nearest–neighbor sites; (c) 60° and 90° rotations without translation. The meanings of the symbols are the same as those in Fig. 1.
We adopted the atomistic kMC method in which the time evolution of atomic states is carried out by the residence-time procedure12) shown in Fig. 4. This is so that it can simulate atomic diffusion without restriction of the time-scale problem such as in the molecular dynamics method. In this procedure, $1/k_{tot}$ is calculated as the expectation value for the probability distribution function, $p(t) = k_{tot}{\rm exp} (- k_{tot}t)$, which represents the probability of the time at which the system first escapes from the state. Thus, the time step is independent of which state is selected next.
Time evolution of states by the residence-time procedure r is a random number.
For the kMC simulation of SIA or mixed dumbbell migration, we employed a cubic system of a bcc lattice with 2,000 Fe atoms that has periodic boundary conditions in all three directions. The state in which an Fe atom chosen randomly in the system was replaced with an SIA or a mixed dumbbell was used as the initial state. The transition rate $W_{i}$ in Figs. 2 and 3 was set to $k_{i}$ in Fig. 4. For the jump frequency and the barrier energy, we used $\nu_{I} = 5.16 \times 10^{12}/{\rm s}$5) and the value in Figs. 2 and 3, respectively. In the simulation, the position of dumbbell, ${\bf R}_{i}$, where i = SIA or MID (the mixed dumbbell), was tracked, and the diffusion coefficient was calculated from the following equation8):
| \[ D_{\rm i} = \frac{\left\langle {\bf R}_{i}^{2}\right\rangle}{6t},\] | (2) |
Our rate-theory model employed previously in Refs. 6, 7) consists of the following three equations:
| \[ \begin{split} & \partial N_{P}/\partial t = - \nabla J_{P},\\ & \partial N_{V}/\partial t = - \nabla J_{V} + K - \alpha N_{I} N_{V} - L_{V},\\ & \partial N_{I}/\partial t = - \nabla J_{I} + K - \alpha N_{I} N_{V} - L_{I}. \end{split} \] | (3) |
| \[ \begin{split} & J_{P} = - (D_{VP}^{P} n_{V} + D_{IP}^{P} n_{I}) \nabla n_{P} - D_{PV}^{P} n_{P} \nabla n_{V} - D_{PI}^{P} n_{P} \nabla n_{I},\\ & J_{V} = - D_{V} \nabla n_{V} - D_{PV}^{V} n_{P} \nabla n_{V} - D_{VP}^{V} n_{V} \nabla n_{P},\\ & J_{I} = - D_{I} \nabla n_{I} - D_{PI}^{I} n_{P} \nabla n_{I} - D_{IP}^{I} n_{I} \nabla n_{P}. \end{split} \] | (4) |
The octahedral P atom is not considered in eqs. (3) and (4), so we modified them to include it. The following eq. (5) contains the modified equations corresponding to eq. (3):
| \[ \begin{split} & \partial N_{sub}/\partial t = - \nabla J_{sub} - \beta N_{I} N_{sub} + \gamma N_{V} N_{int},\\ & \partial N_{int}/\partial t = - \nabla J_{int} - \gamma N_{V} N_{int} + \beta N_{I} N_{sub},\\ & \partial N_{V}/\partial t = - \nabla J_{V} - \gamma N_{V} N_{int} + K - \alpha N_{I} N_{V} - L_{V},\\ & \partial N_{I}/\partial t = - \nabla J_{I} - \beta N_{V} N_{sub} + K - \alpha N_{I} N_{V} - L_{I}. \end{split} \] | (5) |
| \[ \begin{split} & J_{sub} = - D_{V\ sub}^{sub} n_{V}\nabla n_{sub} - D_{sub\ V}^{sub} n_{sub} \nabla n_{V},\\ & J_{int} = - D_{oct} \nabla n_{oct} - D_{MID} \nabla n_{MID} \\ & \hphantom{J_{int}} = - (D_{oct} + R_{IM} D_{MID}) \nabla n_{oct},\\ & J_{V} = - D_V \nabla n_{V} - D_{sub\ V}^V n_{sub} \nabla n_{V} - D_{V\ sub}^{V} n_{V} \nabla n_{sub},\\ & J_{I} = - D_{I} \nabla n_{I}, \end{split} \] | (6) |
The modified rate-theory model represented by eq. (5) was applied to the irradiation-induced GB P segregation simulation. The simulation condition was almost the same as in Ref. 7). The calculation region was a one-dimensional region with length L = 1,500 nm. One edge of the region was the inside of grain and the other corresponded to the GB. The boundary condition for GB for each defect was $J_{sub} = 0$, $J_{int} = - (D_{oct} + R_{IM} D_{MID}) n_{oct}/8a_{0}$, $J_{V} = - (D_{V} + D_{sub\ V}^{V} n_{sub}) n_{V}/8a_{0}$, and $J_{I} = - D_{I} n_{I}/8a_{0}$. The Neumann condition for $N_{i}$ of each defect was set at the boundary for the inside of grain. For the initial condition, the values that satisfy the following equations were given homogeneously to the calculation region: $N_{V} N_{int} + \beta N_{I} N_{sub} = 0$, $\gamma N_{V} N_{int} - \alpha N_{I} N_{V} = 0$, and $\beta N_{V} N_{sub} - \alpha N_{I} N_{V} = 0$. The initial values were obtained by solving the three equations iteratively from $n_{sub}$ = 0.013 at%, $n_{oct}$ = 0, $n_{V} = {\rm exp} (E_{f}^{V}/k_{B}T)$, and $n_{I} = {\rm exp} (E_{f}^{I}/k_{B}T)$. Here, $E_{f}^{V} = 2.12$ eV and $E_{f}^{I} = 3.64$ eV are the formation energies of V and SIA, respectively. We assumed that the generated SIAs become self-interstitial dumbbells immediately. The sink terms in the third and fourth equations of eq. (5), namely $L_{V} = - (D_{V} + D_{sub\ V}^{V} n_{sub}) n_{V} b_{V}^{2}$ and $L_{I} = - D_{I} n_{I} b_{I}^{2}$, respectively, and the sink strengths $b_{V,I}$ were chosen so that the simulation results reproduced the experimental result described in Ref. 7). In the simulation, GB P segregation was evaluated using the following equation that is based on the one proposed in Ref. 13):
| \[ \begin{split} C_{P}^{GB} =& \frac{2}{a_{0}} \int_{a_0/2}^{\infty} dx \left\{ N_{sub}^{init} + N_{int}^{init} - N_{sub} (x,t) - N_{int}(x,t) \right\} \\ & + C_{P}^{GB,init}. \end{split} \] | (7) |
Figure 5 shows the temperature dependence of diffusion coefficients for the mixed dumbbell and octahedral P atom. The diffusion coefficient for P migration by the V mechanism is also shown in Fig. 5. Since in the V mechanism, substitutional P atoms migrate by interacting with Vs as described in the Appendix, P diffusion depends on V and P concentrations. Therefore, $D_{V\ sub}^{sub} n_{V}$ and $D_{sub\ V}^{sub} n_{sub}$ in the first equation of eq. (6) are shown in Fig. 5. The second one results from the V drag effect. Here, we use $D_{V\ sub}^{sub}$ and $D_{sub\ V}^{sub}$ of the case that includes the second-nearest-neighbor binding state and their detail is described in the Appendix. For $n_{V}$, we used the equilibrium value evaluated from $E_{f}^{V}$ and the value estimated by the simulation using eq. (3), and $n_{P}$ is the value of the simulation. According to Fig. 5, the diffusion coefficients of the mixed dumbbell and octahedral P atom are almost identical to each other. They are larger than that for the V mechanism, although that value depends on the P and V concentrations. The diffusion of P atoms by the interstitial mechanism is much faster than by the V mechanism. However, since mixed dumbbells and octahedral P atoms are generated by the processes shown in Figs. 1(b) and 1(c), P transport by mixed dumbbells and octahedral P atoms depends on the concentration of SIA as P transport by the V mechanism depends on the V concentration. Hence, all mechanisms need to be considered for P transport in the simulation for irradiation-induced GB P segregation.
Temperature dependence of diffusion coefficient. The comparison between the interstitial mechanism and the V mechanism is shown. All the coefficients in the figure are evaluated by kMC simulation.
Figure 6 shows the comparison of the SIA diffusion coefficient between the kMC simulation and the analytical method. In the latter, only the rigid translation mode shown in Fig. 2(a) is considered5). In Fig. 6, the SIA diffusion coefficient evaluated by the kMC simulation is slightly smaller than that by the analytical method. This result shows that the suppression by the rotation modes is more influential on diffusion than is the acceleration by the jump between second-nearest-neighbor sites.
Temperature dependence of diffusion coefficient of SIA. The comparison between the kMC simulation and the analytical method is shown.
Figure 7 shows the GB P coverage against dose that was simulated under the following conditions7): a dose rate of 1.7 × 10−10 dpa/s (1.0 × 1011 n/cm2/s) and a temperature of 563 K (290℃). The experimental results shown in Fig. 7 are the same as those in Refs. 6, 7). The sink strength was $b_{V}^{2} = b_{I}^{2} = 2.3$ × 1018/m2, which was chosen so that the simulated GB P coverage reproduced one of the experimental results. In Fig. 7, although the sink strength is larger than the value used in the previous simulations that were carried out using the model of eq. (3)6,7), the simulation seems to reproduce the tendency of several experimental results. According to the distribution of $N_{sub}$ and $N_{int}$ in Fig. 8, the increase in GB P coverage which was calculated using eq. (7) results from the homogeneous decrease in $N_{sub}$ in the calculation region. This homogeneous decrease is brought about by the disappearance of octahedral P atoms at the GB as with V and SIA, which results from the boundary condition at the GB for the octahedral P atom. As described in Section 2.3, since the octahedral P atoms can migrate by themselves, we used $J_{int} = - (D_{oct} + R_{IM} D_{MID}) n_{oct}/8a_{0}$ as the boundary condition at the GB which is similar to that of V and SIA that also migrates by themselves14). Since substitutional P atoms become octahedral P atoms by the process shown in Fig. 1(c), substitutional P atoms gradually decrease due to the decrease of octahedral P atoms, and eventually both substitutional and octahedral P atoms disappear from the calculation region. Therefore it is considered that the GB P coverage calculated by eq. (7) in Fig. 7 continues to increase unrealistically until P atoms disappear from the calculation region.
GB P coverage versus dose simulated using the model of eq. (5) the dose rate is 1.7 × 10−10 dpa/s, temperature is 290℃, and sink strength is $b_{V}^{2} = b_{I}^{2} = 2.3$ × 1018/m2.
Distribution of $N_{sub}$ and $N_{int}$ at 0.11 dpa in the calculation region.
In the simulation using the previous model of eq. (3) without octahedral P atoms, substitutional P atoms transported to the GB by V and SIA accumulate at the GB due to the disappearance of the V and SIA at the GB14). This is realized by the following boundary condition at the GB: $J_{P} = 0$, $J_{V} = - (D_{V} + D_{P\ V}^{V} n_{P}) n_{V}/8a_{0}$, and $J_{I} = - (D_{I} + D_{P\ I}^{I} n_{P}) n_{I}/8a_{0}$. The accumulation of P atoms makes the concentration of P atoms at the GB larger than that in the vicinity of the GB, so that $\nabla n_{P}$ in the first equation of eq. (4) becomes large and then the flow of P atoms from the GB to the inside of grain is formed. In other words, P atoms do not disappear at the GB and can return to the inside of grain. Thus, the GB P coverage evaluated by eq. (7), which calculates the decrease of P atoms in the inside of grain, does not continue to increase unrealistically and can become constant when the P transport to the GB by V and SIA and the flow of P atoms to the inside grain become steady.
From the abovementioned consideration, in the model including the octahedral P atom, it is necessary to make octahedral P atoms stay at the GB instead of disappearance, namely we had better model the GB which can trap octahedral P atoms. In addition, it is also necessary that the GB can detrap octahedral P atoms. As for the substitutional P atom, strictly speaking, the accumulation at the GB is not the trapping at the GB. Thus, substitutional P atoms had better be also trapped and detrapped at the GB as with the octahedral P atom. In the case incorporating the processes of trapping and detarapping at the GB into the model, we consider that the GB P coverage should be evaluated by counting the trapped P atoms at the GB instead of using eq. (7).
We evaluated the diffusion coefficients of SIA and mixed dumbbell in α-iron using kMC simulation based on the first-principles calculation. The evaluated diffusion coefficient of mixed dumbbells was compared with that for octahedral P atoms, which was obtained by the analytical method, and that for the V mechanism, which was evaluated previously by kMC simulation. The diffusion coefficient of SIA was compared with that obtained using the analytical method. In addition, the rate-theory model was modified by incorporating the octahedral P atom and the evaluated diffusion coefficients, and the modified model was applied to the simulation of irradiation-induced GB P segregation. We obtained the following results.
This work was supported by JSPS Grant-in-Aid for Scientific Research (C) Grant Number 15K06429.
Here, we describe the estimation of the partial diffusion coefficient (PDC) of the V mechanism and the cause of the V drag effect, neither of which are clarified in Ref. 7).
Figure A1 shows the model of that V mechanism that is determined from the first-principles calculation9). In the model, a substitutional P atom can migrate by exchanging a site with a V at a nearest-neighbor site. There are two kinds of binding states between a P atom and a V: one is the first-nearest-neighbor (1NN) binding state whose association and dissociation are shown in Figs. A1(c)–(e), and the other is the second-nearest-neighbor (2NN) binding state whose association and dissociation are shown in Fig. A1(f), which cannot be treated by the analytical method10,11). The 2NN binding state is considered as the cause of the V drag effect. However, even without the 2NN binding state, the V drag effect occurs7).
Model of the V mechanism (a) V migration, (b) P migration, (c)–(e) association and dissociation of P and V in the 1NN distance, (f) association and dissociation of P and V in the 2NN distance. The meanings of the symbols are the same as those in Fig. 1.
The PDC of the V mechanism is estimated by the kMC simulation, which is almost the same as for the interstitial dumbbell. However, the PDC is calculated from the transport coefficient $L_{ij} = \frac{1}{V_{s}k_{B}T} \frac{\langle {\bf R}_{i}{\bf R}_{j}\rangle}{6t}$, where ${\bf R}_i$ (i = P or V) is the displacement of P or V and $V_{s}$ is the system volume, using the following equations; $D_{PV}^{P} = \frac{k_{B}TL_{VP}}{(N_{l}/V_{s})n_{P} n_{V}}$; $D_{VP}^P = \frac{k_{B} TL_{PP}}{(N_{l}/V_{s})n_{P}n_{V}}$; $D_{PV}^V = \frac{k_{B} TL_{VV}}{(N_{l}/V_{s})n_{P}n_{V}} - \frac{D_{V}}{n_{P}}$; and $D_{VP}^{V} = (1 + z) D_{V} + D_{PV}^{P}$, where $N_{l}$ is the number of perfect lattice sites in the system, and z is the number of binding states of P and V. We have z = 8 when only the 1NN binding state is considered, and z = 14 when the 2NN binding state is included. Here, $D_{V}$ is calculated analytically. The estimated PDC is given in Table A1. The PDC of the case including the 2NN binding state was used for the rate-theory model in Section 2.3.
| PDC, $D_{jk}^{i}$/m2/s | ||
|---|---|---|
| With 2NN binding state | Without 2NN binding state | |
| $D_{sub\ V}^{sub}$ | 85.3 (±1.2) | 3.68 (±0.12) |
| $D_{V\ sub}^{sub}$ | 85.6 (±1.0) | 3.95 (±0.049) |
| $D_{V\ sub}^V$ | 85.7 (±1.2) | 3.92 (±0.12) |
| $D_{sub\ V}^V$ | 84.9 (±1.7) | 5.69 (±0.68) |
According to Table A1, $D_{sub\ V}^{sub}$ is a positive value. Hence, from $J_{sub}$ in eq. (6), P atoms are transported in the direction of V diffusion caused by the concentration difference. Specifically, the V drag effect can be described by the model. Meanwhile, $D_{sub\ V}^{sub}$, obtained using the analytical method10,11), is a negative value7). It was considered previously that the V drag effect requires the 2NN binding state in Fig. A1(f)8). However, in Table A1, $D_{sub\ V}^{sub}$ for the case excluding the 2NN binding state is a positive value. Therefore, the V drag effect also occurs in the case excluding the 2NN binding state. The cause of this inconsistency is not clarified in Ref. 7), but it can be explained by considering the energy configuration for a V around a P atom shown in Fig. A2.
Energy configuration of V around a P atom (a) the case including the 2NN binding state, (b) the case excluding the 2NN binding state. The encircled number means the distance from the P atom. In (b), the 2NN binding state is excluded by $E_{5}^{V} = E_{6}^{V} = E_{0}^{V}$. In these figures, although $E_{3}^{V}$ and $E_{4}^{V}$ are used representatively, the situation is the same as the case using ${E'}_{3}^{V}$, ${E'}_{4}^{V}$, ${E''}_{3}^{V}$, and ${E''}_{4}^{V}$.
As seen in Fig. A2(a), in the case including the 2NN binding state, both the 1NN and 2NN sites of a P atom become the binding state because the energy level of V at these sites is lower than that at a perfect lattice site. Meanwhile, in the case excluding the 2NN binding state, the 2NN site of a P atom is not the binding state as shown in Fig. A2(b). However, in both cases, the energy barrier in the direction in which a V approaches a P atom is lower than that in the direction in which a V separates from a P atom. This situation raises the possibility that a V stays at the 1NN and 2NN sites for the P atom, so that the V drag effect occurs in both cases. Hence, the V drag effect does not result from the 2NN binding state but rather from the barrier-energy height for V migration around a P atom.
