Incorporating hydrogen in mesoscale models

Hydrogen embrittlement is a major and longstanding challenge in materials science. Mesoscale mechanical models are currently being developed which are providing new insight into hydrogen-dislocation interactions and the subsequent detrimental effect on ductility and ultimate tensile strength. Here, a review of some recent approaches to incorporating hydrogen within crystal plasticity, discrete dislocation plasticity and cohesive zone models are reviewed. Modelling at this scale bridges the gap between atomistic models and continuum en- gineering models which will be essential to fully understand hydrogen embrittlement.


Introduction
In future energy and propulsion systems, engineering materials will be exposed to ever more extreme environments. To accurately simulate the mechanical properties of engineering materials in service will require new approaches for incorporating the influence of the environment on the mechanical properties. Hydrogen is detrimental to structural integrity as it reduces the toughness and ultimate tensile strength of engineering alloys as, shown in Fig. 1. For over a century this has been well documented [1] but is still not well understood. Crucially, hydrogen can promote failure by reducing ductility, this complex and material dependent process is known as hydrogen embrittlement (HE). Engineering applications where understanding HE is critical include high strength steels; used in piping in the oil and gas industry [2]. In hydride forming materials such as zirconium alloys, hydrogen can diffuse to a notch and form brittle hydrides which then fracture, propagating the crack, providing a failure mechanism known as delayed hydride cracking. This is a challenge in the nuclear industry where zirconium alloys, used as fuel cladding due to their low neutron cross section, are exposed to coolant water in certain nuclear reactor designs (CANDU and RBMK) [3]. This review highlights some developments which have taken place over the last two years to simulate HE. Reviews of (HE) [4,5], and mesoscale modelling techniques, discrete dislocation plasticity (DDP) [6,7] and crystal plasticity (CP) [8] already exist. The focus is on coupled mechanical/diffusion modelling of hydrogen at the continuum scale and on incorporating hydrogen within discrete dislocation plasticity. Finally some interesting parallels between hydrogen and radiation embrittlement are briefly discussed.
There are several proposed HE mechanisms and it is not yet clear which dominates in which regime. This is where modelling and simulation are crucial. Hydrogen enhanced decohesion (HEDE) (often called hydrogen induced decohesion HID) and hydrogen enhanced localised plasticity (HELP) are arguably the two main mechanisms. HEDE assumes HE is due to hydrogen reducing the cohesive strength and is supported by atomistic simulations [9,10]. HELP assumes that hydrogen increases dislocation mobility [11] and nucleation [12] as well as reducing the interaction stress between dislocations. It is supported by in situ microscopy [13,14], deformation experiments [15], macroscopic flow stress measurements and fractography [16]. HEDE and HELP can also operate together [17,18], this is not surprising given the complex microstructure typical of high strength engineering alloys. Hydrogen weakens the strength of interfaces such as grain boundaries which act as trap sites for hydrogen [10] but fracture is accompanied by enhanced dislocation activity when hydrogen is present [11]. High dislocation densities generate stress and act as trap sites producing concentration and stress hot spots close to materials interfaces, which can also leads to decohesion. Exciting recent developments allow hydrogen to be incorporated within large scale discrete dislocation plasticity simulations for the first time [19,20] closing the gap between atomistic models, TEM observations and micromechanical tests. Discrete dislocation models can be used to calibrate crystal plasticity and continuum plasticity models. These plasticity models can then be used in combination with cohesive elements or phase field methods to simulate the interaction between hydrogen, plasticity and fracture.
Over the last decade, the availability of powerful graphics processing units (GPU) allow larger samples to be simulated [22] more easily without the need for a supercomputer. Whereas new micromechanical techniques such as focused ion beam milling and nano-indentation, allow smaller samples to be tested. This convergence means it is now possible to perform virtual experiments and make direct comparison between discrete dislocation plasticity simulations and micromechanical tests. Popular geometries include micro pillar compression [23], microcantilever bending [24,20], and nano-indentation [7]. These new techniques are now being adapted and used to study the one hundred year old problem of hydrogen embrittlement. Microcantilever bend tests in a hydrogen environment [25,26] provide definitive observation of a reduced flow stress and enhanced slip localisation and crack propagation. Simulating these tests using discrete dislocation plasticity which is informed by atomistic simulations will provide essential new insight into the mechanisms controlling hydrogen embrittlement. Discrete dislocation dynamics simulations can then be used to calibrate crystal plasticity models to allow accurate simulation of complex microstructures.

Coupled mechanical diffusion modelling
Hydrogen is abundant and, being the smallest atom, readily diffuses through metals by occupying the normal interstitial lattice sites (often refereed to as NILS in the hydrogen literature). Conservation of mass can be expressed as where the first term is the rate of change of the hydrogen concentration, c, (atoms per volume), the second term is the divergence of the flux, and r is a sink term (if > r 0) or source term (if < r 0). The hydrogen flux, j i , along a direction, x i , is down the gradient in concentration, c, and up the gradient in hydrostatic stress, = /3 h kk . Hydrogen therefore flows to regions of low concentration or high tensile stress, so as to minimise the chemical potential. D is the diffusion coefficient and is the volume change due to a hydrogen atom in the lattice where m is the partial molar volume, and N A is the Avogadro constant and k B is the Boltzmann constant. If the temperature is nonuniform then the temperature gradient would also contribute to the hydrogen flux, however this term is usually ignored. For many engineering problems this is a valid assumption and has a major practical advantage. Commercial finite element software (such as Abaqus) can solve coupled thermal-mechanical problems, and because the governing diffusion equation for heat and mass are equivalent (see Table 1), finite element software optimised to solve coupled mechanical/thermal problems can also be used for solving the stress driven mass diffusion equation [27]. There are some subtleties in the equivalence which have been overlooked, which become apparent when solving for the lattice occupancy, rather than the lattice concentration. Temperature, T in a thermal problem is equivalent to lattice occupancy, , in a mass diffusion problem. Specific heat capacity, c m , is equivalent to solubility, s, and mass density, m is equivalent to hydrogen lattice site density (see Table 1 for further details). Methods for solving coupled thermal-mechanical problems can then be readily deployed to solve stress driven diffusion problems, for example by writing a user thermal material (UMATHT) in Abaqus. In order to accurately evaluate the pressure gradient which appears in the hydrogen flux (the second term in Eq. (2)), quadratic elements should be used (although this is not always the case in the literature). In fact, Abaqus imposes a constant pressure at the integration points in coupled linear thermal-mechanical elements (CPE4T) meaning the pressure gradient would be zero so these elements should be avoided when simulating mass diffusion.
Metals and alloys contain microstructural defects which are either static, such as grain boundaries and precipitates, or dynamic, such as dislocations. These defects act as hydrogen traps, as it is energetically favourable for a hydrogen atom to reside at a defect site rather than a lattice, resulting in a higher concentration at traps. Oriani showed that, provided the lattice and trapped hydrogen are in chemical equilibrium, the trapped hydrogen concentration is fully determined by the lattice concentration and the number of trap sites [28]. Therefore, even when a multitude of different trap types are present there is still only one unknown field, the lattice occupancy, which then fully determines the occupancy in the traps and the total hydrogen concentration. When the  Table 1 The equivalence between thermal and mass diffusion. = C c/ is the hydrogen concentration, c, scaled by the hydrogen lattice site density, . The lattice occupancy is and the density and occupancy of a trap of type i are and i i respectively.

Thermal Diffusion Mass diffusion
Energy conservation Mass conservation Thermal energy per unit mass (J kg −1 ) Concentration, H atoms per H lattice site Rate of flow of heat into sinks Rate of increase of trapped H atoms per unit volume Lattice occupancy gradient (m −1 ) Heat flux derivatives Hydrogen flux derivatives trap density is time dependent then an extra term, r, appears in the diffusion equation which is a sink or source depending on if the number of traps is increasing or decreasing. Note that an incorrect implementation can inadvertently violate conservation of mass and it is essential that suitable checks are performed. The coupling between the mechanical problem and mass diffusion is not simply one way. While the pressure gradient drives hydrogen diffusion, the concentration of hydrogen also influences the mechanical stress. Trapped hydrogen can significantly alter the properties of the traps, atomistic simulations quantify how hydrogen reduces the grain boundary cohesive energy [9] and can increase the mobility of screw dislocations in iron [11]. Constitutive laws can then be implemented within larger scale models such as discrete dislocation plasticity [20], crystal plasticity [29] or continuum plasticity [27]. Cohesive zone modelling allows a hydrogen dependent cohesive law to be implemented within continuum plasticity, we initially did this with a steady state hydrogen distribution [30] using an image based mesh generation scheme [31] to simulate a real microstructure. Microcracks formed at several different matrix/carbide interfaces in the region with highest hydrogen concentration. These microcracks were then connected by localised plastic flow in the bulk. More recently researchers have implemented a hydrogen dependent cohesive law [32] in conjunction with the UMATHT method outlined above [27] to solve the stress driven diffusion equation to simulate environmentally assisted fracture [33] and fatigue [34]. These papers are notable for making direct comparison with experimental data. In the future it is likely that models will move away from solving for hydrogen occupancy completely and instead solve for the chemical potential [35] eliminating the need to compute the hydrostatic stress gradient when evaluating the flux.
Many papers in the literature share stress and plastic strain data with the diffusion subroutine by writing to disk which is very slow. Data can be easily shared between different Abaqus subroutines by using a common block in fortran. In order to use multiple cores the Abaqus utility function MutexLock must be used when writing to the common block to prevent multiple threads from attempting to write to it simultaneously. These simple techniques allow much larger problems to be simulated and enable strain gradient crystal plasticity to be implemented within a user material (UMAT) [36]; which is generally less involved than a user element (UEL).

Hydrogen and cohesive zone modelling
Currently the cohesive elements [37] used in the hydrogen literature are formulated so that the traction separation law has a reduced damage initiation stress with increasing hydrogen coverage, a typical example is shown in Fig. 2 but hydrogen flux across the cohesive element is not implemented. This means the cohesive zone acts as an insulating boundary which limits their use to symmetric problems. However this limitation can be overcome by including the additional degrees of freedom for hydrogen occupancy at the nodes of the cohesive element. The cohesive law in 2D then has the form where T i and ũ i are the traction and opening displacement across the cohesive element in direction i J , n is the hydrogen flux and is the occupancy gradient across the interface. So that the mechanical damage, d, degrades the conductivity across the cohesive zone in the same way as it reduces the stiffness k. This allows hydrogen flux to be implemented with only a minor modification to a conventional cohesive element [39]. To simulate the reduction in cohesive strength due to hydrogen one or more cohesive law parameters are reduced. If the cohesive law parameters vary significantly within an increment then this should be accounted for in the tangent stiffness matrix to ensure convergence. Formally the required cohesive contribution (Abaqus AMATRX) to the global Jacobian is then and the residual vector for the cohesive element is = B c u . This extension of a traditional cohesive element implementation should enable complex microstructures representative of high strength engineering alloys to be simulated. Cohesive elments are simple to implement but have the disadvantage that the possible crack paths must be predefined. Alternatively, phase field methods for hydrogen assisted cracking [40] look promising and have the potential to simulate complex crack paths which will be essential to simulate hydrogen embrittlement in industrial components.
In the future hydrogen diffusion and damage need to be modelled in real microstructures in order to understand HE in materials of interest to industry. This can be achieved by combining the approaches outlined above with image based mesh generation methods [41] or by using representative virtual microstructures [42]. We performed an initial study with an image based mesh and a hydrogen dependent cohesive law and plasticity law but used a static hydrogen distribution [30]. This is illustrated in Fig. 3. The model showed how failure occured in interfaces with a high hydrogen concentration linked by regions of high plastic strain.

Hydrogen and discrete dislocation plasticity
Discrete dislocation plasticity (DDP) is crucial for simulating hydrogen as it bridges the gap between molecular dynamics (MD) and crystal plasticity (CP) simulations. DDP is based on elasticity theory so is efficient and can be used to study the emergent behaviour or large numbers of dislocations. The elastic interactions between dislocation segments are computed analytically [43] and supplemented with constitutive rules such as the dislocation mobility law and core force obtained from MD. The dislocation network is usually discretised into straight finite segments bounded by two nodes [43]. The driving force on a dislocation node k is then, where contributions arise from the elastic interactions between dislocation segments, F , the corrective force due to the applied and image stress computed using FEM, F , and the core force, F c , derived from atomistic simulations. To account for the hydrogen elastic shielding effect an additional force can be added, F , which is proportional (but with opposite sign) to the force generated by an edge dislocation segment [19]. Importantly this additional hydrogen force is then determined by the dislocation structure, assuming fast diffusion or chemical equilibrium, or by the initial dislocation structure assuming slow diffusion. That is to say the hydrogen concentration is fully determined by the dislocation stress field, rather than by solving the diffusion equation. This driving force is balanced by a drag force which is proportional to the velocity. Implementing a non-linear mobility law in 3D, requires an iterative scheme to solve for the nodal velocity so is rarely used. The velocity of every node can then be obtained as where the sum is over the l nodes connected to node k, with segment lengths L kl . The drag matrix B kl is dependent on the temperature, material and orientation of the segment. Atomistic simulations show the drag coefficient is also hydrogen dependent as evidenced by the increased mobility of screw dislocations at low concentrations in bcc metals [11]. A simulated load displacement curve for an end loaded microcantilever is shown in Fig. 4. At the low hydrogen concentrations typically found in bcc metals we found the hydrogen elastic stress was negligible [20]. We also implemented a hydrogen dependent mobility law for screw segments based on atomistic simulations [11]. This hydrogen enhanced screw mobility law was found to generate significant softening and a higher dislocation density even at low concentrations, as shown in Figs. 4 and 5. This agrees with the observation of hydrogen softening in microcantilever bend tests [44,45]. In bcc metals edge segments move significantly faster than screw segments, the fast moving edge segments generate long screw dislocations which can cross slip. Hydrogen reduces the anisotropy in the mobility which reduces the proportion of screw segments and in turn the amount of cross slip.

Informing crystal plasticity
Although dislocation dynamics is computationally less expensive than molecular dynamics it is still too costly to simulate systems much larger than around 10 μm even at small strains. Plasticity arises due to the nucleation and motion or large numbers of dislocations. A typical dislocation density in a mechanical test of 10 m 15 2 means that the total dislocation line length in just 1 mm 3 of metal is 1000 km if this were discretised with a 100 nm segment size, a total of 10 13 segments would be needed. Even storing the position of so many segments would be a challenge, requiring around 100 TB let alone computing, even approximately, the 10 26 elastic interactions between every pair of segments. Therefore crystal plasticity (CP) and continuum plasticity models are required at the component scale. The crystallographic slip rate is dis sys (8) where the sum over the rate of area slipped by every dislocation segment with line length l i and velocity v i in a volume of material V is approximated in crystal plasticity by a sum over the slip systems. Each slip system has a density of mobile dislocations, i , moving with an average velocity v i which is a function of the resolved shear stress relative to the critical value required for slip c and b i is the magnitude of the Burgers vector. CP models usually assume either a power or hyperbolic law for the slip rate, implicitly assuming a highly non-linear dislocation velocity v( ). However it is important to distinguish between the velocity which is linear in stress in many fcc and hcp materials with the non-linearity in the slip rate arising due to the dislocation density evolution. Rate equations for the dislocation density are of varying complexity. Usually the total density is divided into mobile and sessile populations with cross terms so that density can switch between mobile and immobile. Terms in the evolution equations can account for dislocation generation due to sources as well as sink terms due to dislocation annihilation and cross slip onto other slip systems [46]. An advantage of these physically based modes is that each term has an unknown coefficient which can be calibrated with discrete dislocation dynamics. The dislocation densities are usually updated at the end of the increment with a simple Euler forward integration scheme although implicit time integration methods where stress and density are solved concurrently are starting to be used [47]. This has the advantage or ensuring accurate convergence although at an increased computational cost. Strain gradients generate geometrically necessary dislocations (GNDs). Typically the GND density is obtained by computing the curl of the elastic rotation tensor e as this can be measured experimentally, and is related to the GND density through some care must be taken as three different definitions of the curl of a second order tensor are used in the literature, and the curl of the elastic or plastic parts can be used and taken in either the reference or deformed configuration; we published a comparison table in [36]. The modelling of GNDs is important because the elevation in stress associated with large dislocation densities in the vicinity of a crack tip could be fundamental in understanding hydrogen embrittlement. The specific stress level attained is not only relevant for decohesion but also because it governs the hydrogen distribution. Thus, dislocation-based models could be pivotal in providing accurate stress and concentration levels that could enable a sound understanding of crack tip conditions. Strain gradient plasticity (SGP) theories have been used in this regard, showing important differences in hydrogen concentration predictions, as well as improved predictive capabilities. SGP predictions of crack tip mechanics resemble those of discrete dislocation calculations. DDP simulations are an effective computational laboratory for guiding interpretation of SGP models [48,49]. We performed a subsurface comparison under a nano-indent between the simulated lattice rotations and GND density predicted by CP and measured using Laue diffraction as shown in Fig. 6. The agreement between the measured and simulated lattice rotations were good and provides some confidence that the GND densities obtained by solving (9) are reasonable. Note that there are only nine known values on the LHS and N sys unknown dislocation densities on the right. Therefore minimisation of either the elastic energy (L1), which is physically correct, or the sum of the squares of the densities (L2), which is easy to do, are typically used. We found the total GND density obtained using a simple L2 minimisation scheme agreed well with the more involved iterative L1 minimisation [36]. L2 minimisation is far easier to implement and can therefore be safely used if the hardening law is based on the total GND density. However if the dislocation density on each slip system must be distinguished then more sophisticated approaches are required.
A crystal plasticity model for hexagonal close packed (hcp) metals based on discrete dislocation simulations was recently developed by Messner et al. [50]. They extended the Kocks-Mecking model by using machine learning to rigorously calibrate CP with discrete dislocation dynamics. The dislocation density evolution equations were sys sys (10) where the coefficients K ijk accounts for the probability of a dislocation reaction i j k and accounts for all energetically favourable binary junctions. The hardening law was of the form where H ij is the hardening coefficient, the increase in the slip resistance on slip system i due to dislocation density on slip system j. By explicitly separating the density evolution and the hardening the model can be calibrated with discrete dislocation dynamics. Messner et al. [50] used machine learning to obtain the optimal H and K by calibration with large scale discrete dislocation dynamics simulations of beryllium. By simulating hydrogen accurately at the discrete dislocation level, it will therefore be possible to calibrate an accurate crystal plasticity model using this type of approach. To accurately include hydrogen in DDP will require calibration of the dislocation mobility law and hydrogen dependent core force with atomistic simulations. The increase in dislocation density, mobility and weakening of junctions due to hydrogen could then naturally be captured within CP based on the underlying physics.

From atoms to components
There are many similarities between hydrogen and irradiation, they both increase the initial yield stress and then lead to strain localisation and softening, they both cause swelling, dislocation formation, embrittlement and premature failure of materials, and both are inherently multiscale in nature. Simulating these processes is challenging as it requires accurately passing information from atomistic simulations up to mesoscale and macroscale models. However this also means that methods developed for one problem might also be applied to the other. Irradiation by high energy neutrons generates point defects. A detailed review of radiation damage models can be found in [51] and is not discussed here. Well established methods exist to compute elastic dipole tensors of point defects [52,53] including a hydrogen atom in zirconium [54]. Crucially the dipole tensor can be used to acurately characterise the stress field produced by a defect in a crystal lattice.
We recently showed how elastic dipole tensors can be defined and evaluated for defect clusters produced by the collapse of entire highenergy collision cascades [55]. It is possible to use the concept of the elastic dipole tensor to calculate stresses because the spatial scale of a defect is always small compared to the scale of a component. The stresses produced by radiation defects have the same microscopic origin as for interstitial hydrogen atoms and stem from the fact that radiation defects have significant elastic relaxation volumes [52], which produces significant local deformation of the lattice. For example the elastic relaxation volume of a self-interstitial atom in tungsten, predicted by density functional theory, is = 1.67 rel 0 , where 0 is the volume of an atom. We defined the density of relaxation volumes of defects, a dimensionless quantity equal to the total relaxation volume of all defect configurations in a unit volume which can be evaluated using atomistic simulations. We showed that the gradient in the density of relaxation volumes generates a body force on the material which is where B is the Bulk modulus. The gradient can be computed either analytically or numerically and implemented within the finite element method (FEM). Mechanical equilibrium means that this can be solved with FEM, where the nodal forces are where N a are the shape function of FE node a and T the applied tractions on the surface S T . We implemented this in Abaqus using a user (DLOAD) subroutine allowing a direct link between atomistic and continuum component simulations. An example simulation of a reactor pressure vessel and the swelling resulting from this body force induced by point defects is shown in Fig. 7. A cap at the top of the sphere was insulated, = r ( ) 0 rel within°20 of the upper pole, to demonstrate that this method can be used with arbitrary boundary conditions. This is powerful as it provides a direct link between atomistic modelling and finite element analysis of engineering components. We studied a static problem with a fixed distribution of point defects. The method could possibly be extended to study the evolution of hydrogen perhaps by first computing the stress and then solving the stress driven diffusion equation to update r ( ) rel at each time increment.

Conclusions
Now is an exciting time as the pieces are in place to develop physics based material models of hydrogen in metals. The dislocation mobility law and core force can be extracted from atomistic simulations and implemented within DDP. Large scale DDP simulations can then be used to calibrate crystal plasticity using machine learning. Initial DDP simulations indicate that hydrogen will reduces the values of H (the coupling between slip resistance and dislocation density) and increase the values of K (the coupling between slip system interactions and dislocation density). The first step would be to assume a linear dependence on occupancy H c   on the same materials with and without hydrogen. Although in its infancy, DDP with hydrogen has already demonstrated that the hydrogen elastic stress field is negligible at the low concentrations found in bcc metals; although it may be significant in fcc. The hydrogen elastic stress field could also be computed using the methods developed to simulate irradiation as discussed briefly in this review. In either case it is essential that the influence of hydrogen on the dislocation core force, F c , in Eq. (6) be fully investigated as this will likely have a significant effect which can not be captured by linear elasticity theory. Another exciting development will be on DDP and CP models which incorporate cohesive elements with a hydrogen dependent traction-separation law. It is important to incorporate hydrogen flux across the cohesive zone in CP and dislocation flux across the cohesive zone in DDP to allow non-symmetric problems to be simulated. Going forward this capability will allow the interaction between plasticity, hydrogen and fracture to be simulated with unprecedented accuracy. The final challenge is to simulate the complex multiphase, polycrystalline microstrucure typical of high strength engineering alloys. It is only at this stage that models will be able to guide the development of new materials as experiments indicate that multiple HE mechanisms operate.