Deformation evolution depending on grain boundary type and hydrogen concentration in nickel investigated by molecular dynamics simulation

The impacts of hydrogen concentration on tensile deformation in the nickel bicrystals with different typical grain boundaries (GB) were investigated by molecular dynamics simulation. The deformation behavior was dependent on the GB type and hydrogen concentration. A critical hydrogen concentration was obtained from a drastic change of the theoretical yield strength. Below critical concentration, hydrogen increased the yield strength and the atomic rearrangement was effectively hindered due to uniform distribution of hydrogen. Above critical concentration, the nickel-hydride formed and caused a sharp decrease in yield strength, which was independent of the GB type.


Introduction
Metallic material working in high pressure and high purity hydrogen environment for long time usually leads to ductile and plastic reduction, crack growth rate acceleration and hysteresis fracture, which is called hydrogen embrittlement (HE) [1][2][3]. So far, researchers have conducted enormous work to study the mechanism of HE, including hydrogen enhanced decohesion (HEDE), hydrogen enhanced local plasticity (HELP), hydride formation and cleavage [4][5][6][7]. However, due to the complexity caused by the interaction between solute hydrogen and dislocation, grain boundary (GB) and crack tip, the mechanism of HE still has great controversy and disagreements.
GB is one of the most important micro-defects in polycrystalline materials, which directly affects the yield strength, brittleness, edge shape and creep of material [8]. In recent years, many attempts have been made to improve the properties of structural material by optimizing grain boundaries. Hydrogen atom and vacancy usually congregate to form defect complexes near GB reducing the cohesive strength of the grain boundary such that it became the weakest link in the system. Whatever the mechanism is, the crystallographic characteristics of the GB are crucial as it is here that the solute hydrogen cause the embrittlement. Bechtle et al [9] experimentally revealed that in the presence of hydrogen concentrations between 1200 and 3400 appm, the lower-energy grain boundary (LEGB) microstructure showed almost double the tensile ductility; also, the proportion of intergranular fracture was significantly lower and the fracture toughness values were some 20%-30% higher in comparison with the higher-energy grain boundary (HEGB) microstructure. The LEGB is characterized by a particular misorientation, low excess free volume and high degree of atomic matching. Martin and Robertson [10] investigated the basic mechanisms of hydrogen-induced intergranular fracture in nickel by a combination of high-resolution scanning and transmission electron microscopy. They proposed that it is necessary to establish a critical hydrogen concentration at the GB to cause intergranular fracture. It can be seen GB structure and hydrogen concentration are two main factors of intergranular failure. However, it is still a challenge to research the deformation mechanisms of materials systematically by means of experimental methods. Due to the complexity of GB structure and the limitations of experimental condition, experimental methods cannot fully obtain all atomic-scale information about the mechanical behavior of GB. By contrast, molecular dynamics (MD) simulation is an effective means to reveal the microscopic behavior. Kuhr et al [11] performed MD simulations of the mechanical response of identical samples with and without hydrogen in the grain boundaries. They emphasized that hydrogen would segregate at the vacancy of GB and weaken the binding force between GB through chemical interaction with metal ions, leading to embrittlement. However, Tehranchi [12] calculated the reduction of theoretical strength on the various symmetric tilt GBs, and concluded that the theoretical strength was not significantly reduced by the presence of hydrogen atoms for all studied GBs. In addition, different deformation mechanisms are reported depending on GB types and bulk hydrogen concentrations [13][14][15].
Based on the facts, it is necessary to deepen understand the influence of hydrogen on the plastic deformation process of GBs to reveal the mechanism of HE. The work aims to reveal the effect of hydrogen concentration on the deformation behavior of the bicrystal with different symmetric tilt grain boundaries (STGB). It provides theoretical support for regulating GB to restrain the initiation and propagation of hydrogen-induced intergranular cracks, reduces the HE sensitivity of the alloy and improves the resistance to hydrogen damage.

Interatomic potential
The accuracy of the interatomic potential will directly affect the accuracy of the simulation results. The widelyused embedded-atom-method (EAM) interatomic potential can be used in the simulation of simple metals and precious metals. Considering the accuracy of simulation in the study of the interaction of pure Ni and the Ni-H binary system, we chose the EAM potential proposed by Angelo [16][17][18]. Therefore, the EAM potential described above is used in the simulation of this work. Therefore, the EAM potential with the r cut of 4.92Å described above is used in the simulation of this work [19].

Simulation details
Two typical types of 〈1 1 0〉symmetric tilt GBs were constructed to study the diffusion behavior of hydrogen and the mechanical response of GB structure to tensile loading, respectively. One was the ∑3 coherent twin boundary (CTB) representing the LEGB (∑3/(1 1 1)/109.47°). The other was the general large-angle grain boundary representing the HEGB (∑9/(2 2 1)/141.06°). The modeling process has been described in detail in previous literature [20,21]. A parallelepiped model containing about 200,000 Ni atoms was established for uniaxial tensile simulation to weaken the interference between GBs. Periodic boundary condition was applied in all directions, while a time-step of 2 fs in NPT (Nose-Hoover thermostat-barostat) ensemble under 300 K and zero pressure were used. After the system reaches an equilibrium state, the virtual strain-controlled tensile deformations were conducted with a constant strain rate of 1×10 8 /s along the y direction. Then, the plastic deformation processes of the two GB structures on the atomic scale can be obtained.
Hydrogen atoms were first introduced into the bicrystal, followed by a thermal relaxation at 300 K. Most studies used artificial initial hydrogen distributions and unrealistically high hydrogen concentrations to research the effect of hydrogen to avoid long time simulations. Therefore, the Monte Carlo (MC) method in grand canonical ensemble (GVT) at 300 K was used to calculate the hydrogen charging, and the reasonable distribution and equilibrium concentration of hydrogen in the system were guaranteed by randomly inserting hydrogen atoms into the matrix [22]. The bicrystal models with different equilibrium hydrogen concentrations were obtained by MC simulation at 300K. Here, the bulk hydrogen atoms concentrations range from 0.0001 at% to 0.08 at%. The initial configurations of the bicrystal with and without hydrogen atoms are shown in figure 1. The local stress distribution was calculated by the Virial theorem. We used the centro-symmetry parameter (CSP) [23] and the common neighbour analysis (CNA) parameter [24] to analyze the structure of the configuration and determine the crystal defect, including the determination of the defect position and type.

Results
The effect of different hydrogen concentrations on the plastic deformation behavior of Ni bicrystals with LEGB or HEGB are studied, respectively. Here, the ∑3/(1 1 1)/109.47°GB representes the LEGB and the ∑9/(2 2 1)/ 141.06°GB representes the LEGB. The bulk hydrogen concentrations for the bicrystal with ∑3/(1 1 1) GB are 0.001 at%, 0.01 at%, 0.03 at%, 0.05 at% and 0.08 at%. The bulk hydrogen concentrations for the bicrystal with ∑9/(2 2 1) GB are 0.0001 at%, 0.001 at%, 0.01 at% and 0.03 at%. Figure 1 shows a detailed view of the ∑3/(1 1 1) and ∑9/(2 2 1) GB structures in Ni. After MC and MD relaxation simulations, the hydrogen atoms diffused along the GB to occupy low energy sites. The viewing direction is along the 〈1 1 0〉crystallographic direction (Z-axis) and atom positions are projected into the X-Y plane for clarity. The ∑3/(1 1 1) GB in figure 1(a) is composed of D type structural units along the entire length of the interface. Similarly, the ∑9/(2 2 1) GB interface is composed entirely of E structural units. The presence of hydrogen in the GBs results in significant changes in the atomic structure of the grain boundary, as shown in figure 1. We note that when the local hydrogen concentration is high enough, the initial structure of the GB is severely damaged. However, the maximum solubility of hydrogen varies significantly for different structures of the bicrystals. Previous studies have shown that the lattice hydrogen content is approximately 2000 appm by thermally charging the specimens in high-pressure hydrogen atmosphere [10]. When the bulk hydrogen concentration is 800 appm in MC simulation, the local hydrogen concentration at the GB is much higher than 2000 appm. Thus, the MD simulations of the bulk hydrogen concentration higher than 800 appm are impractical and meaningless. The maximum hydrogen concentration we calculated is set to 0.08 at% (800 appm). It can be seen from figure 1(g) that when the bulk hydrogen concentration is 800 appm, the ∑3/(1 1 1) CTB is extremely disordered.

Hydrogen trapping maps
The segregation energy criterion reflects the condition in which the GB trapping sites are saturated. Figure 2 shows the hydrogen trapping energy maps for the two typical GBs. The trapping sites of hydrogen atoms are identified by smaller spheres and these are colored according to the segregation energy. All the maps are viewed along a direction paralleled to the tilt axis 〈1 1 0〉. Some sites actually on different layers appear to overlap due to the observation direction. Because of the unique GB structural units, the distribution of dissolved energy for each GB type is also unique. It can be seen from figure 2 that, for the bicrystal models with different structures, the highest segregation energy of hydrogen at the HEGB is much larger than that at the LEGB. To reflect the gap more directly, we plot the segregation energy distribution curves for the two GBs, as shown in figure 3. Figure 3 presents that for the ∑3/(1 1 1) GB, the highest segregation energy is about 0.075 eV, which is far less than that at the ∑9/(2 2 1) GB with the segregation energy of 0.25 eV. In addition to, the possible hydrogen trapping sites with high segregation energy are generally distributed within a region ±5Å from the boundary planes for both the GBs, and the segregation energy at sites located far away from GB planes approaches 0 eV. It indicates that the trapping ability of ∑3/(1 1 1) GB is almost negligible when compared to the HEGB. The trapping ability is strongly depending on the GB type.

Stress-strain response
After investigating the ability of hydrogen segregation at different GB types, we studied the influence of hydrogen segregation on tensile deformation behavior of Ni-H system. The stress-strain curves for the samples with various hydrogen concentrations are shown in figure 4. The stress-strain curve describes the mechanical response of materials under loading [18]. We define the peak value as the theoretical yield strength.
The figures 4(c) and (d) present the yield strength of the bicrystal with ∑3/(1 1 1) LEGB and ∑9/(2 2 1) HEGB as a function of the hydrogen concentration, respectively. On the whole, with the increase of hydrogen concentration, the theoretical yield strength of bicrystals increases to the maximum then decreases rapidly, finally to decrease slowly and finally level off, regardless of the GB structure. However, the maximum yield strength of the bicrystal with ∑3 CTB is much higher than that with ∑9 HEGB. The maximum yield strength for the ∑3 CTB sample with a bulk hydrogen concentration of 0.001 at% is 14.27GPa, while the maximum yield strength for the ∑9/(2 2 1) sample with a bulk hydrogen concentration of 0.001 at% is 6.35GPa. It is worth noting that, the yield strength for different samples reaches the maximum value at the same bulk hydrogen concentration (0.001 at%). On the other hand, when the bulk hydrogen concentration is high enough, the theoretical strength of the samples reduces to around 5GPa, regardless of the GB structure. The decreases of the yield strengths for the bicrystals with CTB and HEGB are 64.96% and 20.82%, respectively. The exceptional but regular deformation behaviour mentioned above indicates that different plastic deformation mechanisms may exist for different samples.

Discussion
So far, the experimental methods cannot directly calculate the decrease in the cohesive energy under a specific hydrogen pressure. However, in order to assess the relationship between the yield strength of the bicrystals with different GBs and the hydrogen concentrations, it is necessary to determine the correlation. Considering even the same hydrogen concentration has very different or even opposite effect on the GB with different structure units, it is absolutely vital that this assessment must be specific to different GB type. Previous studies have shown that the ∑3/(1 1 1) CTB is obviously different from other general GBs. The grain boundary energy of the ∑3/(1 1 1) CTB is about 82.45 mJ m −2 , while the ∑9/(2 2 1) GB composed of E structural units have an energy of 1190.27 mJ m −2 . Thus, we research the properties of the two representative GBs in present work.
Based on the above results, we analyzes the critical bulk hydrogen concentration at different GBs to cause plastic deformation and illustrates the relationship between critical hydrogen concentration for deformation and stress. It is very necessary for revealing the effect of GB types on mechanical properties of materials and the hydrogen embrittlement mechanism.
The hydrogen segregation energy maps reflect the saturation state of the GB trapping sites. Rice et al [25] emphasized that the concept of the segregation energy is a macroscopic average based on the Langmuir-McLean law. By the hydrogen segregation energy maps in the GB, we found that the energy was dependent on the energy of the highest energy sites and the type of GBs. Compared with the interior of the grain, there are more possible positions in the GB and the average dissolution energy is lower. Thus, the local hydrogen concentration at the GB is higher. However, it is practically impossible to fill all available sites. In order to obtain the trapping ability of a GB, we calculated the maximum excess by dividing the total number of probable sites by the area [26]. The results show that the ∑9/(2 2 1) GB has higher maximum excesses of hydrogen than the ∑3/(1 1 1) CTB. In other words, the local hydrogen concentration at the HEGB with E structural unit is higher than that at the LEGB with D structural unit at the same bulk hydrogen concentration, as shown in figures 1(c)-(f).  Figure 4 presents the stress-strain curves for two typical GBs with different bulk hydrogen concentrations and the dependency of the theoretical yield stress on hydrogen content. The results of these simulations show the same trend. The yield strength of the bicrystals increased similarly in the range of lower hydrogen concentration. Obviously, the concentration range is quite different for the bicrystals with different GB structures due to different trapping ability of hydrogen and structural stability. In the low hydrogen concentration ∑3/(1 1 1) CTB, the distribution of hydrogen is uniform, and he imperfections and solute hydrogen could act as alloying elements for tuning properties [11]. For ∑9/(2 2 1) GB, it is clear that the segregation of hydrogen is more obvious at the GB. Since the bulk hydrogen concentration is low, there is sufficient space to trap hydrogen at the GB and the initial E structure remains intact. As the tensile strain increases, Ni atom attempts to slip into the E structure, but this movement is effectively hindered by the uniformly distributed hydrogen atoms in the E structure. Thus, compared to the hydrogen-free case, the yield strength with low hydrogen concentration increases. Of course, when the tensile stress is large enough, the atom eventually be squeezed into the gap, and the GB structure changes irreversibly, as shown in figure 5 [26].
Within a fairly large concentration range, the presence of hydrogen can increase the theoretical yield strength of the bicrystal with ∑3/(1 1 1) CTB. However, the corresponding concentration range of the bicrystal with ∑9/(2 2 1) GB is much smaller. The D structural units at the CTB are compact, which is consistent with Rittner's result [20]. This structure tend to strong local atomic relaxation and low excess volume, showing poor hydrogen trapping ability. The highest segregation energy is about 0.075 eV, only a little higher than the grain interior. Therefore, the hydrogen segregation at the CTB is relatively obvious only when the concentration reaches a certain value. On the contrary, the trapping ability of ∑9/(2 2 1) GB is 3.5 times when compared to the CTB. The E structure is characterized by excess free volume. These provide favorable conditions for hydrogen segregation at GB. Thus, hydrogen segregation and clustering lead to instability even destruction at relatively low concentrations.
When the bulk concentration exceeds a certain value, the yield strength drops rapidly in both the ∑3/(1 1 1) CTB and the ∑9/(2 2 1) GB, as shown in figures 4(c) and (d). In addition, when the bulk hydrogen concentration is high enough, the theoretical strength of the samples reduces to around 5GPa, regardless of the GB structure. It indicates that the GB structure loses its stability and the reasons for the yield decline of both models are probably same. Hinotani et al [27] reported that nickel-hydride is responsible for the intergranular fracture and the lattice parameter of nickel-hydride is about 6% larger than of pure nickel because of a high hydrogen concentration. Therefore, we believe that the primary reason for the decrease in theoretical yield strength of the bicrystals is the formation of hydrides. Specifically, the diffusion of hydrogen to the GB leads to a very high local hydrogen concentration that corresponds, essentially, to the formation of a 'hydride-seed', where the hydrogen concentrations approaching the accepted bulk hydride phase [3]. Figure 6 clearly shows that significant hydride formation occurs at the GB due to energetically favorable hydrogen sites along GB as compared to bulk Ni. The 'hydride-seed' is the weakest link in the samples, and provides convenient condition for dislocation nucleation. Increased loading then allows for further accumulation of hydrogen atoms near the GB, the GB becomes completely surrounded by a growing hydride and finally culminates in brittle cleavage through most of the hydride phase. The results reveal the yield strength of materials with high hydrogen concentration does not depend strongly on the GB structures.
The study seems employed impractically high hydrogen concentrations. However, previous studies have shown that grain-boundary hydrogen concentrations could be magnified in higher-strength materials. Hydrogen accumulation at the crack tip is even more obvious than that at the GB [28]. The simulations at high concentrations make it possible to reveal the relationship between critical hydrogen concentrations and stress.

Conclusions
The impacts of hydrogen concentration on tensile deformation in the nickel bicrystals with different typical GB were investigated by molecular dynamics simulation. The following conclusions were drawn: 1) The relationship between critical hydrogen concentrations and stress with two different GB structures was established. With the increase of hydrogen concentration, the theoretical yield strength increases to the maximum then decreases rapidly and finally leveled off, regardless of the GB structure.
2) Below the critical concentration, the variation of yield strength is strongly sensitive to the GB type at the same lower bulk hydrogen concentration. The difference in yield strength originates from the fact that the expansion or destruction of HEGB due to hydrogen is much stronger than that of LEGB.
3) When the bulk hydrogen concentration exceeds the critical concentration of the corresponding GB, the formation of nickel-hydride is the primary reason for the sharp decrease in yield strength, which is independent of the GB type.