Sage Journals: Discover world-class research

Abstract

Hydrogen diffusion and trapping in ferrite is evaluated by quantum mechanically informed kinetic Monte Carlo simulations in defective microstructures. We find that the lattice diffusivity is attenuated by two to four orders of magnitude due to the presence of dislocations. We also find that pipe diffusivity is vanishingly small along screw dislocations and demonstrate that dislocations do not provide fast diffusion pathways for hydrogen as is sometimes supposed. We make contact between our simulations and the predictions of Oriani's theory of ‘effective diffusivity’, and find that local equilibrium is maintained between lattice and trap sites. We also find that the predicted effective diffusivity is in agreement with our simulated results in cases where the distribution of traps is spatially homogeneous; in the trapping of hydrogen by dislocations where this condition is not met, the Oriani effective diffusivity is in agreement with the simulations to within a factor of two.

This paper is part of a thematic issue on Hydrogen in Metallic Alloys

Keywords

Simulation kinetic Monte Carlo hydrogen diffusion trapping local equilibrium

1 Introduction

Hydrogen can be dissolved in most metals and alloys [1], although it has remarkably low solubility in body-centred cubic (bcc) iron. The interaction of hydrogen with the microstructure is one reason for embrittlement and premature failure in high-strength steels and other engineering metals including nickel, titanium and zirconium alloys [2]. The deleterious effects of hydrogen are frequently connected to its rapid diffusion through the crystal lattice, in particular in bcc and martensitic steels. The interactions of hydrogen with crystal defects and their consequences are fundamentally less well understood than diffusion of hydrogen in the perfect crystal lattice, despite generally dominating the influence of hydrogen in metals.

Hydrogen atoms in the lattice are frequently divided into two categories: diffusing hydrogen, which can move freely through the normal interstitial sites, and trapped or non-diffusing hydrogen, residing around various crystal imperfections (vacancies, dislocations, grain boundaries and so on) [3 –5]. The driving force for diffusion is the chemical potential gradient generated by hydrogen gradients combined with hydrostatic stress gradients. If fast trapping and detrapping kinetics are assumed, the freely diffusing and trapped hydrogen concentrations need to be in equilibrium. Based on the McLean isotherm, Oriani's [6] fractional occupancy of a trap site is connected to the occupancy of normal interstitial sites by (1) where is the corresponding trap binding energy, k is the Boltzmann constant and T is the absolute temperature. If and represent the number of normal interstitial sites and number of traps per unit volume, then the volume concentrations of hydrogen in normal interstitial sites are and in traps . Using Equation (1), Oriani [6] finds an effective diffusivity of hydrogen given by (2) in which is the hydrogen diffusivity in the perfect bcc lattice.

In this work the validity of Equations (1) and (2) is investigated by comparison to quantum mechanically informed kinetic Monte Carlo (kMC) simulations. Using the kMC, we may calculate the diffusivity of hydrogen in the perfect lattice, . We can then create blocks of crystal containing defects which act as trap sites, and from the evolution of the simulation we can extract an effective diffusivity, . From the kMC simulations, the evolution of the occupancy of normal interstitial and trapping sites in pure bcc Fe is extracted and used as inputs to Equations (1) and (2) to calculate corresponding trap barriers and effective diffusivity, , in steady state, which are in turn compared for consistency with the input trap depths and output effective diffusivity, , to the kMC. Simulations of a perfect lattice are used to benchmark the resulting effective diffusivity in the lattice with defects (vacancy, edge or screw dislocation).

2 Diffusion model and parametrisation

Diffusion of hydrogen atoms through normal interstitial sites is governed by a random walk between occupied and empty sites. Under the framework of a kinetic Monte Carlo (kMC) method (cf. [7]), the rate probability for a jump from site A to B is given as

(3) where

is the attempt frequency and

is the energy barrier between metastable position A through the saddle point to the position B, dependent on the local atomic configuration. For a given local configuration, a process k is selected according to

(4) where

is a uniform random number between 0 and 1, and n is the number of all possible events, i, during one transition. The kMC simulates Poisson processes and the simulation time

of each event is evaluated as [8 –10]

(5) where

is a second random variable in

, and directly corresponds to physical time t when only interstitial diffusion is considered. The time-dependent diffusion coefficient is given by the Einstein expression, [11, 12] as

(6) where d is the system dimensionality (d=1,2,3) and

is the average squared displacement of particles in time t, [13]. In this study, all simulations were done at a constant temperature of 300 K, using a box with

bcc Fe unit cells, unless stated otherwise. Simulations were performed with full periodic boundary conditions in all directions. Fe atoms were not moved throughout the simulations and the permitted events were jumps between tetrahedral (T) sites, from T sites to trapping sites, vice versa and between trapping sites. The parameters used in the kMC model were for a vacancy obtained by using the quantum mechanical tight-binding approximation and molecular dynamics methods based around the Feynman path integral formulation of the quantum partition function [14, 15]; or taken from the literature for the hydrogen interaction with

edge and screw [16] dislocations on

slip planes, also determined within the framework of path integral molecular dynamics at 300 K. The use of

bcc Fe unit cells results in very large defect concentrations in the simulation box; therefore, hydrogen concentrations were increased to obtain a realistic ratio between defect and interstitial concentrations. Simulations were performed with 5–30 hydrogen atoms and a vacancy, edge or screw dislocation, and a combination of edge and screw dislocations in the simulation box. The corresponding concentrations of interstitial hydrogen in atomic and weight parts per million are given in Table 1. To demonstrate the effect of simulation box size, test runs were also performed using

and

bcc Fe unit cells and compared with results obtained with

bcc Fe unit cells.

Table 1

Concentrations of interstitial hydrogen used in simulations.

Number of hydrogen atoms	Concentration (appm)	Concentration (wppm)
5	312.4	5.6
10	624.6	11.2
15	936.6	16.8
20	1248.4	22.4
25	1560.1	27.9
30	1871.5	33.5

2.1 Parameters for diffusion around vacancy

For simulations on the movement of interstitial hydrogen, the energy landscape around a defect is needed. Energies are then used as inputs in Equation (3) to calculate the rate of all possible events. In Table 2, energies and calculated rate constants needed for interstitial movement around a vacancy at 300 K are given. Several events are considered: TT is a jump between adjacent tetrahedral sites, TTV1 is a jump from tetrahedral site to an unoccupied vacancy trap, TTV2 is a jump from a tetrahedral site to the vacancy trap if there is already a hydrogen atom occupying one trapping site around the vacancy, TVV is a jump between adjacent trap positions around the vacancy and TVT is used for detrapping from the vacancy. According to quantum calculations up to six hydrogen atoms can be trapped around a vacancy; however, only two have high probability of being found there [14, 15, 17]. Therefore, in the model only up to two hydrogen atoms are allowed to be trapped around the vacancy to increase the computational efficiency of the model without sacrificing significant reality. Figure 1 shows a schematic representation of the events we consider in the kMC simulations of the vacancy microstructures.

Figure 1

Schematic representation of the events considered in our kMC model shown in the bcc lattice, where the vacancy (V) is presented as a dotted cube with six possible binding points around its centre, T is a tetrahedral interstitial site and Fe is an iron position.

Table 2

Energies and rate coefficients used to construct the kMC model for trapping around a vacancy at 300 K and resulting rate coefficients.

Parameter	Energy (eV)	Rate coefficient (Hz)
TT	−0.045
TTV1	−0.033
TTV2	−0.013
TVV	−0.23
TVT	−0.65	63

Notes: TT is a jump between adjacent tetrahedral sites, TTV1 is a jump from tetrahedral site to an unoccupied vacancy binding point or a vacancy trap, TTV2 is a jump from tetrahedral site to a vacancy trap if there is already a hydrogen atom occupying one trapping site around the vacancy, TVV is a jump between adjacent traps around a vacancy and TVT is used for detrapping from a vacancy.

In all the calculations performed at 300 K, we use an activation barrier for lattice diffusion of 0.045 eV, calculated using zero-point energy-corrected density functional theory [18]; and an Arhenius prefactor of.

2.2 Parameters for diffusion around edge and screw dislocation

The energetic landscape of an interstitial hydrogen around edge or screw dislocations has many more local minima in comparison to the vacancy. As a remote hydrogen atom approaches the core, it encounters several metastable binding sites distributed between two parallel adjacent crystal planes which are perpendicular to the dislocation line sense [16]. In particular, the cores of edge and screw dislocations are reached through seven and six deep metastable binding sites, respectively. Movements from one metastable binding site to another and energy barriers for that transition were determined in [16] at 300 K using path integral molecular dynamics. However, these authors did not address transitions from bulk tetrahedral to metastable binding sites. The saddle point energy between a tetrahedral interstitial (T) site and a metastable binding site is assumed here to be similar to values obtained for transitions from a tetrahedral interstitial site to metastable binding site next to a vacancy [19]. In the case of edge dislocations, the present kMC model does not account for the change of binding energies and energy barriers between the metastable sites in compressive and tensile areas on both sides of the planar core. Instead, we assume that each metastable binding site is surrounded by four tetrahedral sites with energies for transition between T site to metastable binding points given in Tables 3 and 4 for edge and screw dislocation, respectively. TB and TZ represent transitions from four nearest T sites to metastable binding sites close to the dislocation, while BT and ZT represent the reverse. Transitions are grouped into TB/BT and TZ/ZT based on the distance between the T sites and metastable binding sites close to the dislocation (cf. Figure 2). A larger number labelling of a binding site corresponds to a position closer to the dislocation core, where, for an edge dislocation, metastable point 7 is the closest to the dislocation line and, in the case of a screw dislocation, metastable point 6 is the closest to the core. As in the case of vacancy, at 300 K we use an activation barrier for lattice diffusion of 0.045 eV and an Arhenius prefactor of Hz.

Figure 2

Schematic representation of the TB/BT and TZ/ZT transitions based on the distance between four nearest tetrahedral interstitial sites T and metastable binding site B.

Table 3

Energies for transitions to and from the metastable binding site used in constructed kMC model for edge dislocation at 300 K.

	Energy (eV)
Binding site	TB1	TB2	BT1	BT2	TZ1	TZ2	ZT1	ZT2
1	−0.41	−0.02	−0.52	−0.17	−0.41	−0.02	−0.52	−0.17
2	−0.41	−0.01	−0.43	−0.2	−0.41	−0.01	−0.43	−0.2
3	−0.22	−0.015	−0.25	−0.125	−0.22	−0.15	−0.26	−0.125
4	−0.12	−0.02	−0.19	−0.09	−0.12	−0.02	−0.19	−0.09
5	−0.085	−0.025	−0.14	−0.07	−0.085	−0.025	−0.14	−0.07
6	−0.07	−0.045	−0.08	−0.045	−0.07	−0.045	−0.08	−0.045
7	−0.02	−0.02	−0.35	−0.065	−0.02	−0.02	−0.35	−0.065

Table 4

Energies for transitions to and from a metastable binding site used in constructed kMC model for screw dislocation at 300 K.

	Energy (eV)
Binding site	TB1	TB2	BT1	BT2	TZ1	TZ2	ZT1	ZT2
1	−0.02	−0.025	−0.27	−0.36	−0.02	−0.13	−0.43	−0.32
2	−0.035	−0.01	−0.18	−0.14	−0.03	−0.01	−0.32	−0.14
3	−0.06	−0.015	−0.13	−0.075	−0.06	−0.15	−0.22	−0.075
4	−0.05	−0.02	−0.1	−0.065	−0.05	−0.02	−0.18	−0.065
5	−0.045	−0.045	−0.08	−0.07	−0.045	−0.045	−0.1	−0.07

Schematic representations of edge and screw dislocation energetic landscapes are shown in Figure 3. In the case of an edge dislocation, the core is depicted as an inverted T, while numbers represent metastable binding points for a hydrogen atom diffusing along the slip plane of the edge dislocation (Figure 3(a)). Figure 3(b) shows schematically energetic landscape and metastable points around a screw dislocation. The dotted lines serve as a guide to view threefold symmetry of the slip planes, and solid squares indicate the three deepest traps closest to the dislocation core. The arrows indicate those transitions of hydrogen atoms along the slip plane in the direction of the core that are included in our kMC simulations.

Figure 3

Schematic representation of an energetic landscape around (a) edge dislocation and (b) screw dislocation.

The activation energies for a proton to jump from a metastable binding site close to the core to another metastable binding site close to the core in the direction of the dislocation line are given in Table 5.

Table 5

Activation energies for pipe diffusion along dislocation core lines for edge and screw dislocations.

Dislocation	Site	Energy/eV
Edge		−0.32
		−0.13
Screw		−0.43
		–0.03

Note: Referring to Figure 3, and are barriers for the two-hop jump along the direction at the centre of the core indicated by the inverted T; and are barriers for the two-hop jump along the direction between sites labelled by the filled square. For details, see [16].

3 Results and discussion

3.1 Diffusivity in ideal lattice

To begin with, we show in Figure 4 the calculated lattice diffusivity, , in a perfect lattice at 300 K as a function of simulation time. (In this and subsequent figures we use ‘diffusivity’ loosely to mean the simulated diffusivity as a function of simulation time. But our calculated diffusivities are of course the steady-state limit of this quantity.) It is vital to note that simulation times of at least one microsecond are required before a convergent result is reached. This is an important observation since it emphasises the need to go beyond timescales accessible to molecular dynamics. We have used an activation barrier for lattice diffusion of 0.045 eV and an Arhenius prefactor of Hz at 300 K, which reproduces the diffusivity measured by Nagano et al. [20], namely . As shown, there is negligible effect of hydrogen concentration on the interactions between interstitial atoms despite the small simulation box, as our results are in the range from to.

Figure 4

Influence of hydrogen concentration (cf. Table 1) on kMC determined lattice diffusivity in an ideal lattice.

3.2 Diffusivity around point defect

In order to determine appropriate simulation box size in the presence of a trap, we performed simulations with different box sizes of , , , , and unit cells and a single point defect (vacancy) with 5 hydrogen atoms in the simulation box. Our results depicted in Figure 5(a) show that box size only has an effect on diffusivity at the outset of the simulation. Since very similar final results for diffusivity was obtained for all tested box sizes (cf. Figure 5(b)), the decision was made to use bcc Fe unit cells (5.74 nm with 16000 Fe atoms and 100800 tetrahedral sites) to keep the simulation runtime achievable and reasonable.

Figure 5

Results of simulations with single vacancy and five hydrogen atoms, (a) effect of simulations box size and (b) details showing final results of (a).

The effects of different initial hydrogen positions were evaluated. It was found that different initial positions of interstitial hydrogen have negligible effects on the diffusivity and time when the steady state is reached. Figure 6 shows evolution of effective diffusivity from kMC simulations for a single vacancy and different hydrogen concentrations. As shown, the time needed to reach steady state for a vacancy is about s, indicating that hydrogen diffusivity around point defects achieves steady state almost instantly. For point defects this is expected since the activation barrier for trapping is lower than that for lattice diffusion, while the barrier for escape from the trap is so large that escape is a very rare event on the timescale of our simulations.

Figure 6

Influence of hydrogen concentration on effective diffusivity in lattice with a single point defect.

A comparison between calculated diffusivity in the perfect lattice and calculated effective diffusivity with different vacancy or hydrogen concentrations determined from kMC simulations is shown in Figure 7. It can be seen that similar ratios between defect and hydrogen concentrations lead to comparable effective diffusivities.

Figure 7

Influence of number of point defects or hydrogen concentrations on effective diffusivity and comparison to ideal lattice diffusivity.

3.3 Diffusivity around dislocations

Simulations with dislocations were performed mainly using a box size of Fe bcc unit cells and different hydrogen concentrations. Results of effective diffusivity for a single screw, a single edge and a combination of screw and edge dislocations are shown in Figure 8. Depending on the defect type and hydrogen concentration, the time needed to reach steady state is for both dislocation types, around s. Results shown in Figure 8 indicate that hydrogen diffusivity around traps, simulated with kMC achieves steady state almost instantly. This discovery is important for mesoscopic and continuum simulations of crack propagation where the redistribution of hydrogen during crack growth plays an important role.

Figure 8

Influence of hydrogen concentration on kMC determined effective diffusivity in lattice with (a) single edge dislocation; (b) single screw dislocation; (c) one edge and one screw dislocation.

During the simulations, we also followed diffusion along dislocation cores (pipe diffusion) and results show that pipe diffusivity is several orders of magnitude lower compared to the effective lattice diffusivity. The resulting diffusivity in the case when hydrogen atoms are trapped inside the most stable points close to the dislocation core and move along the

direction of the dislocation line, clearly indicate that the diffusivity resulting from hydrogen transport along the dislocation line is several orders of magnitude lower compared to the effective lattice diffusivity. It follows that dislocation lines do not act as fast pathways for hydrogen pipe diffusion in bcc iron. In Table 6 are given results obtained for various hydrogen concentrations in a simulation box with a single edge or screw dislocations, where

is the effective diffusivity and

is the diffusivity along the

direction within the dislocation core. Figure 9 shows the evolution of the diffusivity along a screw dislocation core, showing that the diffusivity along the screw dislocation line is more than order of magnitude smaller than the effective diffusivity of the dislocated lattice. Pipe diffusion along the screw dislocation

direction is several orders of magnitude smaller than the lattice diffusivity.

Table 6

Results of kMC simulations for single edge and screw dislocation.

Defect	(wppm)	(m² s⁻¹)	(m²s⁻¹)
1 Edge	5.6	0.86	10.7
	11.2	0.90	11.3
	16.8	0.94	10.4
	22.4	0.98	18.3
	27.9	1.02	29.1
	33.5	1.05	20.9
1 Screw	5.6	0.030	0.0445
	11.2	0.031	0.0415
	16.8	0.031	0.0441
	22.4	0.032	0.0509
	27.9	0.032	0.05
	33.5	0.034	0.0531

Note: Comparison of effective diffusivity, , pipe diffusion, , for various hydrogen concentrations, .

Figure 9

Pipe diffusivity along screw dislocation core () determined from kMC simulations.

From this we can conclude that the screw dislocation core is not a pathway for fast diffusion of hydrogen. The origin of this is twofold: first the activation barrier is large, Table 5, and second the trap occupancy is relatively large leading to blocking of the diffusion path.

In Figure 10 a comparison between calculated effective diffusivity for several cases with different dislocation or hydrogen concentrations is shown. It can be seen that similar ratios between defects and hydrogen concentration lead to comparable effective diffusivities.

Figure 10

Influence of defect and hydrogen concentration on effective diffusivity for dislocations.

3.4 Comparison between kMC and Oriani's trapping theory

To test the validity of Oriani's assumption (cf. Equation (1), [6]) of equilibrium between freely diffusing and trapped hydrogen concentrations, the kMC model was developed in a way that trap and lattice occupancies could be followed. The total time hydrogen atoms spent at each site can be calculated at any given time during the simulation as the time spent in the trap or the lattice. Fractional occupancy of sites is obtained by dividing the time hydrogen atoms spent in a site by the total time. Since a kMC simulation aims to reach thermodynamic steady state, values at the end of simulations were used to determine trap binding energies of defects from resulting data for occupancies between lattice hydrogen and trapped hydrogen using equation (1) and to test the validity of Oriani's theory of local equilibrium by comparing these binding energies with those used as input. After suitable averaging of the trap occupancies, the diffusivity obtained from kMC is compared to the effective diffusivity calculated with Equation (2) using occupancies from kMC simulations. Detailed results and calculated trap binding energies for different defect and hydrogen concentrations are presented in Tables 7–9. Results for the vacancy in Table 7 show a minimal effect of concentration in the cases where the number of hydrogen atoms is smaller than the number of trapping sites. The average value for vacancy trap binding energy calculated using equation (1) is determined to be eV. In the case of a single type of dislocation (cf. Table 8), the number of possible trapping sites far exceeds the number of hydrogen atoms in the box. The trap binding energy determined for an edge dislocation is eV and is independent of the number of dislocations in the box. The trap binding energy of the screw dislocation is slightly dependent on the number of dislocations in the box and is eV for a single screw dislocation in the box and eV in the case of simulations with two screw dislocations. This increase occurs due to overlap between both dislocations. Oriani's effective diffusivity, , calculated using equation (2) for vacancies using occupancies from kMC simulations is presented in Table 7. Comparison between and shown in Figure 11 reveals that the diffusivity obtained from kMC simulations is slightly lower than that obtained using equation (2). The largest differences are in the case of lowest hydrogen concentrations, and the least differences when the hydrogen concentration is highest. This can be explained by the small difference between the number of hydrogen atoms that can freely move compared to possible trap sites. In the cases where the number of traps exceeds the number of hydrogen atoms, we get very similar results when comparing both approaches. From our results, we may show that Oriani's effective diffusivity is a very good approximation for the description of effective diffusivity of hydrogen in the case of microstructures with a single type of trap uniformly distributed, for example, the case of vacancies.

Figure 11

Results of obtained diffusivities for point defect from kMC simulations, (closed symbols) and calculated using Oriani's equation, (open symbols).

Table 7

Results of kMC simulations for various vacancy and hydrogen concentrations.

Vac. /num.	(wppm)	(m² s⁻¹)	(−)	(−)	(eV)	Avg (eV)	(m² s⁻¹)
1	5.6	5.26	2.98	0.333	0.252	0.21±0.02	6.07
	11.2	7.01	7.94	0.333	0.226		7.51
	16.8	7.60	12.90	0.333	0.214		7.95
	22.4	7.89	17.86	0.333	0.205		8.16
	27.9	8.06	22.82	0.333	0.199		8.28
	33.5	8.18	27.78	0.333	0.194		8.37
2	5.6	1.75	0.99	0.333	0.280	0.22±0.03	2.39
	11.2	5.26	5.95	0.333	0.234		6.07
	16.8	6.43	10.91	0.333	0.218		7.05
	22.4	7.01	15.87	0.333	0.208		7.51
	27.9	7.36	20.83	0.333	0.201		7.78
	33.5	7.59	25.79	0.333	0.196		7.95
3	5.6	0.00	0.99	0.278	0.273	0.23±0.03	0.00006
	11.2	3.50	3.97	0.333	0.244		4.38
	16.8	5.26	8.93	0.333	0.223		6.07
	22.4	6.13	13.89	0.333	0.212		6.82
	27.9	6.66	18.85	0.333	0.204		7.24
	33.5	7.01	23.81	0.333	0.198		7.51
5	5.6	0.00045	0.99	0.167	0.256	0.23±0.02	0.00
	11.2	0.00002	1.98	0.300	0.258		1.2
	16.8	2.92	4.96	0.333	0.238		3.76
	22.4	4.38	9.92	0.333	0.220		5.26
	27.9	5.26	14.88	0.333	0.210		6.07
	33.5	5.84	19.84	0.333	0.202		6.57

Note: Vac. is the number of vacancies, is hydrogen concentration, is the effective diffusivity, is the lattice occupancy, is the trap occupancy and is the trap binding energy determined from assumption of local equilibrium using Equation (1) and is the effective diffusivity calculated from Equation (2).

Table 8

Results of kMC simulations for dislocations and hydrogen concentrations.

Defect	(wppm)	()	(−)	(−)	(eV)	Avg (eV)	()
1 Edge	5.6	0.86	6.89	29.75	0.200	0.21±0.01	0.61
	11.2	0.90	4.44	15.05	0.207		0.53
	16.8	0.94	3.13	8.27	0.213		0.42
	22.4	0.98	2.17	7.4	0.206		0.53
	27.9	1.02	1.65	3.82	0.216		0.37
	33.5	1.05	1.11	2.82	0.214		0.42
1 Screw	5.6	0.030	4.1	17.7	0.278	0.28±0.002	0.01
	11.2	0.031	4.63	27.4	0.282		0.086
	16.8	0.031	4.7	33.5	0.281		0.088
	22.4	0.032	4.78	40.9	0.282		0.085
	27.9	0.032	4.58	40.2	0.279		0.097
	33.5	0.034	4.92	43.5	0.279		0.098
2 Edge	5.6	0.96	673.5	290.62	0.202	0.21±0.007	0.37
	11.2	0.98	421.3	147.34	0.207		0.29
	16.8	1.0	306	78.36	0.215		0.21
	22.4	1.02	197.8	66.68	0.207		0.28
	27.9	1.04	160.3	37.21	0.217		0.19
	33.5	1.07	105.3	26.9	0.214		0.21
2 Screw	5.6	0.0038	3699	24.13	0.321	0.34±0.01	0.0038
	11.2	0.0035	6467	43.50	0.335		0.007
	16.8	0.0037	7676	50.70	0.346		0.01
	22.4	0.0039	8231	53.21	0.354		0.013
	27.9	0.0042	8012	55.63	0.349		0.013
	33.5	0.0046	7500	56.63	0.341		0.011

Note: is the hydrogen concentration, is the effective diffusivity, is lattice occupancy, is the trap occupancy and is the trap binding energy determined from assumption of local equilibrium using Equation (1).

Table 9

Results of kMC simulations for dislocations and hydrogen concentrations.

Defect	(wppm)		(−)	(−)	(eV)	Avg (eV)
1 Edge + 1 Screw	5.6	0.31	1746	3.34	0.221	0.22 ± 0.003
	11.2	0.31	815	1.83	0.217
	16.8	0.31	598	1.47	0.215
	22.4	0.31	548	1.17	0.218
	27.9	0.32	449	0.81	0.223
	33.5	0.32	365	0.76	0.219

Note: is the hydrogen concentration, is the effective diffusivity, is the lattice occupancy, is the trap occupancy and is the trap binding energy determined from assumption of local equilibrium using Equation (1).

Results from kMC simulations with edge and screw dislocations in the box are presented in Table 9. The trap binding energy for combination of both types of dislocations in kMC simulation is eV and is close to the simulations where only an edge dislocation is present. Comparison of effective diffusivity shows that though vacancies can bind hydrogen to their nearest tetrahedral sites, the number of these sites is too small to have any large effect unless the number of possible trapping sites exceeds the number of hydrogen atoms. In the case of dislocations effective diffusivity is significantly lower than that in the ideal lattice or even in a lattice with only vacancies. Moreover, results show that the effective diffusivity in the presence of screw dislocations is a couple of decades smaller compared to simulations with edge dislocations or combinations of defects. This is due to deeper trap sites around screw dislocations and interactions between them if several screw dislocations are present in the box. Resulting values for the binding energies for vacancy and edge or screw dislocations are in good agreement with published results from the literature [21, 22]. Figure 12 shows comparison between calculated diffusivity from kMC simulations and determined effective diffusivity using Oriani's theory for dislocations. In the case of a screw dislocation, Oriani's effective diffusivity calculated with inputs for occupancies at the end of kMC simulations is slightly higher compared to diffusivity from kMC simulation (cf. Figure 12(a)). In the case of an edge dislocation and combination of edge and screw dislocations depicted in Figure 12(b), diffusivity obtained from kMC simulations is higher than calculated using equation (2). Larger discrepancies between diffusivities obtained from kMC simulations and calculated Oriani effective diffusivity, , using data for occupancies from kMC simulations can be explained by the spatial inhomogeneity of the energy landscape due to trap sites.

Figure 12

Results of obtained diffusivities from kMC simulations, , (closed symbols) and calculated using Oriani's equation, , (open symbols) for (a) single screw dislocation and (b) single edge dislocation and one edge and one screw dislocation.

4 Conclusions

A kMC model for diffusion of hydrogen has been developed. This is able to simulate diffusion in ideal lattices and with vacancies and edge or screw dislocations with Burgers vector. The following conclusions can be reached from the present work.

Diffusivity in the ideal lattice is comparable to experimental results.

Direct comparison can be made of explicit simulations and the theory of Oriani which assumes local equilibrium and one type of trap.

The kMC confirms the establishment of local equilibrium in all the cases we have looked at.

The Oriani effective diffusivity, , forms a good approximation to that calculated by kMC, . The best agreement is in the case of a homogeneous distribution of vacancies. For the case of dislocations, because these are spatially inhomogeneously distributed and there is a range of trap types with different binding energies, the Oriani effective diffusivity may differ by a factor of two from the kMC value.

As a result of kMC simulations of hydrogen transport in the proximity of a screw dislocation, we can conclude as follows:

Hydrogen atoms are confined to the region around the dislocation line, where they jump between the adjacent binding and low-energy metastable sites.

Hydrogen atoms occupy predominantly the deepest trap sites close to the dislocation line.

Hydrogen diffusivity significantly decreases with increasing dislocation density.

Dislocation pipe diffusion of hydrogen results in significantly lower diffusivity compared to lattice diffusion. We calculated the jump diffusion coefficient for hydrogen diffusion in the cases when hydrogen atoms are trapped inside the core and move in the direction along the dislocation line. Our kMC simulations clearly indicate that the diffusivity resulting from hydrogen transport along the dislocation line is several orders of magnitude lower than that resulting from the lattice diffusion. Thus, the dislocation lines do not act as fast pathways for hydrogen diffusion in bcc Fe. The two reasons for this are the large activation barrier and the near saturation of trap sites leading to a blocking of the diffusion path.

Similar simulations of the hydrogen diffusivity in proximity of edge dislocations show the following.

The diffusivity near edge dislocations is higher compared to that in the vicinity of screw dislocations. Apart from the higher number of trap sites, diffusivity near a screw dislocation is also affected by its non-planar core structure. The planar core of the edge dislocation allows hydrogen atoms to move in directions parallel to the lines of trap sites, which increases their pipe diffusivity. We repeat that, in these preliminary simulations, we have neglected possible attraction and repulsion of hydrogen by the tensile and compressive strain fields above and below the slip plane of an edge dislocation.

There are no dominant trap sites. Hydrogen atoms are confined more or less evenly in all binding sites and jump between trap sites and the adjacent low-energy metastable sites.

There is no low-energy pathway for hydrogen diffusion in the direction along the dislocation line. We used our kMC model for screw and edge dislocations to simulate diffusion of ten hydrogen atoms. Hydrogen diffusivity in this case is higher than the diffusivity around the same number of screw dislocations and lower than the diffusivity around edge dislocations. The reason for this is the higher hydrogen diffusivity in the proximity of edge dislocations.

Footnotes

Disclosure statement

No potential conflict of interest was reported by the author(s).

ORCID

D. Bombac

A. T. Paxton

References

Graham

. On the occlusion of hydrogen gas by metals. Proc R Soc London. 1867–1868;16:422–427. doi: 10.1098/rspl.1867.0090

Archakov

, Grebeshkova

. Nature of hydrogen embrittlement of steel. Metal Sci Heat Treat. 1985;27:555–562. doi: 10.1007/BF00699349

Pressouyre

. A classification of hydrogen traps in steel. Metall Trans A. 1979;10:1571–1573. doi: 10.1007/BF02812023

Pressouyre

. Hydrogen traps, repellers, and obstacles in steel; consequences on hydrogen diffusion, solubility, and embrittlement. Metall Trans A. 1983;14:2189–2193. doi: 10.1007/BF02662391

Myers

, Baskes

, Birnbaum

, . Hydrogen interactions with defects in crystalline solids. Rev Modern Phys. 1992;64:559–617. doi: 10.1103/RevModPhys.64.559

Oriani

. The diffusion and trapping of hydrogen in steel. Acta Metall. 1970;18:147–157. doi: 10.1016/0001-6160(70)90078-7

Bombac

, Kugler

. Influence of diffusion asymmetry on kinetic pathways in binary fe-cu alloy: a kinetic Monte Carlo study. J Mater Eng Perform. 2015;24:2382–2389. doi: 10.1007/s11665-015-1493-2

Young

, Elcock

. Monte Carlo studies of vacancy migration in binary ordered alloys: I. Proc Phys Soc. 1966;89:735–746. doi: 10.1088/0370-1328/89/3/329

Bortz

, Kalos

, Lebowitz

. A new algorithm for Monte Carlo simulation of Ising spin systems. J Comput Phys. 1975;17:10–18. doi: 10.1016/0021-9991(75)90060-1

10.

Gillespie

. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976;22:403–434. doi: 10.1016/0021-9991(76)90041-3

11.

Einstein

. On the movement of small particles suspended in a stationary liquid demanded by the molecular-kinetic theory of heat. New York: Dover Publications; 1956.

12.

Einstein

. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der Physik. 1905;322:549–560. doi: 10.1002/andp.19053220806

13.

Limoge

, Bocquet

. Temperature behavior of tracer diffusion in amorphous materials: A random-walk approach. Phys Rev Lett. 1990;65:60–63. doi: 10.1103/PhysRevLett.65.60

14.

Paxton

, Katzarov

. Quantum and isotope effects on hydrogen diffusion, trapping and escape in iron. Acta Mater. 2016;103:71–76. Available from: http://www.sciencedirect.com/science/article/pii/S1359645415007430 doi: 10.1016/j.actamat.2015.09.054

15.

Katzarov

, Pashov

, Paxton

. Fully quantum mechanical calculation of the diffusivity of hydrogen in iron using the tight-binding approximation and path integral theory. Phys Rev B. 2013;88:054107-1–054107-10. doi: 10.1103/PhysRevB.88.054107

16.

Kimizuka

, Ogata

. Slow diffusion of hydrogen at a screw dislocation core in

-iron. Phys Rev B. 2011;84:024116-1–024116-6. doi: 10.1103/PhysRevB.84.024116

17.

Tateyama

, Ohno

. Stability and clusterization of hydrogen-vacancy complexes in

an ab initio study. Phys Rev B. 2003;67:174105. Available from: http://link.aps.org/doi/10.1103/PhysRevB.67.174105 doi: 10.1103/PhysRevB.67.174105

18.

Ramasubramaniam

, Itakura

, Ortiz

, et al. Effect of atomic scale plasticity on hydrogen diffusion in iron: Quantum mechanically informed and on-the-fly kinetic Monte Carlo simulations. J Mat Res. 2008;23:2757–2773. Available from: https://www.cambridge.org/core/article/effect-of-atomic-scale-plasticity-on-hydrogen-diffusion-in-iron-quantum-mechanically-informed-and-on-the-fly-kinetic-monte-carlo-simulations/E1187E3F2179A5B47C828B4BAB99A522https://www.cambridge.org/core/article/effect-of-atomic-scale-plasticity-on-hydrogen-diffusion-in-iron-quantum-mechanically-informed-and-on-the-fly-kinetic-monte-carlo-simulations/E1187E3F2179A5B47C828B4BAB99A522 doi: 10.1557/JMR.2008.0340

19.

Paxton

. From quantum mechanics to physical metallurgy of steels. Mat Sci Technol. 2014;30:1063–1070. Available from: http://www.tandfonline.com/doi/full/10.1179/1743284714Y.0000000521 doi: 10.1179/1743284714Y.0000000521

20.

Nagano

, Hayashi

, Ohtani

, et al. Hydrogen diffusivity in high purity alpha iron. Scr Metall. 1982;16:973–976. Available from: http://www.sciencedirect.com/science/article/pii/0036974882901363 doi: 10.1016/0036-9748(82)90136-3

21.

Kirchheim

. Hydrogen solubility and diffusivity in defective and amorphous metals. Progress Mater Sci. 1988;32:261–325. Available from: http://www.sciencedirect.com/science/article/pii/0079642588900102 doi: 10.1016/0079-6425(88)90010-2

22.

Stopher

, Rivera-Diaz-delCastillo

PEJ

. Hydrogen embrittlement in bearing steels. Mat Sci Technol. 2016;32:1–10.