Modeling hydrogen diffusion in precipitation hardened nickel-based alloy 718 by microstructural modeling

1 Introduction

Nickel-based alloy UNS N07718 is a high strength, precipitation hardened alloy extensively used in marine environments due to its combination of high general corrosion resistance, good strength, and high resistance to chloride and sulfide stress corrosion cracking (Chen et al. 2014, 2015). The strength of the material is given primarily by precipitation hardening of three intermetallic phases, namely γ′, Ni3(Al,Ti) with L1₂ structure, γ″, a metastable Ni3Nb phase with DO₂₂ structure, and δ phase, the stable counterpart of the γ″ phase with orthorhombic structure. Because of the material heterogenous microstructure and the role of the precipitates in both deformation and environmental assisted cracking behavior, controlled chemistry, solution annealing, and heat treatments are required to obtain a material resistant to hydrogen producing marine environment. Variations in material microstructure are a function of the manufacturing process (Saleem et al. 2020) and minimizing the presence of δ phase in the microstructure is shown to be important in reducing the susceptibility to intergranular cracking. The most commonly used condition in marine applications, particularly in the Oil and Gas sector, is produced to meet the American Petroleum Institute (API) standard 6A 718 specification which, requires a solution annealing temperature between 1870 and 1925 °F, above the δ solvus, and includes a single heat treatment step at a temperature between 1425 and 1475 °F for 6–8 h to achieve a minimum yield strength of 120 ksi (2019). The API specification also requires compliance with NACE MR0175 for all cases involving exposure to H₂S, which limits the maximum hardness to Rockwell Hardness (HRC) 40. The resulting microstructure was shown to have minimal δ precipitation at the grain boundaries and the environmentally assisted cracking (EAC) behavior to be among the best of similar PH Ni-based alloys (Thodla 2018; Thodla et al. 2020).

The use of UNS N07718 (API 6A) in naval applications is increasing, particularly as a fastener material, due its resistance to EAC and hydrogen embrittlement (HE) under cathodic protection systems that maintain an applied potential equal to or more electronegative than −850 mV versus Ag/AgCl (Esaklul and Ahmed 2009; Thodla 2018). The API specification includes two additional grades for 718 with increased yield strength up to a maximum of 150 ksi. EAC data on the latter grades are being investigated with interest as the potential impact of an increase in yield by 25 % is significant (Li et al. 2019). There is considerable interest in understanding the role of microstructure in the EAC performance of the material, with the potential of further refining the heat treatment to achieve even better performance at higher strengths and because of the significant use of 718 in additive manufacturing (Hosseini and Popovich 2019). We have shown in previous investigation that the matrix-precipitate interactions are important in determining the material elastic-plastic response causing increased shear localization and partial dislocation activity at the δ-phase and matrix interface, and activation of specific slip systems and dislocation-density interactions (Arcari et al. 2023). Mechanistically, shear strain localization is considered a strong driving force for hydrogen assisted cracking (Obasi et al. 2018; Thodla 2018). Trapped hydrogen can be more easily transported by dislocations and therefore might accumulate at locations with highly localized strain and at grain boundaries causing progressive degradation of the interfacial strength, leading to decohesion. Additionally, previous microstructural modeling for 718 showed that different slip systems are activated at the precipitate-matrix interface with respect to the main activated slip systems at the grain boundary interior (Arcari et al. 2023). The interaction of dislocations on different slip systems near the interface of the precipitate with the matrix is another potential driving force for increased hydrogen uptake in the vicinity of the precipitates.

In this study we developed a microstructural model representative of UNS N07718 alloy and define the local stress-strain conditions at the tip of a notch of a fracture mechanics sample to calculate the driving forces for hydrogen diffusion within the material microstructure. The driving forces coming from hydrostatic stress and mobile dislocation densities were combined to provide an assessment of hydrogen concentration ahead of an initiated crack at the root of the notch. The role of immobile dislocations as irreversible traps and the changes in hydrogen concentration trends were also identified.

2 Modeling methods and boundary conditions

A global model representative of typical fracture mechanics geometry was developed using the finite element (FE) method to best approximate the boundary conditions of a detailed microstructural model at the crack tip. The model is a 2D plane-strain simulation of a wedge-opening load sample with a sharp notch with a root radius, ρ = 0.5 mm. The load is applied as a surface pressure at the bottom of the sample as indicated in Figure 1a and at a normalized distance of 0.25W from the edge, where W is the width of the sample. The notch length is 0.52 × a/W, where a is crack depth, the total length of the sample is 1.25W and the height of the sample is 1.2W. The geometry is consistent with a typical compact tension C(T) sample, and used for environmental assisted cracking (EAC) assessment of structural materials (ASTM International 2023). The model is subjected to an applied load such that the stress intensity factor (SIF) at the notch root K is 88 MPa√m, assuming the sharp notch has a microstructurally small crack at its root.

Figure 1:

Finite element model of fracture mechanics sample and resulting maximum principal stress.

The model uses a Ramberg–Osgood formulation to describe the material stress-strain behavior:

where σ is Cauchy stress, ε is mechanical strain, E is Young’s modulus, α is the “yield” offset defined as: εyEσy−1, at the condition of yielding, σy is yield strength, and n is the hardening exponent, fitted to the material stress-strain behavior. The API specification gives different minimum strengths depending on the grade, but in this work we refer to the specification with the minimum yield of 120 ksi. Commercial forms produced according to this specification generally have higher yield strength than the minimum. In this work σy=135 ksi, α = 0.5, and n = 9.0, obtained from fitting experimental stress-strain data from the literature (Thodla et al. 2020). The material response according to the Ramberg–Osgood formulation with the linear elastic portion and the 0.2 % offset is shown in Figure 1b.

The contour of the maximum principal stress is shown in Figure 1a. The plastic zone defined as the area with stress above 2σy , reflecting the plain strain constraints, is shown in red. Figure 2a shows the resulting normalized stress gradient σyyσy as a function of distance from the sample notch. The peak is at approximately 0.5 mm from the notch tip, and corresponds to 2.25 times the yield strength for the material. The stress rapidly decreases as a function of distance from the notch root. Figure 2b shows the distribution of strain, εyy, as a function of distance from the notch tip. Strain exceeds 4 % at the notch tip and rapidly decreases to less than 0.2 %.

Figure 2:

Stress, σyy, and strain, εyy, as a function of distance from the notch tip.

The results of the analysis are considered the initial mechanical condition for the sample before a crack initiates from the notch (Vasudevan 2013). This condition also holds when the crack is short when compared to the notch root radius (Suresh 1998). As the crack grows beyond the influence zone of the notch, its mechanical response tends to the ideal sharp crack solution. The results are used as boundary conditions for a local model, representing the microstructure of the material.

3 Local model and assumptions

A local model was developed to better represent the stress-strain condition at the notch when a crack initiates and then propagates as driven by environmental driving forces. The model developed in this work is intended to provide a more in-depth understanding of the effects of local microstructure in promoting the diffusion of hydrogen towards a fracture process zone (FPZ). The mechanical driving force distribution within the material microstructure is used as direct input for a simulation of diffusion of hydrogen. In this work we focus on hydrogen assisted cracking in a regime where the crack growth rates are higher for increasingly electronegative potentials, such as −1.0 V versus Ag/AgCl (Nagao et al. 2012; Thodla et al. 2020). In this regime, the hydrogen at the wake and tip of an initiated crack diffuses towards the fracture process zone, leading to crack propagation. While there is still significant debate over the role of hydrogen in EAC for this material, experimental data for UNS N07718-API alloy points towards the significance of the diffusion of hydrogen through the material microstructure as the rate limiting step for crack growth progression (Gangloff 2003). More recent studies looked at the role of the material microstructure in hydrogen assisted cracking of 718 alloy, identifying the importance of grain size, grain boundary characteristics, and dislocation-precipitates interactions (Jothi et al. 2016; Khalid et al. 2024; Ogawa et al. 2022).

The grain structure was reproduced by using Voronoi tessellation: a set of 18 grains in a 100 × 100 µm microstructure in a two-dimensional plane-strain finite element analysis. The number of grains in the representative area was defined based on an ASTM grain size finer than 2.0, as required by the API specification and literature on the material. The major precipitates are explicitly modeled within the material microstructure. The δ phase was modeled at the grain boundaries as an ellipsoidal phase with a major axis of 1.0 µm, and an aspect ratio of 2.0. The number of δ phases on each grain boundary is maximized at 2.0, for grain boundaries that were long enough to accommodate them. Their number on different grain boundaries also varied as to avoid intersections with other δ phases or parts of the microstructure.

The γ′ and γ″ phases were modeled at the interior of the grain as circular precipitates with a diameter of 0.1 µm and as elliptical precipitates with major axis of 0.2 µm and aspect ratio of 2. These values are representative of an upper bound limit for the expected size of these phases for materials with comparable thermal history (Obasi et al. 2018; Saleem et al. 2020; Slama and Abdellaoui 2000). The γ′ and γ″ phases were randomly distributed at the interior of each grain with constraints dictated by distance from the grain boundary, to avoid intersection with δ phase or grain boundaries, and by distance to the other γ′ and γ″ phase particles.

A crystal plasticity formulation was implemented into a user material subroutine (UMAT) in Abaqus based on the work by Huang (Huang 1991) to model the mechanical response of the metal grains. Before plastic deformations occur, the FCC γ-matrix of the material is modeled as an elastic crystal with constants: C₁₁ = 233 GPa, C₁₂ = 147 GPa, C₄₄ = 112 GPa (Gupta and Bronkhorst 2021). The formulation accounts for shear strain deformation on all 12 slip systems when plastic deformation occurs. The constitutive equations were modified to account for the distribution of mobile and immobile dislocations within the microstructure. Following the work by Zikry (Wu and Zikry 2015; Zikry and Kao 1996), the strength of each slip system was modeled as a function of the mobile and immobile dislocation densities developed upon plastic deformations. The main step of the implementation of the theory is in the following evolution equations and their coupling to the shear strain rate and strength equations as:

where

and γ˙(α) is the shear strain rate and γ˙ref(α) is a reference shear strain rate, τ(α) is the resolved shear stress, τy(α) is the initial strength, m is a rate sensitivity parameter, G is the shear modulus, b(β) is the Burgers vector, and all refer to the slip system α. A value of τy(α) = 330 MPa was used based on previous work with this alloy system (Arcari et al. 2023), and the γ˙ref(α) and m values were determined based on an optimization approach to give the best match with the material yield strength. The selected values were respectively γ˙ref(α) = 0.01 and m = 10. The variable nss is the number of slip systems, aαβ represents the influence of dislocation accumulation along a β slip system, 0.2 for cross-coupling and 1.0 for self-coupling, ρim(β) is the immobile dislocation density for the β slip system. The evolutions of the mobile and immobile dislocation densities are coupled through a separate set of differential equations (Wu and Zikry 2015; Zikry and Kao 1996).

where gsourceα is a coefficient pertaining to the increase in mobile dislocation density due to dislocation sources, gmnterα is related to the trapping of mobile dislocations due to forest intersections or cross-slip around obstacles or dislocation interactions, gimmobα is related to the immobilization of mobile dislocations, grecovα is related to the rearrangement and annihilation of immobile dislocations. The complete formulation of all of the coefficients was developed by Zikry et al. (Wu and Zikry 2015; Zikry and Kao 1996). The dislocation-based crystal plasticity parameters controlling the evolution of mobile and immobile dislocation densities for Nickel alloy were assumed to be consistent with the ones obtained for a 718 material model (Arcari et al. 2023). The calibration of the several model parameters is a significant challenge for the development of a crystal plasticity formulation that can best represent the stress-strain material behavior. In this work we implemented an optimization methodology to identify the τy(α), γ˙ref(α), and m parameters that maximize the correlation with the yield strength of the material. The simulation is performed for a model of grain aggregates subjected to plane stress conditions.

The δ phase was modeled as an ordered tetragonal crystal with elastic constants C₁₁ = 264 GPa, C₁₂ = 175 GPa, C₁₃ = 155 GPa, C₃₃ = 261 GPa, C₄₄ = 142 GPa, C₆₆ = 160 GPa (Xiong et al. 2022). The material orientation of the δ phase was simplified as the average of the Euler angles of the two adjacent grains. The γ′ phase was modeled as an elastic crystalline material with elastic constants C₁₁ = 224 GPa, C₁₂ = 147 GPa, C₁₃ = 147 GPa, C₃₃ = 224 GPa, C₄₄ = 124 GPa, C₆₆ = 124 GPa (Xiong et al. 2022), with orientation inherited from the grain it resides in. The γ″ phase was also modeled as an elastic crystalline material with the same constants as the δ phase. A characteristic set of orientations was observed for the γ″ phase within the γ microstructure (Zhang et al. 2019), this model randomizes instead the orientation with respect to the orientation of the metal grain as a simplification.

The material properties for the crystal plasticity model were developed by simulating the tensile response of a grain aggregate and comparing with the stress-strain response of the material. The mesh was refined around all the δ and γ″ phases, for a minimum element size of 0.0125 µm. The model includes 18 grains, 68 δ phase particles, 5000 γ′ and 5000 γ″ phase particles; the number of nodes was 1,262,155. The model was developed as parametric, allowing all the inputs to be changed systematically and to observe the variations in mechanical response.

The microstructural model is loaded in tension by applying a constant displacement on the top boundary with symmetry conditions at the bottom. The displacement boundary condition for the microstructural model was set as to represent the condition shown in Figure 4 with the model being subjected to the applied strain calculated from the global model. The condition of strain, εyy, is simplified as an effective strain of 3.5 %, which corresponds in Figure 2 to an average value of strain obtained within a distance of 100 μm from the notch tip. The model also included an initiated crack at one side of the grain aggregate corresponding to the notch root of the sample. The crack is modeled as a microstructurally sharp notch with root radius of 1.25 μm.

4 Formulation for hydrogen diffusion

The diffusion of hydrogen through the microstructure was modeled by using Fick’s law of diffusion (Smith 2009). The assumption is that hydrogen is generated at the mouth and surface of the initiated crack and then diffuses within the material microstructure. The model uses a stress-assisted diffusion formulation to account for the role of hydrostatic stress as a contributor to the chemical driving force given by the concentration gradients throughout the microstructure. The diffusion simulation is only weakly coupled to the mechanical analysis as it is performed with input from it, but not at the same time as the mechanical simulation and it does not influence the mechanical response.

The flux of diffusible hydrogen is given by the chemical and hydrostatic stress gradients within the material. The model implementation accounts for the effects of chemical and hydrostatic stress gradients directly. The additional driving force is given by the velocity of mobile dislocations, obtained from the mobile dislocation densities calculated from the mechanical loading of the sample (Dadfarnia et al. 2015). The formulation for hydrogen transport can be summarized as:

where CL and CT are the trapped and diffusible concentrations of hydrogen in the microstructure, JT and J are the corresponding fluxes, with α being the number of sites per trap, and θT the number of traps per solute volume, NTm is the dislocation trap density, VH is partial molar volume of hydrogen in solid solution, σkk is hydrostatic stress. The hydrogen diffusivity, D, and solubility are 2.9 × 10⁻¹¹ cm²/s and 20.85 molH2/m3MPa respectively based on literature data (Jebaraj et al. 2014). Hydrogen trapped by mobile dislocations is transported with a velocity Vρm calculated as:

where ε˙p is obtained by assuming a final plastic strain of 3.5 %, a total loading time for the sample of 1800 s, Burgers vector, b, is 0.2 nm, and ρm is obtained from the mechanical analysis. The flux of trapped hydrogen is determined by the number of available trap sites per volume of solute, αθT, and by the dislocation trap density, expressed as:

where ρm is mobile dislocation density calculated by the model, a is lattice parameter, and ρ=3 for FCC materials. The model does not calculate separate concentrations for diffusible, CL, and trapped, CT, hydrogen and uses the above formulation to calculate a single concentration of hydrogen at steady-state. The concentration is the summation of CL and CT.

The model accounts indirectly for the effects of immobile dislocations as irreversible traps for hydrogen by including a dependence of hydrogen diffusivity on immobile dislocation densities. The diffusivity is modified as a function of immobile dislocation density, changing linearly and becoming 1 2 of its initial value when the immobile dislocation density is 10⁵ 1/μm². The latter value does not have at this stage a supporting experimental basis, and it is only used as variable to estimate the effects of trapping by immobile dislocations. The results will be shown for the case of no trapping by immobile dislocations, considering the effects of the driving forces in Eqs. (7) and (8) alone, and for the case of medium trapping and strong trapping effects, arbitrarily defining a diffusivity reduced to 1 2, and to 1 10 respectively of its initial value when the immobile dislocation density is 10⁵ 1/μm².

The solution domain for the diffusion analysis is the grain aggregate shown in Figure 3. The hydrogen flux is constant at the beginning of the analysis for the “external surfaces” of the representative microstructure, i.e. the surfaces on the left side of the domain in Figure 3, and for the surfaces forming the crack faces and crack tip. These areas correspond to the notch root of the global model, shown in Figure 4; note that the solution is not obtained for the entire global model domain and only limited to the local representative microstructure domain in Figure 3. The rest of the domain is set to have a concentration of 0.0 ppm at the beginning of the simulation, and the diffusion analysis proceed in two steps. In the first step the flux is constant at the boundary until the concentration of hydrogen at the boundary reaches a level of 5.13 ppm, corresponding to the total hydrogen concentration, including trapped hydrogen, for the material as obtained from the literature when the material is polarized to −1 V_SCE (Lillard 1998). In the second step of the analysis, the concentration at the boundary remains at 5.13 ppm and the diffusion equations are continuously solved throughout the domain.

Figure 3:

Local model with representative microstructure for 718-API material.

Figure 4:

Representation of local model in reference to the global model and the resulting plastic strain at the notch.

The boundary condition change is justified, based on electrochemical hydrogen charging experiments at different cathodic potentials performed for 718 by Lillard, followed by measurements of hydrogen concentration by inert gas fusion analysis. A relation between cathodic current and hydrogen concentration showed in Lillard’s work is used in this work along with the polarization curve for 718 from (Thodla et al. 2020) to find the current density at −1.0 V. The flux is set to yield a concentration of 5.13 ppm at the external boundary in 1 h of exposure time. The analysis is continued with a constant concentration boundary condition at the external surfaces for a total of 10⁷ s or approximately 2700 h.

5 Microstructural model results

5.1 Polycrystalline stress-strain response

This section describes the results of the stress-strain simulation performed with the model of the microstructure without the initiated crack to validate the choice of material parameters in the crystal-plasticity model described previously. The polycrystal is subjected to a prescribed strain and the average stress in the direction of straining calculated by the model is used as a metric for comparison with the stress-strain behavior of the material. The analysis assumes plane stress constraints to represent a typical tensile test standard condition. A constant strain is applied as boundary condition at the top portion of the model with the bottom being symmetric. An average measure of stress is obtained from the model as a representation of the engineering stress value.

Figure 5 shows the resulting stress-strain response of the grain aggregate. The stereographic projections representing the orientations of the grains within the microstructure are shown on the right (Arcari et al. 2023); this is representative of a uniformly oriented grain structure as it would be obtained for example after forming into a solid bar. The Euler angles are not from direct measurements, and randomly generated for this set of grains. The stress-strain results are consistent with the elastic portion of the curve and for the initial phase of yielding of the material. The simulation results are shown to deviate for large strain values, with fluctuating stresses for stains higher than 1 %. This is due to the low number of grains simulated and the texture of the grain aggregate, showing the progressive activation of more slip directions as the material is strained beyond yielding (Rees 2006). The engineering stress-strain curve is an upper limit for the analysis for strains within 1.5 %, while being lower than the stress predicted by the crystal plasticity model. The further increase in stress as a function of strain, rather than a progressively lower strain hardening is a function of the directionality of the metal matrix.

Figure 5:

Stress-strain response of the metal grain aggregate subjected to tension (no notch present).

5.2 Mechanical response in the presence of a crack

The results from the mechanical loading of the microstructural model with an initiated crack are presented in this section. Boundary conditions representative of the strain state at a sharp notch of a fracture mechanics samples, as presented in Figure 2, are applied at the top surface (see orientation in Figures 1 and 4). A condition of symmetry is specified at the bottom boundary. A state of plane strain is imposed as solution condition within the model to simulate the constraint effects at the center of a fracture mechanics sample.

Figure 6a shows the total cumulative shear strain at the end of the deformation step, with 5 % being the maximum value shown by the scale. The total cumulative shear strain is defined as the sum of the absolute values of shear strains in all slip systems, as γ=∑(α)∫0t|γ˙(α)|dt. Shear strain is highly localized at the initiated crack, as expected, exceeding 20 % at the tip of the crack, but also localized at grain boundaries near the crack. The slip gradient across some grain boundaries is high, mostly driven by the relative grain misorientation, and specifically for grain boundaries with a 45° orientation with respect to the loading direction.

Figure 6:

Results of the numerical analysis for the application of a constant displacement resulting in 4 % strain of the metal grains aggregate representing the state of deformation at the notch tip; the images show: (a) total slip accumulation, (b) hydrostatic stress, (c) mobile dislocation density, (d) immobile dislocation density, and (e) calculated plastic strain.

Figure 6b shows the distribution of hydrostatic stress within the microstructure. A large area of high hydrostatic stress ahead of the crack tip is determined by the model. Note that the finite element solver used for this problem defines pressure stress as −13σii, hence the results in Figure 6b indicate positive hydrostatic pressure. The average value of hydrostatic stress in this area is above two times the yield strength of the material. The locations throughout the microstructure with increased hydrostatic stress are mostly in the vicinity of grain boundaries, forming a path in the center of the microstructure. The interface between matrix and precipitates also shows increased hydrostatic pressure. This results from the difference in the material properties between the matrix and precipitates, also increased by a different modeling strategy used in this work: the matrix is modeled as an elastic-plastic material, while the particles as elastic, leading to local increase in hydrostatic stress to accommodate the difference in lateral constraints at the interface during inelastic deformations.

Figure 6c shows the mobile dislocation density obtained by the model for the applied 3.5 % strain. The maximum value shown in the image is 10,000 1/μm², which corresponds to a density of 10¹⁶ 1/m². Mobile dislocation densities are highest at the root of the crack; however areas of increased dislocation activity are visible, specifically along grain boundaries towards the interior of the microstructure. The mobile dislocation density path along the grain boundaries results from the higher resistance to slip as determined by the different relative grain orientation and the grain boundary orientation with respect to the applied load. This causes a significant gradient in both mobile dislocation densities across the grain boundary.

Figure 6d shows the immobile dislocation density obtained by the model for the applied 3.5 % strain. The maximum value shown in the image is 20,000 1/μm², which corresponds to a density of 2 × 10¹⁶ 1/m². The immobile dislocation densities calculated by the model are on average higher than the mobile densities. The highest values of immobile dislocation densities are shown around the crack and along grain boundaries leading towards the interior of the microstructure. A significantly large area ahead of the crack tip shows a higher density of immobile dislocations, followed by an area of low density immediately adjacent to it. The immobile dislocation density gradients are less severe than the gradients shown for mobile dislocation densities. This is particularly evident by looking at the immobile dislocation density distribution at or near grain boundaries in comparison with the distribution at the grain interior.

Figure 6e shows the plastic strain calculated by the model for the domain. The plastic strain value at the crack tip is higher than the upper bound value indicated in the figure, with strain exceeding 10 % at the tip of the crack. The limit values of the legend only represent useful bounds to appreciate the location of higher calculated plastic deformations. The areas where highest shear strain was identified earlier clearly show higher plastic deformations in the model.

6 Microstructural model diffusion results

The results of hydrogen diffusion from the surface of the material to the fracture process zone of an initiated microstructurally small crack are presented in this section. The results are obtained first for a model including the effects of hydrostatic stress, mobile dislocation velocity, and trapping by immobile dislocations. For this case the trapping by immobile dislocations is simulated by decreasing the diffusivity of hydrogen linearly with increasing immobile dislocation density, becoming 1⁄2 its initial value when the immobile dislocation density is 10⁵ 1/μm².

Figure 7 shows the resulting hydrogen concentration after a simulation time of 10⁷ s. The maximum simulation time is identified as the time that minimizes the maximum change in concentration as a function of time calculated in the model. The final state is considered achieved when the maximum concentration change in the model is 0.01 ppm per 86,400 s (one day). Figure 7a shows the concentration of hydrogen is highest at a small distance ahead of the crack tip. The concentration is higher than the surface concentration by a factor of 1.2×. The normalized profile of hydrogen concentration as a function of distance from the crack tip is shown in Figure 7b. The data are extracted from the path shown in Figure 7a, and the normalized profiles of hydrostatic stress, mobile and immobile dislocation densities, and accumulated slip are also shown for reference. The hydrogen concentration is normalized by 5.13 ppm, the boundary condition at the sample surface at steady state. The area of highest hydrogen concentration generally corresponds to the area with the highest hydrostatic stress calculated by the model. The location of the highest concentration does not however correspond the location of maximum stress, as it is affected by the presence of a high density of mobile and immobile dislocations. The high immobile dislocation density right at the tip of the crack causes a reduction in the flux and a lower accumulation of hydrogen, resulting in a reduction of hydrogen concentration with respect to the condition at the boundary, followed by a significant increase.

Figure 7:

Diffusion results for microstructural-representative model of UNS 07718 material: (a) distribution of hydrogen through the material microstructure, and (b) normalized concentration as a function of distance from the crack tip, along the path shown in (a).

Other areas of higher concentration of hydrogen can be identified throughout the microstructure: precipitates, such as the δ phase on some grain boundaries show higher concentration than average, γ′ phases on the other end show lower concentration than the surrounding matrix, while γ″ does not seem to play a strong role in either increasing or decreasing the nearby hydrogen concentration.

To assess the contribution of irreversible trapping of hydrogen by immobile dislocations, two additional models are developed, and the results compared to Figure 7. In the first model the trapping by immobile dislocations is not included in the analysis. This is simulated by removing the dependency of diffusivity on immobile dislocation densities. In the second model a more severe change in hydrogen diffusivity as a function of immobile dislocation densities is implemented. This is achieved by modeling a sharper decrease in diffusivity with increasing immobile dislocation densities: the diffusivity decreases to 1/10 of its original value when the immobile dislocation density reaches 10⁵ 1/μm².

The map of hydrogen concentration for all three cases is shown in Figure 8. The scale is changed with respect to the results in Figure 7 to better appreciate the relative differences.

Figure 8:

Distribution of hydrogen within the material microstructure based on the diffusion modeling analysis for three cases of influence of trapping by immobile dislocations.

Figure 9 shows the distribution of hydrogen as a function of distance from the crack tip and compares the results for the case of diffusivity changed by a factor of ½ (baseline case), and by a factor of 1/10 and when trapping by immobile dislocations is suppressed. The normalized profiles of hydrostatic stress, σ_H, immobile dislocation density, ρ_Im, mobile dislocation density, ρm, and accumulated slip, ∑γ, are shown for reference. Trapping by immobile dislocations changes both the peak hydrogen concentration and the distribution of hydrogen near the crack tip. The concentration is highest for the analysis case where no trapping by immobile dislocations. The peak of hydrogen concentration moves closer to the tip of the crack when trapping is not present. An increased rate of trapping on the other hand causes the peak to shift down and towards the interior of the microstructure. To better identify the role of mobile dislocation velocity, its contribution is removed from the analysis. The results are shown in Figure 10 where the normalized concentration of hydrogen as a function of distance from the crack tip is shown. The effect of transport by mobile dislocation is small, with the peak changing by less than 1 %.

Figure 9:

Profile of hydrogen concentration along the paths shown in Figure 8 as a function of distance from the tip of the simulated crack. The normalized profiles of hydrostatic stress, σ_H, immobile dislocation density, ρ_Im, mobile dislocation density, ρ_m, and accumulated slip, ∑γ, are shown for reference.

Figure 10:

Comparative results of hydrogen concentration for the analysis with mobile dislocation velocity set equal to 0, as compared to the baseline case; normalized concentration as a function of distance from the crack tip, along the path shown in Figure 7a.

7 Discussion

In this work we identified the main characteristics of a PH nickel alloy microstructure that play a significant role in the material susceptibility to hydrogen embrittlement and developed a model to represent the response of the material to the applied loads. The microstructure of UNS 07718 is known to be driving the response to hydrogen producing environments; however the role of different constituents in the microstructure is not understood from a perspective of combined mechanical-chemical driving force.

From a mechanical standpoint the local conditions at a critical location, such as a notch, are needed. Despite the condition at the notch tip being nominally of yielding, it is important to advance the qualification of this state with the help of numerical tools currently available. Crystal plasticity formulations are suitable for this purpose and within the limits of the modeling assumptions can help qualify the differences between dissimilar microstructures. We hypothesized that the mechanical response of the material at the level of grains, precipitates, and dislocation densities would inform our understanding of hydrogen diffusion from the tip of an initiated crack to the fracture process zone, determining its advancing rate.

The results show that within an area of the microstructure subjected to significant plastic deformations mechanical driving forces vary significantly, both in terms of localized deformations (slip), and in terms of accumulated dislocations. Our formulation implements a hardening law as a function of generated dislocations within the microstructure. The nucleation, emission, and annihilation of dislocations are empirically accounted for by the evolution equations (Equations (4) and (5)) and associated coefficients (Arcari et al. 2023; Zikry and Kao 1996). Within these assumptions, grain orientation and grain boundary characteristics are observed to influence the spatial gradients of hydrostatic stress, mobile and immobile dislocations.

All of these are considered important in determining the driving force for hydrogen diffusion ahead of an initiated crack, as described by Equations (7) through (8), and discussed in the literature (Djukic et al. 2019). A high localization of slip, mobile and immobile dislocations at the tip of the initiated crack is expected, however the balance of the three main driving forces ahead of the crack tip considered in this work is important. The relative comparisons of hydrogen diffusion within the microstructure show that the hydrostatic pressure gradient ahead of the advancing crack is the strongest driving force for hydrogen accumulation. Hydrogen is drawn towards the fracture process area; however it is not the area with the highest calculated hydrostatic stress, but an area adjacent to it to be the area with the highest concentration. This is due to the combination of high hydrostatic stress and the presence of immobile dislocation density, acting as traps for hydrogen in the material. As shown by comparing the results of Figures 8 and 9, the effects of immobile dislocation density in modifying the location of the highest concentration of hydrogen is significant. Interestingly, the area with the highest concentration of hydrogen is an area that includes part of the microstructure with the lowest values of accumulated slip, within the highest hydrostatic stress region, corresponding to the regions of the microstructure with the lowest density of traps. Diffusion by mobile dislocation velocity is not, in this formulation, as strong of a driving force as hydrostatic stress; results obtained by excluding the effects of mobile dislocation velocity do not differ significantly. The work from Dadfarnia et al. (2015) showed the contribution of dislocation velocity and transport on hydrogen diffusion ahead of the tip of the crack was significant, which is contrary to what we found in our analysis. Qualitatively both models agree that the peak of hydrogen concentration is ahead of the crack tip with the distribution further into the material being a function of hydrogen diffusivity through the lattice. Our work indicates that the local microstructure can influence the contribution of hydrogen transport by dislocations within the material, reducing its effectiveness in some cases. Our work however did not account for a full coupling at this stage, and the mobility of dislocations is not being affected by the increasing concentration of hydrogen, which could increase the transport of hydrogen and balance the trapping by immobile dislocations.

Areas at the interfaces of matrix and precipitates are locations where hydrogen accumulation is influenced by the driving forces in this model. As expressed above, while it is good to note that these areas correspond to the precipitates within the microstructure, this only implies a tendency of hydrogen to redistribute at sites of gradients and discontinuities in deformation response of the materials. In this model matrix and precipitates have been modeled with significantly different behavior, elastic-plastic and only elastic, respectively, to exacerbate the inherent differences in mechanical response. The results certainly stress the importance of the combination of the driving forces modeled in this work for hydrogen diffusion and the consideration of including their effects into a modeling strategy. The model developed in this work also simplifies the deformation behavior of the microstructure of the material, in particular simplifies the interaction of the matrix and the precipitates without considering the shearing of precipitates (Zenk et al. 2021). Because of the choice in constitutive equations modeled the precipitates are considered fully coherent with the matrix during elastic deformations, while they act like obstacles to deformations when the matrix plastically deforms. The intent was to idealize the behavior of the interaction to reduce the computational burden for a numerical simulation of a larger portion of the microstructure ahead of the crack tip. At the level of volume fractions modeled in this work, the assumption seems reasonable.

The use of a computational model is advantageous as separate effects can be quantified independently. While a validation data set would be difficult to obtain for this type of results, making a quantitative evaluation of material behavior based on the results of modeling difficult, relative comparisons can be beneficial to inform our understanding of the relation among different driving forces, such as hydrostatic stresses, mobile and immobile dislocations, and hydrogen diffusion within the material microstructure. Specifically, the local material texture, grain size and orientation, along with grain boundary characteristics could be systematically varied to understand what presents the weakest link in the progression of an initiated crack.

8 Experimental validation efforts

Experimental characterization of the material microstructure is currently underway to inform the microstructural model presented in this work. The goal is to characterize the precipitates within the metal grains, their sizes and characteristics to develop a unique model for variants of 718 material and representing their composition, thermal-mechanical processing and resulting microstructure. Measurements were made using Atom Probe Tomography (APT) to identify at the atomic scale the composition, shape, and conformation of the precipitates for UNS N07718-API alloy produced following the API specification and a thermal-mechanical processing to yield a minimum yield strength of 150 ksi. APT samples were prepared by standard liftout and milling techniques in a ThermoFisher Nova 600 Dual Beam focused ion beam-scanning electron microscope (FIB-SEM), and APT was conducted using a CAMECA 4000X Si local electron atom-probe (LEAP). A sequence of sample preparation steps to yield a sample that can be analyzed with APT are shown in Figure 11a and b. A protective Pt strap was deposited on the 718 sample before lifting out material with a micromanipulator needle, Figure 11a. The liftout material is a triangular prism that’s been shaped using the FIB and then deposited onto a grid of posts on a microtip array. The array is made of pre-shaped posts for placing the materials of interest, which is then sharpened using a series of annular mills with an ion beam into a very sharp needle shape, Figure 11b, so that ions can be field evaporated from it during the APT measurement. The results in Figure 11c show the microstructure of the material visualized by iso-surface contours with a threshold is an APT reconstruction of one of these tips. The interface surfaces indicated correspond to 57 % Ni concentration, and highlight the interfaces between the matrix, colored in pink, and Ni-, Nb-rich precipitates, shown in green. Figure 11d shows a proximity histogram, which gives the atomic concentration measured in the APT sample as a function of distance from the precipitate-matrix surfaces shown in Figure 11c. The data show the quantitative partitioning of elements to the matrix and the precipitates.

Figure 11:

Ongoing experimental validation effort by atom probe tomography: (a) and (b) showing sample preparation, (c) showing the iso-surface reconstruction of precipitates with limiting surface being 57 % Ni in content, and (d) showing a proximity histogram collected during measurement and highlighting the precipitates boundaries and compositions.

The experimental technique is planned to be adapted for post-charging hydrogen measurements to identify at the atomic level the distribution of hydrogen within the material microstructure. This would provide a unique validation data set for the modeling approach in this work. The goal is to electrochemically charge a needle-shaped sample, such as the one shown in Figure 11b, in a polarization test to allow hydrogen to diffuse through the material microstructure. The APT measurement post-charging will enable mapping the distribution of hydrogen within the material, providing a reference window and validation of the microstructural material modeling. Extending the model shown in this work to a 3D representative volume element to represent the material microstructure such as the one in Figure 11c, would also allow a more effective validation methodology. This technique combined with microstructural characterization by electron backscatter diffraction (EBSD) and transmission electron microscopy (TEM) will enable characterization across multiple length scales and of larger volumes, and distributions of microstructural descriptors could be used as inputs for developing a microstructural model of the local conditions at the crack tip (Zhao et al. 2022).

9 Conclusions

In this study we modeled the stress-strain behavior at the notch of a standard fracture mechanics sample by modeling. The material was PH nickel-based superalloy 718, and the relevant material properties were used to identify the local deformations at a sharp notch for a CT sample. A meso-scale microstructural model was developed to study the role of the microstructure in the deformation behavior of a small representative volume of material ahead of the notch. The representative model of the material microstructure included the metal grains, the material texture, and the major precipitates along the grain boundaries and within the grains. The model was subjected to the plastic strain conditions identified at the notch to simulate the metal grain aggregate response from a mechanical standpoint when an initiated crack is introduced. The modeling results were used as input for the diffusion of hydrogen from the surface of the notch, assuming the material had been introduced to a hydrogen producing environment. The results show that hydrostatic stress is a strong driving force for hydrogen, increasing its concentration near the crack tip. Trapping by immobile dislocations, simulated as a decrease in diffusivity corresponding to locations within the microstructure with the highest calculated density of dislocations, can strongly influence the location of the highest hydrogen concentration, moving it towards the interior of the material.