Introduction
In spite of the successful development of new alloys with outstanding strength and toughness, hydrogen (H) embrittlement stays one of the most severe and most controversial types of failure that affect almost all types of metals.[Reference Barnoush and Vehoff1] This can be accredited to the fact that H atoms can easily infiltrate the metal lattice. Hydrogen embrittlement typically leads to H-induced fracture even at low stress intensities.[Reference McMahon2] One of the main hypotheses proposed to explain H embrittlement is the H-enhanced decohesion. In this hypothesis, the increased solubility of H results in a decrease in the atom cohesive forces at the grain boundary (GB). Thus, at adequately high H concentrations, the cohesive strength of the material decreases leading to an early brittle intergranular separation due to this diminishing in the cohesive forces along the GBs.[Reference Song and Curtin3]
Coherent twin boundaries (CTBs) have been observed to show superior corrosion resistance in plethora of engineering alloys.[Reference Seita, Hanson, Gradečak and Demkowicz4] However, in practice, GBs are rarely “ideal” and their structures do not always coincide with the proper coincidence site lattice (CSL). According to the Brandon criterion, the CSL-GB structure can show small angular deviation from its ideal symmetry plane and the proximity of this deviation to a CSL structure can determine the GB traits. Additionally, this deviation is accommodated by misfit GB dislocations. Thus, such a deviation could significantly impact the segregation and diffusivity behavior of solute atoms along the GB. Indeed, it has been shown through experimental observations that many GBs in engineering materials diverge from the ideal symmetry plane.[Reference Herbig, Raabe, Li, Choi, Zaefferer and Goto5, Reference Wright and Larsen6] Herbig et al.[Reference Herbig, Raabe, Li, Choi, Zaefferer and Goto5] has also shown through transmission electron microscopy and atom probe tomography study that small deviations from CSL-GB orientations in ferritic steels cause significant high peaks of segregation at those GBs.
Impurity atoms and solutes tend to segregate toward GBs and other defects due to their atomic mismatch in order to reduce the free energy of the system.[Reference Herbig, Raabe, Li, Choi, Zaefferer and Goto5] Segregation energy is an important factor, particularly in determining the vulnerability of a specific GB to H embrittlement. Molecular statics (MS) and ab initio simulations have been used broadly to calculate H segregation energies in the lattice or at different defects in metals.[Reference Song and Curtin7–Reference Rajagopalan, Tschopp and Solanki9] On the other hand, molecular dynamics (MD) studies were also conducted to predict the H-diffusivity in the lattice of α-iron as well as on special symmetrical tilt GBs.[Reference Kimizuka, Mori and Ogata10, Reference Liu, Xie, Chen and Zhang11]
Nevertheless, a comprehensive study of the impact of small angular deviation from the ideal CSL configuration on H diffusion and segregation is still missing. Consequently, the effect of the deviation angle from a perfect
$\Sigma 3(111)[1\bar 10]$ GB in α-iron on H segregation energy and diffusion along the GB is analyzed. First MS simulations are utilized to assess H segregation energies at the GB and free surface. Moreover, MD simulations are executed to predict the effect of the deviation angle on H diffusivity along the GB. The results are then analyzed with regard to the impact of angular deviation on the tendency for H embrittlement. The research article is organized as follows, the simulation methodology are discussed in details in “Methodology” section. The simulation results are presented in “Results” section and discussed in “Discussion” section. Finally, conclusions and summary are given in “Conclusions” section.
Methodology
All simulations were conducted using the large-scale atomic/molecular massively parallel simulator,[Reference Plimpton12] with the embedded atom method potential for H-αFe developed by Ramasubramaniam et al.[Reference Ramasubramaniam, Itakura and Carter13] and subsequently modified to account for the repulsive H–H interactions by Song and Curtin.[Reference Song and Curtin3] The simulation cells utilized here were rectangular with a simulation cell size of
$20\,\hbox{nm} \times 20\,\hbox{nm} \times 20\,{\kern 1pt} \hbox{nm}$ unless otherwise stated.
For bulk diffusivity simulations, periodic boundary conditions (PBCs) were employed along the x, y, and z coordinates, which are parallel to the [100], [010], and [001] crystallographic directions, respectively. Hydrogen atoms were introduced randomly at tetrahedral sites, which are the most energetically favorable sites for H occupation in BCC iron, with a concentration of
${\rm \sim} 1.0\% $ of all possible occupation sites in the simulation cell.
For free surface diffusivity simulations, free surface boundary conditions were imposed along the z direction, while PBCs are imposed along the other two directions. The H atoms were randomly introduced at tetrahedral sites within a 1.0 nm layer from the top and bottom free surfaces of the simulation cell with a concentration of
${\rm \sim} 1.0\% $ based on all occupation sites available in those regions.
For GB diffusivity simulations, a symmetrical tilt
$\Sigma 3(111)[1\bar 10]$ GB with misorientation angle θ = 109.5° was introduced at the middle of the height of the simulation cell. PBCs were employed along all three directions. In the GB plane, this periodicity is attained by maintaining the corresponding CSL along the crystallographic orientations. In addition, atom deletion with multiple configurations was used to remove overlapping Fe atoms at the boundary to achieve an optimum GB energy and avoid the formation of residual stresses inside the simulation cell.[Reference Tschopp and McDowell14] The H-atoms were then added at the tetrahedral sites within a 0.5 nm-thick layer surrounding the GB plane along its plane normal.
After the simulation cell was initialized, energy minimization was performed at 0 K using the Polak–Ribiere conjugate gradient algorithm[Reference Yuan, Wei and Li15] with energy tolerance of 10−12 and force tolerance of 10−13 eV/Å. The system temperature was then increased to the desired temperature in the range of 50–700 K and then the system was relaxed for a period of 2 ns. The time step was set to be 0.001 ps for all the segregation and diffusion simulations. The isobaric–isothermal NPT ensemble (constant number of atoms, pressure, and temperature)[Reference Mendelev, Han, Srolovitz, Ackland, Sun and Asta16] was used for the bulk and GB simulations, while the canonical NVT ensemble (constant number of atoms, volume, and temperature) was used for the free surface simulations.
As it was mentioned previously, by applying Brandon criterion for conserving CSL structures along the special GB planes, the allowable deviation angle threshold is[Reference Brandon17]:
where the constant ϕ 0 ≈ 15° and Σ is the reciprocal value of the CSL density. Equation (1) gives ϕ limit ≈ 8.66° for the Σ3 boundary. Hence, to measure the effect of a small deviation from the symmetrical tilt plane of an ideal
$\Sigma 3(111)[1\bar 10]$ GB on H-energetics and dynamics, the GB plane is deviated from the (111) ideal symmetry plane by changing the misorientation angle θ between the two grains by a small deviation angle 0° ≤ ϕ ≤ 9.0° measured from the (111) plane, as shown in Fig. 1. In these simulations, free surface boundary conditions are employed in all three directions since the periodicity of the simulation cell can no longer be maintained. Furthermore, the same atom deletion and energy minimization technique discussed above are also utilized for these simulation cells.
The average GB and free surface H-segregation energy within a 2Å layer encompassing the aforementioned defects as a function of the deviation angle from the
$\Sigma 3\,(111)[1\bar 10]$ ideal symmetry plane are presented. The dash-dot parabolic fit and the solid spline interpolation fit represents the H behavior at GB and free surface, respectively (statistical analysis can be found in the Supplementary document). A schematic diagram of the simulation cell for the
$\Sigma 3(111)[1\bar 10]$ GB structures is embedded as well; where the misorientation angle, θ, between the
$(1\bar 1\bar 1)$ plane normal and ideal (111) reference plane normal as well as the deviation angle, ϕ, are shown.

In these simulations, the hydrogen GB segregation energy was computed by varying the H-atom position in the range of 0–15 Å from the GB plane and at a minimum distance of 85 Å from the simulation cell side surfaces to minimize the effect of free surface boundary conditions. Each simulation was repeated at least five times with different atom positions.
In addition, the H-diffusivity is also calculated at 300 K as a function of ϕ by randomly introducing H-atoms in the simulation cell at a minimum distance of 10 nm from the simulation cell side surfaces such that the effects of the free surface boundary conditions are minimized. In addition, the H-atoms are initially placed within a 0.5 nm-thick layer around the GB with a concentration equal to 1% of the available tetrahedral sites in that layer. Each simulation was performed with different random H-atom initialization at least four times. The same energy minimization techniques discussed above were also used here.
GB energy and hydrogen segregation energy calculations
The GB energy, γ GB, is computed as follows[Reference Bhattacharya, Tanaka, Shiihara and Kohyama18]:
where E GB is the simulation cell total energy having both a GB and free surfaces and is calculated from the summation of the potential and kinetic energies of all atoms in the simulation cell, E SC is the single-crystal simulation cell total energy having the same dimensions with free surfaces only, and A GB is the GB area.
In addition, the hydrogen GB-segregation energy,
$E_{{\rm seg,GB}}^\alpha $, is computed as follows[Reference Solanki, Tschopp, Bhatia and Rhodes8, Reference Rajagopalan, Tschopp and Solanki9]:
where
$E_{{\rm GB}}^\alpha $ is the simulation cell total energy with a GB, free surfaces, and H-atoms, while
$E_{{\rm SC}}^\alpha $ is the corresponding single-crystal simulation cell total energy with free surfaces and H-atoms.
Similarly, the hydrogen-free surface segregation energy,
$E_{{\rm seg,FS}}^\alpha $, is calculated as follows[Reference Hamza, Hatem, Raabe and El-Awady19]:
where
$E_{{\rm FS}}^\alpha $ is the single-crystal simulation cell total energy with free surface boundary conditions and charged with H-atoms, E FS is the corresponding pure cell total energy with free surface boundary conditions,
$E_{{\rm SCP}}^\alpha $ is the corresponding hydrogen charged cell total energy with PBCs and finally E SCP is the corresponding pure cell total energy with PBCs.
Hydrogen diffusivity calculations
The mean square displacement (MSD) of all H-atoms during an MD simulation can be computed by[Reference Liu, Xie, Chen and Zhang11]:
$$\hbox{MSD} = \displaystyle{1 \over {2dN}}\sum\limits_{i = 1}^N (r_i(t_0 + \Delta t) - r_i(t_0))^2,$$where d is a dimensionality parameter that is equal to 1.0 for diffusivity calculations along a certain direction, 2.0 for free surface and GB, and 3.0 for bulk, N is the total number of H-atoms, r i(t 0) is the initial position of a H-atom i at time t 0 and r i(t 0 + Δt) is its position after a time step of Δt. The H-diffusivity is thus computed as the slope of MSD versus time. The bulk H-diffusivity is computed as the average between the diffusivity along the x, y, and z directions. For free surface and GB simulations, the H-diffusivity is computed as the average of the diffusion coefficients within the free surface (i.e., x–y plane) and the GB (i.e., x–z plane), respectively.
Results
Effect of small deviations from the Σ3 (111) ideal symmetry plane on H-segregation energy
From the current simulations, the GB energy of the ideal
$\Sigma 3(111)[1\bar 10]$ GB structures is 1.3 J/m2, which is in good agreement with previous published estimates.[Reference Solanki, Tschopp, Bhatia and Rhodes8] In addition, the average H-segregation energy as a function of position from the GB as computed from Eq. (3) is shown in Fig. 1 (Supplementary Fig. S1) for ϕ = 0°, 3°, 5°, and 9°. For all the simulated cases, the segregation energy at the GB is the highest, which indicates that these GBs are all more preferential for hydrogen segregation as compared with the bulk lattice. Moreover, as the angular deviations increase, the GB characteristic H-absorption length increases. This length indicates the influence zone of the GB on H-atoms located at tetrahedral sites away from the GB plane. The H-absorption length for the ideal
$\Sigma 3(111)[1\bar 10]$ GB is
${\rm \sim} 3{\kern 1pt} \,{{\AA}}$, and reaches
${\rm \sim} 10{\kern 1pt} \,{{\AA}}$ for ϕ = 9°. This increase in the characteristic H-absorption length is mainly attributed to the increase in atomic disorder in the GB with increasing deviation angle.
The impact of the deviation angle on the H-segregation energy within a 2 Å layer encompassing the GB is summarized in Fig. 1 (Supplementary Fig. S2). The H-segregation energy is the highest in the ideal GB with a mean value of
${\rm \sim} \hskip-2.5pt- \hskip-2.5pt 0.62{\kern 1pt} \,\hbox{eV}$, decreases to
${\rm \sim} \hskip-2.5pt- \hskip-2.5pt 0.45{\kern 1pt} \,\hbox{eV}$ and remains relatively constant for 3° ≤ ϕ ≤ 7°, then increases to
${\rm \sim} \hskip-2.5pt- \hskip-2.5pt 0.6{\kern 1pt} \,\hbox{eV}$ for ϕ ≤ 9°. Note all the segregation energies aforementioned are mean values calculated using statistical analysis with a 95% confidence interval.
As shown in Fig. 1, the free surface in the
$[\bar 1\bar 11]$ grain, termed hereafter “grain 1”, always has the same orientation in all simulations. However, the free surface in the other grain, termed hereafter “grain 2”, will change for the different cases with a deviation from the ideal GB. Thus, a different H-segregation energy would be anticipated for that later free surface. The average H-segregation energy, as computed from Eq. (4), for each free surface orientation is summarized in Fig. 1 (Supplementary Fig. S3). There are no considerable differences between the free surface H-segregation energy in “grain 1” and “grain 2” for ϕ = 0°, 3°, and 9°, in which the mean value with a 95% confidence level is
${\rm \sim} \hskip-2.5pt- \hskip-2.5pt 0.71{\kern 1pt} \,\hbox{eV}$. On the other hand, the H-segregation energy mean value for the free surface in “grain 2” with a deviation angle of 5° is
${\rm \sim} \hskip-2.5pt- \hskip-2.5pt 0.44{\kern 1pt} \,\hbox{eV}$, which is the lowest value, in absolute terms, in all the simulated cases. It is also interesting to note that this energy is smaller, in absolute terms, than those for the predicted GB H-segregation energies.
Hydrogen diffusivity
The predicted H lattice diffusion coefficient, d, from the current simulation in a bulk single-crystal α-Fe as a function of temperature is shown in Fig. 2a compared with previously published experimental results,[Reference Beck, Bockris, McBreen and Nanis20–Reference Nagano, Hayashi, Ohtani, Isshiki and Igaki22] MD simulations,[Reference Kimizuka, Mori and Ogata10, Reference Liu, Xie, Chen and Zhang11, Reference Zhu and Oda23] centroid molecular dynamics (CMD) simulations,[Reference Kimizuka, Mori and Ogata10] and path integral quantum transition state theory (PI-QTST) calculations.[Reference Katzarov, Pashov and Paxton24] In general, the currently predicted H diffusion coefficient as well as previously published MD simulation predictions fall in the range of experimental scatter in the temperature range of 250–1200 K. This large scatter in the experimental predictions is typically attributed to trapping of H at different defects (dislocation, GBs, precipitates, etc.) among other effects.[Reference Kimizuka, Mori and Ogata10, Reference Kiuchi and McLellan25] Thus, the comparisons between the MD and experimental predictions should be made on a qualitative basis rather than a quantitative one.
(a) The lattice H diffusion coefficient in α-Fe as predicted from the current MD simulations. Other predictions from published experiments,[Reference McMahon2, Reference Yuan, Wei and Li15, Reference Brandon17] MD simulations,[Reference Rajagopalan, Tschopp and Solanki9, Reference Plimpton12, Reference Krom and Bakker30] CMD simulations,[Reference Rajagopalan, Tschopp and Solanki9] and PI-QTST calculations.[Reference Solanki, Tschopp, Bhatia and Rhodes8] The shaded region represents a large rage of experimental predictions as summarized by Kiuchi and McLellan.[Reference Kimizuka, Mori and Ogata10] (b) A semi-log graph of a subset of data in (a) showing the non-Arrhenius relationship for the lattice H diffusion coefficient as a function of reciprocal temperature. The solid line represents the best curve fit.

In addition, as shown in Fig. 2b, the relationship between the diffusion coefficient and the reciprocal temperature as predicted from the current MD simulations follows a nearly non-Arrhenius relationship. There is a distinct deviation from a linear relationship between the logarithmic diffusion coefficient and the reciprocal temperatures for temperatures below 150 K. This is qualitatively in agreement with CMD,[Reference Kimizuka, Mori and Ogata10] PI-QTST,[Reference Katzarov, Pashov and Paxton24] and density function theory[Reference Di Stefano, Mrovec and Elsässer26] predictions. This non-Arrhenius response is attributed to that H diffusion is dominated by the classical jump over the barrier mechanisms above a critical temperature, while below it quantum tunneling becomes dominant.[Reference Kimizuka, Mori and Ogata10, Reference Katzarov, Pashov and Paxton24, Reference Di Stefano, Mrovec and Elsässer26] Nevertheless, it should be noted that quantum tunneling effects could not be captured accurately in classical MD simulations. This explains the large deviation between the predicted diffusion coefficient below room temperature from the current MD simulations and those predicted from other methods that account for quantum tunneling. However, for T ≥ 300, the current MD simulations are in excellent agreement with the PI-QTST and CMD predictions.
The predicted H diffusion coefficient as predicted from the current simulations for lattice, (001) free surfaces and ideal
$\Sigma 3(111)[1\bar 10]$ GB as a function of reciprocal temperature are shown on a semi-log scale in Fig. 3. It is clear that in all three cases, the diffusivity follows a non-Arrhenius-like behavior. Nevertheless, the results shown in Fig. 3 below room temperature should only be viewed in a qualitative manner since quantum tunneling is not accounted for in the current simulations, which if accounted for would most likely lead to an order of magnitude higher diffusivity values at temperatures below room temperature.
A semi-log graph showing the predicted H diffusion coefficient as predicted from the current MD simulations for lattice, (001) free surfaces, and ideal
$\Sigma 3(111)[1\bar 10]$ GB as a function of reciprocal temperature. The solid line represents the best power-law curve fit.

The effect of the deviation from the Σ3 (111) ideal symmetry plane on the GB H-diffusion coefficient at 300° is shown in Fig. 4. The ideal GB structure is shown to result in the highest H-diffusion coefficient D ≈ 3.5 × 10−7 cm2/s. However, the extra atomic disorder induced by the increase in the deviation angle leads to a continuous decrease in the H-diffusivity with increasing deviation angle. In particular, the H-diffusivity for the 9° case shows a decrease by an order of magnitude as compared with the ideal GB case. The only exception for the continuous decrease in H-diffusivity is a sudden increase for the 7° case. This is most likely a special case in which the GB structure enhances H-diffusion as compared with the other GBs. Further quantifications of the GB structure is however needed to quantify this effect further.
A semi-log graph showing the predicted H diffusion coefficient at room temperature as predicted from the current MD simulations as a function of deviation angle from the Σ3 (111) ideal symmetry plane. The solid line represents the best curve fit.

Discussion
GB embrittlement
Generally, solute atoms segregation at a GB can act as either GB-cohesive enhancer or GB-cohesive reducer depending on the resulting solute atom–GB interaction. As an example, Rajagopalan et al.[Reference Rajagopalan, Tschopp and Solanki9] showed that phosphorous segregation at Σ5(210)θ = 53.13° GBs in α-Fe reduces the GB-cohesion strength leading to a weakened GB structure, while vanadium acts as a GB-cohesive enhancer and strengthens this same GB. On the other hand, H typically leads to an increase in the likelihood of GB separation and leads to the formation of micro cracks at the GB, and consequently intergranular fracture in almost all metals.
In order to quantify the effect of the deviation from the Σ3 (111) ideal symmetry plane on the susceptibility to H-embrittlement, the GB-cohesive energy can be computed by the Rice–Wang model[Reference Rice and Wang27, Reference Zhong, Wu, Freeman and Olson28]:
where 2γ int is the GB-cohesion energy with the existence of H atom, (2γ int)0 is the cohesion energy of the pure GB, and τ is the GB and free surface coverage constant. Based on Eq. (6), as the difference between the GB and free surface segregation energies increases, the GB-cohesion energy will be less than that of the pure GB, leading to an enhanced embrittlement effect.[Reference Yamaguchi, Shiga and Kaburaki29]
As shown in Fig. 1, the hydrogen-free surface segregation energy is always more negative than that of the GB, and hence the GB cohesion energy as predicted from Eq. (6) will be reduced for the ideal and deviated GB cases. This implies that H is more energetically stable at the free surfaces than GB, leading to decrease in the Griffith work for brittle fracture, hence intergranular fracture is observed. Moreover, the embrittlement potency can be viewed as a linear function of the difference between the free surfaces and GB segregation energies. Accordingly, the ideal GB seems to be more resistant to H embrittlement as compared with the deviated cases.[Reference Yamaguchi, Shiga and Kaburaki29] Compared with all the simulated cases here, the deviated GB with ϕ = 3° seems to show the greatest susceptibility to decohesion and subsequently to boundary separation as compared with all other cases considered here.
In the case of GB saturation, where all trapping sites in the GB are occupied, the GB segregation energy will be effectively zero.[Reference Solanki, Tschopp, Bhatia and Rhodes8] Thus, the free surface segregation energy will be the key factor in determining the GB potency to embrittlement. As shown in Fig. 1 (Supplementary Fig. S1), the 5° deviation case has the lowest free surface H-segregation energy, hence this case is expected to show the highest resistance to embrittlement, followed by the 7° deviation case. However, for the 0°, 3°, and 9° cases, there are no noticeable differences between their free surface segregation energies, thus they will all show similar embrittlement effects, which is relatively higher than the other GB cases simulated here.
Effect of deviation angle on hydrogen diffusion in the GB
The deviation from the ideal GB symmetry plane will significantly affect the diffusivity of H within the GB structure. Angular deviations can be classified as either deviation from the ideal symmetry plane or deviation from the CSL misorientation angle. Brandon[Reference Brandon17] indicated that with a small deviation from the ideal symmetry plane, the GB will be accommodated with a forest of screw dislocations superimposed on the coincidence boundary. This occurs when the sub-boundary axis of rotation is perpendicular to the CSL boundary ideal symmetry plane. On the other hand, for the case of deviation from the misorientation angle, Brandon's model suggests a superimposed boundary of mixed dislocations formed within the CSL boundary with an axis of rotation being parallel to the ideal symmetry plane.[Reference Brandon17] In practice, the axis of rotation of this sub-boundary structure lies at an arbitrary angle from the ideal reference plane, and thus both sub-boundary suggested models will occur simultaneously. It can be deduced from Fig. 1 that in our current simulations the angular deviation would lead to the superposition of both types of sub-boundary models.
It can thus be inferred that with the increase of dislocation densities within the GB, more trapping sites will be available,[Reference Krom and Bakker30] which subsequently impede hydrogen motion. This can also be correlated with the diffusivity trend noticed in Fig. 4. As with extra atomic mismatch and incoherence, the diffusivity along those boundaries decreases and the boundary looses its hydrogen transport ability.
Conclusions
Hydrogen segregation and diffusion along the
$\Sigma 3(111)[1\bar 10]$ α-iron GB structure were comprehensively analyzed using MS and dynamics simulations. The results strongly recommend that angular deviations from the proper CSL significantly affect the hydrogen energetics and kinetics within the GB. The simulation results indicate that the ideal GB structure acts as a preferable hydrogen sink and has the highest binding energy to the segregated atoms. However, with increasing deviation angle from the ideal symmetry plane within the limits defined by the Brandon's criterion for the classification of GBs, it is shown that the hydrogen segregation energy at the GB decreases, in absolute terms, with increase in deviation angle up to 5°, and subsequently decreases. The hydrogen diffusivity along the ideal GB is observed to be the highest and decreases continuously with increasing deviation angle due to the extra atomic disorder in the deviated GBs. Finally, by utilizing the Rice–Wang model, it is shown that the 3° deviated GB has the highest susceptibility to hydrogen embrittlement, while the ideal GB shows the highest resistance to decohesion below the hydrogen saturation limit. These results suggest the importance of considering all crystallographic aspects of GBs in real materials when quantifying segregation and embrittlement. The results can also help explain the experimentally observed susceptibility of CTBs to crack initiation.[Reference Seita, Hanson, Gradečak and Demkowicz4]
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1557/mrc.2018.186.
Acknowledgments
This work was supported by the National Science Foundation CAREER Award #CMMI-1454072 and Academy of Scientific Research and Technology JESOR grant #17.

