Numerical Analysis of the Coupling between Hydrogen Diffusion and Mechanical Behavior near the Crack Tip of Titanium

Yang, Fu-qiang; Zhan, Wen-juan; Yan, Tao; Zhang, Hai-bing; Fang, Xiu-rong

doi:https://doi.org/10.1155/2020/3618589

Mathematical Problems in Engineering

On this page

Abstract Introduction Discussion Conclusion Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article Corrigendum

!

A Corrigendum for this article has been published. To view the article details, please click the ‘Corrigendum’ tab above.

Research Article | Open Access

Volume 2020 | Article ID 3618589 | https://doi.org/10.1155/2020/3618589

Numerical Analysis of the Coupling between Hydrogen Diffusion and Mechanical Behavior near the Crack Tip of Titanium

Fu-qiang Yang,¹Wen-juan Zhan,²Tao Yan,²Hai-bing Zhang,³and Xiu-rong Fang²

Academic Editor: Mijia Yang

Received04 Jan 2020

Revised18 Mar 2020

Accepted26 Mar 2020

Published25 Apr 2020

Abstract

Hydrogen plays a detrimental effect on the degeneration of titanium and its alloys, and it is very important to quantify the hydrogen concentration when estimating the microstructure evaluation of titanium and its alloys in a hydrogen environment. In this paper, the hydrogen atoms are assumed to reside in interstitial sites and in trapping sites such as dislocations, and a mechanic-diffusion coupled model was proposed to describe the stress effects on the diffusion of hydrogen in titanium. A titanium plate with a central crack was modeled to verify the mechanic-diffusion model, and it is solved by the finite element method in commercial software COMSOL. The results indicate that hydrogen diffusion near the crack is determined by the stress state. When the stress state of the crack tip is elastic, the hydrogen will diffuse from both sides of the crack towards the tip and lead to the highest hydrogen concentration in the crack tip. When a plastic zone exists in front of the crack tip, the highest hydrogen concentration at crack surface deviates to the side near the crack tip; a hydrogen concentration peak exists at a characterized distance in front of crack tip initially and then diminishes with the diffusion process. The proposed model is expected to solve the hydrogen concentration under stress in more complex structures.

1. Introduction

Titanium and its alloys are attractive materials for many structural applications in aerospace, industrial, and marine, because of their excellent specific strength, stiffness, and corrosion resistance. Although titanium is considered to be immune to corrosion, the hydrogen embrittlement, which manifests as a reduction in mechanical properties of titanium and its alloys, becomes a common phenomenon when they come in contact with a hydrogen-containing environment [1]. Several mechanisms have been summarized to reveal the reason for hydrogen embrittlement, including the stress-induced hydride formation and cleavage mechanism [2, 3], decohesion [4, 5], and the hydrogen-enhanced localized plasticity mechanism [6, 7]. Of all the mechanisms mentioned, it is very important to accurately model the hydrogen transportation process in metals when determining the likelihood of cracking.

It was found that the hydrogen transport process is strongly influenced by the hydrostatic stress and the plastic strain [8, 9], so Sofronis and McMeeking carried out coupled diffusion elastic-plastic stress finite element analyses incorporating the effect of hydrostatic stress and trapping [10]. And then Abramov and coworkers [11, 12] used a model accounting for diffusion and trapping based on McNabb and Foster’s work [13] and Oriani’s extension [14] for local equilibrium between lattice and trap sites. Following Sofronis and Mcmeeking’s work [10], Cui et al. developed a fully implicit numerical algorithm to superimpose the hydrogen-induced strain onto the equivalent plastic strain during the stress update process within each iteration on each integration point in the implicit backward Euler algorithm, by which the interactive plastic strain and the hydrogen concentration could be solved simultaneously [15]. Considering that sufficient time is needed for the diffusion process to reach a steady state of hydrogen distribution between the trapping sites and the normal interstitial lattice sites [16]. Luo et al. firstly determined the loading speed so that a steady-state hydrogen distribution is achieved and then studied the conjoint effect of hydrogen and stress state on ductile fracture by unit cell model [17]. Gobbi et al. proposed a weakly coupled approach with three steps to present the couples between hydrogen diffusion and stress-strain analysis in Abaqus [18]. However, the hydrogen concentration is imported as a predefined field and it is not updated during crack propagation, so a fully coupled model was developed to simulate the coupling of thermal stress, mass diffusion, and heat transfer [19]. Díaz et al. [20] repeatedly revisited the coupled hydrogen diffusion simulation proposed by Sofronis et al., but they considered the plastic strain rate, the coupled diffusion, and stress-state dependent boundary conditions in the model. In the model proposed by Stopher et al. [21], the hydrogen-trapping approach is combined with physical modeling of nucleation, growth, and coarsening of second-phase particles within the matrix, and the hydrogen-trap binding enthalpies and mechanical loading conditions are of relevance too and are considered. Thus, the role of microstructure on hydrogen diffusion is associated with hydrogen trapping in multi-precipitate distributions simultaneously.

Similar hydrogen diffusion models coupling with stress were also developed by other researchers to study different problems [22–26]. However, when dealing with the coupled mechanical and hydrogen transport problems, a so-called Oriani’s assumption and McNabb and Foster equation are widely used in the models. Many analytical solutions have been developed for the McNabb and Foster equation, such as approaches proposed by Javanmardi and Bashiri [27] or Azizian [28] to model the adsorption process and by Toribio et al. [29, 30] for hydrogen-trapping estimation. These formulations, however, though quite exacting, might be complex to compute numerically, especially when both hydrogen concentration and mechanical loading fields are unknown. So, an approximation of the analytical solution is proposed by Benannoune et al., in which the Macroscopic Rate Equation is used [31].

Other methods were also mentioned to solve the hydrogen diffusion and trappings, such as molecular dynamics simulation [32, 33], first-principles calculation [34], and phase-field models [35]. However, there is no coupling between diffusion and stress in these models.

In this paper, a hydrogen transport equation is proposed to determine how hydrogen is distributed in titanium and how the stress state influences the hydrogen diffusion. This paper aims to explain in detail how to solve the coupled hydrogen diffusion-mechanics problems by the finite element method in commercial software COMSOL. By using the diffusion model built in this paper, the hydrogen diffusion in titanium structure with crack could be solved, and the quantitative prediction of hydrogen embrittlement could be achieved.

2. Hydrogen Transport Equation

2.1. Hydrogen Position in Metals

Hydrogen is an interstitial atom that exists in the crystal lattice of metals, and the vibration causes their diffusion by a random jump to a neighbor site. The diffusion only model was built based on Fick’s laws to describe the hydrogen diffusion in metals, but the results significantly deviated from experiment [36]. One reason is that the metals were assumed to have ideal crystal lattices in the model; in fact, there are defects such as inclusions, vacancies, dislocations, grain boundaries, and second phases existing in the metal lattice. These defects were called hydrogen traps, and they could retain the hydrogen, thus increasing the apparent solubility and decreasing the apparent diffusivity of hydrogen [36]. To incorporate the hydrogen traps into the hydrogen transport equation, we follow the equilibrium theory presented by Oriani [14] to consider a lattice consisting of two kinds of sites for occupancy by hydrogen: a vast majority of sites called “normal interstitial” site and a minor fraction of sites called “trapping” sites. The total hydrogen concentration is given bywhere C_tot is the total hydrogen concentration in metals, mol·m⁻³; C_L denotes the hydrogen concentration in the normal lattice, mol·m⁻³; and C_T denotes the concentration associated with the hydrogen in the traps, mol·m⁻³. The two kinds of hydrogen concentration C_L and C_T are related to the number of sites:where α is the number of interstitial sites per solvent atom; β = 1 denotes the number of trapping sites; θ_L and θ_T are the fraction of lattice sites and trapping sites occupied by H atom, respectively, 0 ≤ θ_L, θ_T ≤ 1; N_A is Avogadro’s number, 6.023 × 10²³ atom·mol⁻¹; and N_L is the number of atoms of solvent per volume, mol·m⁻³, which could be calculated aswhere ρ (g·m⁻³) and M (g·mol⁻¹) are the density and relative atomic mass of the metal. And N_T is the number of trapping sites per volume.

The number of trapping sites N_T is not a characteristic value for metals, and it is very difficult to determine due to the large variety of microstructural defects that can be included. However, various investigations have shown that the trap population is associated with dislocation density and depends on the level of plastic strain (ε_P). Here we only consider the hydrogen trapped at dislocations as the only kind of saturable and reversible trap. To assume a trap for each atomic plane intersected by a dislocation, the number of trapping sites N_T could be calculated as [36, 37]where b is the Burgers vector of unit dislocation, m, and it is the magnitude of a/2<111> in BCC and a/3 <1120> in HCP, which equals and a, respectively; a is lattice constant, m; ρ_d is the dislocation density measured in dislocation line length per cubic meter, m·m⁻³, and it depends on the accumulated plastic strain as follows [36]:where ρ₀ denotes the dislocation density for the annealed material, m⁻², and γ is defined as the dislocation strengthening factor caused by plastic deformation.

Although there are also other empirical relationships for N_T − ε_P [10, 37, 38], it should be noted that this is an important simplification; the number of traps depends not only on dislocations but also on other types of traps.

2.2. Diffusion Driving Force and Hydrogen Flux

The hydrogen diffusion in metals is a matter exchange process, and it will continue until thermodynamic equilibrium is reached. However, this equilibrium does not always mean an absence of concentration gradients [39]; other potential factors could also cause the flux of hydrogen. Usually, the diffusion could be caused by many factors such as chemical potential μ, temperature T, pressure P, and other external potentials, so the diffusion driving force could be the gradients of one or more of these potential forces. When several of these processes occur at the same time, the problem becomes too complex to solve, so it is considered that all driving forces can be reduced into gradients of chemical potential to simplify the problem [40]. And this simplification allows including both random effects and the deviation from them in the chemical potential as the only driving force. In this case, the hydrogen flux iswhere J is the hydrogen flux, mol·m⁻²·s⁻¹; D_L is the diffusion coefficient for hydrogen, m²·s⁻¹; μ_L is the chemical potential of the hydrogen in the lattice sites; R is the universal gas constant, 8.3144 J·mol⁻¹·K⁻¹; and T is the absolute temperature, 273 K.

In equation (7) the diffusion is governed by chemical potential gradient; the equation is traditionally expressed in terms of concentration. For a dilute solution with low hydrogen occupancy of θ_L << 1, the relationship between chemical potential of hydrogen (disregarding the stress state) and hydrogen could be achieved by replacing the activity with hydrogen concentration aswhere μ_L0 is the chemical potential at a reference temperature and pressure.

For tetrahedral or octahedral sites of the crystal, both the sites are too small for the hydrogen atom, so a volume increment will be induced when hydrogen dissolved into the lattice. Studies have shown that, in a lattice subjected to a traction stress state, the sites will be wider and hence the chemical potential will be lower. For stressed solid in different thermodynamics [41], the following has been found:where σ_h is the hydrostatic stress, which is calculated by the on the diagonal terms of the stress tensor, σ_h = ∑σ_ii/3; V_h is the partial molar volume of hydrogen, m³·mol⁻¹; μ(σ_h0) and μ(σ_h) are the chemical potential without stress and stress level of σ_h, respectively.

From equation (9) it could be seen that the energy of interaction of hydrogen atoms with the stress field only depends on the diagonal terms of the stress tensor σ_ii, and the hydrogen partial molar volume V_h.

By incorporating equations (8) and (9), the chemical potential of hydrogen in a stressed metal could be described as

Substituting equation (10) into equation (7), we can obtain the hydrogen flux expressed by the concentration of hydrogen as

2.3. Mass Balance Equation

The diffusion of hydrogen in titanium is defined from the requirement of mass conservation for the total hydrogen concentration:where is the partial derivative with respect to time; V is any volume whose surface is S, n denotes the outward normal to ; J is the hydrogen flux through the surface , mol·m⁻²·s⁻¹; and J·n is the hydrogen concentration flux leaving S.

By substituting equations (1) and (11) into equation (13), we obtain

Applying the divergence theorem, equation (13) is expressed as

This constitutive equation represents a modified form of the second Fick’s law; however, there are two unknows concentration C_L and C_T. So, the relationship between C_L and C_T should be assumed before solving the equation with a numerical method.

During hydrogen diffusion in metals, it jumps between lattices interstitial sites and traps; Oriani’s theory states that the variation of concentration in traps can be expressed as the total jumps from L to T minus the total jumps from T to L [14]; when the total number of interstitial sites is much greater than the total number of traps N_L >> N_T, and there is a low occupancy θ_L << 1, the relationship between θ_L and θ_T could be acquired in the case of equilibrium:where K_T is the trap equilibrium constant, aswhere E_b is the trap binding energy concerning the lattice site; it is inherently negative energy [10], J·mol⁻¹.

By incorporating equations (2) and (15) into equation (3), the hydrogen concentration in trap sites could be expressed aswhere and are defined as

Since the number of lattice sites N_L is a constant and the temperature will be kept constant, the hydrogen concentration in trap sites C_T is a function of the hydrogen concentration in lattice sites C_L, and the number of trap sites N_T depends on the plastic deformation level described by equations (5) and (6).

So the partial derivative of the hydrogen concentration in trap sites with respect to time in equation (14) becomeswhere the partial derivative items and could be obtained from equation (17), asand the item could be induced by combining equations (5), (6), and (9), as

Finally, by substituting equation (20) into equation (14), we obtain the mass constitutive equation as

3. Analogy with the Diffusion Problem in COMSOL

In COMSOL, the mass conservation equation for one species transport in diluted solution is given as [42]where c is the concentration of the species, mol·m⁻³; D denotes the diffusion coefficient m²·s⁻¹; R_a is a reaction rate expression for the species, mol·m⁻³·s⁻¹; and u is the velocity vector, m·s⁻¹.

In equation (24), the mass transport mechanisms contain diffusion and convection. The second term accounts for the diffusive transport, accounting for the interaction between the dilute species and the solvent. The third term on the left side describes the convective transport due to a velocity field u. On the right-hand side of the mass balance equation, R_a represents a source or sink term, due to a reaction or desorption.

To implement the coupled hydrogen transport equation given in equation (23), it is analogized with equation (24) and changed to the following form:where d_a is mass coefficient defined as

Define effective diffusion coefficient D_eff, convection velocity u, and reaction rate R_a as

It is found in equation (25) that the hydrogen transport mechanism in stressed metals also contains diffusion and convection. The effective diffusion coefficient D_eff is a variable, which is a function of the hydrogen concentration in lattice sites C_L and the number of trapping sites associated with plastic deformation. The hydrogen convection velocity u is determined by the hydrogen concentration in lattice sites, the degree of plastic deformation, and also the hydrostatic pressure gradient. The right item of equation (25) corresponds to the source or sink of hydrogen; as this equation solves the hydrogen concentration in lattice sites, if there is no hydrogen injected or flowing out of the metal, the reaction rate R_a in equation (25) could be considered as the total rate for hydrogen adsorption or dispersion in the trapping sites.

4. Physical Model and Coupling Simulation Process in COMSOL

4.1. Geometric Model, Mesh Model, and Material Model

To analyze the coupling between stress and hydrogen diffusion in a crack tip, a plate with a hole crack in the center is modeled in 2D as shown in Figure 1. The length and width of the plate are 2H and 2W, respectively; the radius of the hole is r₁. To simplify the calculation, only a quarter of the model was considered. In this simulation, W and H are equal to 60 mm, r₁ = 4 mm.

(a)

(b)

(c)

To strike the balance between calculation speed and calculation accuracy, the quarter model is divided into three regions to acquire finer mesh near the crack and coarse mesh away from the crack. As shown in Figure 1, the three regions are divided by two arcs with a radius of r₂ = 10 mm and r₃ = 20 mm. The whole model was meshed to 6412 free triangular elements, and region 1 has 4838 triangular elements. With a thickness of 0.6 mm, the plane strain condition was applied to the model.

For a hardening material, the elastoplastic properties could be given by the Swift power-law equation [42]:where σ_Y (MPa) is the local flow stress with the plastic strain of ε_p; σ₀ (MPa), ε₀, and n are the yield stress, yield strain, and hardening exponent measured in the absence of hydrogen, respectively, and σ₀· = ε₀E, E is Young’s modulus, MPa.

4.2. Boundary Conditions and Loads

As an initial condition, the hydrogen is assumed to have already diffused into the stress-free solid and a uniform concentration C_L = C₀, and it is assumed that there is no hydrogen flux through the edges of the specimen during the simulation. A uniform stretch load was applied on the upper edge, and the geometric symmetry conditions are considered for both diffusion processes and the structural mechanics problem.

4.3. Simulation Parameters

The material of simulated plate in Figure 1 was considered as Ti-6Al-4V titanium alloy, and its physical and mechanical properties are listed in Table 1, the other parameters used in the simulation are also listed in Table 1, and they will be discussed in the following:(1)Binding energy (E_b) Perssouyer et al. found that the interaction energy of Ti-H is in the range of 0.22 eV∼0.27 eV [43, 44], and well below the value necessary to form an actual hydride, as expected. However, it is found that the binding energy E_b between hydrogen atoms and dislocations is only 0.035 eV in the β-phase titanium alloy Ti-10V-2Fe-3Al, which is one order of magnitude smaller in comparison with that in austenitic steels [45]. In this study, the binding energy is selected as 0.25 eV (about 26.05 kJ·mol⁻¹) for dislocation trapping according to Rath and Bernstein [46].(2)Dislocation density With the deformation at 40%, 60%, and 80% reduction levels in commercially pure titanium specimens, the dislocations densities calculated are 5 × 10¹⁴ ± 5% m⁻², 8 × 10¹⁴ ± 10% m⁻², and 10¹⁵ ± 7% m⁻², respectively [47], while the dislocation densities are 2.1 × 10¹⁴ m⁻², 3.58 × 10¹⁴ m⁻², and 4 × 10¹⁴ m⁻² for vertical rolling direction and 1.15 × 10¹⁴ m⁻², 2 × 10¹⁴ m⁻², and 3 × 10¹⁴ m⁻² for horizontal rolling direction with 40%, 60%, and 80% cold rolling [48]. With a strain rate of 0.01 s⁻¹, the average dislocation densities of near-α titanium alloy Ti-5.4Al-3.7Sn-3.3Zr-0.5Mo-0.4Si alloy is measured as 4.18 × 10¹³ m⁻² at 1233 K [49]. For the α + β type dual-phase TC6 titanium alloy, α and β phases have the pre-existing dislocations of 1.17 × 10¹⁴ m⁻² and 2.45 × 10¹⁴ m⁻², respectively. During the deformation, the dislocation density of the α phase increases significantly from 1.88 × 10¹⁴ m⁻² (ε = 0.73%) to 3.36 × 10¹⁴ m⁻² (ε = 1.02%) and then plateaus until ε = 1.30%. The dislocation density of the β phase increases slightly from 3.36 × 10¹⁴ m⁻² to 3.79 × 10¹⁴ m⁻² at first and then increases sharply to 4.87 × 10¹⁴ m⁻² [50]. Even though the dislocation density is related to plastic degree, material type, and strain rate, the magnitude of the dislocation density is 10¹⁴ m⁻², so the initial dislocation density, ρ₀, is considered as 10¹⁴ m⁻², and the dislocation density keeps constant as 10¹⁵ m⁻²; when plastic strain ε_p is larger than 0.5, then the dislocation strengthening factor γ is 1.8 × 10¹⁵ m⁻².(3)Partial molar volume of hydrogen (V_h) When hydrogen dissolved in metals, the partial molar volume of hydrogen is related to the stress state, hydrogen concentration, and temperature. It is proposed that the data for strain vs H concentration can be used to calculate V_h, which gives an independent determination of the partial molar volume by the following equation [51]: where V_s is the molar volume of the solution at the composition in question (not the partial molar volume, and not the molar volume of the pure titanium) and N_H and ΔN_H are the atom percent of hydrogen and its change with a strain Δe₁₁. Waisman et al. compared the partial molar volume results obtained by different materials and methods; the values were acquired at different temperatures and different H atom fraction [51], so the partial molar volume of hydrogen is arbitrarily selected 1.6 × 10⁻⁶ m³·mol⁻¹.(4)Density The elemental titanium’s density is nominally computed on the basis of 4.512 g·cm⁻³, and then the densities of α, α + β, and β titanium alloys are determined approximately by applying the factors of 1.019, 1.025, and 1.094 [52]; the calculated densities of the three types titanium are 4.597 g·cm⁻³, 4.624 g·cm⁻³, and 4.936 g·cm⁻³, respectively. The relative atomic mass of Ti is 47.9 g·mol⁻¹, and it is used instead of the values of different titanium alloys.(5)Number of interstitial sites The crystal structure of α titanium is HCP with 6 atoms included in a crystal cell, which has 6 octahedral interstices and 12 tetrahedral interstices, so the number of interstitial sites per solvent (Ti) atom of α titanium is 3. Even though the β titanium with BCC structure also has 6 octahedral interstices and 12 tetrahedral interstices, only 2 atoms are involved in a crystal cell, so the number of interstitial sites per solvent (Ti) atom of β titanium is 9. The α + β titanium has both HCP and BCC structures, so the values are in the range of 3∼9, and it is determined as 6 arbitrarily.(6)Hydrogen diffusivity Many hydrogen diffusivities are obtained at elevated temperatures; however, it is not suggested to extrapolate the data obtained at elevated temperatures to room temperature. The most important reason for not extrapolating is that, in polycrystalline materials, the activation energy is radically lower at lower temperatures. This change is caused by differences in diffusion paths at high and low temperatures primarily bulk diffusion at the elevated temperatures, with both bulk and grain boundary diffusion at the lower temperatures. This effect would make the room temperature value of diffusivity orders of magnitude higher than that predicted from extrapolation from the elevated temperature value [53]. Even though the information of hydrogen diffusion at room temperature is extremely scarce, the diffusion coefficient at room temperature collected by Madina et al. is 3.2 × 10⁻¹⁶ m²·s⁻¹ for α phase titanium and 3.3 × 10⁻¹¹ m²·s⁻¹ for the β-phase titanium [54]. Waisman et al. have summarized that a diffusivity of 4.95 × 10⁻¹³ m²·s⁻¹ in Ti-6Al-4V alloy has a good fit with the pressure vessel measurements [53] at 27°C, and this meets the key point that the diffusivities for β-Ti alloys, with or without minor alloying elements, are comparable to that for palladium, of the order of 10⁻¹¹ m²·s⁻¹ [55]. It can be seen that the diffusivity of hydrogen in α phase titanium is much lower than others, but Malyshev et al. obtained much higher data varying from 2.7 × 10⁻¹⁰ m²·s⁻¹ to 2.7 × 10⁻⁹ m² s⁻¹ [56] for pure VT1-00 titanium at room temperature, and it is not used in this study. It should be noted that the diffusivities acquired above are not lattice diffusivity but apparent diffusivity; they are different for different alloy [55]; however, it is not distinguished in this study.

4.4. Coupling Simulation Process

In COMSOL, there are many built-in modules to solve different physical problems. To solve the interaction effects between stress and hydrogen diffusion, three modules should be involved. One is the solid mechanics module, which solves the stress and strain state of the plate. The second is the transport of diluted species module, which solves the hydrogen concentration in lattices C_L with equation (25). The modules of solid mechanics and transport of diluted species are built-in modules in COMSOL, so they were used directly. However, the hydrogen concentration in traps C_T, which is described by a partial differential equation (PDE) of equation (20), should also be solved; thus a general form PDE was involved as the third module. The coupling of the three modules is by transferring the same simulation parameters.

According to equation (25), the hydrogen diffusion is affected by stress and strain, so the solid mechanics module should be solved independently at first to obtain the mechanically related variables used in the transport of diluted species module. Then the three modules could be solved simultaneously, and the coupling between hydrogen diffusion and stress field near the crack tip could be achieved.

5. Numerical Results and Discussion

5.1. Model Verification with Pure Diffusion

To verify the proposed hydrogen transport equation, the boundary conditions of the physical model were modified, so that the model can be simplified to pure diffusion, and to compare the results with the analytic solution. The load applied on the upper edge DE is set to 0 to avoid the stress effects. The hydrogen of the right edge CD is 40 mol/m³ during the diffusion, which is considered as the hydrogen source. The initial hydrogen in the plate is still uniformly distributed of 20 mol/m³. Without the effects of stress, it is assumed that the hydrogen only diffuses in the lattice, and there is no interaction between C_L and C_X, which indicates that

Without the effects of stress and trap, equation (25) could be simplified as

For the diffusion in a semi-infinite plate, the analytic solution readswhere C_b and C₀ are the hydrogen concentration in boundary and the plate, mol/m³. x is the distance away from the boundary, mm. C_L(x, t) is the hydrogen concentration at position x at time t. erf is error function with the formation as

Figure 2 shows the comparison between analytic solution and numerical results of the transient hydrogen concentration distribution along path CA, where x = 0 is point C. It can be seen that the numerical results almost perfectly agree with the analytic solution of equation (31).

5.2. Coupled Diffusion-Mechanics under Constant Load

Figures 3 and 4 demonstrate the redistribution of hydrogen in normal lattice and traps under different constant load. Even though Figures 3 and 4 look the same, they denote different variables, which should be concerning when considering the hydride formation in titanium alloy, so we show them independently. When the load is in the range of 200 MPa to 300 MPa, it can be seen that the hydrogen left the free end of the crack and diffusing into the crack tip; this causes the reduction of hydrogen in the free end and the increase of hydrogen in the crack tip. Moreover, this tendency was strengthened with the increase of load. With the remote load equal to 200 MPa, the highest hydrogen concentration in lattice increased to 23 mol/m³ at the crack tip, and the value increased to 25 mol/m³ when remote load increases to 300 MPa. The hydrogen concentration in traps increases from 9 × 10⁻³ mol/m³ to 1.13 × 10⁻² mol/m³ when remote load increases from 200 MPa to 300 MPa. However, when the remote load continues to increase, the hydrogen distribution characteristic changes. The free end of the tip still maintained as the source of hydrogen flowing out; the highest hydrogen concentration changed from the crack tip to a location between the crack tip and free end around the crack contour. The crack tip hydrogen is somewhat lower than the initial value. Another notable change is the reduced redistribution of hydrogen. Even though the remote load increased to 350 MPa and 400 MPa, the highest hydrogen concentration in the lattice is 20.65 mol/m³, and the hydrogen in traps is 5.81 × 10⁻³ mol/m³, and the values decrease with the increased remote load.

(a)

(b)

(c)

(d)

(e)

(a)

(b)

(c)

(d)

(e)

A more detailed distribution of hydrogen in normal lattice and traps is shown in Figures 5 and 6. Focusing on Figure 5, it is found that the hydrogen concentration decreases from crack tip to the free end of crack when load increases from 200 MPa to 300 MPa, and their differences increase with the increased remote load. However, the highest hydrogen concentration deviates to the sides of the crack tip, and deviation was strengthened with the increased remote load. Figure 6 shows the hydrogen concentration in front of the crack tip; it can be seen that when the remote load is in the range of 200 MPa∼300 MPa, the hydrogen concentration decreases with the increases of the distance from the crack tip, and the hydrogen concentration at crack tip changes dramatically when remote load changes from 200 MPa to 300 MPa, but the hydrogen concentration changes a little away from the crack tip. When the remote load increases to 350 MPa, the lowest hydrogen concentration appears at the crack tip, and the value increases when away from the crack tip. When the remote load continually increases to 400 MPa, the crack tip hydrogen concentration reduces to a lower level.

(a)

(b)

(a)

(b)

According to equation (25), the diffusion is driven by the concentration gradient and Mises stress gradient, where the hydrogen will diffuse from high concentration region to low concentration and diffuse from low hydrostatic stress region towards high hydrostatic stress. As the hydrostatic stress distribution demonstrated in Figure 7(a), when the remote stress is under 300 MPa, the largest hydrostatic stress locates at the crack tip and decreases as the distance increases in front of the crack tip, which leads to negative hydrostatic gradient stress, and causes the hydrogen diffusion towards the crack tip. The numerical results demonstrated in Figures 5 and 6 are as expected as the hydrogen concentration decreases in front of the crack tip with the increased distance. However, the hydrogen concentrations C_L and C_T are higher than the initial value of 20 mol/m³ and 5.81 × 10⁻³ mol/m³. This indicates that the area in front of the crack tip is the region of hydrogen inflow, even a negative value of . The source of hydrogen to outflow is the region furthest away from the crack tip along the crack contour; this is because a more negative hydrostatic stress gradient exists along the crack contour, as shown in Figure 7(b).

(a)

(b)

When the remote load increases, the strengthened hydrostatic stress leads to a more intense crack tip hydrogen concentration. But the situation is quite different when the remote stress increases above 300 MPa; the highest hydrogen concentration is not located at the crack tip, and the highest hydrogen is only a little over the initial value. We attribute this change to the stress state variation. As shown in Figure 8, when the remote stress is above 300 MPa, the stress in the tiny zone of the crack tip is above 800 MPa, which has exceeded the yield stress of the material, and a plastic zone exists at the crack tip. The boundaries of plastic and elastic locate about 0.278 mm and 0.491 mm in front of the crack tip with a remote load equal to 350 MPa and 400 MPa, respectively. By comparing Figures 7(a) and 8(a) it is found that the largest hydrostatic stress also locates at these positions. The same boundary could also be found around the crack contour, as shown in Figure 8(b). The Mises stress change is not smooth at 18.62 and 24.23 with corresponding remote load equal to 350 MPa and 400 MPa; as the Mises stress is above 800 MPa when the angle is smaller than 18.62 or 24.23, which indicates the yield of material, the two positions are the boundaries of plastic and elastic. The existing boundary of elastic and plastic also causes the unsmooth hydrostatic stress change shown in Figure 7(b) and finally leads to the different hydrogen diffusion behavior.

(a)

(b)

The hydrogen concentration distribution is much different when remote load is above 300 MPa, so a much detailed hydrogen concentration change during the diffusion process is shown in Figures 9–12. It can be seen from Figure 9(a) that the distribution tendency of hydrogen in lattice along the crack contour is consistent; however, when the diffusion time increases from 1 × 10⁴ s to 2 × 10⁵ s, the hydrogen concentration increases slightly within ∼24° of the crack tip, but as hydrogen continues to diffuse to 4 × 10⁵ s and 1 × 10⁶ s, the hydrogen concentration decreases. The hydrogen concentration at the position larger than ∼24° always diminishes during the diffusion process, and the decrease is higher with the increased degree value. When the remote load increases to 400 MPa, it can be seen that the hydrogen concentration near the crack tip keeps decreasing from 1 × 10⁴ s to 1 × 10⁶ s, and a much lower value was obtained when diffusion time is 1 × 10⁶ s. The largest value appears at 2 × 10⁵ s, which takes short diffusion time, and the position is at ∼30°, which is farther from the crack tip than the remote load of 350 MPa. As the hydrogen concentrations in the lattice and traps are related to each other according to the models built, the same hydrogen distribution tendency in traps during the diffusion could be found in Figure 10, so no more detailed discussion here.

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

Figure 11 demonstrates the hydrogen in lattice in front of the crack tip at different diffusion stage. When the remote load is 350 MPa, it can be seen that, at the initial diffusion stage, the hydrogen concentration increases to the highest value with a very short distance of 0.001 m and then also quickly decreases to the initial concentration level within 0.01 m. With the continuous diffusion process, the concentration peak increases the highest value of 20.28 mol/m³ until 2 × 10⁵ s; the concentration peak will decrease when the diffusion goes on. When the diffusion time increases to 1 × 10⁶ s, even though there is also a crest at about 0.001 m, the concentration is smaller than the initial value, so the value decreases from the peak to trough and then increases to the initial hydrogen concentration with the increasing distance. As shown in Figure 11(b) with remote load increasing to 400 MPa, it can be seen that the load effects on the hydrogen diffusion process are dramatic. When the diffusion time is 1 × 10⁴ s, the highest hydrogen concentration of the peak is only 20.03 mol/m³; with the increasing diffusion time, the concentration peak disappears, and the hydrogen concentration monotonically increases to the initial level with the increasing distance in front of the crack tip. And the lower hydrogen concentration region at crack tip expands with the increasing diffusion time. The diffusion tendency of hydrogen in traps is the same as hydrogen in the lattice, and it is shown in Figure 12; the two variables are related to each other, so detailed discussion is ignored here.

6. Conclusion

A coupled mechanic-diffusion model is proposed to describe the stress effects on the diffusion of hydrogen in titanium. By incorporating a plate model with a central crack at different remote load conditions, it is found that the hydrogen diffusion near the crack is determined by the stress state. When the deformation is elastic at the crack tip, the hydrogen will diffuse from both sides of the crack towards the tip, and hydrogen concentration decreases with the increasing distance away from the crack tip. The hydrogen accumulation will be enhanced with the increased stress. When a plastic zone exists in front of the crack tip, the highest hydrogen concentration at crack surface deviates to the side near the crack tip, and the hydrogen concentration is much smaller than the elastic stage. Hydrogen will concentrate initially at a characterized distance in front of the crack tip, and the hydrogen concentration peak will diminish with the increasing diffusion time.

Data Availability

The data used to support the findings of this study were calculated according to the finite element method, and they are included in the article. The parameters used in the calculation model were cited from the references listed.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 51775427) and the Research Fund of State Key Laboratory for Marine Corrosion and Protection of Luoyang Ship Material Research Institute (LSMRI) under the contract no. 6142901180101.

References

D. Eliezer and T. H. Böllinghaus, “Hydrogen effects in titanium alloys,” in Gaseous Hydrogen Embrittlement of Materials in Energy Technologies, R. P. Gangloff and B. P. Somerday, Eds., pp. 668–706, Woodhead Publishing Ltd., Cambridge, UK, 2012.
View at: Google Scholar
D. A. Meyn, “Cleavage in Ti-8Al-1Mo-1V caused by hydrogen gas,” Metallurgical Transactions, vol. 3, no. 8, pp. 2302–2305, 1972.
View at: Publisher Site | Google Scholar
D. S. Shih, I. M. Robertson, and H. K. Birnbaum, “Hydrogen embrittlement of α titanium: in situ tem studies,” Acta Metallurgica, vol. 36, no. 1, pp. 111–124, 1988.
View at: Publisher Site | Google Scholar
R. A. Oriani and P. H. Josephic, “Equilibrium and kinetic studies of the hydrogen-assisted cracking of steel,” Acta Metallurgica, vol. 25, no. 9, pp. 979–988, 1977.
View at: Publisher Site | Google Scholar
D. F. Teter, I. M. Robertson, and H. K. Birnbaum, “The effects of hydrogen on the deformation and fracture of β-titanium,” Acta Materialia, vol. 49, no. 20, pp. 4313–4323, 2001.
View at: Publisher Site | Google Scholar
C. D. Beachem, “A new model for hydrogen-assisted cracking (hydrogen “embrittlement”),” Metallurgical and Materials Transactions B, vol. 3, no. 2, pp. 441–455, 1972.
View at: Publisher Site | Google Scholar
G. M. Bond, I. M. Robertson, and H. K. Birnbaum, “On the mechanisms of hydrogen embrittlement of Ni₃Al alloys,” Acta Metallurgica, vol. 37, no. 5, pp. 1407–1413, 1989.
View at: Publisher Site | Google Scholar
C. E. Johnson Jr. and F. A. Bovey, “Calculation of nuclear magnetic resonance spectra of aromatic hydrocarbons,” The Journal of Chemical Physics, vol. 29, no. 5, pp. 1012–1014, 1958.
View at: Publisher Site | Google Scholar
A. R. Troiano, “The role of hydrogen and other interstitials on the mechanical behavior of metals,” Trans ASM, vol. 52, pp. 54–80, 1960.
View at: Google Scholar
P. Sofronis and R. M. McMeeking, “Numerical analysis of hydrogen transport near a blunting crack tip,” Journal of the Mechanics and Physics of Solids, vol. 37, no. 3, pp. 317–350, 1989.
View at: Publisher Site | Google Scholar
E. Abramov and D. Eliezer, “Gas trapping and release in polycrystalline nickel preimplanted with helium,” Metallurgical and Materials Transactions A, vol. 25, no. 5, pp. 949–959, 1994.
View at: Publisher Site | Google Scholar
R. G. Macaulay-newcombe, M. Riehm, D. A. Thompson, W. W. Smeltzer, and E. Abramov, “Deuterium diffusion, trapping and release in ion-implanted nickel,” Radiation Effects and Defects in Solids, vol. 117, no. 4, pp. 285–297, 2006.
View at: Publisher Site | Google Scholar
A. McNabb and P. K. Foster, “A new analysis of the diffusion of hydrogen in iron and ferritic steels,” Transactions of the Metallurgical Society of AIME, vol. 227, pp. 618–626, 1963.
View at: Google Scholar
R. A. Oriani, “The diffusion and trapping of hydrogen in steel,” Acta Metallurgica, vol. 18, no. 1, pp. 147–157, 1970.
View at: Publisher Site | Google Scholar
T. C. Cui, P. F. Liu, and C. H. Gu, “Finite element analysis of hydrogen diffusion/plasticity coupled behaviors of low-alloy ferritic steel at large strain,” International Journal of Hydrogen Energy, vol. 42, no. 31, pp. 20324–20335, 2017.
View at: Publisher Site | Google Scholar
A. Taha and P. Sofronis, “A micromechanics approach to the study of hydrogen transport and embrittlement,” Engineering Fracture Mechanics, vol. 68, no. 6, pp. 803–837, 2001.
View at: Publisher Site | Google Scholar
T. Luo, C. Huang, and X. Gao, “An investigation of the effect of hydrogen on ductile fracture using a unit cell model,” International Journal of Hydrogen Energy, vol. 44, no. 16, pp. 8627–8640, 2019.
View at: Publisher Site | Google Scholar
G. Gobbi, C. Colombo, S. Miccoli, and L. Vergani, “A weakly coupled implementation of hydrogen embrittlement in FE analysis,” Finite Elements in Analysis and Design, vol. 141, pp. 17–25, 2018.
View at: Publisher Site | Google Scholar
G. Gobbi, C. Colombo, S. Miccoli, and L. Vergani, “A fully coupled implementation of hydrogen embrittlement in FE analysis,” Advances in Engineering Software, vol. 135, Article ID 102673, 2019.
View at: Publisher Site | Google Scholar
A. Díaz, J. M. Alegre, and I. I. Cuesta, “Coupled hydrogen diffusion simulation using a heat transfer analogy,” International Journal of Mechanical Sciences, vol. 115-116, pp. 360–369, 2016.
View at: Publisher Site | Google Scholar
M. A. Stopher, P. Lang, E. Kozeschnik, and P. E. J. Rivera-Diaz-del-Castillo, “Modelling hydrogen migration and trapping in steels,” Materials & Design, vol. 106, pp. 205–215, 2016.
View at: Publisher Site | Google Scholar
C. Huang, T. Luo, X. Gao, and S. M. Graham, “Modeling the effect of hydrogen on ductile fracture,” Materials Performance and Characterization, vol. 7, no. 2, pp. 101–118, 2018.
View at: Publisher Site | Google Scholar
T. Das, E. Legrand, S. V. Brahimi, J. Song, and S. Yue, “Evaluation of material susceptibility to hydrogen embrittlement (HE): an approach based on experimental and finite element (FE) analyses,” Engineering Fracture Mechanics, vol. 224, Article ID 106714, 2020.
View at: Publisher Site | Google Scholar
Y. Wang, X. Wu, Z. Zhou, and X. Li, “Numerical analysis of hydrogen transport into a steel after shot peening,” Results in Physics, vol. 11, pp. 5–16, 2018.
View at: Publisher Site | Google Scholar
V. Olden, A. Alvaro, and O. M. Akselsen, “Hydrogen diffusion and hydrogen influenced critical stress intensity in an API X70 pipeline steel welded joint—experiments and FE simulations,” International Journal of Hydrogen Energy, vol. 37, no. 15, pp. 11474–11486, 2012.
View at: Publisher Site | Google Scholar
S. Benannoune, Y. Charles, J. Mougenot, M. Gaspérini, and G. De Temmerman, “Numerical simulation by finite element modelling of diffusion and transient hydrogen trapping processes in plasma facing components,” Nuclear Materials and Energy, vol. 19, pp. 42–46, 2019.
View at: Publisher Site | Google Scholar
A. H. Javanmardi and H. Bashiri, “Analytical solution of Langmuir behavior of statistical rate theory: adsorption at solid/solution interface,” Journal of Environmental Chemical Engineering, vol. 5, no. 4, pp. 4024–4030, 2017.
View at: Publisher Site | Google Scholar
S. Azizian, “Kinetic models of sorption: a theoretical analysis,” Journal of Colloid and Interface Science, vol. 276, no. 11, pp. 47–52, 2004.
View at: Publisher Site | Google Scholar
W. Rudzinski and W. Plazinski, “Kinetics of solute adsorption at solid/solution interfaces: on the special features of the initial adsorption kinetics,” Langmuir, vol. 24, no. 13, pp. 6738–6744, 2008.
View at: Publisher Site | Google Scholar
J. Toribio and V. Kharin, “Review of hydrogen diffusion models for the analysis of hydrogen embrittlement of materials,” in Proceedings of the 13th International Conference on Fracture, Beijing, China, June 2013.
View at: Google Scholar
S. Benannoune, Y. Charles, J. Mougenot, and M. Gaspérini, “Numerical simulation of the transient hydrogen trapping process using an analytical approximation of the McNabb and Foster equation,” International Journal of Hydrogen Energy, vol. 43, no. 18, pp. 9083–9093, 2018.
View at: Publisher Site | Google Scholar
L. Chen, J. L. Fan, and H. R. Gong, “Atomistic simulation of mechanical properties of tungsten-hydrogen system and hydrogen diffusion in tungsten,” Solid State Communications, vol. 306, Article ID 113772, 2020.
View at: Publisher Site | Google Scholar
E. Torres, J. Pencer, and D. D. Radford, “Atomistic simulation study of the hydrogen diffusion in nickel,” Computational Materials Science, vol. 152, pp. 374–380, 2018.
View at: Publisher Site | Google Scholar
X. Wang, C. Xu, S. Hu, H. Xiao, X. Zhou, and S. Peng, “First-principles calculation on the structure stability, hydrogen trapping behaviour, and adhesion properties of the Zr(0001)ZrC(100) interface,” Applied Surface Science, vol. 508, Article ID 144825, 2020.
View at: Publisher Site | Google Scholar
X. Q. Ma, S. Q. Shi, C. H. Woo, and L. Q. Chen, “The phase field model for hydrogen diffusion and γ-hydride precipitation in zirconium under non-uniformly applied stress,” Mechanics of Materials, vol. 38, no. 1-2, pp. 3–10, 2006.
View at: Publisher Site | Google Scholar
A. H. M. Krom and A. Bakker, “Hydrogen trapping models in steel,” Metallurgical and Materials Transactions B, vol. 31, no. 6, pp. 1475–1482, 2000.
View at: Publisher Site | Google Scholar
P. Sofronis, Y. Liang, and N. Aravas, “Hydrogen induced shear localization of the plastic flow in metals and alloys,” European Journal of Mechanics—A/Solids, vol. 20, no. 6, pp. 857–872, 2001.
View at: Publisher Site | Google Scholar
W. Dietzel, M. Pfuff, and G. G. Juilfs, “Hydrogen permeation in plastically deformed steel membranes,” Materials Science, vol. 42, no. 1, pp. 78–84, 2006.
View at: Publisher Site | Google Scholar
A. J. Kumnick and H. H. Johnson, “Deep trapping states for hydrogen in deformed iron,” Acta Metallurgica, vol. 28, no. 1, pp. 33–39, 1980.
View at: Publisher Site | Google Scholar
A. N. Gorban, H. P. Sargsyan, and H. A. Wahab, “Quasichemical models of multicomponent nonlinear diffusion,” Mathematical Modelling of Natural Phenomena, vol. 6, no. 5, pp. 184–262, 2011.
View at: Publisher Site | Google Scholar
J. G. Wijmans and R. W. Baker, “The solution-diffusion model: a review,” Journal of Membrane Science, vol. 107, no. 1-2, pp. 1–21, 1995.
View at: Publisher Site | Google Scholar
A. G. McLellan, “Non-hydrostatic thermodynamics of chemical systems,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 314, no. 1518, pp. 443–455, 1970, https://search.crossref.org/?q=Non-hydrostatic+thermodynamics+of+chemical+systems.
View at: Publisher Site | Google Scholar
COMSOL, Corrosion Module User’s Guide v5.2.
G. M. Pressouyre and I. M. Bernstein, “A quantitative analysis of hydrogen trapping,” Metallurgical Transactions A, vol. 9, no. 11, pp. 1571–1580, 1978.
View at: Publisher Site | Google Scholar
G. M. Pressouyre and I. M. Bernstein, “An electrical analog model of hydrogen trapping in iron alloys,” Corrosion Science, vol. 18, no. 9, pp. 819–833, 1978.
View at: Publisher Site | Google Scholar
S. M. Teus, D. G. Savvakin, O. M. Ivasishin, and V. G. Gavriljuk, “Hydrogen migration and hydrogen-dislocation interaction in austenitic steels and titanium alloy in relation to hydrogen embrittlement,” International Journal of Hydrogen Energy, vol. 42, no. 4, pp. 2424–2433, 2017.
View at: Publisher Site | Google Scholar
B. B. Rath and I. M. Bernstein, “The relation between grain-boundary orientation and intergranular cracking,” Metallurgical Transactions, vol. 2, no. 10, pp. 2845–2851, 1971.
View at: Publisher Site | Google Scholar
I. C. Dragomir, D. S. Li, G. A. Castello-Branco et al., “Evolution of dislocation density and character in hot rolled titanium determined by X-ray diffraction,” Materials Characterization, vol. 55, no. 1, pp. 66–74, 2005.
View at: Publisher Site | Google Scholar
H. ALkhazraji, M. Z. Salih, Z. Zhong et al., “Estimation of dislocation density in cold-rolled commercially pure titanium by using synchrotron diffraction,” Metallurgical and Materials Transactions B, vol. 45, no. 4, pp. 1557–1564, 2014.
View at: Publisher Site | Google Scholar
H. Li, Z. Zhao, H. Guo, Y. Ning, and Z. Yao, “Dislocation density-based model for flow behavior of a near-α titanium alloy considering effects of initial lamellar thickness,” Journal of Materials Engineering and Performance, vol. 28, no. 4, pp. 2477–2487, 2019.
View at: Publisher Site | Google Scholar
R. Shi, Z. Nie, Q. Fan, F. Wang, Y. Zhou, and X. Liu, “Correlation between dislocation-density-based strain hardening and microstructural evolution in dual phase TC6 titanium alloy,” Materials Science and Engineering: A, vol. 715, pp. 101–107, 2018.
View at: Publisher Site | Google Scholar
J. L. Waisman, G. Sines, and L. B. Robinson, “Diffusion of hydrogen in titanium alloys due to composition, temperature, and stress graduates,” Metallurgical Transactions, vol. 4, no. 4, p. 1190, 1973.
View at: Publisher Site | Google Scholar
M. J. Donachie, Titanium: A Technical Guide, ASM International, Metals Park, OH, USA, 2nd edition, 2000.
J. L. Waisman, R. Toosky, and G. Sines, “Uphill diffusion and progressive embrittlement: hydrogen in titanium,” Metallurgical Transactions A, vol. 8, no. 8, pp. 1249–1256, 1977.
View at: Publisher Site | Google Scholar
V. Madina and I. Azkarate, “Compatibility of materials with hydrogen. Particular case: hydrogen embrittlement of titanium alloys,” International Journal of Hydrogen Energy, vol. 34, no. 14, pp. 5976–5980, 2009.
View at: Publisher Site | Google Scholar
B. G. Pound, “Hydrogen trapping in aged β-titanium alloys,” Acta Materialia, vol. 45, no. 5, pp. 2059–2068, 1997.
View at: Publisher Site | Google Scholar
L. G. Malyshev, S. A. Lebedev, R. A. Ryabov, P. V. Gel’d, and S. A. Kuznetsov, “Hydrogen diffusion in titanium at room temperature,” Soviet Physics Journal, vol. 25, no. 3, pp. 232–235, 1982.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2020 Fu-qiang Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1897

Downloads

1137

Citations