Molecular Dynamics Analysis of Hydrogen Diffusion Behavior in Alpha-Fe Bi-Crystal Under Bending Deformation

Abstract

The hydrogen embrittlement (HE) phenomenon occurring in drawn pearlitic steel wires sometimes results in dangerous delayed fracture and has been an important issue for a long time. HE is very sensitive to the amount of plastic deformation applied in the drawing process. Hydrogen (H) atom diffusion is affected by ambient thermal and mechanical conditions such as stress, pressure, and temperature. In addition, the influence of stress gradient (SG) on atomic diffusion is supposed to be crucial but is still unclear. Metallic materials undergoing plastic deformation naturally have SG, such as residual stresses, especially in inhomogeneous regions (e.g., surface or grain boundary). In this study, we performed molecular dynamics (MD) simulation using EAM potentials for Fe and H atoms and investigated the behavior of H atoms diffusing in pure iron (α-Fe) with the SG condition. Two types of SG conditions were investigated: an overall gradient established by a bending deformation of the specimen and an atomic-scale local gradient caused by the grain boundary (GB) structure. A bi-crystal model with H atoms and a GB structure was subjected to bending deformation. For a moderate flexure, bending stress is distributed linearly along the thickness of the specimen. The diffusion coefficient of H atoms in the bulk region increased with an increase in the SG value. In addition, it was clearly observed that the direction of diffusion was affected by the existence of the SG. It was found that diffusivity of the H atom is promoted by the reduction in its cohesive energy. From these MD results, we recognize an exponential relationship between the amount of H atom diffusion and the intensity of the SG in nano-sized bending deformation.

1. Introduction

In recent years, to manage many environmental issues, fuel-efficient vehicles have been pursued, and lightweight metallic materials with high strength are extensively sought in industries. Metal forming using the technology of plasticity has been a good solution to provide materials with good mechanical properties for industrial use. Examples are found in the rolling of high-tensile steel sheets for automobiles, the drawing (squeezing) of pearlitic steel wires for suspension bridges, the automobile tire code, and so on. In plastic processing, as the material is continuously provided a plastic strain, it generally gains higher strength. However, as seen in drawn pearlitic steel wires, the material is brittle and frequently breaks, resulting in so-called delayed fracture, which is quite dangerous and must be avoided in the materials in service. The feature of delayed fracture in steel has been extensively studied and is thought to be possibly caused by “hydrogen (H atom) diffusion” inside the material. Thus, the hydrogen embrittlement (HE) of steel has been well recognized [1,2]. There have been many theories concerning HE. One theory notes that the existence of H atoms promotes nucleation and glide of dislocations, and plastic deformation is eventually enhanced there (hydrogen-enhanced localized plasticity: HELP) [3]. Another theory posits that H atoms chemically reduce the interatomic interaction around them, and the formation energy of a fractured surface becomes lower (hydrogen-enhanced decohesion: HEDE) [4]. In spite of the progress of these theories of HE, an experimental approach to capture the dynamic motion of H atoms is still difficult, and therefore any unified view of HE has not yet been reached.

When observing the fracture behavior in high-strength steel due to HE, grain boundaries (GBs) often play an important role in making a fracture path, i.e., many GB fractures are confirmed in HE. As the composition of manganese (Mn) increases, the strength of steel, including the H atoms, largely decreases [5], and its fracture surface becomes more smooth. In that case, the large amount of defects is confirmed by thermal desorption analysis (TDA), which means that H atoms near the GB region are strongly correlated with plastic deformation [6]. As for the deformation mode, H atoms are experimentally visualized using the “hydrogen micro-print” (HMT) technique [7]. In such experiments, when a pearlitic steel wire is forcibly bended, it is confirmed that the segregation of H atoms into GBs is enhanced by the bending stress there.

The behavior of H atoms depending on deformation modes is important, and it can be studied by atomistic simulations with quantitative evaluation as well as qualitative insight. Atomistic simulation has been applied in many previous studies concerning H atoms in metallic materials. They focus on various topics, e.g., the amount of H atoms that are trapped in the GB region and their effect on cohesive energies [8] and the relationship between the GB structure and the diffusion constant of H atoms [9]. These studies have applied molecular dynamics (MD) simulations in cases where the time evolution of individual atoms and atomic systems can clearly and directly be obtained by just a numerical calculation.

Generally speaking, it is well known that the diffusivity of H atoms in metal is greatly affected by ambient pressure and possibly stress fields and temperature. A metallic material that has experienced some plastic working often contains a very high level of residual stress. The residual stress in the outer region of the material, which is nearly in contact with the roll or draw die, exhibits a tensile nature, while the stress near the center region is compressive so as to balance the total stress of the material [10,11]. This means that after severe plastic processing, materials inevitably include a gradient of residual stress, like bended state. The more plastic deformation introduced into a material, the larger the stress gradient (SG) observed in it. In newly developed process techniques for severe plastic deformation, the resulting material usually shows high sensitivity to the existence of H atoms. Thus, the correlation between the SG provided during plastic deformation and the diffusivity of H atoms is certain and should be well recognized. Generally speaking, for the diffusion phenomenon in metals, the diffusion coefficient increases as the hydrostatic stress (isotropic component of the stress tensor) increases in an expanding (positive) direction. This tendency has already been discussed for H atoms, where the diffusion coefficient is higher on the tensile side than on the compression side [9]. However, to the authors’ knowledge, there are few reports on the diffusion behavior of H atoms under the existence of SG. Crucial mechanisms of H atom diffusion that are related to the amount and directionality of atomic movement still remain unknown so far, both in experiments and simulations. Therefore, in this study, we examined this topic using the atomistic simulation method in order to reveal the mechanism.

In this study, we analyzed the diffusion behavior of H atoms near the GB of pure iron (α-Fe) crystals using the MD method. The reason why we focused on the vicinity of the GB region is that HE is remarkably generated there due to a local stress concentration. In MD modeling, several intensities of SG are applied to the bi-crystal specimen of α-Fe in relation to bending deformation of the overall region. This is because drawn steel wires often show a wide distribution of residual stress from the center to the surface region, and the situation is very similar to the free bending deformation of the specimen [11]. In this study, bending deformation was adequately realized by a preliminary deformation calculation. In the bended specimen, the diffusion coefficient of H atoms and its directionality in the bulk region were examined in detail. Finally, we tried to formulate a mathematical relationship between the intensity of SG and the diffusion coefficient of H atoms.

2. Theory and Methods

2.1. MD Theory and Calculation

In molecular dynamics (MD), Newtonian equations of the motion for each H atom, including interatomic interaction potential energy, are numerically solved, and the dynamic behavior and properties of the total atomic system are evaluated. The interatomic potential function used here is an embedded atom method (EAM) formulation that has relative high accuracy needed for H atoms as well as the iron (Fe) crystalline structure. In addition to the Fe-Fe interaction function, which was derived by Mendelev et al. [12], Fe-H and H-H interactions were formulated by Ramasubramaniam et al. [13].

The EAM potential by Mendelev et al. [12] was determined not only by fitting to a perfect α-Fe crystal but also by considering the energetics of self-interstitial atoms; therefore, it will accurately express the energy, especially when defects are included. In the EAM potential function proposed by Ramasubramaniam et al. [13], potential parameters were determined using first principle (ab initio) calculation so as to match the Fe-Fe energy expressed by Mendelev et al. Ramasubramaniam et al. checked the effect of H atom as impurity using a larger calculation cell than that for an ordinary ab initio calculation. Thus, the properties of H atoms in the vicinity of the defect can be reproduced with high accuracy by this potential function, which is expected to be suitable for the purpose of this research.

The accuracy of the interatomic potential used here was first verified. As shown in Figure 1, one H atom is located near the center of a cubic α-Fe single crystal model, which is surrounded by free surfaces in all directions and has a cell length of $a=a_0\times 10$ = 2.867 nm (where the lattice constant $a_0$ of α-Fe is supposed to be 0.287 nm). First, an MD calculation of the structural relaxation was conducted with a constant-temperature statistical ensemble. Figure 2 shows the time average of the radial distribution function (RDF) for the averaged distance between Fe-Fe and Fe-H (orange and blue solid lines are for Fe-Fe and Fe-H, respectively). From their first peak, the stable distances between Fe-Fe and Fe-H are approximately 0.246 and 0.170 nm.

Table 1. Comparison of interatomic distances in α-Fe crystals using MD, ab initio, and theoretical calculations.

Method	The Nearest Neighbor Interatomic Distance [nm] (Error [%])
	Distance Between Fe and Fe Atoms	Distance Between Fe and H Atoms
Molecular dynamics calculation (a)	0.246 (−0.806)	0.170 (6.25)
Ab initio calculation (b),(c)	0.236 (−4.84)	0.158 (−1.25)
Theoretical value (d)	0.248	0.160

^(a) This work, ^(b) Mendelev et al. [12], ^(c) Ramasubramaniam et al. [13], ^(d) calculated from lattice constant $(a_0=0.287 \mathrm{n}\mathrm{m})$ and ordinary theory of interstitial [14,15].

shows the error of the present MD value compared with theoretical and ab initio calculated values. The error of the MD value is −0.806% for the Fe-Fe distance, but it is 6.25% for the Fe-H distance, which is relatively large. However, they are better than ab initio calculations (−4.84% between Fe and Fe). Therefore, this interatomic potential gives a reasonable calculated value regarding the behavior of H atoms.

2.2. MD Modeling

The calculation model used in MD calculation is shown in Figure 3, where the orange and blue spheres indicate the position of the Fe and H atoms, respectively, and the pink spheres show the atoms with velocity constraints. Periodic boundary condition is applied only in the z direction. Because α-Fe crystal is assumed, Fe atoms are arranged in body centered cubic (bcc) lattices so that the $xy$ plane corresponds to the $\left(001\right)$ plane, as shown in Figure 4. Initial positions of H atoms are tetrahedral sites (hereafter abbreviated as T-sites), which are selected at random from all T-sites (as many as needed). They are represented by a red solid circle in Figure 4.

The four-point bending test is adopted here. In this material testing configuration, the maximum bending moment occurs in the wide region between two loading points by indenters. Therefore, near the center of the specimen, the bending stress measured along the thickness should exhibit a linear distribution of normal stress with a stress gradient (SG) according to the linear elastic theory. The initial positions of the H atoms stated above are restricted inside the red rectangular frame in Figure 3, where they are subjected to almost maximum bending moment and are less affected by local stress concentration produced in the vicinity of the two indenters. Additionally, one grain boundary (GB) gets created at the center of the specimen along the x direction. Two crystals separated by the GB plane are rotated oppositely by the same angle around the $z$ axis. The reason why the model includes the GB structure is that both diffusive and non-diffusive H atoms are supposed to exist preferentially in the vicinity of the intergranular boundary. Although mechanical properties strongly depend on the type of GB structure and its energy, the GB structure used here is a well-defined symmetrical tilt boundary, which is specified by a parameter: $\Sigma =$ 5 (the misorientation angle is $\theta =$ 36.5 deg.). This GB can be constructed by relatively small periodic units along the boundary and is therefore supposed to be energetically the most stable [16]. Two orthogonal crystal orientations, $\left[310\right]$ and $\left[3\overline{1}0\right]$, correspond to the $x$ and $y$ directions in Figure 3, respectively. In order to scrutinize the effect of GB, we should compare this bi-crystal model including GB with a single-crystal model. However, unfortunately, when we preliminarily performed the bending simulation of the single crystal model, uniform distribution of SG along the thickness was broken by a strong curvature of the specimen due to the concentrated deformation around the contact with indenters. Technically speaking, by including the GB structure as a kind of buffer of inelastic distortion of the perfect lattice, such unwanted mechanical concentration could be eliminated.

The calculation conditions are summarized in

Table 2. Calculation conditions.

Property		Value
Lattice Constant of α-Fe $a_0$ [nm]		0.287
The number of atoms [−]	Iron (Fe): $N_{\mathrm{F}\mathrm{e}}$	51,095
	Hydrogen (H): $N_{\mathrm{H}}$	10
Hydrogen content $C_{\mathrm{H}}$ [at.%]		0.196 × 10−3
Cell size in x, y, z [nm]		45.00, 12.00, 1.433
Misorientation angle of GB $\theta$ [deg.]		36.9
Temperature T [K]		300.0
Stress gradient (SG) [GPa/nm]		0.000, 0.222, 0.552, 0.935

. The number of H atoms inserted in the lattice is $N_{\mathrm{H}}=$ 10, which seems low, while the number of Fe atoms is $N_{\mathrm{F}\mathrm{e}}$ = 51,095. However, the concentration of H atoms in the present MD model is much larger than the experimental estimation, $C_{\mathrm{H}}=5.602\times {10}^{-8}$ at.%, which is estimated in iron at room temperature via an empirical formula of Kubashewski [17]. The concentration of H atoms near a lattice defect, such as the GB structure, should be elevated, as in the present MD model, though it is still unclear in practice. In any case, if we assume that most H atoms are located near the lattice defect, a concentration of H atoms in the MD model larger than the experimental value will be quite acceptable.

Figure 5 shows the calculation flow. First, a calculation of the structural relaxation at constant temperature is performed with no applied deformation. Then, in order to form SG, a preliminary calculation with deformation by four-point bending is performed. After that, maintaining the deformed state, the diffusivity of H atoms is evaluated as the main calculation. In the preliminary deformation of the four-point bending, two indenters (the constrained atoms colored in pink beneath the MD specimen in Figure 3) are provided with constant velocity in the y direction, $v_y$, and the constrained atoms at the two ends in the x direction (the constrained atoms colored in pink above the MD specimen in Figure 3) are completely fixed in the space, i.e., $v_y$ = 0; bending deformation then takes place on the $xy$ plane. The calculated temperature $T$ is assumed to be in the standard state, and the velocity scaling method is utilized in all calculations to keep a constant temperature of $T=$ 300 K. The SG values after applying bending deformation can be theoretically estimated as shown in the last row of Table 2 (using the elementary bending theory of a straight beam). The SG values obtained by MD calculation will be explained and discussed in the next section.

2.3. Stress Gradient (SG) Obtained by Bending

In the calculation of bending deformation, both of the indenters move with a constant speed at $v_y=$ 0.0, 5.0, 10, and 20 m/s each for the duration of 90 ps, and the displacements ${\delta }_{ymax}$ occur at the maximum value of 0.00 (the case without SG) and 0.45, 0.90, 1.80 nm (the cases with SG), respectively. Hereafter, those four conditions will be called “case 1” (${\delta }_{ymax}=$0.00), “case 2” (${\delta }_{ymax}=$ 0.45 nm), “case 3” (${\delta }_{ymax}=$ 0.90 nm) and “case 4” (${\delta }_{ymax}=$ 1.80 nm), respectively. For reference, the experimental value of actual SG comes in the range of around ${10}^9$ Pa/m, which was estimated from the measurement of residual stress distribution in a drawn steel wire [11]. The SG values obtained in the present MD calculation were in the range of around ${10}^{18}$ Pa/m and were much larger than the experimental values. It is understood that the MD specimen is free from defects and is almost a single crystal, so the elastic range will be much larger than an actual polycrystalline material. However, in “case 4”, which was applied the largest bending, it was observed that mobile dislocations inevitably occurred in the specimen from the contact region to the indenter. It should be noted that if bending deformation exceeds the level as realized in “case 4”, the condition for the elasticity mechanism is no longer satisfied and a different plasticity-dominated mechanism has to be considered regarding the effect of SG.

The distribution of normal stress ${\sigma }_x$ (bending stress) obtained in the four cases is shown in Figure 6. These distributions are produced from atomic stresses of limited atoms located along the measuring plane that shifted in the $x$ direction by −0.5 nm from the GB plane. Basically, the stress and y position have a linear correlation. However, as shown in Figure 6b,c for “case 3” and “case 4”, the increase in compressive stress is saturated in the outer region (minus y position and minus ${\sigma }_x$ stress) of the specimen. This is supposedly due to dislocation emissions or strong constraints near the indenters. Therefore, either outer region $y>$ −2.3 nm in “case 3” or $y>$ −2.0 nm in “case 4” should be technically excluded in estimating the SG values. The individual plots of each atom in Figure 6 show large deviations from the mean stress value. This is naturally caused by the thermal fluctuation of atomic positions and forces at a finite temperature. Nevertheless, the averaged stress values (solid lines in each figure) show a mainly linear relationship with regard to the direction toward the thickness (in the y direction). It also seems that the averaged stresses match the theoretical bending stresses (those estimated by linear elasticity).

Table 3. Stress gradients (SG) values obtained by bending simulation.

Calc. Case	Maximum Displacement of Indenter ${\mathit{\delta }}_{\mathbf{y}\mathbf{m}\mathbf{a}\mathbf{x}}$ [nm]	Stress Gradient $\partial \mathit{\sigma }/\partial \mathit{y}$ [GPa/nm]
1	0.000	0.000
2	0.450	0.222
3	0.900	0.552
4	1.80	0.935

summarizes the calculated values of SG. The resulting state of those specimens with bending will be effectively put to use in the next procedure of this study, in which calculation models will be accompanied by a deformation field of linear SG.

3. Results and Discussion

3.1. Diffusion Behavior of H Atoms

The trajectories of H atoms under the stress state with SG are shown in Figure 7. Here, orange-colored plots show Fe atoms, while blue-colored plots and lines show H atoms and their trajectories (other indications are as follows: pink-colored plots show the atoms with velocity constraint, a red-colored broken line represents the GB plane, and the green-colored solid line represents the position of the free surface). As indicated by the black arrows, the moving direction of each H atom in the bulk region has some directionality (which may be called diffusion direction). The diffusion direction of H atoms for “case 1” shows no common feature, whereas in “case 3” and “case 4”, it can be observed that many H atoms are moving toward the positive y direction. In “case 4”, the region of Fe atoms more condensed near the end of the figure presents mobile dislocations.

Table 4. Averaged diffusion coefficient of hydrogen atoms.

Calc. Case	Maximum Bending Stress ${\mathit{\sigma }}_{\mathit{m}\mathit{a}\mathit{x}}$ [GPa]	Diffusion Coefficient of H Atoms ${\mathit{D}}_{\mathit{H}}$ [m2/s]	Diffusion Coefficient of Fe Atoms ${\mathit{D}}_{\mathit{F}\mathit{e}}$ [m2/s]
1	0.000	3.898 × 10−10	9.590 × 10−12
2	1.810	$5.202\times {10}^{-10}$	$1.385\times {10}^{-12}$
3	2.698	$7.925\times {10}^{-10}$	$1.378\times {10}^{-12}$
4	5.158	$16.44\times {10}^{-10}$	$3.048\times {10}^{-12}$
-	Experimental value	$D_{\mathrm{H}}{,}_{\mathrm{e}\mathrm{x}\mathrm{p}}=$ $95.48\times {10}^{-10}$(a)	$D_{\mathrm{F}\mathrm{e}}{,}_{\mathrm{e}\mathrm{x}\mathrm{p}}$ = $5.41\times {10}^{-23}$ (b)

^(a) Hagi [18], ^(b) Iijima, et al. [19].

shows the diffusion coefficients $D_H$ and $D_{Fe}$ of H and Fe atoms obtained under each condition, compared with the experimental values $D_{\mathrm{H}}{,}_{\mathrm{e}\mathrm{x}\mathrm{p}}$ [18] and $D_{Fe}{,}_{exp}$ [19]. Here, $D_{\mathrm{F}\mathrm{e}}{,}_{\mathrm{e}\mathrm{x}\mathrm{p}}$ is the self-diffusion coefficient of α-Fe at 766 K. The diffusion coefficient $D_i$ of the objective element i (Fe or H) is estimated from its mean square displacement (MSD) during MD calculation using the Einstein equation shown in Equation (1).

$$D_i={\displaystyle \frac{1}{N_i}}{\sum }_j^{N_i}{\displaystyle \frac{{\left\{{\mathit{r}}_j\left(t\right)-{\mathit{r}}_j\left(0\right)\right\}}^2}{6t}}$$

(1)

where $N_i$ is the number of $i$-type atoms, $t$ is the duration of diffusion measurement ($t=$ 0.5 ns here), and $r_j \left(t\right)$ is the coordinate of atom $j$ at time $t$. From Table 4, $D_{\mathrm{H}}$ tends to increase as the value of SG increases. $D_H$ is approximately 40 times larger than $D_{\mathrm{F}\mathrm{e}}$ for “case 1”, 375 times for “case 2”, 575 times for “case 3”, and 540 times for “case 4”, that is, the diffusivity of the H atom becomes significantly larger as SG increases, unlike $D_{\mathrm{F}\mathrm{e}}$. This means that the diffusion behavior is largely different between H and Fe atoms. It can also be seen that the calculated diffusion coefficient $D_{\mathrm{H}}$ is much smaller than the experimental value ${\mathrm{D}}_{\mathrm{H}}{,}_{\mathrm{e}\mathrm{x}\mathrm{p}}$. In the experiment by Hagi, non-trapped H atoms were found to be responsible for the observed diffusion coefficient [18]. However, in MD calculations, as confirmed in Figure 7, H atoms tend to stagnate in the GB region, and so $D_{\mathrm{H}}$ becomes low value. As the GB structure works as a strong trap site of H atoms and their atoms are bonded with stronger energy than other trap sites such as vacancies and precipitates [20], it is reasonable that H atoms stay there for a long time. The concentration of H atoms to a strong binding site (e.g., lattice defects such as dislocation or GB) has already been discussed [21]. For instance, using Mclean’s equation described in [21], we can estimate a theoretical value that is hundreds of times larger than the present averaged value, $C_H=0.196\times {10}^{-3}$. The present MD results reflect this tendency of concentration qualitatively well.

These results indicate that directional diffusion of H atoms occurs due to the variation of strain or stress provided by bending. It also means that the amount of diffusion increases due to the reduction in the cohesive energy around the H atom. These factors will be discussed in detail in Section 3.2 and Section 3.3, respectively. Because it is supposed that diffusion of H atoms near the surface region is not largely influenced by SG, it is omitted here. Therefore, we will focus on the behavior of H atoms in the bulk near the GB region in this study.

3.2. Influence of SG on Diffusion Directions

Figure 8 shows the relationship between diffusion coefficients in the x and y directions, i.e., $D_{\mathrm{H}}\left(x\right)$ and $D_{\mathrm{H}}\left(y\right)$, where the broken line serves as a reference when $D_{\mathrm{H}}\left(x\right)$ and $D_{\mathrm{H}}\left(y\right)$ are equal. In “case 1”, $D_{\mathrm{H}}\left(x\right)$ and $D_{\mathrm{H}}\left(y\right)$ show very close values. On the other hand, in “case 3” and “case 4”, where the SG is relatively intense, $D_{\mathrm{H}}\left(y\right)$ dominates $D_{\mathrm{H}}\left(x\right)$, that is, the diffusivity of H atoms has a certain anisotropy so that it increases in the direction in which bending stress changes (y-axis direction here).

Figure 9 shows the time transition of the y position of each H atom located in the bulk region. Broken lines represent their positions as the MD calculation starts, while solid lines are those at final time. The moving direction of the H atom is indicated by a horizontal arrow. For “case 3” and “case 4” with larger SG, it is confirmed that a lot of H atoms move to the positive y direction, i.e., toward the tensile side. Because of the bending deformation, the lattice structure of the α-Fe crystal is distorted so as to gain free volume for tensile area but to reduce it for compressive area; consequently, H atoms diffuse more likely toward the tensile side.

3.3. Influence of SG on the Amount of Diffusion

Figure 10 shows the relationship between the cohesive energy $E_{\mathrm{c}\mathrm{o}\mathrm{h}}$ of H atoms calculated by Equation (2) and the SG value, which is denoted by $\partial \sigma /\partial y$.

$$E_{\mathrm{c}\mathrm{o}\mathrm{h}}=E_{\mathrm{H}}-E_{\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{k}},$$

(2)

where $E_{\mathrm{H}}$ and $E_{\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{k}}$ are averaged potential energies of Fe atoms just around the objective H atom and in the whole crystal. Averaging for E_H is performed within a spherical region with a radius $r$ from the H atom, where $r$ is chosen as the arithmetic mean of the first and second nearest interatomic distances $r_1$, $r_2$ as to the T-site, i.e., $r=(r_1+r_2)/2=$ 0.209 nm. While cohesive energy is strictly defined as the energy difference between the whole system with and without the H atom, the cohesive energy discussed here is an approximation of the local vicinity of the targeted H atom. However, such evaluation using Equation (2) will provide us with a good insight about diffusion behavior. As shown in Figure 10, the cohesive energy of H atoms tends to decrease as SG becomes larger. This means that the SG reduces the energy barrier required for diffusive movement of H atoms.

As usually explained in textbooks [22,23], the diffusion coefficient of jumping atoms is largely enhanced by the increase in the activation volume around it. Intuitively, in the tensile region, the activation volume will increase and the diffusivity will increase. The dependency of $E_{\mathrm{c}\mathrm{o}\mathrm{h}}$ on the SG exhibited in Figure 10 should be taken into account considering the effect of the activation volume and the stress state itself.

3.4. Relationship Between SG and the Diffusion Coefficient

As found in the textbook of diffusion [22,23] and explained above, in general, the diffusion coefficient of atoms enormously depends on the local volumetric environment, that is, the stress or strain state inside the material as well as ambient pressure [24]. On the other hand, it has been pointed out that the hydrostatic stress hardly influences H atom diffusion because these atoms need a negligible migration volume for diffusion barriers. However, it has also been indicated that hydrostatic stress evidently alters the binding free energy of H atoms to the vacancy or defect and thus changes the escape rate [25]. There have been lots of investigations regarding the dependency of H atom diffusion in the Fe system. The direct approach by MD simulation can capture the dependency of stress as well as atomic mechanisms [26]. Discussion of our atomistic results presented here can also be organized based on analysis of local stress state and its formulation, where a function of diffusivity is made simply in terms of stress or pressure. In fact, experimental results have clearly shown the dependence on pressure (e.g., the case of self-diffusion of gold [22]). On the other hand, the interesting point of our model is that it has nano-sized (very local) stress concentrations accompanying strong SG. Keeping in mind the fact that there have been lots of formulations of atomic diffusion using the stress or pressure state, at this point, we would like to first look at how much such microscopic and strong SG affects the diffusivity of H atoms.

Figure 11 shows the relationship between the averaged diffusion coefficient of H atoms and the SG supplied by bending. The graph employs a logarithmic scale just on the horizontal axis for the value of SG (i.e., semi-logarithmic representation). It is confirmed that there is a linear relationship, so the diffusion coefficient increases exponentially with regard to the SG value. Using an exponential-type least squares method, the mathematical expression may be obtained as Equation (3).

$$D_{\mathrm{H}}=A\hspace{0.17em}\mathrm{e}\mathrm{x}\mathrm{p}\left(B{\displaystyle \frac{\partial \sigma }{\partial y}}\right),$$

(3)

where A and B are parameters, and for the present results, they are $A=3.385\times {10}^{-10}$ ${\mathrm{m}}^2/\mathrm{s}$ and $B=1.633\times {10}^{-18}$ m/Pa, respectively. The approximation by Equation (3) is shown as a dashed line in Figure 11, and it matches the calculated value very well. Some coincidence can be found between Equation (3) and Arrhenius’s relationship, which is expressed by an exponential function. A general Arrhenius’s expression for atomic diffusion as thermal activation process is

$$D=D_0\hspace{0.17em}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{E_a}{k_{B}T}}\right),$$

(4)

where $D$, $D_0$, and $E_a$ are the diffusion coefficient, frequency factor, and activation energy, and $k_B$ and $T$ are Boltzmann constant and temperature. Furthermore, $D_0$ will be decomposed as in Equation (5) [22]:

$$D_0=gf\hspace{0.17em}{v}_0\hspace{0.17em}{a}_0^2\hspace{0.17em}\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{∆S}{k_B}}\right),$$

(5)

where $g$ and $f$ are called the geometric factor and the correlation factor, ${\nu }_0$ is the trial frequency that physically means thermal oscillation, $a_0$ is a lattice constant, and $\Delta S$ is an entropy of diffusion. The change in lattice constant directly affects $D_0$ as in Equation (5) and therefore $D$ as in Equation (4). As described above in this paper, intuitively, it is due to the change in activation volume for diffusion. However, if, for simplicity, $a_0$ is assumed to be constant and therefore the stress does not affect the value, then another factor, i.e., SG, will play a role. To confirm this, we would like to extract SG as a key factor in changing diffusivity after temperature, which must be the main factor.

The exponential argument of Equation (3) does not depend on the temperature, so it ought to be included in the prefactor $D_0$ in Equation (5). Therefore, parameter $A$ in Equation (3) will change according to the temperature and the applied mechanical deformation. However, parameter $B$ inside the exponential argument of Equation (3) must be constant (it can be validated from Figure 11), so the variable of SG, $\partial \sigma /\partial y$, is depend on the diffusion entropy $\Delta S$. In the present MD results, as the SG increases, diffusion of H atoms occurs not in a uniform fashion but with the direction from compressive to tensile fields. As the SG or strain gradient becomes strong in the crystal, there exist compressive and tensile (with respect to mean stress field) regions at the same time. Therefore, such a stress field promotes the diffusion of impurity atoms in a more directional manner. At present, in our MD study, it is true that ΔS has not been well understood and not been directly evaluated. However, the results clearly suggest that $\Delta S$ will increase with the increase in non-uniformity (directionality) in a deformed state (i.e., bending) and stress (i.e., SG), and such a situation increases the H atom diffusion.

Thus, we have performed the formulation of hydrogen diffusion in terms of SG. The system temperature and stress distribution almost completely govern the diffusion of hydrogen atoms in iron. Our results confirm that this dependency can also be applied to atomistic (very tiny-scale) stress concentration by bending of the specimen and by the structure of GBs in it. SG is surely included in such a heavily deformed structure, and it may play an important role for hydrogen to diffuse there.

4. Conclusions

Using MD simulation, the diffusion behavior of H atoms in α-Fe bi-crystals and its response to the applied stress gradient (SG) provided by bending deformation were investigated, and the following findings were obtained.

• Directionality for diffusion of H atoms occurs as the value of SG increases.
• The amount of diffusion of H atoms increases as the SG increases. The diffusion occurs in the direction from compressive to tensile stress field. This is caused by the slight distortion of the crystal lattice by bending deformation, where the moving route of H atoms toward the tensile side is generally expanded.
• The approximated cohesive energy around H atoms becomes smaller as the SG increases. This means that the energy barrier for an H atom to escape from the trap site (T-site) will be reduced by the SG, and the diffusion becomes energetically easier.
• The MD simulations show that there is a certain relationship between the diffusion coefficient of the H atom, DH, and the stress gradient, $\partial \sigma /\partial y$, as expressed by the following equation (A and B are parameters):

  $$D_{\mathrm{H}}=A\hspace{0.17em}\mathrm{e}\mathrm{x}\mathrm{p}\left(B{\displaystyle \frac{\partial \sigma }{\partial y}}\right)$$

Because parameter $B$ in the exponent argument is almost a constant value in constant temperature, it is inferred that the SG is expected to be correlated to entropy of diffusion.

However, in reality, system temperature and stress distribution almost completely govern the diffusional behavior of hydrogen atoms. In addition, we clarified that the atomistic stress concentration by bending of the specimen and by the structure of grain boundaries also play a certain role.