Theoretical framework for predicting solute concentrations and solute-induced stresses in finite volumes with arbitrary elastic fields

A theoretical model for computing the interstitial solute concentration and the interstitial solute-induced stress field in a three-dimensional finite medium with any arbitrary elastic fields was developed. This model can be directly incorporated into two-dimensional or three-dimensional discrete dislocation dynamics simulations, continuum dislocation dynamics simulations, or crystal plasticity simulations. Using this model, it is shown that a nano-hydride can form in the tensile region below a dissociated edge dislocation at hydrogen concentration as low as χ0=5×10−5, and its formation induces a localized hydrogen elastic shielding effect that leads to a lower stacking fault width for the edge dislocation. Additionally, the model also predicts the segregation of hydrogen at Σ109(13 7 0)/33.4∘ symmetric tilt grain boundary dislocations. This segregation strongly alters the magnitude of the shear stresses at the grain boundary, which can subsequently alter dislocation-grain boundary interactions and dislocation slip transmissions across the grain boundary. Moreover, the model also predicts that the hydrogen concentration at a mode-I central crack tip increases with increasing external loading, higher intrinsic hydrogen concentration, and/or larger crack lengths. Finally, linearized approximate closed-form solutions for the solute concentration and the interstitial solute-induced stress field were also developed. These approximate solutions can effectively reduce the computation cost to assess the concentration and stress field in the presence of solutes. These approximate solutions are also shown to be a good approximation when the positions of interest are several nanometers away (i.e. long-ranged elastic interactions) from stress singularities (e.g. dislocation core and crack tip), for low solute concentrations, and/or at high temperatures.


Introduction
Interstitial solute atoms (such as carbon (C), hydrogen (H), oxygen (O), lithium (Li), etc.) play an important role in controlling the physical and mechanical properties (e.g. yield strength (Gavriljuk et al. 1998; Barrera et al. 2016;Cui et al. 2018), and phase composition (Schwarz and Khachaturyan 2006)) of different metals and alloys. Despite many studies over the years, understanding the effects of interstitial solute atoms on the mechanical response of metals remains challenging. In particular, accurately predicting the solute concentration distribution and the solute-induced stress field in real defected threedimensional (3D) finite volumes remains an unsolved problem. Developing a theoretical framework to quantify the role of solid solution on the mechanical properties of metals is thus a necessity in the process of understanding and predicting the response of metals and alloys.
In 1949, Cottrell and Bilby developed a theory of yielding and strain ageing in iron in which they explained that the segregation of C around dislocations lead to their pinning (Cottrell and Bilby 1949). This was latter termed the "Cottrell atmosphere". In 1982, Barnett et al. studied the pinning force exerted by the solute atmosphere around edge dislocations (Barnett et al. 1982a; Barnett et al. 1982b). They used the Fermi-Dirac statistics over the Boltzmann distributions to correctly calculate atomic solute fractions.
Wolfer and Baskes incorporated the effect of lattice dilatations caused by solute atoms to derive the general solute distribution and the coherency stress field associated with the solute atoms trapped around an infinitely long edge dislocation (Wolfer and Baskes 1985). Furthermore, Sofronis and Birnbaum accounted for the effect of the elastic moduli changes due to solute atoms and used iterative finite element analysis to calculate the stress field around edge dislocations under plane strain conditions (Sofronis and Birnbaum 1995). Chateau et al. used the same technique to calculate the stress field in the vicinity of a crack to study the solute-crack-dislocation elastic interaction (Chateau et al. 2002).
On the other hand, Cai and his co-workers pointed out that it is not correct to include the homogenized self-stresses of the solute atoms to calculate the equilibrium solute distribution. They reformulated the solute concentration by excluding the self-stress inside existing inclusions (Cai et al. 2014). In their work, the coherency stress field induced by the solute atoms around an infinitely long edge dislocation were solved using the Papkovich-Neuber scalar potential. More recently, Gu and El-Awady adopted the solute concentration formulation in Ref. Cai et al. (2014) and developed a three-dimensional (3D) formulation for the solute-induced stresses in an infinite medium to quantify the effect of hydrogen on dislocation dynamics in 3D (Gu and El-Awady 2018). This was then extended by Yu et al. to a finite boundary value problem via the superposition principle (Yu et al. 2019a), which enabled the study of the influence of H on junction formation (Yu et al. 2019;Yu et al. 2019b).
The above studies are all based on deriving the solute concentration in an infinite mediums, where the image interaction between solutes are not considered. Some studies do solve for the image effect but for special geometries (e.g. cylindrical volumes Hirth et al. (2017); Song et al. (2018)).
As discussed above, the applicability of the formulations developed in the existing literature to calculate the solute concentration and solute-induced stress fields are limited to simple dislocation configurations in infinite media or finite media with specific geometries. Moreover, the evaluation of the solute-induced stress usually requires numerical integrations, that typically have a high computation cost. In this paper, a numerical scheme is presented to calculate the solute concentrations and solute-induced stress fields in arbitrary three-dimensional finite volumes with arbitrary pre-existing stress fields. A linearized closed-form solution is also developed to enable the efficient investigation on the role of interstitial solute atoms on deformation mechanisms and subsequently the mechanical properties of the material.

Interstitial solute concentration in a finite medium
Consider the case of interstitial solute atoms in a finite isotropic material, B. Let the purely dilatational volume change induced by a solute atom be denoted by V , while χ, c, and μ 0 be the solute mole fraction, solute volume concentration, and the reference chemical potential, respectively. The solute volume concentration is related to the mole fraction by c = c max χ, where c max is the maximum solute volume concentration. In addition, the chemical potential of solute atoms can be expressed as Sills and Cai (2018): where x is the position vector of any point in B, k is the Boltzmann constant, T is the temperature, σ kk (x) is the trace of the stress tensor σ (x), G is the shear modulus, ν is the Poisson's ratio,( ) indicates the volume average over the entire volume (the volume average of a quantity f is defined asf = B f dV B dV ), Z is the number of nearest neighbors for a solute site, and E ss is the interaction energy of two solute atoms residing at neighboring sites. The five terms on the right hand side of Eq. (1) represent in order of appearance: (1) the reference chemical potential; (2) the entropy associated with the solute concentration; (3) the work against the hydrostatic stress in the volume; (4) the solute image interaction energy; and (5) the solute-solute interaction energy arising from chemical and electronic effects that are short ranged (Wagner 1971;Udyansky et al. 2009;Khachaturyan;Barouh et al. 2014). Here, as a first order approximation, the solute image interaction energy between solutes accounts only for the spatially homogeneous part (Cai et al. 2014).
First, consider the special case of the finite volume, B, being subject only to the imposed traction, T 0 , and displacement, U 0 , boundary conditions, and in the absence of any crystalline defects (dislocations, twin boundaries, grain boundaries, etc.). Here, cracks and notches are considered parts of the boundary conditions. In this case, the elastic field and the stress tensor in B are denoted by ( ) A and σ A , respectively. Thus, according to Eq. (1), the chemical potential in this case is The solute mole fraction in the defect-free state is taken as the reference mole fraction, χ A 0 , and is evaluated from Eq.
(2) at chemical equilibrium (i.e. μ A (x) = 0) to be: where η 1 = 4G(1+ν) . Additionally, the solute mole fraction in a perfect crystal in the stress-free state, χ 0 , is an intrinsic property of the material and is assumed uniform throughout the entire medium. Thus, χ 0 is defined as the intrinsic reference mole fraction and is evaluated from Eq. (3) in the absence of any hydrostatic stress in the volume (i.e. σ A kk (x) = 0) to be: For dilute solute concentrations and in the limit when σ A ij (x)/T 1, both χ 0 and χ A 0 (x) can be related by linearizing and combining Eqs. (3) and (4), as follows Integrating both sides of Eq. (5) over B, and recalling the definition ofχ A 0 andσ A kk , one obtains: Secondly, consider the same finite volume B subject to the same boundary conditions, T 0 and U 0 , but now in the presence of crystalline defects. The total elastic field in this case is denoted by( ). The solute equilibrium concentration in this case can be evaluated by combining Eqs. (1) and (3), such that (Cai et al. 2014): For dilute solute concentrations and in the limit when σ A ij /T 1, the following linearized relationships may be derived for χ andχ: andχ Combining Eqs. (5) and (8), and using Eqs. (6) and (9), the excess solute concentration can be shown to be: The first term on the right-hand side of Eq. (10) is the contribution from the work against the pre-existing elastic field and the solute-solute interaction, while the second term corresponds to the contribution from the solute image interaction energy.
Accordingly, in the absence of any defects, the solute concentration in a finite volume B subject to boundary conditions, T 0 and U 0 , can be computed numerically from Eq. (3) or approximated analytically by Eq. (5) for dilute solute concentrations and in the limit when σ A ij /T 1. Similarly, in the presence of defects, the solute concentration can be computed numerically from Eq. (7) or approximated analytically by Eq. (10) for dilute solute concentrations and σ A ij /T 1.

The solute-induced stress field in a finite medium
The excess solute concentration, which is the difference between the true and the reference solute concentrations, induces an extra lattice distortion in the volume. This is accompanied by a solute-induced elastic field, denoted by ( ) I , in order to preserve the coherency of the crystal lattice. The solute-induced elastic strain field, ε I,el ij , is the difference between the solute-induced strain, ε I ij , and the eigenstrain, e * ij , and is given locally by: Thus, the solute-induced stresses are and the equilibrium equations for this solute-induced stress field are Substituting Eqs. (11) and (12) into Eq. (13), the equilibrium equations with respect to the displacement field, u I i (x), become Using the Papkovich-Neuber scalar potential B 0 (x) (Papkovich 1932;Neuber 1934), the solute-induced filed can be shown from Eq. (14) to be ("Derivation of the solute-induced stress field" section) It should be noted that B 0 (x) is derived without the constraints of boundary conditions. Therefore, a correction field is needed to satisfy the boundary conditions on B. This will be discussed in the next section.
Using the solute concentrations derived in Eqs.
(3), (4) and (7), the solute-induced stress field, σ I ij (x), can be evaluated numerically from Eq. (15). Since the linear elasticity theory usually breaks down at the core of a defect, the solute atoms residing at those trapping sites typically contribute to the solute-induced core effects rather than the elastic interactions. For the sake of clarity, we limit our analysis to the solute absorbed in normal interstitial lattice sites in the small deformation theory.
In the limit of σ A ij /T 1 and for dilute solute concentrations, the linearized solution derived in Eq. (10) can be used to simplify the calculation of the solute-induced stress field and to develop an analytical solution for σ I ij (x). This analytical solution is derived for plane strain problems and 3D problems in the following.

Approximate analytical solution for plane strain problems
For plane strain problems, the Airy stress function φ(x) is utilized ("Airy stress function" section). Here, can be shown to be harmonic functions with respect to the Airy stress function, wherẽ Using Eq. (16), the scalar potential B 0 (x) can be shown to be: Combining Eqs. (17) and (37), the solute-induced stress field in plane strain problems can be shown to be: It should be noted that the components of the solute-induced stress tensor can be related to the corresponding components of the stress tensor resulting from the imposed boundary conditions and the pre-existing defects through the following relationships:

Approximate analytical solution for three-dimensional problems
For general 3D problems with prescribed boundary conditions,σ kk (x) − σ A kk (x) and σ kk −σ A kk are harmonic functions with respect to the Galerkin vector representations ("Galerkin vector representation" section): Accordingly, B 0 (x) can be determined from Eq. (36): The solute-induced stress field in 3D problems can thus be computed from Eq. (37), for a given Galerkin representation of the pre-existing defects: and the corresponding stress components are

The boundary value problem formulation
Given the imposed traction, T 0 , and displacement, U 0 , boundary conditions on the finite volume, B, the boundary value problem (BVP) formulation can be expressed as follows: where n is the outer normal unit vector to the boundary ∂B.
As schematically shown in Fig. 1, by imposing the superposition method (Van der Giessen and Needleman 1995), the elastic field in the finite volume can be decomposed as follows: σ =σ + σ I,tot .
In Eq. (25),( ) represents the elastic field induced by both the imposed boundary conditions and the presence of crystalline defects in B. It satisfies the following BVP: The( ) field is the superposition of the infinite-body field induced by the crystalline defects and the corrective field accounting for the boundary conditions. For some types of defects, the infinite-body fields are known analytically (e.g. dislocations (Hirth and Lothe 1982) and grain boundaries (Li 1972)). However, this is not the case for general defects and those solutions should be obtained numerically.
Additionally, ( ) I,tot in Eq. (25) represents the solute-induced field, which is related to the excess solute concentration based on Eq. (12), and is determined by the following BVP: It should be noted that the( ) field is used in Eq. (7) to evaluate the solute mole fraction distribution, χ(x). On the other hand, the ( ) A field in the defect-free state subject to the same boundary conditions is needed in Eq. (3) to calculate the reference mole fraction, χ A 0 , and is defined by the following BVP: Furthermore, as illustrated in Fig. 2, the ( ) I,tot field can be decomposed into two fields: the ( ) I field in an infinite volume generated by the excess solute concentration, and thê ( ) field to enforce the actual boundary conditions on B. The ( ) I filed can be numerically The total elastic field is thus the superposition of the following three fields: the( ) field, the ( ) I field, and the( ) field. In this case, the BVP for the( ) field is where the boundary conditions are determined by the ( ) I field: The BVP described by Eqs. (26), (28) and (27) can be solved using the finite difference method (FDM) (Sadd 2009), finite element method (FEM) (Zienkiewicz and Taylor 2005), or boundary element method (BEM) (Becker 1992). The solution can also be approximated in the case of dilute solute concentrations and when σ A ij (x)/T 1 by applying the linearized results in "The solute-induced stress field in a finite medium" section. The complete numerical and the approximated solution schemes are summarized schematically in Fig. 3.
In particular, for plane strain problems with traction boundary conditions only (i.e. ∂B = ∂B f ), the total solute-induced stress, σ I,tot ij , and total stress, σ , are directly calculated using Eq. (19), and:

Numerical examples
In this section, the formulation developed above is utilized to solve three demonstrative numerical problems that are of general interest for understanding the effect of hydrogen (H) on the response of metals. In all cases, the material properties of FCC Ni are adopted: G = 134GPa, ν = 0.201, V = 1.4Å 3 , a 0 = 3.524 Å, Z = 12, E ss = −0.026 eV, and atomic volume = 10.77Å 3 . Additionally, H solutes are assumed to occupy octahedral interstitial sites. The transient diffusion process is not modeled here and the H concentration is assumed to be in equilibrium with the pre-existing stress field. All FEM calculations are conducted using the MATLAB PDE Toolbox™(a 2018) and the MATLAB FEM code (Alberty et al. 2002).

Hydrogen segregation to dislocations and its effect on the dislocation stacking fault width
Dislocations in FCC metals typically split into two partial dislocations separated by a stacking fault. The separation distance between the two partial dislocations, known as the stacking fault width, has strong effects on dislocation mechanisms (e.g. glide and crossslip) and the mechanical properties of the crystal (Aubry and Hughes 2006). Atomistic simulations show that H atoms affect not only the stress states, but also the stacking fault energy (SFE) landscapes and subsequently the stacking fault width (Lu et al. 2002;Wen et al. 2007;Tang and El-Awady 2012). This is the proposed mechanism for H-induced slip planarity (Wen et al. 2007;Lu et al. 2001;Barnoush and Vehoff 2008) and H-assisted fracture process (Song et al. 2010;Tarzimoghadam et al. 2017). In this section, the numerical framework developed in "The boundary value problem formulation" section is used to evaluate the effect of H on the stacking fault width of an infinitely long dissociated edge dislocation.
Consider an infinitely long straight edge dislocation positioned at the center of a 2D square simulation cell, with its Burgers vector being parallel to the simulation cell horizontal direction. The edge dislocation line direction in this case is perpendicular to the 2D simulation cell plane. The perfect dislocation is assumed to decompose into two partial dislocations, p 1 and p 2 , that are separated by a stacking fault, as shown schematically in Fig. 4. The repulsive elastic force between the two partials is denoted by f p , and the H-induced elastic force between the partials is denoted by f H . These two forces are calculated using the Peach-Koehler equations: f = (σ · b) × ζ , where σ is the corresponding stress field, b is the Burgers vector, and ζ is the line direction of the dislocation. The stress field of a partial dislocation was derived using the Peierls-Nabarro (PN) model   (Hirth and Lothe 1982) and the H-induced stress field can be computed as discussed in "The boundary value problem formulation" section. It should be noted that both f p and f H vary as a function of the separation distance. Additionally, the stacking fault area gives rise to a constant attractive force (per unit length) between the partials, f SF = γ ζ × (0, 0, 1), where γ is the SFE. The partials are in equilibrium with each other when the elastic forces on each partial balance the attractive SF force, i.e. f p + f H + f SF = 0 (Cai and Nix 2016). The equilibrium separation distance between the partials (i.e. the stacking fault width) is denoted by l sep . In the current analysis we assume that the elastic constants do not change due to the presence of H. Additionally, the SFE is taken as γ = 0.16 J/m 2 (von Pezold et al. 2011).
The predicted separation distance at T = 300K as a function of χ 0 is shown in Fig. 5a. In the absence of H the predicted stacking fault width is 1.85nm, which is in good agreement with that predicted from atomic relaxations of an ideal dislocation configuration obtained from linear elasticity theory using embedded atom method (EAM) potential fitted to the same SFE used in the current analysis (von Pezold et al. 2011). Furthermore, when H is introduced in the simulation cell, the predicted stacking fault width does not change for χ 0 ≤ 2.5 × 10 −5 . In Fig. 5b, the equilibrium H concentration contours around the dissociated dislocation for χ 0 = 2.5×10 −5 is shown and indicates that while there is some segregation of H at each partial dislocation, the H-induced elastic force by this segregation is negligible and has no effect on the dislocation stacking fault width.
On the other hand, when χ 0 is increased further to 5×10 −5 , the H-induced elastic forces increase, which leads to decreasing the stacking fault width to .55nm. Nevertheless, this width remains unchanged with further increase in χ 0 . The equilibrium H concentration contours around the dissociated dislocation for χ 0 = 5 × 10 −5 is shown in Fig. 5c, and it shows the formation of a nano-hydride (the region corresponding to χ ≈ 1) in the tensile region below the dissociated dislocation. The elastic interaction between the hydride and the partials, f H , is comparable to f p and f SF . Therefore, the formation of nano-hydride induces a localized H elastic shielding effect even at such a low χ 0 . These predictions are in excellent agreement with those from combining ab initio electronic structure calculations, semi-empirical embedded atom method potentials, and a lattice-gas Hamiltonian (von Pezold et al. 2011).
It should be noted that the effect of H on the elastic moduli change is not considered in the current analysis. Nevertheless, the elastic constants at high H concentrations can also strongly influence the separation distance, and thus, should be accounted for in the future to make more accurate predictions of the effect of H on the separation distance. This is, however, beyond the scope of this paper.

Hydrogen segregation to grain boundaries
Grain boundaries (GBs) are planar defects that play an important role during hydrogenenhanced plasticity (e.g. slip transfer across grain boundaries that cause intergranular fracture (Lee et al. 1990;Wang et al. 2016) and enhanced dislocation emission from grain boundary sources at high deformation levels (Kuhr et al. 2016)). Additionally, GB engineering can also lead to improved hydrogen embrittlement (HE) resistance in metals (Kwon et al. 2018). In this section, we use the framework developed in "The boundary value problem formulation" section to gain further insights on the role of GBs on the HE resistance.
The GB structure can be described by the disclination structural unit model (DSUM) (Li 1972). A GB is represented by a contiguous and alternating sequence of majority and minority structural units with associated misorientation angles (Wang and Vitek 1986;Hurtado et al. 1995). The stress field induced by the GB in an infinite medium is therefore calculated as the summation of the majority and minority structure units. Consider a square bicryastal containing a centered 109(13 7 0)/33.4 • symmetric tilt GB in the vertical direction, as shown schematically in Fig. 6. This GB can be represented by a periodic 6 majority unit with a characteristic length 0.3938nm and 1 minority unit with a characteristic length 0.2387nm. In all the following analysis the surfaces of the bicrystal are Page 13 of 22 assumed traction-free. The problem is solved using the superposition method described by Eq. (25).
The H concentration contour plots for different H mole fractions, different temperatures, and different bicrystal edge lengths, are shown in Fig. 7. The H concentration as predicted by the direct numerical solution of Eq. (7) is shown in the left column of Fig. 7 and is designated as the benchmark solution, while that predicted using the linearized solution given by Eq. (10) is shown in the middle column of Fig. 7. The contour plots of the relative errors between both cases are shown in the right column of Fig. 7.
These results indicate that H segregation to the GB dislocations is in good qualitative agreement with molecular dynamics simulation results (Chandler et al. 2008). While the benchmark solution and the linearized solution both predict similar H concentration contours, the relative errors between both solutions are highest within a 1nm thick slab around the GB plane. However, outside this local region both solutions converge with the relative error being below 0.01.
Additionally, the predicted total H-induced shear stresses (i.e. σ I,tot xy ) calculated numerically from Eq. (27) (the benchmark solution) as well as that calculated using the linearized solution in Eq. (31) for different H mole fractions are shown in the left and middle columns of Fig. 8, respectively. The relative errors between both solutions along the perpendicular direction to the GB and at different heights are shown in the right column of Fig. 8. It is observed that the relative errors between both solutions are highest within a 2nm thick slab around the GB plane, while outside this local region both solutions converge.
Moreover, the results in Figs. 7 and 8 indicate that the relative errors for both the H concentration and H-induced stress increase with decreasing temperatures, increasing intrinsic reference H mole fractions, or with decreasing bicrystal sizes. Nevertheless, the linearized solutions for H concentrations and H-induced stresses provide very reasonable approximations in large simulation volumes, with low intrinsic reference mole fractions, Finally, it should be noted that the local stress state and any compositional and structural modifications in the GB can play a key role in characterizing slip transmissions across grain boundaries (Lee et al. 1990;Wang et al. 2016) and establishing the conditions for environment-induced intergranular fracture (Robertson et al. 2015). Along those lines, Fig. 9 shows the contours of the GB dislocations-induced stresses, the H-induced shear stresses and the total shear stresses in the bicrystal with χ 0 = 0.005, L x = L y = 26nm, and T = 300K. It is observed that the magnitudes of the shear stresses are strongly altered along the GB in the presence of H. The different changes of shear stresses due to the existence of H at different GB positions lead to varying capabilities to promote dislocation slip transmissions across the GB. These results indicate that the current numerical framework could be used to systematically study the H effect on slip transmissions, which is beyond the scope of the current paper. Here, the H distribution is predicted in a rectangular plate having a mode-I central crack under plain strain conditions as shown schematically in Fig. 10. The plate edge lengths are L x = 1μm and L y = 2μm, and the crack length is L c . The plate is loaded uniformly in the y-direction under a constant axial stress σ 0 , and the temperature is fixed at 300K. The overall stress field in the H-charged medium in the absence of pre-existing defects is expected to be the same as that for the H-free case (i.e. σ ij =σ ij = σ A ij ) since both cases are modeled by the same governing equation, Eq. (24), and the elastic constants are assumed to be unchanged throughout the whole plate.
The contours of the H concentrations in the plate for different intrinsic reference H mole fractions, crack lengths, and imposed axial stresses are shown in Fig. 11 as predicted from the benchmark numerical solution of Eq. (7) (left column) and from the linearized approximation given by Eq. (10) (middle column). The relative errors, defined as |χ−χ benchmark | χ benchmark , between both solutions are shown in the right column. It is observed that the maximum relative error when σ 0 = 10MPa is well below 5 × 10 −3 , which implies that Eq. (10) is a good approximation for small stress cases. However, the relative errors increase when the imposed axial stress increases. Furthermore, it is observed that the relative errors near the crack tip (i.e. within a distance of 200nm) are typically higher than those away from the crack tip. These relatively large errors suggest that the linearized approximation may lose some degree of accuracy around the regions of large stress levels. Nevertheless, the errors throughout the whole plate are small enough to be neglected in this case.
Furthermore, Fig. 11 shows that H concentration at the crack tip increases with increasing external loadings, higher intrinsic H concentration, or larger crack lengths. The larger H concentration at the crack tip for larger crack lengths will promote further acceleration of the crack. This is consistent with experimental reports that in the presence of H the crack growth rate increases under cyclic loading (Matsuda et al. 2016).

Conclusions and final remarks
Here a theoretical model for computing the solute concentration and solute induced stress field in a three-dimensional finite medium with any arbitrary elastic fields was developed. This method expands on a previous three-dimensional model for solutes in an infinite medium with dislocation fields only (Gu and El-Awady 2018). The developed model can be incorporated into two-dimensional or three-dimensional discrete dislocation dynamics simulations (Gu and El-Awady 2018;Yu et al. 2019a;Yu et al. 2019b;El-Awady et al. 2016;Lavenstein and El-Awady 2019), continuum dislocation dynamics simulations (Zhu and Xiang 2015), or crystal plasticity simulations (Castelluccio and McDowell 2017;Castelluccio et al. 2018). For instance, in the crystal plasticity model developed by Castelluccio et al. (2018), the flow rule can be improved by introducing a weighted average of edge and screw dislocation densities since the solute effects are dependent on the dislocation character.
The model was used to quantify the concentration of hydrogen solutes around an edge dislocation in Ni and the effect of hydrogen on the stacking fault width, hydrogen segregation at a 109(13 7 0)/33.4 • symmetric tilt grain boundary and the effect of hydrogen-induced stresses on the local stress concentration near the grain boundary, and the hydrogen segregation at a mode-I central crack and the dependence of the concentration on the crack length. The result suggests that a nano-hydride forms in the tensile region below the dissociated edge dislocation at hydrogen concentration as low as χ 0 = 5 × 10 −5 . The formation of such nano-hydrides induces a localized hydrogen elastic shielding effect on the dislocation partials that leads to the formation of a lower stacking In all cases, L x = 1μm and L y = 2μm. The horizontal black line represents the crack fault width. These results are in excellent agreement with those from the combination of ab initio electronic structure calculations, semi-empirical embedded atom method potentials, and a lattice-gas Hamiltonian (von Pezold et al. 2011). Additionally, the analysis show the segregation of hydrogen at the grain boundary dislocations. This segregation strongly alters the magnitude of the shear stresses at the grain boundary, which can alter dislocation-grain boundary interactions and subsequently dislocation slip transmissions across the grain boundary. Finally, the hydrogen concentration at a mode-I central crack tip is shown to increases with increasing external loadings, higher intrinsic hydrogen concentrations, or larger crack lengths. The larger hydrogen concentration can promote further acceleration of the crack, which is consistent with experimental reports that in the presence of H the crack growth rate increases under cyclic loading (Matsuda et al. 2016).
A linearized approximated closed-form solution was also derived. The linearized closed-form solution can also effectively reduces the computation cost to assess the concentration and stress field in the presence of solutes, which enables large scale simulations on the effect of solutes. The linearized closed-form solution is also shown to be a good approximation of the numerically solved solute concentration and solute-induced stress field when the positions of interest are several nanometers away (i.e. long-ranged elastic interactions) from stress singularities (e.g. dislocation cores and crack tips), for low solute concentrations, and/or at high temperatures. On the other hand, the linearized solution fails when dealing with short-range interactions (e.g. dislocation junction formation and destruction, dislocation jog-formation, etc.) or when dealing with materials at high solute concentrations.
In addition, it is important to note that the volume change induced by a solute atom is not always purely dilatational as has been assumed in the current formulation. Thus, a more general expression of the work required to insert a solute atom, W sol , is needed. In this case, to rigorously solve the solute-induced stress field, numerical methods (e.g. anisotropic truncated kernel method (Jiang et al. 2017)) are required to solve the following Papkovich-Neuber scalar potential (for an infinite medium): While the solute effect on elastic moduli is neglected in the current analysis, the presence of high local solute concentration will have an influence on the elastic moduli (Wen et al. 2011). To take this effect into account, the elastic interactions arising from the local modulus changes per unit volume 1 2 (C ijkl − C ijkl )ε ijεkl should be included in the equation of chemical potentials, where C ijkl and C ijkl are locally-changed elastic constants in the absence and presence of solute atoms, respectively (Sofronis and Birnbaum 1995). The elastic moduli change due to the solute concentration induces additional inheterogeneity in the stress field, which in turn affects the solute concentration. To find the equilibrium solute concentration, this problem needs to be solved iteratively.
Finally, another simplification adopted in the current model is that the image effect on the solute atom is spatially homogenized. The complete evaluation of the piece-wise image effects requires the solution on specific boundary geometries (Hirth et al. 2017;Song et al. 2018).
Substituting Eq. (33) into Eq. (14), the following Poisson's equation is obtained When c(x) ≡ c A 0 , Eq. (34) is valid if the following relation holds : Substituting Eq. (10) into Eq. (35), one obtains Ifσ kk and¯σ kk are harmonic functions, Eq. (36) can be directly solved by canceling out the Laplace operators in both sides. The solute-induced stress field is then expressed as

Airy stress function
It was first proposed by Maxwell (1870) and then rigorously proved by Rostamian (1979) that all elasticity solutions admit the following Maxwell stress function representation: where is the Maxwell stress function being a symmetric second-order tensor whose all off-diagonal elements are set to zero, i.e.