A review of modelling and simulation of hydrogen behaviour in tungsten at different scales

A series of short-ranged two-body potentials, most of which exhibit a Morse-type function, started to be developed at the beginning of the 1960s [14–21]. The main properties/data used to fit these potentials were the lattice constant, the elastic constants, the cohesive energy and the bulk modulus [14–19], as well as the formation energies of the vacancy and the self-interstitial atom (SIA) [20, 21]. In the 1980s, a big step forward was achieved with the birth of many-body potentials such as the embedded atom method (EAM) [22], Finnis–Sinclair (FS) [23], effective medium [24] and other similar types of potentials [25–28]. In the EAM scheme, the total energy of the system V_TOT can be expressed as

$$\begin{equation} V_{{\rm TOT}} =\sum\limits_{i,j} {V_{{\rm pair}} (r_{ij} )+} \sum\limits_i {F_{i}} \left( {\sum\limits_j {\rho_{ij} (r_{ij} )}} \right). \end{equation} \tag{ 1 }$$

By introducing a many-body contribution, i.e. the second term of equation (2), the description of metal cohesion was greatly improved. Furthermore, this method is very convenient to combine other elements, i.e. to build potentials for alloys. In the FS scheme, the F_i term of the EAM formulation is replaced by a square root $(-\sqrt x)$, which seems to be appropriate for bcc or half-filled d-band transition metals. The original FS potential was recently successfully applied by Fu et al to study the thermal conductivities of W at different temperatures [29]. In 1987, Ackland and Thetford developed a modified FS potential [30], which well describes the properties of both bulk W and the vacancy. Another modified EAM potential with second-nearest-neighbour extension was constructed in 2001 [31], and the calculated properties of W were very close to those of the FS–Ackland potential. In 2007, Derlet et al [32] developed yet another FS potential to study point defects in W and the mobility of SIAs. This potential was used to investigate the behaviour of SIAs formed in displacement cascades by Park et al [33]. In 2013, Marinica et al [34] developed a series of EAM potentials, which can well describe the physics of screw dislocations in W in comparison with other potentials. Recently, one of these EAM potentials has been used to investigate the finite temperature properties of screw dislocations in W by Cereceda et al [35].

Potentials based on the bond order formalism have been developed recently, which better describe the different bonding states of the atoms. In 2007, Mrovec et al [36] introduced a bond order potential (BOP) to further study defects such as GB and dislocations in W. Whereas, in 2010, Ahlgren et al [37] introduced a BOP predicting accurate formation and migration energies for point defects in W, as these properties are crucial when one wants to model high-energy collisions such as the ones created by neutrons or high-energy particles impinging a metallic matrix. Finally, an electron-temperature-dependent empirical interatomic potential was developed for W to address the electronic excitations that can take place during ion bombardment [38]. Using their potential to study sputtering erosion, Khaskhouri and Duffy [39] demonstrated the sensitivity of surface damage to the model used to describe the electronic thermal transport. A comparison of some of the potentials described above can be found in [40]. The sensitivity of the primary damage to the interatomic potential is now clearly established [41], and the predictions of these potentials in radiation damage studies can be found in the recent work [42].

Recently, Sand et al [43] investigated the dynamics of defect production under 150 keV collision cascades in W using two different EAM potentials and one BOP, respectively. The EAM potentials show evidence of the formation of interstitial-type dislocation loops, in good agreement with experimental observations. However, the BOP did not reproduce experimental findings. This should be attributed to the difference of the interatomic bonds description between the EAM potential and the BOP.

The first attempt to study the behaviour of H in W by MD can be ascended to the 1970s, in which Creery et al investigated the recombination and adsorption dynamics of H on the W surface with a modified two-body London–Eyring–Polyani–Sato (LEPS) potential [44]. After that, no more potential was developed to study the W–H system for a long period, until 2005, where a valuable progress in the simulation of H radiation on W was made by Juslin et al [45]. They derived a BOP based on a Tersoff potential for carbon (C) [46] for modelling non-equilibrium processes in the W–C–H system. Compared with the EAM potentials, the Tersoff potential has the advantage of being well suited for modelling compounds involving metals (W) and non-metals (C, H), despite the lower efficiency in comparison with the EAM potential. The BOP has been used to simulate low-energy bombardment of D impinging onto crystalline as well as amorphous WC surfaces [47, 48], which demonstrates that the C-terminated surfaces exhibit larger sputtering yields than the W-terminated surfaces.

In 2011, Li et al [49] developed another BOP-type W–H potential, focusing more on the interaction between H and point defects. The calculated formation and migration energies of vacancy and interstitials in W are closer to the experimental data for Li's potential than those for Juslin's potential, as can be seen in table 1. Better consistency with the first-principles calculations was also achieved for the W–H interaction with Li's potential as regards the formation energies of H at different sites: tetrahedral interstitial site (TIS), octahedral interstitial site (OIS) and substitutional site (SS) (table 2). However, Juslin's potential can better describe the W–H_n complex than Li's potential, and it was used very recently to model D bombardment of monocrystalline W [50]. The results obtained from Juslin's potential indicate that the formation of gas bubbles is caused by the near-surface D super-saturation region and the subsequent plastic deformation induced by the local high gas pressure.

Finally, to better understand the synergistic effects of H and He, Li et al [60] fitted a W–H–He BOP in 2012 based on the W–H potential they developed. With these potentials, the properties of W–H and W–He are in agreement with the DFT calculations, and so far it is the only potential that can be used to study the behaviour of W containing He and H atoms simultaneously.

It should be noted that the accuracy and the reliability of interatomic potentials can be assessed by their specific applications. For the potential of W, the EAM potential developed by Derlet et al [32] can well describe the properties of vacancy, SIAs and the 〈1 1 1〉 crowdion. The BOP developed by Mrovec et al [36] can predict the reconstruction of W (1 0 0) surface successfully, making it suitable for both the bulk and the surface of W. The EAM potential by Marinica et al can well describe the physics of screw dislocations in W [34]. For the W–H potential, the BOP developed by Li et al [49] is suitable to describe the interaction between H and W at the defects.

A single H atom can occupy three possible sites in W: the TIS, the OIS and the SS. The DFT calculations demonstrate that the TIS is the most favourable site for H in W, as is also the case in Fe. Table 3 shows the solution energies of H in W calculated by different groups along with the experimental data.

The necessity of considering the ZPE for the light element of H has been demonstrated in the dissolution of H in W [61, 63]. The zero point energy (ZPE) corrections can be obtained as $E_{{\rm ZPE}} =\textstyle{1 \over 2}\sum {hv_{j}}$, in which v_j are the vibration frequencies. The ZPE of H is shown to have a significant effect on the solution energy, but the relative stability of H at the TIS and the OIS remains unchanged. When one takes the ZPE into account, the calculated solution energies are in good agreement with the experimental ones [65].

To shed light on the relative stability of H at the TIS and the OIS, Liu et al [59, 66] explored the underlying physical mechanisms. First, distribution maps of valence charge density were drawn to reveal the features of the atomic bonds. The strongest directional bond forms between H and W for the TIS case, and gives rise to the lowest formation energy. Secondly, local density of states (LDOS) for H and W as the first nearest neighbour (1NN) of H shows a stronger hybridization between H-1s and W-5d for the TIS case as compared with the OIS one. This also directly explains the fact that the TIS is the most stable site for H. Thirdly, the electron localization function (ELF) was found to be higher for the TIS than that for the OIS [66]. Generally, for H dissolution in metals, sites with more localized electrons are favoured to achieve a more stable 1s² valence electronic structure. This explains the higher stability of H at the TIS than at the OIS.

To explore the possibility of forming an H₂ molecule, the interactions between two atoms in W have been extensively studied [59, 67–69]. Two H atoms are shown to have a very low binding energy with no more than 0.1 eV along the 〈1 1 0〉 direction at an optimum distance of ∼2.2 Å. Since the H–H bond length for an H₂ molecule is 0.75 Å, such a low binding energy clearly indicates that self-trapping of H in W is almost impossible, as confirmed by the KMC simulations [67].

Normally, H diffuses through two possible pathways in bcc metals, i.e. TIS → TIS and TIS → OIS → TIS, as suggested by several modelling works [59, 64, 70, 71]. For W, the DFT results [59, 63, 64] show that the H diffusion barrier of the later path is ∼0.18 eV higher than that of the former, indicating that the TIS → TIS path is energetically more favourable for H diffusion. The diffusion barriers of H obtained by different methods are listed in table 4.

Despite similar system size, kinetic energy cutoff, and k-meshes adopted in nudge elastic band (NEB) and climbing image nudged elastic band (CI-NEB) calculations, large discrepancies in the diffusion barriers for H in W can be found (table 4). The reason is unknown, but may come from the fact that the CI-NEB method constitutes a slight but nontrivial modification to the regular NEB, and provides an improved estimate of the reaction coordinate close to the saddle point. As for the discrepancy between experimental and computational results, it is well known that real materials always contain various defects, which can act as traps for the H atoms and make the diffusion barrier appear higher than that in the perfect bulk.

According to the harmonic transition state theory (TST), the Wert and Zener's formulation [72] of diffusion can be applied to calculate the diffusivity of H in W, which can be expressed by

$$\begin{equation} D=D_{0} {\rm e}^{-\frac{E_{{\rm a}}}{k_{{\rm B}} T}}=\frac{{\rm 1}}{6}\lambda^{2}v_{0} {\rm e}^{-\frac{E_{{\rm a}}}{k_{{\rm B}} T}}, \end{equation} \tag{ 2 }$$

where λ is the jump length, v₀ is the vibration frequency, D₀ is the pre-exponential factor and E_a is the diffusion barrier. Using equation (2) and the data of table 4, Johnson and Carter [64] obtain a diffusivity D = 8.93 × 10⁻⁷ exp(−0.38 eV/k_BT) m² s⁻¹. Figure 1 shows the calculation results of H diffusivity in W together with the experimental ones.

Zoom In Zoom Out Reset image size

In addition, Heinola and Ahlgren [63] reported that the quantum effects are significant for estimating the diffusivity of H at low temperatures, similar to the quantum effects in simulating the motion of dislocations observed by Proville et al [73]. For example, at 29 K, the H diffusion coefficient is calculated to be 6.3 × 10⁻¹⁴ m² s⁻¹ with the quantum effect correction, in much better agreement with the experimental result ((1–10) × 10⁻²² m² s⁻¹) than the extrapolated experimental D from high-temperature measurements (6.9 × 10⁻⁷⁵ m² s⁻¹).

Experimentally, H diffusivity in W has been determined based on degassing and permeation measurements. However, the temperature ranges used in these experiments are limited above room temperature, mainly because quantum effects as mentioned above are estimated to play a significant role at lower temperatures. The experimental value of D = 4.1 × 10⁻⁷exp(−0.391 eV/k_BT) m² s⁻¹ obtained by Frauenfelder [65] in the temperature range 1200–2400 K is considered to be the best estimate of H diffusivity in W [8, 9]. Other experimental results [74, 75] are also shown in figure 1 for comparison.

The equilibrium concentration of H in W can be determined by

$$\begin{equation} C=\sqrt {P/P_{0}} \exp \left( {\frac{\Delta ST-E_{{\rm H}}^{{\rm S}}}{kT}} \right) \end{equation} \tag{ 3 }$$

where P is the background pressure and P₀ is the reference pressure (taken as the standard pressure in comparison with the experimental results), k and T are the Boltzmann constant and the absolute temperature, respectively. ΔS is the solution (formation) entropy in reference to H₂ gas of standard pressure, which is chosen as −5.5 k for H in W [76]. Liu et al [77] calculated H equilibrium concentration at typical temperatures ranging from 600 to 2000 K, as listed in table 5. For example, the solution concentration is 2.3 × 10⁻¹⁰ at 600 K and 1.8 × 10⁻⁷ at 1000 K, respectively.

The internal gas pressure in H bubbles is estimated to be up to a few tens of GPa, and H bubbles will thus exert a large strain to the surrounding lattice. This is believed to be an important factor affecting H bubble formation. The effects of strain on H dissolution and diffusion behaviours have been investigated employing DFT [78, 79].

Under an isotropic (hydrostatic) strain [78, 79], H solution energy decreases 'monotonically' with increasing tensile strain but it increases with increasing compressive strain, as shown in figure 2(a). This is in good agreement with the linear elasticity theory expectations. The TIS remains to be the preferred site under all strains. In addition, the H diffusion energy barrier under isotropic strain has also been demonstrated to decrease (increase) with increasing tensile (compressive) strain [78]. Consequently, the isotropic tensile strain makes the dissolution and diffusion of H easier in W, while the isotropic compressive strain makes them more difficult.

Zoom In Zoom Out Reset image size

When an anisotropic strain is applied, the lattice symmetry is broken and two different TIS and OIS configurations need to be considered, denoted here as TIS-I, TIS-II, OIS-I and OIS-II, respectively [79]. As shown in figure 2(b), the results at the TIS-II and the two OISs again show a monotonic strain dependence. Surprisingly, at the TIS-I the solution energy shows a non-monotonic dependence on strain, and most importantly, it 'effectively' decreases with the increase in both signs of strain. It is important to realize that the common expectation of a monotonic dependence of H solution energy on strain has a prerequisite condition that the H atom stays at the same location under strain, which is also a presumed condition in most studies so far in the literature. However, this condition can break down because of the broken symmetry due to the anisotropic strain, which is especially important for H at the TIS-I. Figure 2(c) shows an example of H movement in W under strain when initially placed at the TIS-I, as obtained from our first-principles calculations. It is such unusual strain-induced H motion that is responsible for the observed anomalous non-monotonic dependence of H solution energy on strain, leading to energy decreasing under both signs of strain.

The strain in the W lattice surrounding the H bubble must be highly non-uniform and anisotropic due to the low-symmetry surface environment and the irregular size and shape of H bubble. Since the TIS-I remains to be the preferred site with lower energy than the TIS-II and the two OISs under strain, as shown in figure 2(b), this unusual TIS-I behaviour will dominate H solubility in W. This implies that H solubility can be enhanced by applying either a tensile or compressive anisotropic strain. A strain-triggered cascading effect on H bubble growth in W as well as other bcc metals is proposed [79], which is schematically shown in figure 3. H solution energy will be lower with increasing anisotropic strain, leading to an increase in the H concentration, which facilitates further growth of H bubble. As the bubble grows bigger it induces even larger anisotropic strain around it which in turn attracts more H to make the bubble grow even larger.

Zoom In Zoom Out Reset image size

The accumulation and the retention of H can significantly influence the mechanical properties of W. The ideal strength is the stress required to yield or break a perfect crystal. It forms an upper limit to the strength of the material, which is of both scientific and engineering interest and can be determined by a first-principles computational tensile test (FPCTT). To understand the mechanical properties of materials, one has to investigate the ideal strength of both perfect and defective systems. In the FPCTT, a uniaxial tensile strain is applied to the chosen crystalline direction, and the corresponding stress is calculated according to the Nielsen–Martin scheme [80, 81]. The tensile stress σ is calculated from

$$\begin{equation} \sigma {\rm =}\frac{{\rm 1}}{\Omega (\varepsilon )}\frac{\partial E}{\partial \varepsilon}, \end{equation} \tag{ 4 }$$

where E is the total energy and Ω(ε) is the volume at a given tensile strain of ε.

The ideal tensile strength of W has been calculated using DFT by many groups [82–84]. As shown in figure 4(a), the ideal tensile strength of bulk W in the [0 0 1] direction is 29.1 GPa according to the calculation by Liu et al, which agrees with the values of 28.9 GPa obtained by Sob and Wang [82] and 29.5 GPa by Roundy et al [83], respectively. For the [1 1 0] and [1 1 1] directions, Liu et al found the ideal tensile strength to be 49.2 GPa and 37.6 GPa (figure 4(a)), whereas Sob and Wang [82] obtained 54.3 GPa and 40.1 GPa, respectively. The [0 0 1] direction is shown to be the weakest direction in W.

Zoom In Zoom Out Reset image size

When an H atom is embedded, the ideal tensile strength of W along the [0 0 1] direction decreases to 27.1 GPa at a strain of 12% (figure 4(b)), which is about 6.9% lower than that of pure W, suggesting that H will degrade the ideal tensile strength of W [66]. The degradation can be attributed to the fact that H weakens the binding and cohesion among W atoms.

Norskov et al reported that H isotopes have a strong binding to the vacancy in metals [85, 86]. Effective-medium theory indicates that the energy of H in a metal is a function of the metal electron density [85]. The binding energy of D with a vacancy is then calculated to be a result of the interaction between D and the metal valence electrons, which only depends on the average electron density surrounding D. There exists an optimal electron density for which the D-vacancy binding energy reaches a minimum [86]. However, W was not considered in these early works.

The first full and detailed DFT description of H trapping and bubble formation in metals was reported by Liu et al [87]. By calculating H embedding energy in the vicinity of a vacancy, Liu et al found that the most stable site for a single H atom is a position ∼1.28 Å off the vacancy but close to the OIS, with a binding energy of 1.18 eV. Furthermore, H was suggested to bind onto an isosurface of the same charge density (0.11 electron Å⁻³) with six minimum sites surrounding the vacancy [87]. Johnson and Carter [64] reported similarly that H prefers to sit at an octahedral-like site with a binding energy of 1.24 eV. If the ZPE of H is considered, the binding energy becomes 1.41 eV, 0.17 eV higher than that without the ZPE. Heinola et al [88] also showed that the lowest energy site for H in the vicinity of a vacancy is the distorted octahedral site, where H is bound to its nearest W, and the binding energy is 1.28 eV (1.43 eV with the ZPE correction).

In addition to the vacancy, the SIA as another point defect also plays an important role on the H behaviour in W. However, only limited studies focused on the interaction between H and the SIA in W. The 〈1 1 1〉 SIA geometry has been demonstrated to be the most stable in W [88, 89]. Heinola et al [88] reported that an SIA cannot serve as an effective H trapping centre close to room temperature. However, in the light of the work on the role of strain on H dissolution, it can be expected that the presence of H could facilitate the formation of SIA clusters and loops, which can trap H efficiently.

When many H atoms are trapped by one vacancy, the binding energy varies with the numbers of H atoms. Figure 5 represents the binding energy as a function of H number obtained by different groups all using the first-principles calculations [64, 71, 87, 88, 90–92]. Despite the difference in the absolute energy values, the binding energy exhibits the same trend: it decreases when the number of H in the vacancy increases.

Zoom In Zoom Out Reset image size

Liu et al [87] propose that the additional H binding depends on the isosurface of the optimal charge density for H, and that the H–H interaction at the vacancy is repulsive. As more H atoms are added, the isosurface of optimal charge density of H shrinks to the vacancy centre and disappears finally, leading to a binding energy increase. Johnson and Carter [64] considered this trend as a result of subsequent H–H repulsion related to H–H distances. Two H atoms prefer to occupy opposite OISs separated by a distance of 3.14 Å, rather than sit in neighbouring OISs (separated by 2.22 Å). This suggests that the H–H interaction becomes repulsive when the H atoms get close. It can be noticed, however, that the first two H atoms have the same binding energies (1.24 eV), which implies that the repulsion is weak when the H atoms stay apart from each other. The repulsion between H atoms becomes larger when the number of H atoms increases because the distance between them decreases (e.g. 2.07 Å for five H atoms and 1.95 Å for six H atoms) resulting in a lower binding energy. According to Heinola et al [88], the binding energies of the first two H atoms were found to be identical because H couples to the surrounding W. In addition to the binding energy of H atoms, Johnson and Carter determined the energy barrier to release one H atom from a vacancy trap as 1.79 eV [64].

The binding energies obtained correspond to specific H configurations. Liu et al [87] gave a complete picture of the H atoms trapped onto the isosurface of optimal charge density (see figure 6). The isosurface eventually disappears with the formation of the H₂ molecule with a H–H distance of 0.78 Å at the centre of the vacancy when 10H atoms are added. Heinola et al [88] found that multiple H atoms stay at opposing distorted O sites, forming triangle, tetrahedron, square pyramid and square bipyramid depending on the number of H atoms trapped by the vacancy. Different H configurations are proposed by Ohsawa et al [90, 91], as shown in figure 7, in which H moves close to TIS from OIS when the H number increases and the H₂ molecule forms with 14H atoms at the vacancy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The maximum number of H that can be trapped by the monovacancy is reported to be at least 6 by both Johnson and Carter [64] and Heinola et al [88], 10 by Liu et al [87], 14 by Ohsawa et al [90, 91] and 12 by You et al [92]. Johnson and Carter [64] and Heinola et al [88] found that the lowest energy site for H to sit next to the vacancy is the OIS, thus they only explored the 6 possible OISs up to 6H atoms in W. Note that a monovacancy can trap up to 6H in Fe [91], V [91, 93] and Cr [91], 9H in Mg [94] and 10H in Al [94]. Heinola et al [88] further considered temperature effects, and reported that 5H can be trapped at the monovacancy in W at room temperature.

Sun et al [95], using a newly developed W–H potential [49], found that as many as 20H atoms can be trapped at the monovacancy in W [95]. They also investigated the effect of temperature and the number of H that one vacancy can hold. A statistical average was made to evaluate the number of H atoms trapped at the vacancy at different temperatures. As shown in figure 8, this number decreases when the temperature increases, suggesting that temperature weakens the trapping capacity of vacancy. At 300 K, only ∼6.5H atoms can be trapped by a monovacancy, and it decreases to ∼5H at 900 K. This confirms that the thermodynamic estimation of the number of trapped H atoms at the vacancy in W is basically correct, as it gave 5H at room temperature [88]. Temperature will largely decrease the number of H atoms that a vacancy can hold at 0 K and alter the configurations.

Zoom In Zoom Out Reset image size

All these calculations allowed us to propose a scenario on how H bubbles nucleate from the formation of H₂ molecule in W [87, 95, 96]. As shown in figure 9, when H atoms are trapped by the vacancy, they prefer to occupy the internal surface of the vacancy instead of the vacancy centre because of the H–W interaction. With H atoms continuously diffusing in, they will saturate the internal surface of the vacancy to form a 'screening layer', which can effectively screen the interaction between the further trapped H and the surrounding W atoms, making H stay at the vacancy centre. This will lead to the formation of an H₂ molecule at the vacancy centre finally, which is considered to be the preliminary stage of H bubble nucleation. This mechanism is generic and can be applied to other vacancy-type defects such as GB [96] and dislocations as well as other metals. The mechanism can explain the experimental observations that H bubbles prefer to form at vacancies or vacancy clusters [13, 97], GBs [9] and dislocations [8] in metals.

Zoom In Zoom Out Reset image size

On the other hand, H can induce vacancy formation in W, which has been investigated via MD simulations [98]. It is shown that vacancies are created due to the presence of H in W and it is closely related to the formation of certain self-interstitial structures. In particular, a hexagonal cluster structure with a 〈1 1 1〉 central axis is observed in W, and it is the primary structure coming along with vacancy formation. In addition to the hexagonal cluster, a linear crowdion structure is also observed to be associated with vacancy formation, which can evolve into a part of the hexagonal cluster afterwards. This contributes to the understanding of the early stage of the H bubble formation in W under H plasma irradiation.

According to the mechanism proposed above, it can be predicted that an H bubble can form with the condition that the H concentration is beyond a certain value. Xiao et al [99] reported that the spontaneous formation of mono-vacancies as well as SIAs will occur with the ratio of number of atoms of H to W exceeding 0.5, which can be considered as a preliminary step for the bubble formation. It is thus necessary to determine concentrations of both H atoms as well as H–vacancy complexes. Ning et al [100] examined the nD–vacancy complex concentration as a function of depth in W employing an improved CD model. The dynamic processes as well as different mechanisms of D accumulated into the vacancy in a long depth scale in W were revealed. The concentration of nD–vacancy complexes increases when the number of D atoms increases in the complexes at the near surface (0–0.2 µm) and micrometre depths (3–20 µm) due to the trapping effect of D by the vacancies created by the ion implantation. The nD–vacancy complex concentration decreases at the sub-surface (0.2–3 µm) because the SIAs created by ion implantation recombine with the vacancies in this region.

Using a combination of first-principles data with thermodynamics models, Sun et al derived the equilibrium concentration of H as a function of pressure and temperature in both W and Mo. It was found that when the pressure reaches a critical value, the trapping probability for H in a vacancy exhibits a sharp increase, suggesting a considerable formation and rapid growth of H–vacancy complexes. Such a critical pressure and the corresponding critical H concentration are defined as the values when the concentration of H at one certain mH–V complex first becomes equal to that of H at the interstitial, which are 24 ppm/7.3 GPa and 410 ppm/4.7 GPa at a temperature of 600 K in W and Mo, respectively. The predicted critical H concentrations in W are basically consistent with the experimental observations [101].

When considerable H accumulates into a vacancy, a bubble will form and blistering will occur. An attempt at estimating blistering using FEM has been made by Enomoto et al [102]. Mechanical parameters such as Young's modulus and the critical yield stress of the layers in the bulk continuum W were used to determine the internal gas pressures of the blisters. A thin circular opening was introduced, allowing the internal hydrostatic pressure to increase until the upper layer was plastically deformed to the shape of an H blister. The final internal hydrostatic pressure needed to reproduce the experimentally observed H blister shape was found to be 1.7 GPa. From the internal gas pressure derived, the number of H atoms contained in the internal volume of the H blister was estimated using the van der Waals equation of state and was found to be 3.0 × 10²¹. Such a value is about six times higher than the amount of D retention determined in the experiments (5 × 10²⁰). The reason for this discrepancy may be that the formed H₂ molecules can escape from the blisters.

The solution energies of a single H atom at the GB at all potential interstitial and substitutional sites and the effects of vacancy have been investigated using the DFT [96]. At pure W Σ5 (3 1 0)/[0 0 1] GB, it was found that a single H prefers to sit at the interstitial site with a solution energy of −0.23 eV, and it binds onto an isosurface of charge density (0.14 e Å⁻³), as shown in figure 10(a). Compared with the monovacancy in W, the GB provides a smaller isosurface of optimal charge density, indicating that a GB can hold fewer H atoms than the monovacancy in the bulk. In the presence of a vacancy in the GB, H is stable at the vacant space with a solution energy of −0.63 eV, which is much lower than that of the vacancy in the bulk.

Zoom In Zoom Out Reset image size

Figure 10(b) shows the large isosurface shrinkage that takes place when one H atom is added in the GB (figure 10(a)). Since the optimal density region still exists, one can predict that it is energetically favourable for the GB to trap a second H atom. The largest binding energy for 2H in the GB was found to be −0.13 eV, with an equilibrium distance of 2.15 Å. The second H atom still prefers to bind onto the isosurface of optimal charge density, but the isosurface almost disappears when the second H is introduced (figure 10(c)). Consequently, the GB can trap at the most 2H atoms without the formation of an H₂ molecule. These results indicate that the experimentally observed H bubble formation at GBs in W should take place via a vacancy trapping mechanism due to the much lower formation energy of vacancy at the W GB as compared with that in the bulk W [96].

On the other hand, both GBs and dislocation loops are demonstrated to be able to serve as precursors of H bubble formation based on the DFT calculations [106]. The Σ3(1 1 1) title GB can trap as many as six H atoms before significant sliding occurs at low temperature, and that the trapped H atoms can weaken severely the cohesion across the boundary. At room temperature, a dislocation loop in a (1 0 0) plane can trap 4H per (1 × 1) unit, and the detrimental effect of H was found to be strong enough to break the crystal. Indeed, the decohesion effect introduced by the trapped H is associated with the free volume of the defects, i.e. larger free volume induces stronger decohesion. Furthermore, no H₂ molecule was found to form in the GB and the dislocation loop before fracture of W bonds starts, consistent with the previous study by Zhou et al [96].

Zhou et al [96] examined the diffusion paths for H jumping from the W bulk to the W Σ5(3 1 0)/[0 0 1] GB via the interstitial mechanism. It is shown that the diffusion barrier becomes lower when H approaches the GB. The barrier is in the range 0.13–0.16 eV at one-atomic layer separation from the GB centre. This suggests that it is quite difficult for a segregated H to escape out of the GB due to the much lower segregation energy.

Toussaint et al demonstrated the role of the GB on H diffusion in W by employing the MD simulations [107]. The migration paths with minimum activation energy of H in the GB have been determined, and the network of possible migration paths and the corresponding activation energies are shown to be local and asymmetric. The effective diffusion rates for H in the GB was demonstrated to be enhanced as compared with that in the bulk, suggesting that the GB can act as fast transport channels for H.

The effect of H on the mechanical properties of GBs has also been investigated using the first-principles calculations [108]. The tensile strength of a clean Σ5 GB was calculated to be 36.23 GPa, about 22% lower than that of bulk W along the [3 1 0] direction. When H is introduced, the tensile strength of the GB decreases to 32.85 GPa, about 9% lower than that of the clean GB. The Griffith fracture energy of the H segregated GB is 2.48 J m⁻², lower than that of the clean GB by about 4%. The difference of solution energies for the H atom at the GB and the fracture free surface is calculated to be −0.08 eV, indicating that H is a GB embrittler according to the Rice–Wang thermodynamic theory [109]. The Bader charge analysis suggested that the physical origin for H induced strength reduction is related to the charge transfer induced by the presence of H in the GB.

The geometry, electronic structure and stability of the low-index W surfaces, i.e. (1 0 0), (1 1 0) and (1 1 1), have been investigated both experimentally [110–114] and theoretically [115–123]. Yonehara and Schmidt [110] reported $(\sqrt 2 \times \sqrt 2 )R45^{\circ}$ reconstruction of the W (1 0 0) surface by cooling below ∼370 K. Debe and King [111] found the cooled W (1 0 0) $(\sqrt 2 \times \sqrt 2 )R45^{\circ}$ surface at a temperature below 190 K displaying 2mm point-group symmetry. Weng et al [113] observed three bands of surface resonances based on an extensive study of electronic states of the W(1 0 0) surface above room temperature. Fasolino et al [117] concluded that the interaction between surface atoms are repulsive based on calculations of the W(1 0 0) surface photon spectrum and the known reconstruction of the W(1 0 0) surface. However, Joubert [119] suggested that the attractive interaction between the surface atoms drives the reconstruction of the W(1 00 ) surface. Furthermore, Chen et al [120] reported that the interaction between the ideal W(1 0 0) surface is weakly attractive by investigating the surface strain of W(1 0 0) surface. Biswas et al [123] reported highly dispersive linear bands corresponding to the surface states and signature of Dirac cones by investigating the surface electronic structure of the W(1 1 0) surface. Holzwarth et al [121] found that the W(1 1 1) surface has stronger relaxations due to its more open atomic structure in comparison with other surfaces.

Much effort has been devoted to understanding the interaction between H and the W(1 0 0) and (1 1 0) surfaces. H adsorption on the surface will induce many changes in the structural and electronic properties. Barnes and Willis investigated the vibrational modes of atomic H adsorbed on the W (1 0 0) surface, and found that the reconstruction of the W (1 0 0) surface is induced by adsorption of H [112]. Barker and Estrup reported the transition temperature T_c, at which the surface undergoes reconstruction, was found to increase almost linearly from 250 to 380 K with H coverage θ varying from 0 to 0.3 [114]. The atoms on the reconstructed (1 0 0) surface are laterally displaced along the $[1\,\bar{{1}}\,0]$ direction due to the bonding to their nearest neighbours located in the sub-surface layer, thus forming a zigzag pattern advancing perpendicularly to the bonds with the sub-surface atoms. The phase diagram of adsorbed H on the (1 1 0) surface is well known [124–126]. There exist several ordered phases at low temperatures: a 2 × 1 phase at coverage below 0.5 monolayer (ML) and a 2 × 2 structure above 0.5 ML. Disordered structures will be present at temperatures above 200 K and 250 K, respectively. Further adsorption restores the 1 × 1 periodicity of the pure surface when the H coverage reaches saturation. Blanchet et al reported a reduction in the symmetry of the surface for an H coverage greater than 0.5 ML [127]. Kwak et al found the reconstruction of the W (1 1 0) surface induced by H adsorption, which consists of a shift of the W top layer in the $[1\,\overline {1\,0} ]$ direction [128].

The stability of W surface has been investigated by the MD method. The W (1 1 0) surface has been demonstrated to be the most stable W surface by most of the potentials in comparison with the (1 0 0), (1 1 1) and (2 1 1) surfaces, as listed in table 1 [23, 31, 32, 36, 49], in good agreement with the DFT calculations. However, the value of the surface energies given by these potentials is lower than that obtained from the DFT [23, 31, 32, 36, 49]. For surface reconstructions, employing the MD method, Mrovec et al [36] correctly predict the features of the reconstructed (1 0 0) surface using the BOPs, and Juslin et al [45] reproduce the reconstructed (2 1 1) surface.

6.2.1. The (1 0 0) surface.

On the reconstructed W (1 0 0) surface, a single H atom prefers to be adsorbed at the bridge sites according to the experimental studies [112, 129]. The first-principles calculations based on DFT [64, 130] show that the adsorption site turns out to be two different twofold bridge sites, i.e. a short-bridge (SB) site and a long-bridge (LB) site, where the W–W distances are 2.82 Å and 3.57 Å, respectively (figure 11(a)). The single H atom is predicted to be preferentially adsorbed on the SB site, as shown in table 6. Heinola and Ahlgren [130] calculated the H vibrations on the SB site, and the vibration energies of the different modes agree well with those obtained from the experiments [112]. In addition, Johnson and Carter [64] found that the SB site contracts further when H is adsorbed [64].

Zoom In Zoom Out Reset image size

Table 6. Adsorption energy of H on the W surfaces. The ZPE correction is taken into account for the values in parentheses.

	Experimental		Computation
Surface	(1 0 0)	(1 1 0)	(1 0 0)		(1 1 0)
			SB	LB	TF
$E_{{\rm H}}^{{\rm ads}}$ (eV)	0.7 [131],	0.71 [133]	0.91 [130],	0.47 [130],	0.75(0.71) [64]
	0.82 [132]		0.92(0.88) [64]	0.49(0.46) [64]

6.2.2. The (1 1 0) surface.

6.3.1. The (1 0 0) surface.

6.3.2. The (1 1 0) surface.

6.4.1. The (1 0 0) surface.

6.4.2. The (1 1 0) surface.

Alloying elements on the W surface can have significant effects on the microstructure and its evolution, which will affect the behaviour of H. The Ni, Co and Fe alloying elements were found to wet all the W surfaces for at least one monolayer, which leads to the formation of thermodynamically stable ultrathin coatings of these elements [147]. The adsorption of Be will form a very thin film of Be–W alloy on the W surface [148, 149]. The interaction between low-energy atomic C with W surface indicates that the unscattered C atoms are absorbed on the top W surface layer in the low-energy range below 10 eV, while the mean range increases with increasing incident energy above 10 eV, which is attributed probably to the channelling effect [150].

The behaviour of H with these 'alloyed' surfaces has been investigated so far only for C on the WC surface, which indicates that the adsorption of H on the WC surface differs from that of H on the W surface. Gaston and Hendy [151] found that both the adsorption energies and partial charges of the adsorbed H atoms depend strongly on the surface termination (C or W) and the nature of the binding site. Marinelli et al [152] reported that H is energetically favourable sitting on the top site at 1.0 ML coverage of H, while the surface occupation of both W and C sites becomes possible beyond this coverage. Furthermore, the barrier to desorption of H can be lowered by adding additional layers of C to the WC surface. Träskelin et al [153] investigated the low-energy bombardment of D impinging onto the WC surfaces by means of MD simulations, and it is shown that prolonged bombardment will lead to the formation of an amorphous WC surface layer, regardless of the initial structure. Furthermore, the sputtering yield of the C-terminated surfaces is larger than that of the W-terminated surfaces [153].

7.1.1. Dissolution and diffusion.

7.1.2. Accumulation.

7.2.1. Interaction of H and He in intrinsic W.

7.2.2. Single H atom at a He–vacancy complex.

7.2.3. Multiple H atoms at a He–vacancy complex.

Theoretically, there exist 24 equivalent most stable sites for H atoms surrounding the He–V complex. However, all these sites cannot be occupied simultaneously by H atoms because of the H–H repulsive interaction [59]. Diffusion of H towards a vacancy depends on the H concentration. Qualitatively, under high H concentration, many H atoms can diffuse to the He–vacancy complex at the same time (the simultaneous way), whereas under low H concentration, H atoms will diffuse to the He–vacancy complex sequentially (the sequential way). The trapping energy of additional H atoms segregating to the He–V complex has been calculated for both 'the simultaneous way' [61] and 'the sequential way' [61, 157], as shown in figure 12.

Zoom In Zoom Out Reset image size

It was found that, on the one hand, for 'the simultaneous way', the trapping energy per H atom increases with the number of H atoms. When the number of H atoms is up to 13, the average H trapping energy is still negative (−0.50 eV), as shown in figure 12. It is lower than that of H at the TIS far away from the He–V complex, which indicates that the number of H atom trapped by the He–V complex should be larger than 13. After carefully checking the atomic configurations of each H addition, Zhou et al found that a He–V complex can only accommodate 12 H atoms. The 13th embedded H atom will fall into one of the TISs out of the vacancy, although the average H trapping energy is still lower than that of H in the TIS. On the other hand, for 'the sequential way', the overall trend of the H trapping energy also increases with the number of H atoms up to 12. Finally, the 13th H atom addition gives rise to a positive trapping energy as compared with that of H at the TIS in the bulk (figure 12). This suggests that the 13th H atom will prefer to occupy the TIS in the bulk rather than the He–V complex. Therefore, the maximal number of H atoms that can be trapped by the He–V complex is 12 via both 'the sequential way' and 'the simultaneous way'.

One He–V complex can thus hold up to 12 H atoms, independent of the trapping sequences [61, 157]. The calculations show furthermore that the shortest distance between H atoms is always much longer than that of the H₂ molecule (0.75 Å), which implies that no H₂ molecule forms with He at the vacancy centre. Interestingly, the configuration of H–He–V complex [61] indicates that H forms an icosahedron in the (HeV)H₁₂ complex, consisting of 12 isosceles and 8 equilateral triangles, as shown in figure 13. The distance between H and its 1NN is ∼1.74 Å, and that between H and its 2NN is ∼2.00 Å.

Zoom In Zoom Out Reset image size

7.2.4. He_nH_nV_n complexes.

The first-principles calculations based on DFT performed by Lee et al [62] show that He tends to move into the vacancy while H tends to stay around the TIS when a vacancy appears close to an He–H interstitial pair. Complexes consisting of He, H and vacancies were built, and the number of vacancies was set to be less than the total number of He and H atoms to create competition between the two species. The He and H atoms were initially positioned at the face centres, so that they have an equal probability to move towards the vacancies. The calculations show that all four He atoms move towards the vacancy cluster to form a four-atom cluster of a square shape, while all four H atoms stay away from the vacancies and relax slightly towards the W atoms (figure 14). Zhou et al also demonstrated that it is energetically more favourable that He atoms occupy the vacancy rather than H atoms, and the He atom is more likely to stay at the vacancy centre [61].

Zoom In Zoom Out Reset image size

7.2.5. Suppression of H bubble formation by He.

According to the first-principles calculations based on the DFT, a single C atom is shown to be energetically favourable sitting at the OIS in W with a solution energy of 0.78 eV [174]. It is found that the hybridization between C 2p and 1NN W 5d is more pronounced for C at the OIS based on the LDOS calculations, and a broad energy interval is obviously present near the Fermi level. Also, such a strong hybridized interaction is favourable to increase the formation of the directional C–W covalent-like bond, leading to the stability of C at the OIS. The OIS → TIS → OIS path has been demonstrated to be the optimal diffusion path for interstitial C. The energy barrier is calculated to be 1.46 eV, which suggests that it is very difficult for an interstitial C atom to migrate at a lower temperature.

Two C atoms tend to pair up at two neighbouring OISs along the 〈2 1 0〉 direction with an equilibrium distance of ∼3.57 Å [174]. Such a favourable pairing configuration appears to be mediated by a bended 'bridging' W atom. The binding energy of the C–C pair is ∼0.50 eV, which indicates an attractive interaction between C atoms in W. This should be responsible for the higher local C concentration and further the formation of carbide in W.

When a monovacancy is present, the most stable site for C was found to be at a position of ∼1.31 Å off the vacancy centre close to an OIS, with a solution energy of −1.15 eV, much lower than that in the intrinsic W [175, 176]. Such a strong trapping of C by the vacancy implies that vacancies serve as trapping centres and drive large amount of C atoms to segregate, forming C_nV complexes. Liu et al investigated the trapping of multiple C atoms in the vacancy and found that the trapping energy per C atom has a tendency to increase with the number of C atoms, except for a minor decrease when the second C atom is introduced.

Zhou et al investigated the energetics and site preference of C atoms in a W Σ5(3 1 0)/[0 0 1] GB [177]. The C atom prefers to sit in an interstitial site rather than a substitutional site in the GB. The segregation energy difference of C between the GB and the free surface is positive (+0.19 eV per C atom), indicating that C enhances the GB cohesion according to the Rice–Wang thermodynamic theory [109].

C impurity plays an important role in the blister formation. Experimentally, Ueda et al [178] indicate that a very small amount of C (∼0.3% or more) significantly enhances blister formation in W, which has been explained by the formation of a carbide layer that prevents H from leaving W. The first-principles calculations indicate that at a short distance the interaction between C and H is repulsive [179, 180], and it becomes negligible (∼0.05 eV) at a distance of ∼2.5 Å [179]. Such repulsion increases the diffusion coefficient of C in pure W [180], and enhances the migration barrier of H to travel through the WC layers, which prevents the implanted H from escaping the W matrix, contributing to the blister formation in W [181]. Lasa et al calculated the D retention in W carbides (WC and W₂C) using the empirical potentials [182], from which the D retention is shown to increase linearly with the C content, and the penetration of D is prevented by amorphization. In the W carbides, the D atoms are trapped around the C atoms, which leads to the formation of D₂ as well as different hydrocarbons. This is contrary to pure W, in which atomic D is present mainly.

Kong et al further investigated multiple H accumulation around C in pure W and their influence on the stability of nearby lattice sites [183]. It is found that the vacancy formation energies at the nearest-neighbour sites of C and H are lower than that in the intrinsic bulk, which implies that the introduction of C and H at the interstitial site decreases the stability of their nearby lattice sites. This may be attributed to the lattice distortion and charge redistribution induced by the interstitial atoms.

The multiple H atoms are shown to accumulate around the C atom to form C–H_n complexes [183]. When the C atom is trapped at the vacancy, no more than 9 H atoms can stay within the vacancy and the H atoms cannot bind together to form an H₂ molecule [179, 183]. The strong trapping of a single H by the vacancy remains almost unchanged in the vicinity of a C atom. However, when the number of H atoms increases, the C effect is enhanced due to the C–H interaction. The trapping energy of H increases with the presence of C, making the H_n–V complex less stable than in the C-free case. It will be thus easier for H to escape from the vacancy in the presence of a C atom.

A single O atom prefers to sit at the TIS with a solution energy of −1.75 eV according to the first-principles calculations [184]. Two O atoms tend to be paired up at the two neighbouring TISs along the 〈3 1 0〉 direction with a large binding energy of 1.60 eV. The equilibrium distance of the O–O pair is ∼2.28 Å, much longer than the bond length in the O₂ molecule, indicating that O atoms cannot bind together to form O₂ molecules directly in intrinsic W. However, because of this large binding energy, two O atoms can easily attract each other, leading to the formation of O clusters at a high concentration of O impurities. The preferable diffusion path for O is to migrate directly between two nearest neighbouring TISs with a diffusion barrier of 0.17 eV. The diffusion coefficient of O in W is calculated to be 1.50 × 10⁻⁹ m² s⁻¹ at a temperature of 500 K.

When a monovacancy is present, O will sit 1.28 Å off the vacancy centre close to an OIS with a strong trapping energy of −3.05 eV [183, 184]. This indicates a strong interaction between O and vacancy in W. The extra trapping energy of O with the vacancy suggests that the vacancy formation energy will be lower when more O atoms are trapped in the vacancy. Multi O atoms are repulsive at the vacancy, making the vacancy hold three O atoms without the formation of O₂ or O₃ molecule. It is suggested according to the bond lengths and the LDOS results that O and neighbouring W atoms form stronger W–O bonds similar to the WO₃ compounds.

Recent experimental studies reported that O can act as a trapping site, making large number of H isotope ions to aggregate and further form H blisters [185, 186]. DFT calculations predict a positive binding energy for the O–H pair (∼0.31 eV), suggesting that O has a potential to trap H atoms [183]. The maximum number of accumulated H atoms is 5 surrounding an O atom. The binding energy of H atoms to an O–vacancy complex decreases when the number of H atoms increases, as shown in figure 15. The H binding energy to the O–vacancy complex decreases slightly in comparison with the O-free case. Furthermore, the vacancy formation energies of some lattice sites near the O–H_n complex are found to be negative, suggesting that a vacancy can form spontaneously. Similar case occurs for the C as well (figure 15). The formation energy of the complex is so low that abundant O–H_n/C–H_n complexes may exist and widely distribute in W, and the accumulation of these abundant complexes can serve as the origins for the formation of H bubbles.

Zoom In Zoom Out Reset image size

W is not really a practical engineering material because its machining and welding are generally difficult. To make up for this problem, one well-known solution is the addition of alloying element rhenium (Re), which increases recrystallization temperature and ductility of W. Furthermore, under fusion irradiation conditions, transmutation of W atoms by 14 MeV neutron irradiation can lead to the formation of Re. In certain applications and after long service, as many as 25% of the W atoms can be transmuted [187]. For this reason, the W–Re alloys have been extensively investigated experimentally. Other elements such as Ta and V are, for the same purpose, also under investigation. Muzyk et al [188] indicate that the addition of V and Ta can result in an improved ductility of W and that V atoms are expected to trap SIA defects in dilute W–V alloys. The DFT calculations performed by Li et al [189] indicate that the 1/2 〈1 1 1〉 dislocation slip plane is altered by Re alloying.

Concerning the interactions of solute atoms and H, Golubeva et al [190] investigated D retention in Re-doped W by means of the thermodesorption technique. In the range 1–10% Re, the Re concentration in the W–Re alloy will not influence noticeably the D retention. The DFT calculations by Heinola et al [88] indicate that there is a negligible interaction between Re/Os and H when the solute atoms are in substitutional sites, but such an interaction becomes much more significant when the solute atoms are in interstitial sites. Figure 16 compares the interaction of a series of transition elements (substitution in W) with H, He and a single vacancy predicted by the first principles based on the DFT. For the elements of the first two rows, the interaction of solute with He atoms is always stronger than those with the vacancy or an H atom. The binding energy of the solutes investigated is always positive and thus will favour H solution but by no more than 0.5 eV except for Fe, Ni and Cr which limits the range of temperatures at which they can be effective traps.

Zoom In Zoom Out Reset image size