The effect of hydrogen content and yield strength on the distribution of hydrogen in steel a diffusion coupled micromechanical FEM study

In this study, we investigate the effect of the heterogeneous micromechanical stress ﬁelds resulting from the grain-scale anisotropy on the redistribution of hydrogen using a diffusion coupled crystal plasticity model. A representative volume element with periodic boundary conditions was used to model a syn- thetic microstructure. The effect of tensile loading, initial hydrogen content and yield strength on the redistribution of lattice ( C L ) and dislocation trapped ( C x ) hydrogen was studied. It was found that the heterogeneous micromechanical stress ﬁelds resulted in the accumulation of both populations primarily at the grain boundaries. This shows that in addition to the well-known grain boundary trapping, the interplay of the heterogeneous micromechanical hydrostatic stresses and plastic strains contribute to the accumulation of hydrogen at the grain boundaries. Higher yield strength reduced the amount of C x due to the resulting lower plastic deformation levels. On the other side, the resulting higher hydrostatic stresses increased the depletion of C L from the compressive regions and its diffusion toward the tensile ones. These regions with increased C L are expected to be potential damage initiation zones. This aligns with the observations that high-strength steels are more susceptible to hydrogen embrittlement than those with lower-strength.


Introduction
The use of hydrogen as an energy carrier is a main milestone for achieving widespread green and sustainable energy [1] . One of the challenges to realize this milestone is establishing reliable and safe structures for transporting hydrogen [2] . Metals in contact with hydrogen are well known to suffer from loss of integrity known as hydrogen embrittlement (HE) [3][4][5][6] . Several mechanisms can lead to HE; however, for steels, it is widely accepted that two main mechanisms dominate: hydrogen enhanced decohesion (HEDE) and hydrogen enhanced localized plasticity (HELP) [7] . The underlying principle of HEDE is that a threshold of hydrogen con-centration will reduce the cohesive strength between the metal atoms. When these atoms are interface atoms (grain boundaries, carbides, secondary phases etc.), the characteristic intergranular damage of HE occurs, which was recently termed "IG-HEDE" [8] . This is to distinguish it from the less common transgranular decohesion "T-HEDE" in the bulk atoms. The other mechanism HELP, is based on the assumption that hydrogen interaction with dislocations increases their mobility [9] . Macroscopically, this results in plastic softening and shear localization [10] . Another scenario is that hydrogen shields the interaction among dislocations, allowing their pile-up along interfaces and leading to their decohesion [11,12] .
Although HE has been studied for many decades [13] , a general mechanistic explanation is lacking. This is due to the synergistic operation of multiple mechanisms outlined above promoted by the wide range of steel microstructures [7,[14][15][16] . A step toward discerning the interaction of these mechanisms is to understand the hydrogen distribution and the driving forces behind it. Total hydrogen concentration ( C tot ) can be partitioned into lattice ( C L ) and trapped ( C x ) concentrations [17][18][19] . C L resides at the interstitial lattice sites, while C x can be trapped in imperfections like vacancies, dislocations, interfaces etc. The driving force for hydrogen diffusion, the difference in chemical potential, depends not only on the concentration but also on the hydrostatic stresses. Hydrostatic stress gradient drives hydrogen from compressive to tensile regions. Following Oriani [20] , trapped hydrogen and lattice hydrogen are in equilibrium. Limiting our discussion to dislocation traps, the density of trap sites stems from the dislocation density, which evolves with plastic deformation. This highlights the importance of considering mechanical loading when analysing the hydrogen distribution.
Being the smallest atom in the periodic table, hydrogen has very fast diffusivity and low interaction with electrons making it difficult to measure using electron microscopy [21] . Probing hydrogen with high spatial resolution is a technical challenge [21][22][23] especially for samples under stress. Several studies used diffusion coupled mechanics models to study the effect of stresses at a crack tip [18,[24][25][26][27] . Although they provided very useful insights on the distribution of hydrogen at the crack tip, they were unable to provide detailed insights at the microstructural level. Therefore, more knowledge is needed into how hydrogen is distributed between lattice and trapped hydrogen concentrations, which is critical in understanding the competition between different HE modes.
It is well known from empirical observations that high-strength steels are more sensitive to HE than lower-strength steels [28,29] . However, it is impossible to separate the effects of grain size, composition, yield strength, local mechanical loading, and other factors in experimental studies because these cannot be independently controlled in real steels. Therefore, the different contributions of these effects remain unknown. An appealing way to resolve these questions is to use micromechanics based finite element method (FEM). In this approach, an artificial microstructure is generated using representative volume element (RVE) [30] . Recently, several studies reported the effect of polycrystalline anisotropy on the hydrogen distribution using RVE and crystal plasticity model. Ilin et al. [31] investigated the effect of mechanical loading and grain boundary trapping in 316L stainless steel. They found that hydrogen accumulates at regions with high hydrostatic stresses while trapping had negligible effect. Hassan et al. [32] modeled a RVE under biaxial loading and came to similar conclusions. It is worth noting that in both studies, very large initial hydrogen content was used. Charles et al. [33] studied the diffusivity by permeating hydrogen in a predeformed synthetic polycrystal.
In this study, we use a RVE and diffusion coupled crystal plasticity FEM to investigate how the yield strength, grain-scale anisotropy, and mechanical loading impact the sensitivity to hydrogen embrittlement and the localized distribution of lattice and trapped hydrogen in the absence of grain boundary trap. This is done by simulating the effect of the stresses on the redistribution of hydrogen at the microstructural scale in α-iron. A synthetic microstructure subjected to plane strain tensile load was modeled using a periodic RVE with randomly oriented grains. The heterogeneous deformation resulting from the polycrystalline anisotropy was modeled using crystal plasticity. Dislocations were used as the trapping sites in the hydrogen transport model. Understanding the relative partitioning of C L and C x can aid in identifying the type of hydrogen population, and thus, the possible HE triggering mechanism [6,25,34] . We elucidate the effect of initial hydrogen content on this partitioning. A statistical analysis was performed in order to quantitatively address the partitioning of both populations. Furthermore, it is well known that high strength steels are especially susceptible to HE [28,29] . In order to shed light on this behavior, we investigate the effect of higher yield strength on the redistribution pattern of hydrogen.

Representative volume element and grain morphology
The workflow for artificial microstructure, or RVE, generation involves four steps: geometry generation, generation of synthetic orientations, meshing and pre-processing. The artificial grain morphology was generated using an in-house MATLAB script for periodic Voronoi tessellation. The synthetic orientations and all the subsequent crystallographic analysis was performed using the open source code MTEX [35] . For meshing we used the open source code Gmsh [36] . The preprocessing for ABAQUS was performed using a python script for imposing periodic boundary conditions and assigning grain orientations and properties. The Voronoi tessellation is a widely used method for the geometric representation of polycrystals [30,37] . In 2D, Voronoi tessellation is the partitioning of a plane domain into n convex polygons based on a predetermined distribution of n points or seeds. All the points within a polygon are closer to the polygon's seed than the seeds of other polygons. In this work, the plane domain boundary is a square. To use smaller RVEs with the minimal impact on the accuracy of the results, periodic geometry was used. Together with periodic boundary conditions (PBCs), described in details in section 2.2 , the RVEs correspond to infinitely repeating microstructure. Physically, these geometric processes represent a crystallization process where all seeds concurrently nucleate and grow isotropically in 2D.
The morphology of the grains (relative size and shape) can be controlled by controlling the number and distribution of the seed points within the enclosing domain. In order to improve Voronoi tessellation uniformity and reduce extreme aspect ratios, the Voronoi seed points were sampled from a Halton sequence [38] . A RVE with an edge size of 500 μm is shown in Fig. 1 (a). The corresponding grain area histogram is shown in Fig. 1 (b). Random orientations were assigned to the grains. First, an orientation distribution function (ODF) was generated in MTEX. Then individual orientations were randomly sampled from this ODF and assigned to each grain. It should be noted that periodic grains at the boundary edges were assigned similar orientations in order to maintain microstructure periodicity. The {100}, {110} and {111} pole-figures of these sampled orientations are shown in Fig. 1 (c). To assess the randomness of the orientations, a misorientation distribution function (MDF) was fitted to the sampled orientations. A plot of the MDF in Fig. 1 (d) clearly shows a Mackenzie distribution, confirming the randomness of the orientations. The RVE geometry is then meshed with second-order six-node triangular elements using Gmsh.

Periodic boundary conditions
The application of PBCs results in a compatible deformation of each pair of opposite boundary edges, i.e. have identical shape during deformation. For implementing such boundary conditions, periodic nodes (red circles in Fig. 2 ) on opposite edges of the bounding square were enforced during the meshing process. Finally, the mesh file was pre-processed within ABAQUS using a python script to assign material properties and apply PBCs to the boundary nodes. Here we adopt the implementation described in [39][40][41] . This can be summarized as follows. The vertex nodes V 1 −4 (yellow circles in Fig. 2 ) are defined by the intersection of two edges, and thus, will have a different formulation from the inner edge nodes (red circles in Fig. 2 ) and are retained for load application. V 1 is a fixed node, while in-plane displacements u i in directions 1 and 2 can be applied through V 2 and/or V 4 to result in the desired loading state. V 3 is not independent and its deformation is kinematically tied to the other vertices as In this study, we apply uniaxial far-field strain ε 22 in direction 2 on a RVE with an edge length l. This can be formulated as The displacement degrees of freedom u i (x, y ) of the inner edge nodes on each pair of opposite edges are kinematically linked to the vertex nodes V 1 , V 2 and V 4 as where x and y are the coordinates of the boundary nodes. In ABAQUS, Eq. (1) and (4) are implemented as linear constraint equations on the boundary nodes.

Crystal plasticity model
A crystal plasticity model [42] is used for the mechanical behaviour of individual grains. A finite strain framework where the multiplicative decomposition of the total deformation gradient F into elastic F * and plastic F p parts as A slip system α in the initial configuration is defined by its slip direction s α and normal to the slip plane m α . Here, we assume the 12 slip systems {110} 1 11 in BCC iron as the active slip systems for simplicity. In the deformed configuration they are given by The evolution of the slip rate ˙ γ is modeled using the classical viscoplastic formulation [43] where ˙ a α o is the reference strain on slip system α, τ α is the resolved shear stress, or Schmid stress, on that system and g α is the function describing its current strength. τ α is expressed as (s * α m * α + m * α s * α ) : τ τ τ (8) where τ τ τ is the Kirchhoff stress. The current strengths g α evolve as where the summation is over all the active slip systems β. h αβ is the matrix of the hardening moduli with the diagonal terms

Table 1
Crystal plasticity model parameters for α-iron [44] Cubic elasticity constants Flow and hardening parameters A simple form for self-hardening is given by Pierce, Asaro and Needleman [43] h where h o is the initial hardening modulus, τ s is the saturation stress and τ o is the initial critical resolved shear stress. γ is the cumulative shear strain on all slip systems The latent hardening is calculated as where q is the ratio of latent to self-hardening. The parameters of single crystal α-iron were obtained by calibrating the crystal plasticity model with the results from Franciosi et al . [44] are shown in Table 1 . It was implemented in ABAQUS as a user defined material (UMAT) subroutine developed by Huang [42,45] and further modified to be compatible with the hydrogen transport model as discussed next.

Hydrogen transport model
Here we use the hydrogen transport model of Sofronis and McMeeking [18] to include the effect of trapping and hydrostatic stress, and later improved by Krom et al. [24] to include the effect of plastic strain rate. This model assumes only dislocations as trap sites, which are saturable and reversible. The molar concentration (mol/mm 3 ) of hydrogen in lattice sites C L = N L θ L (13) where N L is the number of solvent atoms per unit lattice volume and θ L is the fractional occupancy of lattice sites. Similarly, the trapped hydrogen concentration (14) where N x is the number of trap sites per unit volume and θ x is the fractional occupancy of trap sites. The trapped C x is not independent of the mobile C L , and indeed, there is a local equilibrium between both populations. Oriani [20] expressed the local equilibrium between both populations as where K T is the trap equilibrium constant, W B is the trap binding energy, R is the universal gas constant and T is the absolute temperature. Krom and Bakker [17] showed that this condition can be assumed for strain rates less than 1 s -1 . Combining Eq. 13 -(15) Assuming that hydrogen diffuses only through lattice sites, the hydrogen mass transport equation can be expressed as where J L is the lattice hydrogen flux. For a system under external loading, the hydrostatic stress term σ H in the chemical potential is not constant. High hydrostatic stresses will reduce the chemical potential and result in an influx of hydrogen from regions with compressive hydrostatic stresses to regions with tensile hydrostatic stresses. This can be formulated in the lattice hydrogen flux as where D L is the lattice-diffusion coefficient, V H is the partial molar volume of hydrogen and ∇σ H is the gradient of the hydrostatic stress. As mentioned previously, in this study, we only consider dislocations as trapping sites. Kumnick and Johnson [46] studied the evolution of trap site density in terms of plastic deformation. While their model is empirical in nature, it has been recently combined with RVE-level simulations (e.g. [32,33] ) to provide insight into microrstructural behavior. They found a sharp increase in trap site density in the initial levels of plastic deformation which saturates gradually with further deformation. A fit of their results as a function of the equivalent plastic strain ε p is The time derivative of the trap site concentration can then be expressed as [17] Using Eq. (18) , (20) and (21) , the mechanics coupled hydrogen mass transport Eq. (17) becomes Material parameters for hydrogen diffusion in α-iron [24] Quantity Value The term D e f f in Eq. (22) is called the effective diffusivity and is given by The material parameters used for the hydrogen transport analysis are summarized in Table 2 . Eq. (22) is a non-standard mass transport equation, and its implementation in commercial FEM software requires a special treatment. Oh et al. [47] took advantage of the analogy between heat and mass transfer to implement Eq. (22) in ABAQUS using user-defined thermal material (UMATHT) subroutine. Several works implemented the same methodology for hydrogen transport analysis both at the macroscale [48][49][50] and the RVE level [32,33,51] . A detailed discussion on the implementation can be found in [47,49] . Furthermore, the solution of Eq. (22) requires the calculation of the hydrostatic stress gradient terms ∇σ H , which are obtained from the UMAT subroutine as discussed in the appendix. For solving this coupled mechanicsdiffusion problem, second order triangular elements CPE6MT were used for the coupled analysis in ABAQUS.

Initial hydrogen concentration and loading
The initial partitioning of both C L and C x depends on the initial hydrogen concentration, initial trap site concentration N x and the trap binding energy W B . At equilibrium and for stress-free samples, the solubility of hydrogen in the lattice sites, which for an initial value problem is evaluated as the initial lattice concentration C Lo [18,24] , can be calculated from Sieverts' law as [52] where p is the hydrogen gas pressure, and the solubility of hydro- kJ/mol respectively for hydrogen in iron. At large hydrogen gas pressures, the pressure in Eq. (24) must be replaced with the fugacity [53] such that where z 1 = 1 . 51 × 10 −6 K / Pa and z 2 = −1 . 04 × 10 −11 Pa −1 are fitting coefficients [53] . A plot of Eq. (24) and (25) is shown in Fig. 3 . It can be seen that hydrogen pressures of up to 100 bar, both values almost coincide. In this study, we chose hydrogen pressures of 1, 10 0, 20 0 and 10 0 0 bar, and their corresponding C Lo was calculated using Eq. (25) . From a practical point of view, pressures of 200 bar are observed in the petrochemical industries, and for future use in energy transport, hydrogen pressure of 10 0 0 bar is being considered to maximise the energy density to a practical level [54] . The corresponding values of C Lo are shown in Table 3 .
Unless otherwise mentioned, we assume that all samples are precharged under equilibrium conditions. The corresponding initial trap site concentration C xo can be calculated from Eq. (16) and (19) with ε p = 0 . The total initial concentration is then C tot o = C Lo + Fig. 3. The effect of pressure and fugacity on the equilibrium initial lattice concentration C Lo of unstressed sample at 300 K. 1 . 112 × 10 −10 C xo . A far-field plane strain tensile loading of ε 22 = 3 × 10 −2 was applied to the RVE over the duration of 100 s, which corresponds to a strain rate of 3 × 10 −4 s −1 . The hydrogen concentration at the boundary of the RVE was maintained constant during the loading. This corresponds to tensile testing under constant hydrogen gas pressure [55] allowing for increasing the hydrogen content in the model.

Effect of hydrogen gas charging pressure
Upon loading, a heterogeneous distribution of stress and strain within the RVE develops due to the randomly distributed misorientations and the grain-scale elastic and plastic anisotropy. This results in the redistribution of the initially uniformly distributed hydrogen as shown in Fig. 4 . This concentration distribution does not change with increasing the period over which the load is applied. Our results (not shown here for brevity) showed that increasing the loading time up to 30 0 0 s did not show a significant change in the hydrogen concentration, indicating a steady state condition. This is due to the very fast diffusion of hydrogen in iron as well as the relatively smaller length scale of RVE computations compared to structural scale models. In order to discuss these effects, a detailed study accounting for the rate dependent parameters like diffusion coefficient, loading rate and viscoplastic effects is under investigation. Fig. 4 (a & b) show that the lattice hydrogen C L diffuses toward regions with tensile hydrostatic stresses and depletes from regions with compressive hydrostatic stresses. These values deviate from the initial concentrations shown in Table 3 . For both cases, these regions are primarily the grain boundary triple points. On the other side, the trapped hydrogen C x accumulates at regions of large plastic strains. Fig. 4 (c & d) shows that these regions form bands of approximately 45 • to the loading direction, which is to be expected for plastic shear bands. The intensity of these bands is maximum along the grain boundaries, which manifests the accumulation of hydrogen due to dislocation pile-up [11] . Furthermore, the effect of geometrically necessary dislocations (GNDs) on the hydrogen distribution at a crack tip has been recently reported [56] . Kumar and Mahajan [51] showed that strain gradients lead to the accumulation of GND trapped hydrogen at the grain boundaries.
Next, we discuss the influence of the hydrogen pressure, i.e. the initial hydrogen concentration C tot o , on the partitioning of lattice and trap site concentrations. This is an important factor in understanding the relative contribution of each population to HE. The distribution of lattice hydrogen within the RVE at 1 bar and 10 0 0 bar are shown in Fig. 4 (a) and (b) respectively. It can be seen that the distribution in both cases is almost identical, except that the quantities are two orders of magnitudes higher for the 10 0 0 bar case. This behaviour is expected, as the loading conditions and the RVE parameters were kept constant, while the only variable is the hydrogen pressure according to Eq. (25) . A quantification of this behaviour can be achieved through the distribution of probability density of the integration point values of C L as shown in the first row of Fig. 5 for different hydrogen pressures. C L shows a normal distribution that is almost identical in shape, i.e. its relative standard deviation (RSD) is constant and is approximately 1.75%. However, the mean of the distribution increases with increasing the charging pressure from 2 . 69 × 10 −12 mol/mm 3 for 1 bar up to 1 . 09 × 10 −10 mol/mm 3 for 10 0 0 bar.
In the case of trapped hydrogen, the distributions are similar and there is only a negligible increase in the magnitudes for both pressures as shown in Fig. 4 (c & d). The probability density distribution of C x is shown in the second row of Fig. 5 . C x has a lognormal distribution which maintains its shape (RSD = 69.23%) and mean of approximately 1 . 31 × 10 −11 mol/mm 3 with increasing hydrogen pressure. As plastic deformation increases, more dislocations are created, and consequently more trap sites according to Eq. (19) . Hydrogen diffuses to fill in these trap sites, while maintaining equilibrium with C L . Due to the strong binding energy of −60 kJ/mol, C x is strongly dependent on ε p and its dependence on C L using Oriani's hypothesis in Eq. (15) is negligible. Weaker dislocation trap binding energies, between −40 and −20 kJ/mol, for different types of steels have been investigated [25,34] . How this might affect C x is a subject of future analysis.
The effect of the loading on the hydrogen uptake for both hydrogen populations can be represented by the mean value of the distributions. The mean values for C L are 2 . 69 × 10 −12 , 2 . 76 × 10 −11 , 4 . 01 × 10 −11 and 1 . 09 × 10 −10 mol/mm 3 for hydrogen pressures of 1, 100, 200 and 1000 bar respectively. These values are approximately equal to the initial lattice concentration C Lo shown in Table 3 . The trap site concentration increases with loading from initially 1 . 55 × 10 −12 to 1 . 31 × 10 −11 mol/mm 3 . This indicates that tensile loading does not affect C L uptake, rather it only leads to its redistribution. On the contrary, tensile loading significantly affects C x uptake. Indeed, it was reported that tensile loading has negligible effect on C L uptake and only affects C x uptake for iron [57,58] . Following Di Leo and Anand [59] , Díaz et al. [48] and Martínez-Pañeda et al. [60] used constant chemical potential boundary conditions instead of constant concentration. The dependence of hydrogen solubility on hydrostatic stress in their formulation showed an increased C L . We believe that this will have a minor effect on C L uptake in our model. Unlike their local evaluation of C L at a blunting crack tip, we evaluate C L uptake throughout the whole RVE with heterogeneous stress distribution. Nevertheless, the effect of boundary conditions on C L uptake within a RVE will be thoroughly investigated in our future studies, especially on hydrogen partitioning.
Another way of modeling the interaction of fracture and hydrogen diffusion is to model a discrete fatigue crack tip. However, the diameter of a fatigue crack tip is the same size or smaller than a single grain in the present simulation. We chose to model an RVE in order to capture a homogenized response of the microstructure and to avoid the possibility that the results depend entirely on the behavior of one or two grains located at the crack tip. It is worth noting that some references that use cracks feature an assumed continuum [59][60][61] as opposed to a microstructural model.
Applied state of stress and total strain are also important loading parameters. Increasing the far-field strain beyond 3% will increase both concentrations, especially the dislocation-trapped hydrogen. Increasing the stress triaxiality from the plane strain condition presented here to that of a notch or crack tip would increase the hydrostatic tension and therefore increase the C L . Cracks are also a critical area for hydrogen embrittlement. Strain, stress, and stress state gradient will also dramatically affect the hydrogen distribution, as it will inherently be drawn to areas of high stress triaxiality and plastic strain.
Having described the effect of the charging pressure on C x and C L , we now discuss the implications on the total concentration C tot shown in Fig. 6 . At 1 bar, the distribution and magnitudes are very similar to the distribution of C x . This can be confirmed from the last row of Fig. 5 . The distribution of C tot has a log-normal distribution with values order of magnitude larger than C L indicating that at low pressures, the hydrogen distribution is dominated by the trapped hydrogen. With increasing hydrogen pressure, the effect of C L becomes more tangible. This can be seen from the change in the distribution and magnitudes of C tot in Fig. 6 (a) through (d). This can also be seen from Fig. 5 where the shape of the distribution changes from log-normal toward more normal distribution, i.e. more relative contribution of C L . The spread of the distribution indicates the level of heterogeneity of hydrogen [31] . Increasing hydrogen pressure increases the heterogeneity of hydrogen, and thus, increases the susceptibility to HE.
The micromechanical stress state results in high local concentration of C L and C x at triple points and grain boundaries, respectively. Since the initiation of hydrogen cracking is triggered by a hydrogen concentration threshold [62,63] , these regions are likely to contribute to fracture given appropriate conditions. Connolly et al . [8] discussed that the trapped hydrogen at the grain boundaries due to dislocation pile up can lead to a HEDE mediated HELP damage mechanism. This behaviour will become more prominent with an increasing hydrogen content. This argument can be supported by the experimental observations reporting an increase in the fraction of intergranular fracture with an increase in the hydrogen content [64,65] . Ayas et al. [34] analysed the effect of lattice and trapped hydrogen at different trap sites on the failure of AISI 4135 steel. They concluded that intergranular fracture doesn't correlate with grain boundary hydrogen, and its physical basis needs to be elucidated. Our results show that the micromechanical stress state, even in the absence of grain boundary trapping, can result in the accumulation of hydrogen at the grain boundaries. To dis- cern the effect of different hydrogen populations at grain boundaries, RVE modeling of boundary trapping [66] with appropriate representation of grain boundary character is being considered for future work.

Effect of yield strength
In this section, we describe the effect of yield strength on the redistribution of hydrogen. Here we used yield strength similar to AISI 4340. The crystal plasticity parameters were calibrated by fitting a stress-strain curve as in [67] . For simplicity, we use the parameters in Table 2 for the mass-transport formulation in section 2.4 . We call this hypothetical material high strength steel (HSS). The stress-strain curve from a RVE with these fitted parameters compared to α-iron is shown in Fig. 7 .
The redistribution of hydrogen due to the heterogeneous micromechanical stresses in HSS with increasing hydrogen pressure follows a very similar trend as that discussed in section 3.1 . In Fig. 8 , we only show the hydrogen distribution at 10 0 0 bar of hydrogen pressure for brevity. Fig. 8 (a) shows that C L accumulates in the same regions as those for α-iron in Fig 4 (b). This is to be expected as the only variable we changed in the model was the constitutive behavior, while the microstructure features were kept identical. In fact, this is an advantage of computational models which provide total control over the parameters used in com- parative studies. The magnitudes of the maximum C L regions are approximately three times larger than in α-iron, while the minimum regions are three times less. This is due to the higher yield strength of HSS, which develops larger tensile and compressive hydrostatic stresses. C x shown in Fig. 8 (b) also shows similar distribution as in Fig. 4 (d), however, the maximum value in HSS is  relatively smaller than that in α-iron. This is due to the relatively smaller values of plastic strains in HSS. Although in both cases the displacement boundary conditions are the same, the plastic strains developing within the individual grains of the RVE are smaller due to the higher yield strength in HSS as shown in Fig. 7 .
To put the numbers in perspective, in Fig. 9 we plot the hydrogen distribution for α-iron and HSS with the same limits in the color map. It can be clearly seen that the amount of C L is much larger for HSS than α-iron. Furthermore, the regions with lower C L are more depleted for HSS. This larger heterogeneity can be seen from superimposing the fitted distributions as shown in the first row of Fig. 10 . HSS has a RSD of 17.3%, i.e. 10 times more than for α-iron. Although we used the same RVE and hydrogen diffusion parameters in both cases, we expected that the lattice hydrogen uptake, i.e. the mean values of C L , will increase due to increased hydrostatic stresses resulting from increasing the yield strength. Our results show that higher yield strength doesn't lead to increasing the lattice hydrogen uptake, rather, it leads to larger heterogeneity in its distribution. On the other side, trapped hydrogen in HSS has less heterogeneity and lower mean value of 8 . 35 × 10 −12 mol/mm 3 compared to 1 . 31 × 10 −11 mol/mm 3 αiron. C tot has lower mean values and larger heterogeneity as can be seen in Fig. 9 (f) and the last row in Fig. 10 . Therefore higher yield strength surprisingly results in lower hydrogen uptake. However, it effectively increases the hydrogen distribution heterogeneity mainly due to lattice hydrogen. In this case, trapped hydrogen has less relative contribution. Although the RVE we used represents a single phase microstructure, high strength steels are characterized by having multiple phases like carbides and/or martensite. These hard phases are expected to behave as stress concentration regions that would facilitate damage initiation and propagation due to hydrogen [29] . This was shown both in modeling [68] and experimental [69] studies. Other experimental studies showed a higher hydrogen concentration at martensite and carbide particles [21] . Thus, our approach for using an idealized microstructure can be thought of as a lowerbound with respect to hydrogen content.

Conclusions
In this work we used a representative volume element (RVE) and a diffusion coupled crystal plasticity model to investigate the effect of micromechanical stresses on the redistribution of lattice ( C L ) and dislocation trapped ( C x ) hydrogen concentrations in αiron. The model represented a precharged sample subjected to a constant hydrogen gas pressure. A uniaxial loading was applied to the RVE under plane strain conditions. The redistribution of the hydrogen after loading is due to the micromechanical anisotropy and the random orientations applied to the grains. Furthermore, we studied the effect of increased yield stress in order to represent a scenario for high strength steel. Our main findings are: • Both C L and C x were found to accumulate at the grain boundaries. This suggests that in the absence of to grain boundary trapping, the heterogeneity of the micromechanical stress fields leads to the accumulation of hydrogen at the grain boundaries. The competition between grain boundary trapping and heterogeneity of micromechanical stress fields will be the subject of future studies. • According to Sieverts law, the initial C L uptake is a function of the hydrogen gas pressure. After applying the load, C L redistributes according to the hydrostatic stress fields, where it depletes from compressive regions and diffuses toward tensile ones. The hydrostatic stress fields had a negligible effect on the uptake of C L . On the other side, the hydrogen gas pressure had a negligible effect on the uptake of C x . Instead, C x was found to be strongly dependent on the plastic deformation levels. In order to generalize this behaviour, further studies on the effect of the trap binding energy, loading state and strain hardening behaviour are required.
• Our results show that at low hydrogen gas pressure, C x dominates C tot . The contribution of C L to C tot increases with increasing the hydrogen gas pressure. • Higher yield strength had a negligible effect on the uptake of C L compared to the lower yield strength case. Rather, it led to an increase in the heterogeneity of the redistribution of C L . Higher yield strength resulted in increasing the localized hydrostatic tensile and compressive stresses. The tensile regions were more enriched in hydrogen, while the compressive zones were more depleted. On the other side, higher yield strength resulted in lower C x . This is due to the lower plastic deformation levels associated with higher yield strength. Therefore, the increased heterogeneity of the redistribution of C L can be one of the reasons behind the increased susceptibility to hydrogen embrittlement observed for materials with high yield strength.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. cedure used to obtain strain gradients for strain gradient plasticity [70][71][72] and other studies implementing the hydrogen transport model used in this study [6,4 8,4 9] . In this procedure, the integration points of a second-order element are assumed to be nodal points of a pseudo-element as shown by the dashed triangle in