Evaluation and Modeling of Anisotropic Stress Effect on Hydrogen Diffusion in Bcc Iron

Abstract

The diffusion behavior of interstitial hydrogen in steel should be clarified to reveal the mechanism of hydrogen embrittlement. In this study, we performed molecular dynamics (MD) simulations to elucidate the relationship between various stress conditions and the diffusion coefficients of interstitial hydrogen in body-centered-cubic (bcc) iron. The results reveal that the diffusion coefficients are independent of isotropic stress and exhibit strong anisotropic stress dependence under uniaxial stress along the ⟨100⟩ direction. The stress-induced deformation of the atomic structure around the octahedral interstitial site (O site) was examined to elucidate the origin of the anisotropy of the stress dependence. Moreover, a model that predicts the activation enthalpy was derived by quantitatively evaluating the deformation around the O site. The activation enthalpy predicted by the model was consistent with the results of MD simulations when the deformation around the O site was not substantial.

 

This Paper was Originally Published in J. Soc. Mater. Sci., Japan 69 (2020) 119–125. All captions of figures and tables were slightly modified.

Trajectory of hydrogen diffusion in body-centered-cubic iron.

1. Introduction

Environmental issues such as global warming are attracting attention and are spurring interest in new energy alternatives to fossil fuels. Among them, hydrogen is expected to be used as an alternative energy because it generates only water, does not generate carbon dioxide, and has little environmental load in combustion, which is the most important reaction in energy utilization.¹⁾

When hydrogen is used as a secondary energy source, it is typically stored at high pressures, and steel materials are expected to be used as structural members for hydrogen storage on the basis of cost. In a hydrogen environment, the strength of not only steels but also various metals is known to decrease because of hydrogen embrittlement;^2,3) however, the embrittlement mechanism is unclear. Clarifying the mechanisms of hydrogen embrittlement in steel materials requires knowledge of the exact distribution of hydrogen. Therefore, understanding the diffusion behavior of hydrogen has become a major issue.

Hydrogen diffusion in metals is described by the diffusion equation (eq. (1)) that includes the isotropic-stress-gradient term (second term on the right side):^4–6)   

\begin{equation} \boldsymbol{{J}} = - D\nabla c + \frac{DcV_{\text{H}}}{RT}\nabla \sigma_{\text{h}} \end{equation} (1)

where J is the diffusion flux of hydrogen, D is the diffusion coefficient, c is the hydrogen concentration, V_H is the molar volume of hydrogen, R is the gas constant, T is temperature, and σ_h is the isotropic stress. Equation (1) shows that hydrogen transport occurs in the direction of higher isotropic stress. However, the dependence of diffusion coefficient D in eq. (1) on the stress has not been sufficiently studied. Taketomi and Matsumoto have observed strong anisotropy in the hydrogen diffusion behavior around dislocations in iron and have speculated that this anisotropy is caused by the contribution of anisotropic stress.⁷⁾ In addition, when the hydrogen microprint method is used for iron, silver particles precipitated at the position of hydrogen release have been reported to not be evenly distributed at grain boundaries or in grains. This uneven distribution may reflect the anisotropy of hydrogen diffusion due to residual stress per grain.^8,9) However, a detailed evaluation of the dependence on anisotropic stress and its reflection in the diffusion equation has not yet been reported.

Regarding the temperature dependence of the diffusion coefficient, quantum effects such as the tunneling effect are known to occur in the low-temperature region, and the diffusion mechanism changes remarkably.¹⁰⁾ However, the classical treatment is effective because hydrogen diffusion via the classical thermal activation process is dominant at temperatures higher than room temperature. The purpose of the present study is to analyze the diffusion behavior of hydrogen in iron by a molecular dynamics (MD) method and to clarify the stress dependence of the classical mechanical diffusion at temperature greater than room temperature, which is comparatively easy to handle. Furthermore, we construct a model that can explain the stress dependence of the activation enthalpy of diffusion by considering the stress-induced change in atomic configuration.

2. Analysis Methods and Analysis Conditions

2.1 Overview

To evaluate the effect of anisotropic stress on the hydrogen diffusion behavior in body-centered-cubic (bcc) iron, an MD method was used to simulate hydrogen atom diffusion. The time history of the coordinates of a hydrogen atom was obtained, and the mean square displacement (MSD) of a hydrogen atom was evaluated. In addition, the diffusion coefficient in the ⟨100⟩ direction under each investigated temperature was calculated from the Einstein relation, and the activation enthalpy and frequency factor of the diffusion were obtained from the corresponding Arrhenius plot. This analysis was performed under various stress conditions to evaluate the stress dependence of the activation enthalpy of diffusion and the frequency factor.

2.2 Analysis model and analysis conditions

The [100], [010], and [001] directions of the bcc iron lattice were taken along the x-, y-, and z-axes, respectively. The initial position of hydrogen atoms was tetrahedral interstitial sites (T sites) of the bcc iron lattice. The temperature was controlled by the velocity scaling method under 11 temperatures: 300, 330, 360, 400, 450, 500, 600, 700, 800, 900, and 1000 K. The stress conditions were isotropic stress σ_h (−3 GPa ≤ σ_h ≤ 3 GPa), shear stress τ_xy (0.2 GPa ≤ |τ_xy| ≤ 2 GPa), and uniaxial stress in the x-axis direction σ_x (−1.5 GPa ≤ σ_x ≤ 1.5 GPa). The positive signs of isotropic and uniaxial stresses indicate tension. In the uniaxial stress state, the cell size was controlled so that the stress in the two axes (y- and z-axes) other than the stress loading direction (x-axis) became 0. When isotropic stress and uniaxial stress were applied, the number of atoms was 2000 iron and 1 hydrogen, totaling 2001. The cell dimensions under no stress were 2.87, 2.87, and 2.87 nm along the x-, y-, and z-axes, respectively. When a shear stress was applied, to give a shear stress τ_xy acting on a plane perpendicular to the x- and y-axes, the $[1\bar{1}0]$, [110], and [001] directions of the bcc lattice were taken as the x′-, y′-, and z′-axes, respectively. In this model, the normal stresses in the x′ and y′ directions (σ_x′ = τ and σ_y′ = −τ) are given to realize a pure shear stress state (τ_xy = −τ). The total number of atoms was 1961, comprising 1960 iron atoms and 1 hydrogen atom. The cell dimensions under no stress were 2.84, 2.84, and 2.87 nm in the x′, y′, and z′ directions, respectively.

In all analyses in the present study, periodic boundary conditions were applied in all directions. For the time integration method, the multiple time integration method by reversible reference system propagator algorithms (rRESPA)¹¹⁾ was used; the time steps were Δt = 2 fs for iron atoms and Δt/24 for the hydrogen atom, and the total calculation time was 3 ns. For the interatomic potential, we used the EAM potential proposed by Wen et al.¹²⁾ Matsumoto et al. performed comparative calculations using this potential and first-principles calculations with respect to the heat of solution of hydrogen into the T site and O site of bcc iron. Their results showed that the heats of solution obtained by the two methods exhibited a similar strain dependence.¹³⁾ Therefore, this potential is suitable for the analysis of hydrogen diffusion behavior under stress.

2.3 Method for evaluating diffusion coefficient and activation enthalpy of diffusion

The diffusion coefficient of hydrogen atom D in the d dimension was obtained from the following Einstein relation on the basis of the time history of the atomic coordinates of hydrogen:   

\begin{equation} D = \lim_{\Delta\tau \to \infty}\frac{1}{2d\Delta \tau} \langle| \boldsymbol{{r}}(t_{0} + \Delta \tau) - \boldsymbol{{r}}(t_{0}) |^{2}\rangle \end{equation} (2)

where r is the position coordinate of the hydrogen atom and ⟨|r(t₀ + Δτ) − r(t₀)|²⟩ is the MSD during Δt (1 ps ≤ Δt ≤ 5 ps). From eq. (2), the diffusion coefficient can be obtained by dividing 2d into the slope of the line fitted to the plot of MSD as a function of the time width of the MSD, Δτ. The slope was determined using the least-squares method. When the MSD was calculated, the first 10 ps in the time history data of the MD calculation was excluded as the relaxation time and the MSD(Δτ) was derived from the average of the square displacements of various initial time t₀.

The activation enthalpy of diffusion E_a was obtained from the following Arrhenius equation:   

\begin{equation} D = D_{0}\exp \left(-\frac{E_{\text{a}}}{k_{\text{B}}T} \right) \end{equation} (3)

where D₀ is the frequency factor, k_B is the Boltzmann constant, and T is the temperature. According to eq. (3), the slope and the intercept of the Arrhenius plot, which displays the logarithm of D (ln D) plotted as a function of the reciprocal of the temperature (1/T), are −E_a/k_B and ln D₀, respectively. That is, the activation enthalpy E_a and frequency factor of diffusion D₀ can be obtained by linear approximation from the least-squares method to the Arrhenius plot of the diffusion coefficient obtained for each temperature condition. Furthermore, if D is replaced with one-dimensional diffusion coefficients D_x, D_y, and D_z, where d = 1 in eq. (2), the activation enthalpies of diffusion in the x, y, and z directions are obtained; these enthalpies are expressed as E_{a_x}, E_{a_y}, and E_{a_z}, respectively.

3. Analysis Results

3.1 Evaluation of diffusion coefficient under no stress

Arrhenius plots of the diffusion coefficients obtained under no-stress conditions are shown in Fig. 1. As discussed in Section 2.3, the activation enthalpy E_a and frequency factor D₀ were obtained from the slope and the intercept of the approximate straight line. Here E_a = 0.034 eV and D₀ = 1.59 × 10⁻⁸ m²/s were obtained.

Fig. 1

Arrhenius plot of diffusion coefficients against temperature: The activation enthalpy E_a and frequency factor D₀ are obtained from the slope and intercept of the fitted line, respectively.

3.2 Evaluation of activation enthalpy of diffusion under stress

The activation enthalpy of diffusion obtained under isotropic stress is shown in Fig. 2. The activation enthalpy was confirmed to be approximately the same value (0.032–0.036 eV) irrespective of the magnitude of the isotropic stress.

Fig. 2

Relationship between activation enthalpy of diffusion E_a and isotropic strain σ_h.

The activation enthalpy of diffusion obtained under shear stress is shown in Fig. 3. Here, activation enthalpies E_{a_x}, E_{a_y}, and E_{a_z} in the x-, y-, and z-axis directions introduced in Section 2.3 are shown. We confirmed that the activation enthalpy is approximately the same value (0.030–0.035 eV) irrespective of the magnitude of stress. In addition, the activation enthalpies in the x- and y-axis directions, where shear stress is applied, and in the z-axis direction, where shear stress is not applied, are also approximately the same. Thus, the activation enthalpy of hydrogen diffusion in bcc iron is not affected by either isotropic stress or shear stress along the {100} plane.

Fig. 3

Relationship between activation enthalpy of diffusion along the x-, y-, and z-axes (E_{a_x}, E_{a_y}, and E_{a_z}) and shear stress τ_xy.

The activation enthalpy of diffusion obtained under uniaxial stress is shown in Fig. 4. Because the activation enthalpies in the y and z directions perpendicular to the x-axis to which uniaxial stress is applied are equal as a consequence of symmetry, they are expressed as $E_{\text{a} \bot }$ and are determined as the average of E_{a_y} and E_{a_z}. The activation enthalpy in the direction parallel to the uniaxial stress $E_{\text{a}\mathrel{/\!/} }$ is defined using the activation enthalpy in the x-axis direction E_{a_x}. Figure 4 shows that, in the direction parallel to the uniaxial stress, the activation enthalpy $E_{\text{a}\mathrel{/\!/} }$ decreases with compression and increases with tension. By contrast, in the direction perpendicular to the uniaxial stress, the activation enthalpy $E_{\text{a} \bot }$ increases under both compression and tension. As previously described, strong anisotropy is observed in the stress dependence of activation enthalpy of diffusion under uniaxial stress along ⟨100⟩.

Fig. 4

Relationship between activation enthalpy of diffusion ($E_{\text{a}\mathrel{/\!/} }$ and $E_{\text{a} \bot }$) and uniaxial stress σ_x: $E_{\text{a}\mathrel{/\!/} }$ and $E_{\text{a} \bot }$ are the activation enthalpy of diffusion parallel to uniaxial stress and that perpendicular to uniaxial stress, respectively.

3.3 Frequency factors under stresses

The frequency factors obtained under isotropic stress and shear stress are shown in Figs. 5 and 6, respectively. In both cases, the value increases and decreases slightly in response to stress but is clearly almost constant.

Fig. 5

Relationship between frequency factor D₀ and isotropic stress σ_h.

Fig. 6

Relationship between the frequency factor along each axis (D_{0_x}, D_{0_y}, and D_{0_z}) and shear stress τ_xy.

The frequency factor obtained under uniaxial stress is shown in Fig. 7. Here, as in the case of the activation enthalpy, the frequency factor in the direction parallel to the uniaxial stress $D_{0\mathrel{/\!/} }$ and the frequency factor in the vertical direction $D_{0 \bot }$ are defined. Although the frequency factor $D_{0\mathrel{/\!/} }$ is almost constant, the frequency factor $D_{0 \bot }$ increases in both compression and tension.

Fig. 7

Relationship between frequency factor ($D_{0\mathrel{/\!/} }$ and $D_{0 \bot }$) and uniaxial stress σ_x.

Ramasubramaniam et al. performed first-principles calculations showing that the energy barriers for diffusion and the frequency factor are independent of isotropic stress and shear stress and reported the uniaxial strain dependence.¹⁴⁾ Our previously described results are consistent with this trend.

3.4 Diffusion path of hydrogen atom

The bcc lattice includes three crystallographically distinct octahedral interstitial sites (O sites) that differ by the direction of the four-fold rotationally symmetric axis of the octahedron. Hereinafter, the sites whose axes of symmetry are parallel to the x [100], y [010], and z [001] axes are referred to as O_x, O_y, and O_z sites, respectively.

Figure 8 shows the trajectory of hydrogen atoms for 11.5 ps at 0 GPa stress and at 300 K. Hydrogen is vibrating at the T site and diffusing along the O site to another T site. In addition, hydrogen diffusing in the x-axis direction passes near O_y and O_z sites, hydrogen diffusing in the y-axis direction passes near O_z and O_x sites, and hydrogen diffusing in the z-axis direction passes near O_x and O_y sites.

Fig. 8

Trajectory of hydrogen under the condition where the temperature is 300 K and the stress is 0 GPa: The hydrogen atom vibrates at one T site and jumps to a neighboring T site via the neighboring O site between them.

4. Discussion

As evident from the results described in Section 3.2, compared with the frequency factor, the activation enthalpy changes remarkably with ⟨100⟩ uniaxial stress. Therefore, in the following discussion, we focus on the activation enthalpy and discuss its stress dependence. As shown in Section 3.4, because hydrogen passes near the O site when it diffuses between two T sites, the activation enthalpy of diffusion is thought to depend approximately on the energy when hydrogen occupies the T site or the O site. Here, to consider the cause of stress dependence and anisotropy of the activation enthalpy shown in Section 3.2, we consider the stability of hydrogen at each site. The uniaxial strain-dependent anisotropy of the heat of solution is known to be negligible when hydrogen is inserted into the T site of bcc iron.¹³⁾ However, the heat of solution of hydrogen at the O site depends on the direction of uniaxial strain.¹³⁾ Thus, the anisotropy of the stress dependence of the activation enthalpy can be qualitatively explained from the effect of strain at the O site.

Compared with a regular octahedron, an octahedron surrounding the O site has a collapsed shape along one axis of symmetry. When the strain in comparison with that of a regular octahedron is expressed as the degree of anisotropy (DOA), a smaller DOA indicates lower energy at the O-site occupancy, and the energy difference between the T and O sites (i.e., the activation enthalpy) becomes smaller. Conversely, the greater the DOA, the greater the activation enthalpy. For example, when a uniaxial tensile strain of 0.5% is applied to reduce the DOA of the octahedron around the O site, the energy difference is 0.027 eV, which is smaller than the energy difference of 0.039 eV in the unstrained state. However, when a uniaxial compressive strain of 0.5% was applied to increase the DOA, the energy difference was 0.053 eV, which is larger than that in the unstrained state.

On the basis of the previous results, we consider the activation enthalpy under uniaxial stress. First, when compressive stress is applied in the x-axis direction, the DOA of the O_x site increases and the DOA of the O_y and O_z sites decreases. As shown in Fig. 8, the hydrogen diffusion path in the x-axis direction (parallel to the stress) passes near the O_y and O_z sites, both of which have a small DOA; thus, the activation enthalpy decreases. This result matches that in Fig. 4. In the diffusion path along the y- [z-] axis direction (normal to the stress), the activation enthalpy is rate-limiting toward the larger path because it passes near the O_x site with the larger DOA and the O_z [O_y] site with a smaller DOA. However, when tensile stress is applied in the x-axis direction, the DOA of the O_x site becomes small and the DOA of the O_y and O_z sites becomes large. The diffusion path in the x-axis direction (parallel to the stress) passes near O_y and O_z sites, both of which have a large DOA, and the diffusion path in the y- [z-] axis direction (perpendicular to the stress) passes near O_x sites with a small DOA and O_z [O_y] sites with a large DOA. That is, the activation enthalpy increases in all axis directions, which matches the results shown in Fig. 4.

5. Modeling of Stress Dependence of Activation Enthalpy

5.1 Formula for predicting activation enthalpy

For the activation enthalpy of diffusion under stress, the activation energy is corrected by (ε^*_xxσ_xx + ε^*_yyσ_yy + ε^*_zzσ_zz) × V using the axial strain of the calculation cell at saddle-point configuration (ε^*_xx, ε^*_yy, ε^*_zz), the applied stress (σ_xx, σ_yy, σ_zz), and the calculation cell volume V. In this method, the activation enthalpy cannot be expressed independently of the isotropic stress. Therefore, a formula is derived to predict the activation enthalpy of diffusion on the basis of the DOA of the octahedron surrounding the O site, as introduced in Section 4. First, we consider the index of the DOA quantitatively. In the case of the O_x site, when the longitudinal strain due to the stress in the x-axis direction and the lateral strain due to the stress in the y- and z-axis directions are considered, the following equation α_x is expressed as the DOA of the O_x site (Fig. 9):   

\begin{equation} \alpha_{x} = \sigma_{x} - \frac{1}{2}(\sigma_{y} + \sigma_{z}) \end{equation} (4)

Similarly, the degrees of anisotropy of the O_y and O_z sites are represented by the expressions   

\begin{equation} \alpha_{y} = \sigma_{y} - \frac{1}{2}(\sigma_{z} + \sigma_{x}) \end{equation} (5)

  

\begin{equation} \alpha_{z} = \sigma_{z} - \frac{1}{2}(\sigma_{x} + \sigma_{y}) \end{equation} (6)

Because the diffusion path in the x-axis direction passes near O_y and O_z sites, the activation enthalpy becomes a function of their DOA α_y and α_z, as shown in the following equation:   

\begin{equation} E_{\text{a}{\_}x} = F(\alpha_{y},\alpha_{z}) \end{equation} (7)

Because diffusion is the same along all three axes, the same function F is used to express E_{a_y} and E_{a_z}:   

\begin{equation} E_{\text{a}{\_}y} = F(\alpha_{z},\alpha_{x}) \end{equation} (8)

  

\begin{equation} E_{\text{a}{\_}z} = F(\alpha_{x},\alpha_{y}) \end{equation} (9)

Because α_x, α_y, and α_z always become 0 in the case of isotropic stress and because the term for shear stress is not included, the activation enthalpy does not depend on isotropic stress and shear stress.

Fig. 9

Stress-induced deformation of the atomic structure around the O site.

When the symmetry of hydrogen diffusion is considered, F(u, v) is symmetric about u and v; thus, F can be expressed by a polynomial in two variables:   

\begin{align} F(u,v)&=a_{00} + a_{01}(u + v) + a_{02}(u^{2} + v^{2}) + a_{11}uv \\ &\quad + a_{12}(u^{2}v + uv^{2}) + a_{22}u^{2}v^{2} \end{align} (10)

The activation enthalpy obtained by MD calculation was plotted against the DOA α and its coefficient was determined by fitting eq. (10) (Fig. 10). When the fitted surface was obtained, in addition to the activation enthalpy for uniaxial stress σ_x (−1.5 GPa ≤ σ_x ≤ 1.5 GPa) shown in Section 3.2, those for biaxial stresses in the x and y directions, σ_x and σ_y (0 GPa ≤ σ_x, σ_y ≤ 2 GPa), were used. The analytical model and analytical conditions for biaxial stress conditions are similar to those for uniaxial stress conditions (Section 2.2). The cell size was controlled such that the stress in the direction where no stress was applied became 0. The obtained coefficients are shown in Table 1.

Fig. 10

Fitted results of F(α₁, α₂) to the MD simulations.

Table 1 Fitted constants in F(α₁, α₂).

5.2 Validating prediction formulas

Equations (7)–(9) for predicting the activation enthalpy were verified under multiaxial stress conditions. To apply shear stress τ_xy in the directions of x [100], y [010], and z [001], the $[2\bar{1}0]$, [120], and [001] directions were set to x′′, y′′, and z′′ axes, respectively (Fig. 11), and normal stresses σ_x′′, σ_y′′, and σ_z′′ were added. Thus, the normal stresses σ_x, σ_y, and σ_z and shear stress τ_xy in the xyz coordinate system become   

\begin{equation} \left. \begin{array}{c} \sigma_{x}\\ \sigma_{y} \end{array} \right\} = \frac{1}{2}(\sigma_{x''} + \sigma_{y''}) \pm \frac{3}{5}\left(\frac{\sigma_{x''} - \sigma_{y''}}{2} \right) \end{equation} (11)

  

\begin{equation} \sigma_{z} = \sigma_{z''} \end{equation} (12)

  

\begin{equation} \tau_{xy} = - \frac{4}{5}\left(\frac{\sigma_{x''} - \sigma_{y''}}{2} \right) \end{equation} (13)

The analytical model comprises 2500 iron and 1 hydrogen atoms, for a total of 2501 atoms. The cell dimensions under no stress are 3.21 nm, 3.21 nm, and 2.87 nm on the x′′, y′′ and z′′ axes, respectively.

Fig. 11

Simulation model used to verify the accuracy of the prediction model.

Table 2 shows the verified stress conditions, the activation enthalpy obtained by MD calculation, the activation enthalpy obtained by prediction formula (10), and the error of the prediction formula relative to the MD calculation. They also show the DOA under each stress condition. Under conditions (a) and (b), where the DOA is small, the activation enthalpy is determined with high accuracy by the prediction formula; however, under conditions (e), (f), and (g), where the DOA is far from 0, such that the equivalent stress σ_eq exceeds 1.5 GPa, a substantial difference is observed between the activation enthalpy obtained by MD calculation and that obtained by prediction formula (10). Because the stress condition of the MD calculation used for fitting the prediction formula is not satisfied in this case, the error of the prediction formula is increased. Such stress conditions occur in very limited cases, such as in the immediate vicinity of dislocations. The accuracy is improved by adopting the sample data of the fitting more widely when greater accuracy is necessary for such stress conditions.

Table 2 Comparison of the activation enthalpy obtained by MD simulations and by the prediction model.

6. Conclusion

The following are the findings from MD simulations of hydrogen diffusion in bcc iron by determining the diffusion coefficients in the [100], [010], and [001] directions and evaluating the stress dependence of the diffusion coefficients.

• •    The value of the diffusion coefficient under isotropic stress and {100} plane shear stress is almost constant, irrespective of the stress.
• •    Diffusion coefficients under ⟨100⟩ uniaxial stress show different stress dependence in the directions parallel and perpendicular to the stress.

In addition, regarding the stress dependence of the obtained diffusion coefficient, we examined the cause of the stress dependence of the activation enthalpy of diffusion by analyzing the effect of stress on the deformation of the octahedron composed of neighboring atoms around the O site in bcc iron. We found that anisotropy occurs in the stress dependence when the diffusion path passes near O sites with different deformation depending on the direction of diffusion. Furthermore, by quantitatively expressing the deformation of the octahedron around the O site, we developed a formula to predict the activation enthalpy and subsequently verified its accuracy. When the deformation of the O site is not large, the prediction of our formula agrees with the result of MD calculations with high accuracy, and its effectiveness can be shown.

REFERENCES

• 1)   Y.  Murakami: Int. J. Fract. 138 (2006) 167–195.
• 2)  M. Nagumo: Fundamentals of Hydrogen Embrittlement, (Springer, Berlin, 2016) pp. 103–133.
• 3)   G.  Han,  J.  He,  S.  Fukuyama and  K.  Yokogawa: Acta Mater. 46 (1998) 4559–4570.
• 4)   P.  Sofronis and  R.M.  McMeeking: J. Mechan. Phys. Solids 37 (1989) 317–350.
• 5)   H.  Kotake,  R.  Matsumoto,  S.  Taketomi and  N.  Miyazaki: Int. J. Press. Vessels Piping 85 (2008) 540–549.
• 6)   T.  Takayama,  R.  Matsumoto,  S.  Taketomi and  N.  Miyazaki: Int. J. Hydrogen Energ. 36 (2011) 1037–1045.
• 7)   S.  Taketomi,  R.  Matsumoto and  N.  Miyazaki: Acta Mater. 56 (2008) 3761–3769.
• 8)   N.  Ishikawa,  H.  Sueyoshi and  A.  Nagao: ISIJ Int. 56 (2016) 413–417.
• 9)   M.  Koyama,  D.  Yamasaki,  T.  Nagashima,  C.C.  Tasan and  K.  Tsuzaki: Scr. Mater. 129 (2017) 48–51.
• 10)   H.  Kimizuka,  H.  Mori and  S.  Ogata: Phys. Rev. B 83 (2011) 094110.
• 11)   M.  Tuckerman,  B.J.  Berne and  G.J.  Martyna: J. Chem. Phys. 97 (1992) 1990–2001.
• 12)   M.  Wen,  X.  Xu,  S.  Fukuyama and  K.  Yokogawa: J. Mater. Res. 16 (2001) 3496–3502.
• 13)   R.  Matsumoto,  Y.  Inoue,  S.  Taketomi and  N.  Miyazaki: Scr. Mater. 60 (2009) 555–558.
• 14)   A.  Ramasubramaniam,  M.  Itakura,  M.  Ortiz and  E.A.  Carter: J. Mater. Res. 23 (2008) 2757–2773.