Online ISSN : 1347-5320
Print ISSN : 1345-9678
ISSN-L : 1345-9678
Mechanics of Materials
Evaluating Solubility and Diffusion Coefficient of Hydrogen in Martensitic Steel Using Computational Mechanics
Kazuhito FukudaAkihiro TojoRyosuke Matsumoto
Author information

2020 Volume 61 Issue 7 Pages 1287-1293


The applications of high-strength martensitic steels in engineering are limited due to concerns regarding delayed fracture. Recent experimental and theoretical studies revealed that the presence of mobile hydrogen moving toward cleavage crack surfaces in these steels plays an important role in delayed fracture. Although the accumulation rate of hydrogen is known to be affected by its concentration and diffusion coefficient, it is still unclear how these properties are influenced by the carbon content in martensitic steel. In this study, we estimated the hydrogen accumulation rate in martensitic steel under different carbon contents (0.395 and 0.787 mass%) and temperatures (300, 700, and 1,000 K) by calculating the hydrogen solubility and diffusion coefficient using computational methods, including density functional theory, molecular dynamics simulation, and finite element method. Results revealed that the hydrogen accumulation rate tends to increase with the carbon content.


This Paper was Originally Published in J. Soc. Mater. Sci., Japan 69 (2020) 134–140. All captions of figures and tables were slightly modified.

Diffusion Coefficient of Hydrogen along c-axis in Martensitic Steel.

1. Introduction

There is a growing demand of high-strength steels such as martensitic steel for reducing the weight of automobiles. In general, it is known that the higher the strength of steel, the stronger the effect of hydrogen embrittlement. Particularly in high-strength steel with yield stress of 1 GPa or more, delayed fracture1) due to hydrogen embrittlement becomes a problem.

The phenomenon of hydrogen embrittlement has been known for a long time, and various theories on its causative mechanisms have also been proposed. In steel materials ruptured by delayed fracture, fracture surfaces at grain boundaries are frequently observed. Therefore, many studies focusing on the interactions between grain boundaries and hydrogen have been conducted.2,3) Recently, Takai et al. reported materials that show intergranular fracture at room temperature by hydrogen and intragranular fracture at liquid nitrogen temperature.4) This strongly suggests the involvement of diffusion phenomena in the fracture process. In addition, Yamaguchi et al. performed density functional theory (DFT) calculations and showed that the reduction of cohesive energy by hydrogen accumulated at grain boundaries only is small and that in order to observe embrittlement, it is necessary to include the contribution of mobile hydrogen moving to the fracture surfaces.5) Thus, it is important to predict the accumulation of hydrogen at grain boundaries and crack tips in delayed fracture.

The accumulation rate of hydrogen depends on the product of the amount of solute hydrogen and the diffusion coefficient of hydrogen. It is well known that the diffusion behavior of hydrogen is affected by the gradient of hydrostatic stress.6) In addition, Taketomi et al. showed that strong anisotropy appears in the diffusion coefficient of hydrogen in the strained lattice around dislocations in pure iron.7) Nagase et al. reported the non-hydrostatic stress dependence of the diffusion behavior in detail.8) Here, martensitic steel has a body-centered tetragonal (bct) structure elongated to one ⟨100⟩ axis of the body-centered cubic (bcc) structure by carbon atoms. Therefore, it has been suggested that the influence of the non-hydrostatic stress mentioned above is strong, but a detailed evaluation has not been carried out. Based on evidence from DFT calculations, Matsumoto et al. showed that in pure iron, the heat of solution of hydrogen changes greatly due to uniaxial strain along ⟨100⟩.9)

In this study, the dependences of solubility and averaged diffusion coefficient of hydrogen on carbon content in perfect martensitic phase are clarified. We estimated the hydrogen solubility based on DFT calculations and determined the diffusion coefficient by a combination of results from DFT calculations, molecular dynamics (MD) method, and finite element method (FEM). From these results, we estimated the effect of carbon on the accumulation rate of hydrogen.

2. Evaluation of Hydrogen Solubility

2.1 Calculation conditions

In this study, we used Vienna Ab initio Simulation Package,1012) which is a first-principle calculation software based on DFT. We also used generalized gradient approximation (GGA)13) in which spin was taken into account for the approximation of exchange correlation energy, and the projected augmented wave method14) for the treatment of inner shell electrons. Here, the cutoff energy was 425 eV. The cell size and the positions of all atoms were relaxed. We described the number of k-points later in Section 2.2. In this study, zero-point vibrational energy for hydrogen atom was not considered, but the effect was small in our calculations of the energy difference between hydrogen atoms in different occupied positions in solute state.

2.2 Atomic arrangement in martensitic steel

The diffusionless transformation15) occurred by quenching the austenite phase of the face-centered cubic (fcc) structure and allowing the martensitic phase of the bct structure to form. Martensitic steel was considered as having a bcc structure at low carbon (0.25 mass% or less) but becomes bct structure as carbon content increases. We defined the extended axial lattice constant as c, otherwise a. The relationship between the axial ratio c/a and the concentration of carbon M (mass%) was obtained experimentally and given as eq. (1):16)   

\begin{equation} \frac{c}{a} = 1.000 + 0.046 \times M \end{equation} (1)
There are three kinds of O sites in the bcc structure, depending on the direction of the fourfold rotational symmetry axis. The carbon atom occupies specific O sites, resulting in the bct structure. The atomic arrangement of the martensitic steels used in this study is shown in Fig. 1(a) to (d). The mass percentage concentrations of carbon in (a) Fe128C1, (b) Fe54C1, (c) Fe54C2, and (d) Fe16C1 are 0.167, 0.395, 0.787, and 1.3 mass%, respectively. The second carbon in (c) Fe54C2 is the farthest from the first carbon.17) The number of k-points was 7 × 7 × 7 for (d) Fe16C1 only and 5 × 5 × 5 for the others.

Fig. 1

Atomic structure of martensitic phase with different carbon content.

The relationship between carbon concentration and the axial ratio c/a obtained by our calculation is shown in Fig. 2. It can be seen that the calculated and experimental values are in good agreement except for Fe16C1 (1.3 mass%). A possible cause of the deviation from the experimental value in Fe16C1 is that in the present calculation, extremely high concentration of carbon is regularly introduced into O sites of the same kind due to the restriction of supercell size. Moreover, it is known experimentally that at 0.25 mass% or less, the axial ratio does not follow eq. (1). We therefore consider only Fe54C1 and Fe54C2 in the following evaluations.

Fig. 2

Relationship between tetragonality and carbon content.

2.3 Calculation of solubility

This section presents a method for evaluating the amount of hydrogen dissolved in the martensitic phase. Firstly, the calculation procedures of the heat of solution, trap energy, and occupancy at the interstitial sites are discussed.

[Heat of solution] The heat of solution of hydrogen EHOS was computed through the following equation:   

\begin{equation} E^{\text{HOS}} = \phi_{\text{XH}} - \left(\phi _{\text{X}} + \frac{1}{2}\phi_{\text{H}_{2}}\right) \end{equation} (2)
where ϕX and ϕXH are the ground-state energy of a system without and with hydrogen, respectively, and $\phi _{\text{H}_{2}}$ is the energy of a hydrogen molecule. Here we use $\phi _{\text{H}_{2}} = - 6.81$ eV.18)

[Trap energy] The trap energy of hydrogen Etrap_x at trap site x is expressed by eq. (3) as follows:   

\begin{equation} E^{\text{trap_${x}$}} = E^{\text{HOS_Tsite}} - E^{\text{HOS_${x}$}} \end{equation} (3)
where EHOS_x and EHOS_Tsite are the heat of solution at trap site x and tetrahedron (T) site in non-deformed bcc lattice composed by Fe atoms, respectively.

Because hydrogen occupies the T site, the interstitial site (trap site) of hydrogen in Fe54C1 considering symmetry is shown in Fig. 3. There are 36 independent positions. When symmetry is taken into account for Fe54C2, the number of independent hydrogen position also becomes 36.

Fig. 3

Position of hydrogen in Fe54C1.

[Hydrogen occupancy at each site] The relationship between the hydrogen occupancy Cx at trap site x and the hydrogen occupancy at the T site in the non-deformed bcc-Fe lattice CTsite can be expressed by eq. (4):   

\begin{equation} \frac{C_{x}}{1 - C_{x}} = \frac{C_{\text{Tsite}}}{1 - C_{\text{Tsite}}}\exp\left(\frac{E^{\text{trap_${x}$}}}{k_{\text{B}}T}\right) \end{equation} (4)
where kB and T are Boltzmann constant and temperature, respectively.

The hydrogen concentration in pure iron (T site in bcc phase) under a gaseous environment is given by Matsumoto et al.19) If the hydrogen gas pressure, p, is fixed, i.e., p = 70 MPa, the occupancies at the T site are CTsite = 4.68 × 10−7, CTsite = 1.45 × 10−4, and CTsite = 5.08 × 10−4 at T = 300 K, T = 700 K, and T = 1,000 K, respectively. The relation between trap energy Etrap_x and occupancy Cx is shown in Fig. 4. When the trap energy is 0, Cx becomes CTsite, and becomes higher at positive Etrap_x and lower at negative Etrap_x.

Fig. 4

Relationship between site occupancy Cx and trap energy Etrap_x under p = 70 MPa.

[Hydrogen concentration in martensitic phase] We define the number of iron atoms in the supercell as NFe (= 54), and the number of T sites equivalent to each individual hydrogen introduction position shown in Fig. 3 as nx. The hydrogen concentration CH (atomic ratio) was calculated in accordance with eq. (5) as follows:   

\begin{equation} C^{\text{H}} = \frac{1}{N_{\text{Fe}}}\sum_{x = 1}^{36}n_{x}C_{x} \end{equation} (5)
where $\sum\nolimits_{x = 1}^{36}n_{x} $ gives the number of T sites in the supercell, and it gives NTsite = 324. The hydrogen concentration in the martensitic phase normalized by the hydrogen concentration in bcc phase (pure Fe) under the same hydrogen environment (hereinafter referred to as normalized hydrogen concentration) is given by CH/(CTsite × NTsite/NFe).

Here, CTsite ≪ 1 is satisfied. Cx ≪ 1 also hold because the maximum value of Etrap_x is about 0.15 eV, which will be discussed in the following (see also the occupancy at this trap energy as shown in Fig. 4). Therefore, the normalized hydrogen concentration is not influenced by CTsite except under extremely high-pressure conditions.

2.4 Results and discussion on the solubility

Two main effects of hydrogen on trap energy are considered: the distance between carbon and hydrogen (C–H distance) and the strain of iron lattice due to interstitial carbon. The relationship between C–H distances and trap energy after relaxation is shown in Fig. 5. The C–H distances of Fe54C2 indicate the distances between the hydrogen and the carbon atom nearest to the hydrogen in the supercell. When the C–H distance was less than 2 Å, the trap energy value became very small (<−0.2 eV), and therefore it was excluded. As shown in Fig. 5, it can be confirmed that the fluctuation of trap energies increases as the C–H distance becomes closer, but the relationship cannot be described as linear.

Fig. 5

Relationship between trap energy and C–H distance.

We then consider the relationship between the strain of iron lattice and trap energy. Here, we consider the volumetric strain of an octahedron, which consists of the plane in which the site locates at and two iron atoms at both sides of the plane, as the strain at T site as shown in Fig. 6. The relationship between the volumetric strain and the trap energy is shown in Fig. 7. The volumetric strain was calculated by using the lattice constant of pure iron as the reference lattice constant. Sites with C–H distances less than 2 Å were excluded. Figure 7 shows that the relationship between the volumetric strain and the trap energy of hydrogen in pure iron can be fitted by a straight line. On this fitted straight line, the larger the volumetric strain, the larger the value of the trap energy. This is because the interstitial sites are expanded by the distortion of the iron lattice, allowing hydrogen to occupy easily. We can see that there are more data points along the straight line in Fe54C1 (Fig. 7(a)) with lower carbon concentration than in Fe54C2 (Fig. 7(b)). As the volumetric strain increases, the data points deviate from the straight line. This deviation from the line is due to the effect of the interaction between carbon atom and hydrogen atom and the effect of anisotropic strain.9) The trap energy value of Fe54C2 is generally higher than that of Fe54C1, and all values are positive.

Fig. 6

Octahedron formed by six Fe atoms which are surrounding four T sites.

Fig. 7

Relationship between trap energy and volumetric strain in (a) Fe54C1 and (b) Fe54C2.

The hydrogen concentration calculated from the obtained trap energy and eqs. (4) and (5) is shown in Fig. 8. The vertical axis represents the atomic ratio of hydrogen to iron, and the values on the bar graph represent the normalized hydrogen concentration. It has been known that the higher the temperature and the carbon concentration, the greater the amount of hydrogen dissolved. The lower the temperature, the greater the carbon concentration dependence of the normalized hydrogen concentration.

Fig. 8

Atomic ratio of H to Fe atom at T = 300, 700, and 1,000 (K): The values on the bars indicate the hydrogen concentrations normalized by the concentration in Fe54.

3. Evaluation of Diffusion Coefficient of Hydrogen

3.1 Calculation procedure

In this study, the martensitic phases were decomposed into bcc lattices subjected to different strains. The diffusion coefficients of these phases were individually evaluated, and the averaged diffusion coefficients of the structures were obtained by continuum approximation. Specifically, we first obtained the strain in each bcc lattice in Fe54C1 and Fe54C2 from the results of DFT calculations in Section 2.2 (Section 3.2). The diffusion coefficient of hydrogen in the bcc iron lattice at the obtained strain was then estimated using the MD method (Section 3.3). Finally, the averaged diffusion coefficients of hydrogen in the periodic structure composed of bcc lattice with different diffusion coefficients were obtained by FEM (Sections 3.4 and 3.5). Although the validity of an approach to obtain firstly the distribution of physical properties and then the average physical properties of atomic-scale structures by continuum approximation is controversial, a similar approach has been, for example, adopted for the evaluation of elastic constants of materials with two-phase structures at the nanoscale, and its validity has been confirmed.20)

3.2 Distortion of bcc lattice

The abc axis is shown in Fig. 9. Because Fe54C1 and Fe54C2 are symmetric about the ac and bc planes, the bcc lattice in the model is classified into eight types. Table 1 shows the strains obtained by DFT along the a, b, and c direction of the bcc lattices (n = 1∼8) in Fig. 9.

Fig. 9

bcc lattice with different strain in Fe54C1.

Table 1 Strain [%] of bcc lattice around carbon atom.

3.3 Diffusion coefficient of hydrogen in strained bcc lattice

We used atomic system composed of 2,000 Fe atoms, where 10 bcc lattices were arranged along each direction and 1 H atom which was initially located at T site. The periodic boundary conditions were applied along all directions. We perform MD calculations of 3 ns to obtain trajectory of H atom.

The first 10 ps was considered as the relaxation time and was not used to evaluate the diffusion coefficient. We employed the embedded atom method potential developed by Wen.21) The calculation was performed for 11 different temperatures in the range from 300 K to 1,000 K. Strain was given by using the equilibrium lattice constant at each temperature as a reference. The evaluation method of diffusion coefficient by MD calculation was shown in Ref. 8). The temperature dependence of the one-dimensional diffusion coefficient along a, b, and c direction in Fe54C1 is shown in Fig. 10. Diffusion coefficient has strong anisotropy, and it is known that the values vary greatly depending on the lattice. Particularly in lattice 1 where carbon exists, the diffusion coefficients ( Point) become very small. In lattice 1, a sufficient number of samples cannot be obtained on the low-temperature side, and the diffusion coefficients were not correctly determined.

Fig. 10

Diffusion coefficient of hydrogen in deformed bcc lattice in Fe54C1.

3.4 Evaluation of averaged diffusion coefficient by FEM

We evaluated the averaged diffusion coefficient of periodic structure composed of lattices with different diffusion coefficients by FEM calculation. In this study, we performed the analysis using MSC Marc by replacing diffusion analysis with thermal conduction analysis because of the similarity between Fourier’s law and the Fick’s law. In this evaluation, the effect of the difference in the hydrogen concentrations between bcc lattices on the apparent diffusion coefficient was ignored. In order to apply the symmetric boundary condition in FEM, the arrangement of each lattice is changed as shown in Fig. 11, and the analytical model is divided into 9 × 9 × 9 elements. The thermal (hydrogen) flux was obtained by giving different temperatures (hydrogen concentration) to both end faces in the direction in which the thermal diffusivity (diffusion coefficient) was evaluated. The symmetric boundary conditions were applied in the direction perpendicular to the evaluation direction. The diffusion coefficients of hydrogen in Fe54C1 and Fe54C2 were determined based on their ratio to Fe54. Here, it is known from the results in Section 2.4 that repulsion acts between carbon and hydrogen atom, and hydrogen does not occupy the nearest neighbor lattice of carbon. Therefore, the diffusion coefficients of lattice 1 for Fe54C1 and lattice 1 and 7 for Fe54C2 are set as 0. The diffusion coefficients of other lattices were acquired according to the result (Fig. 10) obtained by MD calculation. Note that 1 to 8 in Figs. 9, 10, and 11 correspond to the same bcc lattice.

Fig. 11

FEM model for the calculation of averaged diffusion coefficient.

3.5 Results and discussion on the diffusion coefficient

The flux distribution at 300 K is shown in Fig. 12. The color changes from blue to red and yellow as the flux becomes larger. It can be confirmed that there is a large flux away from the carbon and a strong anisotropy between the a (= b) and c axes.

Fig. 12

Flux distribution in Fe54C1 and Fe54C2 models.

Figure 13 shows the temperature dependence of the diffusion coefficients. Table 2 shows the diffusion coefficients at 300 K, 700 K, and 1,000 K. It can be seen that the higher the carbon concentration, the smaller the diffusion coefficient, and that the influence is more pronounced in the c-axis direction.

Fig. 13

Diffusion coefficient of hydrogen in Fe54C1 and Fe54C2 along (a) a-, b-axis and (b) c-axis.

Table 2 Diffusion coefficient D [10−10 m2/s] of hydrogen in Fe, Fe54C1, and Fe54C2 at 300, 700, and 1,000 K.

4. Accumulation Rate of Hydrogen

Table 3 shows the ratios of the product of the hydrogen concentration and the averaged diffusion coefficient of hydrogen for Fe54C1 and Fe54C2 to the value of pure iron. It can be seen that the accumulation rate of hydrogen is smaller along the c-axis with respect to a- (= b-) axis under all conditions. Furthermore, comparison of carbon concentrations reveals that the rate of Fe54C1 is slower than that of pure iron, while the rate of Fe54C2 is faster than that of pure iron, except along the c-axis at the high-temperature region.

Table 3 Hydrogen accumulation rate in Fe54C1 and Fe54C2 at 300, 700, and 1,000 K normalized by it in Fe54.

5. Conclusion

To clarify the carbon concentration dependence of hydrogen accumulation rate in martensitic steel, hydrogen concentration was calculated using DFT calculation, and the diffusion coefficient of hydrogen was estimated by the results from DFT, MD, and FEM calculations. As a result, it was found that the accumulation rate of hydrogen decreases when the carbon concentration is low and increases when the carbon concentration is high. It was also found that diffusion anisotropy of hydrogen occurs in the martensitic phase.

© 2020 The Society of Materials Science, Japan