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

Authors: Fu-qiang Yang; Wen-juan Zhan; Tao Yan; Hai-bing Zhang; Xiu-rong Fang
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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 by(1)Ctot=CL+CT,where Ctot is the total hydrogen concentration in metals, mol·m−3; CL denotes the hydrogen concentration in the normal lattice, mol·m−3; and CT denotes the concentration associated with the hydrogen in the traps, mol·m−3. The two kinds of hydrogen concentration CL and CT are related to the number of sites:(2)CL=αθLNLNA,(3)CT=βθTNTNA,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; NA is Avogadro’s number, 6.023 × 1023 atom·mol−1; and NL is the number of atoms of solvent per volume, mol·m−3, which could be calculated as(4)NL=ρNAM=const,where ρ (g·m−3) and M (g·mol−1) are the density and relative atomic mass of the metal. And NT is the number of trapping sites per volume.The number of trapping sitesNT 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 NT could be calculated as [36, 37](5)NT=ρdb,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 3/2a and a, respectively; a is lattice constant, m; ρd is the dislocation density measured in dislocation line length per cubic meter, m·m−3, and it depends on the accumulated plastic strain as follows [36]:(6)ρd=ρ0+γεp,εp<0.5,1015,εp≥0.5,where ρ0 denotes the dislocation density for the annealed material, m−2, and γ is defined as the dislocation strengthening factor caused by plastic deformation.Although there are also other empirical relationships forNT − ε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 is(7)J=−DLRTCL∇μL,where J is the hydrogen flux, mol·m−2·s−1; DL is the diffusion coefficient for hydrogen, m2·s−1; μL is the chemical potential of the hydrogen in the lattice sites; R is the universal gas constant, 8.3144 J·mol−1·K−1; 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 as(8)μL=μL0+RTlnCL,where μ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:(9)μσh=μσh0−σhVh,where σh is the hydrostatic stress, which is calculated by the on the diagonal terms of the stress tensor, σh = ∑σii/3; Vh is the partial molar volume of hydrogen, m3·mol−1; μ(σ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 Vh.By incorporating equations (8) and (9), the chemical potential of hydrogen in a stressed metal could be described as(10)μL=μL0+RTlnCL−σhVh.Substituting equation (10) into equation (7), we can obtain the hydrogen flux expressed by the concentration of hydrogen as(11)J=−DL∇CL+DLRTCLVh∇σh.


### 2.3. Mass Balance Equation

The diffusion of hydrogen in titanium is defined from the requirement of mass conservation for the total hydrogen concentration:(12)∂∂t∫VCtotdV+∫∂VJ⋅ndS=0,where ∂/∂t is the partial derivative with respect to time; V is any volume whose surface is S, n denotes the outward normal to ∂V; J is the hydrogen flux through the surface ∂V, mol·m−2·s−1; and J·n is the hydrogen concentration flux leaving S.By substituting equations (1) and (11) into equation (13), we obtain(13)∂∂t∫VCL+CTdV+∫∂V−DL∇CL+DLRTCLVh∇σh⋅ndS=0.Applying the divergence theorem, equation (13) is expressed as(14)∂CL∂t+∂CT∂t+∇⋅−DL∇CL+∇⋅DLCLVhRT∇σh=0.This constitutive equation represents a modified form of the second Fick’s law; however, there are two unknows concentrationCL and CT. So, the relationship between CL and CT 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 fromL 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 NL >> NT, and there is a low occupancy θL << 1, the relationship between θL and θT could be acquired in the case of equilibrium:(15)θT1−θT=KTθL,where KT is the trap equilibrium constant, as(16)KT=exp−EbRT,where Eb is the trap binding energy concerning the lattice site; it is inherently negative energy [10], J·mol−1.By incorporating equations (2) and (15) into equation (3), the hydrogen concentration in trap sites could be expressed as(17)CT=NTNAθT=NTNAKTθLKTθL+1=NTNAKTθLKTθL+1=NTKTCLαNL+NAKTCL=N¯TKTCLαN¯L+KTCL,where N¯L and N¯T are defined as(18)N¯L=NLNA=const,N¯T=NTNA.Since the number of lattice sitesNL is a constant and the temperature will be kept constant, the hydrogen concentration in trap sites CT is a function of the hydrogen concentration in lattice sites CL, and the number of trap sites NT depends on the plastic deformation level described by equations (5) and (6).(19)CT=CTCL,N¯Tεp.So the partial derivative of the hydrogen concentration in trap sites with respect to time∂CT/∂t in equation (14) becomes(20)∂CT∂t=∂CT∂CL∂CL∂t+∂CT∂N¯TdN¯Tdεp∂εp∂t,where the partial derivative items ∂CT/∂CL and ∂CT/∂N¯T could be obtained from equation (17), as(21)∂CT∂CL=αKTN¯LN¯TαN¯L+KTCL2,∂CT∂N¯T=KTCLαN¯L+KTCL,and the item dN¯T/dεp could be induced by combining equations (5), (6), and (9), as(22)dN¯Tdεp=2γaNA,εp<0.5,0,εp≥0.5.Finally, by substituting equation (20) into equation (14), we obtain the mass constitutive equation as(23)1+∂CT∂CL∂CL∂t+∇⋅−DL∇CL+∇⋅DLCLVhRT∇σh+∂CT∂N¯TdN¯Tdεp∂εp∂t=0.


## 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 by(1)Ctot=CL+CT,where Ctot is the total hydrogen concentration in metals, mol·m−3; CL denotes the hydrogen concentration in the normal lattice, mol·m−3; and CT denotes the concentration associated with the hydrogen in the traps, mol·m−3. The two kinds of hydrogen concentration CL and CT are related to the number of sites:(2)CL=αθLNLNA,(3)CT=βθTNTNA,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; NA is Avogadro’s number, 6.023 × 1023 atom·mol−1; and NL is the number of atoms of solvent per volume, mol·m−3, which could be calculated as(4)NL=ρNAM=const,where ρ (g·m−3) and M (g·mol−1) are the density and relative atomic mass of the metal. And NT is the number of trapping sites per volume.The number of trapping sitesNT 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 NT could be calculated as [36, 37](5)NT=ρdb,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 3/2a and a, respectively; a is lattice constant, m; ρd is the dislocation density measured in dislocation line length per cubic meter, m·m−3, and it depends on the accumulated plastic strain as follows [36]:(6)ρd=ρ0+γεp,εp<0.5,1015,εp≥0.5,where ρ0 denotes the dislocation density for the annealed material, m−2, and γ is defined as the dislocation strengthening factor caused by plastic deformation.Although there are also other empirical relationships forNT − ε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 is(7)J=−DLRTCL∇μL,where J is the hydrogen flux, mol·m−2·s−1; DL is the diffusion coefficient for hydrogen, m2·s−1; μL is the chemical potential of the hydrogen in the lattice sites; R is the universal gas constant, 8.3144 J·mol−1·K−1; 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 as(8)μL=μL0+RTlnCL,where μ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:(9)μσh=μσh0−σhVh,where σh is the hydrostatic stress, which is calculated by the on the diagonal terms of the stress tensor, σh = ∑σii/3; Vh is the partial molar volume of hydrogen, m3·mol−1; μ(σ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 Vh.By incorporating equations (8) and (9), the chemical potential of hydrogen in a stressed metal could be described as(10)μL=μL0+RTlnCL−σhVh.Substituting equation (10) into equation (7), we can obtain the hydrogen flux expressed by the concentration of hydrogen as(11)J=−DL∇CL+DLRTCLVh∇σh.


## 2.3. Mass Balance Equation

The diffusion of hydrogen in titanium is defined from the requirement of mass conservation for the total hydrogen concentration:(12)∂∂t∫VCtotdV+∫∂VJ⋅ndS=0,where ∂/∂t is the partial derivative with respect to time; V is any volume whose surface is S, n denotes the outward normal to ∂V; J is the hydrogen flux through the surface ∂V, mol·m−2·s−1; and J·n is the hydrogen concentration flux leaving S.By substituting equations (1) and (11) into equation (13), we obtain(13)∂∂t∫VCL+CTdV+∫∂V−DL∇CL+DLRTCLVh∇σh⋅ndS=0.Applying the divergence theorem, equation (13) is expressed as(14)∂CL∂t+∂CT∂t+∇⋅−DL∇CL+∇⋅DLCLVhRT∇σh=0.This constitutive equation represents a modified form of the second Fick’s law; however, there are two unknows concentrationCL and CT. So, the relationship between CL and CT 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 fromL 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 NL >> NT, and there is a low occupancy θL << 1, the relationship between θL and θT could be acquired in the case of equilibrium:(15)θT1−θT=KTθL,where KT is the trap equilibrium constant, as(16)KT=exp−EbRT,where Eb is the trap binding energy concerning the lattice site; it is inherently negative energy [10], J·mol−1.By incorporating equations (2) and (15) into equation (3), the hydrogen concentration in trap sites could be expressed as(17)CT=NTNAθT=NTNAKTθLKTθL+1=NTNAKTθLKTθL+1=NTKTCLαNL+NAKTCL=N¯TKTCLαN¯L+KTCL,where N¯L and N¯T are defined as(18)N¯L=NLNA=const,N¯T=NTNA.Since the number of lattice sitesNL is a constant and the temperature will be kept constant, the hydrogen concentration in trap sites CT is a function of the hydrogen concentration in lattice sites CL, and the number of trap sites NT depends on the plastic deformation level described by equations (5) and (6).(19)CT=CTCL,N¯Tεp.So the partial derivative of the hydrogen concentration in trap sites with respect to time∂CT/∂t in equation (14) becomes(20)∂CT∂t=∂CT∂CL∂CL∂t+∂CT∂N¯TdN¯Tdεp∂εp∂t,where the partial derivative items ∂CT/∂CL and ∂CT/∂N¯T could be obtained from equation (17), as(21)∂CT∂CL=αKTN¯LN¯TαN¯L+KTCL2,∂CT∂N¯T=KTCLαN¯L+KTCL,and the item dN¯T/dεp could be induced by combining equations (5), (6), and (9), as(22)dN¯Tdεp=2γaNA,εp<0.5,0,εp≥0.5.Finally, by substituting equation (20) into equation (14), we obtain the mass constitutive equation as(23)1+∂CT∂CL∂CL∂t+∇⋅−DL∇CL+∇⋅DLCLVhRT∇σh+∂CT∂N¯TdN¯Tdεp∂εp∂t=0.


## 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](24)∂c∂t+∇⋅−D∇c+u⋅∇c=Ra,where c is the concentration of the species, mol·m−3; D denotes the diffusion coefficient m2·s−1; Ra is a reaction rate expression for the species, mol·m−3·s−1; and u is the velocity vector, m·s−1.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, Ra 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:(25)∂CL∂t+∇⋅−DLda∇CL+DLVhRTda∇σh⋅∇CL=−1daDLVhRT∇2σhCL+∂CT∂N¯TdN¯Tdεp∂εp∂t,where da is mass coefficient defined as(26)da=1+∂CT∂CL=1+αKTN¯LN¯TαN¯L+KTCL2.Define effective diffusion coefficientDeff, convection velocity u, and reaction rate Ra as(27)Deff=DLda,u=DLVhRTda∇σh,Ra=−1daDLVhRT∇2σhCL+∂CT∂N¯TdN¯Tdεp∂εp∂t.It is found in equation (25) that the hydrogen transport mechanism in stressed metals also contains diffusion and convection. The effective diffusion coefficient Deff is a variable, which is a function of the hydrogen concentration in lattice sites CL and the number of trapping sites N¯T 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 Ra 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 r1. To simplify the calculation, only a quarter of the model was considered. In this simulation, W and H are equal to 60 mm, r1 = 4 mm.Figure 1
A quarter of the plate model with hole crack: (a) geometric model, (b) mesh of the whole model, and (c) mesh of the crack tip.
(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 Figure1, the three regions are divided by two arcs with a radius of r2 = 10 mm and r3 = 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]:(28)σY=σ01+εpε01/n,where σY (MPa) is the local flow stress with the plastic strain of εp; σ0 (MPa), ε0, and n are the yield stress, yield strain, and hardening exponent measured in the absence of hydrogen, respectively, and σ0· = ε0 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 concentrationCL = C0, 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 Figure1 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 (Eb)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 Eb 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−1) for dislocation trapping according to Rath and Bernstein [46].(2)
Dislocation densityWith the deformation at 40%, 60%, and 80% reduction levels in commercially pure titanium specimens, the dislocations densities calculated are 5 × 1014 ± 5% m−2, 8 × 1014 ± 10% m−2, and 1015 ± 7% m−2, respectively [47], while the dislocation densities are 2.1 × 1014 m−2, 3.58 × 1014 m−2, and 4 × 1014 m−2 for vertical rolling direction and 1.15 × 1014 m−2, 2 × 1014 m−2, and 3 × 1014 m−2 for horizontal rolling direction with 40%, 60%, and 80% cold rolling [48]. With a strain rate of 0.01 s−1, 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 × 1013 m−2 at 1233 K [49]. For the α + β type dual-phase TC6 titanium alloy, α and β phases have the pre-existing dislocations of 1.17 × 1014 m−2 and 2.45 × 1014 m−2, respectively. During the deformation, the dislocation density of the α phase increases significantly from 1.88 × 1014 m−2 (ε = 0.73%) to 3.36 × 1014 m−2 (ε = 1.02%) and then plateaus until ε = 1.30%. The dislocation density of the β phase increases slightly from 3.36 × 1014 m−2 to 3.79 × 1014 m−2 at first and then increases sharply to 4.87 × 1014 m−2 [50]. Even though the dislocation density is related to plastic degree, material type, and strain rate, the magnitude of the dislocation density is 1014 m−2, so the initial dislocation density, ρ0, is considered as 1014 m−2, and the dislocation density keeps constant as 1015 m−2; when plastic strain εp is larger than 0.5, then the dislocation strengthening factor γ is 1.8 × 1015 m−2.(3)
Partial molar volume of hydrogen (Vh)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 calculateVh, which gives an independent determination of the partial molar volume by the following equation [51]:(29)Vh=3Vs1−NHΔe11ΔNH,whereVs 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 NH and ΔNH are the atom percent of hydrogen and its change with a strain Δe11.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−6 m3·mol−1.(4)
DensityThe elemental titanium’s density is nominally computed on the basis of 4.512 g·cm−3, 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−3, 4.624 g·cm−3, and 4.936 g·cm−3, respectively. The relative atomic mass of Ti is 47.9 g·mol−1, and it is used instead of the values of different titanium alloys.(5)
Number of interstitial sitesThe 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 diffusivityMany 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−16 m2·s−1 for α phase titanium and 3.3 × 10−11 m2·s−1 for the β-phase titanium [54]. Waisman et al. have summarized that a diffusivity of 4.95 × 10−13 m2·s−1 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−11 m2·s−1 [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−10 m2·s−1 to 2.7 × 10−9 m2 s−1 [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.Table 1
Physical and mechanical properties of the titanium alloy.
Parameters
Values
Density,ρ (g·cm−3)
4.624
Number of interstitial sites per solvent (Ti) atom
6
Elastic modulus,E (GPa)
100
Poisson’s ratio,υ
0.34
Yield strength,σ0 (MPa)
500
Yield strain,ε0 (%)
0.5
Hardening exponent,n
5
Lattice parameter,a (10−9 m)
0.295
Lattice hydrogen diffusivity,DL (m2·s−1)
4.95 × 10−13
Partial molar volume of hydrogen,Vh (m3·mol−1)
1.8 × 10−6
Initial dislocation density,ρ0 (m−2)
1014
Dislocation strengthen factor,γ (m−2)
1.8 × 1015
Trap binding energy for hydrogen to dislocation,Eb (kJ·mol−1)
26.05
Relative atomic mass,M (g·mol−1)
47.9
Avogadro’s number,NA (atom·mol−1)
6.023 × 1023
Universal gas constant,R (J·mol−1·K−1)
8.3144
Absolute temperature,T (K)
273


### 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 latticesCL 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 CT, 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.


## 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 r1. To simplify the calculation, only a quarter of the model was considered. In this simulation, W and H are equal to 60 mm, r1 = 4 mm.Figure 1
A quarter of the plate model with hole crack: (a) geometric model, (b) mesh of the whole model, and (c) mesh of the crack tip.
(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 Figure1, the three regions are divided by two arcs with a radius of r2 = 10 mm and r3 = 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]:(28)σY=σ01+εpε01/n,where σY (MPa) is the local flow stress with the plastic strain of εp; σ0 (MPa), ε0, and n are the yield stress, yield strain, and hardening exponent measured in the absence of hydrogen, respectively, and σ0· = ε0 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 concentrationCL = C0, 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 Figure1 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 (Eb)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 Eb 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−1) for dislocation trapping according to Rath and Bernstein [46].(2)
Dislocation densityWith the deformation at 40%, 60%, and 80% reduction levels in commercially pure titanium specimens, the dislocations densities calculated are 5 × 1014 ± 5% m−2, 8 × 1014 ± 10% m−2, and 1015 ± 7% m−2, respectively [47], while the dislocation densities are 2.1 × 1014 m−2, 3.58 × 1014 m−2, and 4 × 1014 m−2 for vertical rolling direction and 1.15 × 1014 m−2, 2 × 1014 m−2, and 3 × 1014 m−2 for horizontal rolling direction with 40%, 60%, and 80% cold rolling [48]. With a strain rate of 0.01 s−1, 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 × 1013 m−2 at 1233 K [49]. For the α + β type dual-phase TC6 titanium alloy, α and β phases have the pre-existing dislocations of 1.17 × 1014 m−2 and 2.45 × 1014 m−2, respectively. During the deformation, the dislocation density of the α phase increases significantly from 1.88 × 1014 m−2 (ε = 0.73%) to 3.36 × 1014 m−2 (ε = 1.02%) and then plateaus until ε = 1.30%. The dislocation density of the β phase increases slightly from 3.36 × 1014 m−2 to 3.79 × 1014 m−2 at first and then increases sharply to 4.87 × 1014 m−2 [50]. Even though the dislocation density is related to plastic degree, material type, and strain rate, the magnitude of the dislocation density is 1014 m−2, so the initial dislocation density, ρ0, is considered as 1014 m−2, and the dislocation density keeps constant as 1015 m−2; when plastic strain εp is larger than 0.5, then the dislocation strengthening factor γ is 1.8 × 1015 m−2.(3)
Partial molar volume of hydrogen (Vh)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 calculateVh, which gives an independent determination of the partial molar volume by the following equation [51]:(29)Vh=3Vs1−NHΔe11ΔNH,whereVs 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 NH and ΔNH are the atom percent of hydrogen and its change with a strain Δe11.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−6 m3·mol−1.(4)
DensityThe elemental titanium’s density is nominally computed on the basis of 4.512 g·cm−3, 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−3, 4.624 g·cm−3, and 4.936 g·cm−3, respectively. The relative atomic mass of Ti is 47.9 g·mol−1, and it is used instead of the values of different titanium alloys.(5)
Number of interstitial sitesThe 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 diffusivityMany 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−16 m2·s−1 for α phase titanium and 3.3 × 10−11 m2·s−1 for the β-phase titanium [54]. Waisman et al. have summarized that a diffusivity of 4.95 × 10−13 m2·s−1 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−11 m2·s−1 [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−10 m2·s−1 to 2.7 × 10−9 m2 s−1 [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.Table 1
Physical and mechanical properties of the titanium alloy.
Parameters
Values
Density,ρ (g·cm−3)
4.624
Number of interstitial sites per solvent (Ti) atom
6
Elastic modulus,E (GPa)
100
Poisson’s ratio,υ
0.34
Yield strength,σ0 (MPa)
500
Yield strain,ε0 (%)
0.5
Hardening exponent,n
5
Lattice parameter,a (10−9 m)
0.295
Lattice hydrogen diffusivity,DL (m2·s−1)
4.95 × 10−13
Partial molar volume of hydrogen,Vh (m3·mol−1)
1.8 × 10−6
Initial dislocation density,ρ0 (m−2)
1014
Dislocation strengthen factor,γ (m−2)
1.8 × 1015
Trap binding energy for hydrogen to dislocation,Eb (kJ·mol−1)
26.05
Relative atomic mass,M (g·mol−1)
47.9
Avogadro’s number,NA (atom·mol−1)
6.023 × 1023
Universal gas constant,R (J·mol−1·K−1)
8.3144
Absolute temperature,T (K)
273


## 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 latticesCL 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 CT, 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/m3 during the diffusion, which is considered as the hydrogen source. The initial hydrogen in the plate is still uniformly distributed of 20 mol/m3. Without the effects of stress, it is assumed that the hydrogen only diffuses in the lattice, and there is no interaction between CL and CX, which indicates that(30)∂CT∂CL=0.Without the effects of stress and trap, equation (25) could be simplified as(31)∂CL∂t+∇⋅−DL∇CL=0.For the diffusion in a semi-infinite plate, the analytic solution reads(32)CLx,t=Cb+C0−Cberfx2DLt,where Cb and C0 are the hydrogen concentration in boundary and the plate, mol/m3. x is the distance away from the boundary, mm. CL(x, t) is the hydrogen concentration at position x at time t. erf is error function with the formation as(33)erfΨ=2π∫0Ψe−Ψ2dΨ.Figure2 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).Figure 2
Hydrogen concentration considering without trap and stress effect, comparison between analytic solution and numerical results.


### 5.2. Coupled Diffusion-Mechanics under Constant Load

Figures3 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/m3 at the crack tip, and the value increased to 25 mol/m3 when remote load increases to 300 MPa. The hydrogen concentration in traps increases from 9 × 10−3 mol/m3 to 1.13 × 10−2 mol/m3 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/m3, and the hydrogen in traps is 5.81 × 10−3 mol/m3, and the values decrease with the increased remote load.Figure 3
Distribution of hydrogen adopted in lattice near the tip, mol/m3: (a) σy = 200 MPa, (b) σy = 250 MPa, (c) σy = 300 MPa, (d) σy = 350 MPa, and (e) σy = 400 MPa.
(a)
(b)
(c)
(d)
(e)Figure 4
Distribution of hydrogen adopted in traps near the tip, mol/m3: (a) σy = 200 MPa, (b) σy = 250 MPa, (c) σy = 300 MPa, (d) σy = 350 MPa, and (e) σy = 400 MPa.
(a)
(b)
(c)
(d)
(e)A more detailed distribution of hydrogen in normal lattice and traps is shown in Figures5 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.Figure 5
Distribution of hydrogen in (a) normal lattice and (b) traps around the crack contour.
(a)
(b)Figure 6
Distribution of hydrogen in (a) normal lattice and (b) traps in front of the crack tip.
(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 CL and CT are higher than the initial value of 20 mol/m3 and 5.81 × 10−3 mol/m3. This indicates that the area in front of the crack tip is the region of hydrogen inflow, even a negative value of ∂σh/∂x. 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).Figure 7
Distribution of hydrostatic stress (a) in front of the crack tip and (b) around the crack contour.
(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 Figure8, 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.Figure 8
Distribution of Mises stress (a) in front of the crack tip and (b) around the crack contour.
(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 Figures9–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 × 104 s to 2 × 105 s, the hydrogen concentration increases slightly within ∼24° of the crack tip, but as hydrogen continues to diffuse to 4 × 105 s and 1 × 106 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 × 104 s to 1 × 106 s, and a much lower value was obtained when diffusion time is 1 × 106 s. The largest value appears at 2 × 105 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.Figure 9
Hydrogen in lattice around crack contour at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure 10
Hydrogen in traps around crack contour at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure 11
Hydrogen in lattice in front of the crack tip at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure 12
Hydrogen in traps in front of the crack tip at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure11 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/m3 until 2 × 105 s; the concentration peak will decrease when the diffusion goes on. When the diffusion time increases to 1 × 106 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 × 104 s, the highest hydrogen concentration of the peak is only 20.03 mol/m3; 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.


## 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/m3 during the diffusion, which is considered as the hydrogen source. The initial hydrogen in the plate is still uniformly distributed of 20 mol/m3. Without the effects of stress, it is assumed that the hydrogen only diffuses in the lattice, and there is no interaction between CL and CX, which indicates that(30)∂CT∂CL=0.Without the effects of stress and trap, equation (25) could be simplified as(31)∂CL∂t+∇⋅−DL∇CL=0.For the diffusion in a semi-infinite plate, the analytic solution reads(32)CLx,t=Cb+C0−Cberfx2DLt,where Cb and C0 are the hydrogen concentration in boundary and the plate, mol/m3. x is the distance away from the boundary, mm. CL(x, t) is the hydrogen concentration at position x at time t. erf is error function with the formation as(33)erfΨ=2π∫0Ψe−Ψ2dΨ.Figure2 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).Figure 2
Hydrogen concentration considering without trap and stress effect, comparison between analytic solution and numerical results.


## 5.2. Coupled Diffusion-Mechanics under Constant Load

Figures3 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/m3 at the crack tip, and the value increased to 25 mol/m3 when remote load increases to 300 MPa. The hydrogen concentration in traps increases from 9 × 10−3 mol/m3 to 1.13 × 10−2 mol/m3 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/m3, and the hydrogen in traps is 5.81 × 10−3 mol/m3, and the values decrease with the increased remote load.Figure 3
Distribution of hydrogen adopted in lattice near the tip, mol/m3: (a) σy = 200 MPa, (b) σy = 250 MPa, (c) σy = 300 MPa, (d) σy = 350 MPa, and (e) σy = 400 MPa.
(a)
(b)
(c)
(d)
(e)Figure 4
Distribution of hydrogen adopted in traps near the tip, mol/m3: (a) σy = 200 MPa, (b) σy = 250 MPa, (c) σy = 300 MPa, (d) σy = 350 MPa, and (e) σy = 400 MPa.
(a)
(b)
(c)
(d)
(e)A more detailed distribution of hydrogen in normal lattice and traps is shown in Figures5 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.Figure 5
Distribution of hydrogen in (a) normal lattice and (b) traps around the crack contour.
(a)
(b)Figure 6
Distribution of hydrogen in (a) normal lattice and (b) traps in front of the crack tip.
(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 CL and CT are higher than the initial value of 20 mol/m3 and 5.81 × 10−3 mol/m3. This indicates that the area in front of the crack tip is the region of hydrogen inflow, even a negative value of ∂σh/∂x. 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).Figure 7
Distribution of hydrostatic stress (a) in front of the crack tip and (b) around the crack contour.
(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 Figure8, 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.Figure 8
Distribution of Mises stress (a) in front of the crack tip and (b) around the crack contour.
(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 Figures9–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 × 104 s to 2 × 105 s, the hydrogen concentration increases slightly within ∼24° of the crack tip, but as hydrogen continues to diffuse to 4 × 105 s and 1 × 106 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 × 104 s to 1 × 106 s, and a much lower value was obtained when diffusion time is 1 × 106 s. The largest value appears at 2 × 105 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.Figure 9
Hydrogen in lattice around crack contour at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure 10
Hydrogen in traps around crack contour at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure 11
Hydrogen in lattice in front of the crack tip at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure 12
Hydrogen in traps in front of the crack tip at different diffusion stages: (a)σy = 350 MPa and (b) σy = 400 MPa.
(a)
(b)Figure11 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/m3 until 2 × 105 s; the concentration peak will decrease when the diffusion goes on. When the diffusion time increases to 1 × 106 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 × 104 s, the highest hydrogen concentration of the peak is only 20.03 mol/m3; 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.



---

*Source: 3618589-2020-04-25.xml*