Effect of a Long-Range Dislocation Pileup on the Atomic-Scale Hydrogen Diffusion near a Grain Boundary in Plastically Deformed BCC Iron

1. Introduction

Compared to rail, trucks and tankers, oil pipeline remains as the most economical and efficient means for crude oil transportation in term of the routes flexibility and the transported quantities. As a mixture of solid, liquid and gaseous phases, crude oil contains sediments, water, salts, carbon monoxide, and acid gases such as $H_{2}S$. The presence of $H_{2}S$ facilitates the H diffusion into pipeline steels, which interacts with material defects, such as dislocations and grain boundaries (GBs), resulting in a hydrogen-induced cracking (HIC) and then oil leakage. It was reported that more than 50% of the failure experienced in oil and gas production industry is caused by HIC in oil pipelines [1]. In order to combat with oil spills, a computational framework that can predict the onset of HIC in polycrystalline oil pipeline materials is needed. Up to date, such a predictive computational framework, however, has not been well established yet due to an intrinsic complexity associated with the coupling between dislocation plasticity and HIC, which occurs across a broad range of length scales. That is: (a) the H-enhanced localized plasticity (HELP) [2,3] spanning tens of micrometers when a large number (tens or hundreds) of dislocations are blocked by an atomically structured GB [4,5,6]; (b) the buildup of a high internal stress field ahead of such a microscale pileup not decaying to zero until several micrometers away from pileup tip [7,8]; and (c) the H-enhanced decohesion (HEDE) [9] of the GB at the atomic scale if such a long-range high internal stress cannot be fully relieved by a dislocation cross-slip, transmission, twinning, phase transformation, hydride formation, and so on [10,11,12]. Clearly, a mechanistic understanding of the fundamental mechanisms underlying HIC in plastically deformed polycrystalline materials necessitates computer simulations ranging from the atomistic up to the continuum level.

At the atomic scale, molecular dynamics (MD) has been one most powerful tool for studying H diffusion, dislocation motion, GB structure evolution, as well as their interaction with each other. Through tracking the motion of each atom in the simulation cell, MD retains the atomistic nature of materials. It is thus applicable to link GB and dislocation core structural details with the local H diffusivity, which play a key role in HIC but difficult to resolve in experiments. For example, MD simulations were performed to: (1) provide a direct evidence that H diffuses significantly slower around a screw dislocation core than it does in lattices away from the core [13]; (2) quantify the effects of dislocation types, GB and triple junction structures on the H diffusion barrier [14,15]; (3) discover a new mechanism responsible for the ductile-to-brittle transition in bcc iron referred as the H aggregation-suppressed dislocation emission from the crack tip [16]; (4) help explain the H-induced pop-in load reductions in the force-depth curves of a H-charged sample under indentation [17,18]; (5) identify the possible interstitial H sites at a given GB by dividing it into polyhedral packing units with the H segregation energies at these sites being calculated and then fed into the Rice-Wang thermodynamic theory of interfacial embrittlement [19]; (6) survey the possible nanoscale mechanisms responsible for HELP through quantifying the effects of H on edge dislocation mobility [20]; and (7) analyze the interactions between dislocation and the H-segregated GB [21]. Despite its great popularity, the results from MD simulations need to be understood cautiously because, due to the high computational cost, typical MD simulations employ only one single or a few dislocations very near a GB embedded in a simulation cell with a limited volume [21,22,23,24]. With such a nanoscale simulation cell, the dislocation density becomes unrealistically high. In addition, in many MD simulations, only equilibrium GBs (EGBs) composed of periodic structure units have been considered. As a consequence, the long-range structure heterogeneity of non-equilibrium GBs (NEGBs), especially the GBs charged with H, has been cut off. In such scenario, the forces induced by the interaction between periodic images of dislocations or NEGBs are non-negligible, which may have polluted the results [24,25]. As such, nanoscale MD simulations alone might not be sufficient if one relies on them to develop/calibrate the higher-scale constitutive rules for predicting how the H-charged polycrystalline materials behaves in engineering practices.

By contrast, with a capability of scaling up in length continuum models have been widely used in simulating the deformation behavior of polycrystalline materials exposed to complex environments. In the past several decades, a variety of such models have been developed for characterizing the coupling between plastic flow, H diffusion, and HIC. Here, for the first time, we categorize them into four types as follows.

Type-I: one-way coupling models without considering any underlying grain structures. This type of models follow a pioneering approach by Sofronis and co-workers [26,27,28]. In such models, the H diffusion is considered in conjunction with the material’s elasto-plastic deformation. A two-field boundary value problem is formulated. The displacement field and the H concentration are the two primal kinematic variables. As a first step of such approaches, the traditional elasto-plastic constitutive rules, such as the traditional J2 flow rules, are used in solving for the displacement, strain, and stress fields. Thereafter, H diffusion equation is solved to determine how H is distributed under those stresses and plastic strains. Numerically, this is usually realised via the standard finite element (FE) procedures in a readily available software, such as ABAQUS, through two steps: (a) computing the displacement, strain, stress, and a gradient of them by calling UMAT subroutines; and then (b) solving the H diffusion equation for the H distribution by calling UMATHT subroutines which impose the hydrostatic stress and plastic strains calculated from UMAT on the diffusion. Such models are believed to be "one-way coupling" because only the mechanical effect on H diffusion is considered while the H-induced material property degradation is not. Moreover, such models are always built upon four key assumptions: (i) the material microstructure complexity is smeared out; (ii) H in these models is considered to reside at either normal interstitial lattice sites (NILS) or trapping sites (GBs, dislocation cores, and so on). The H occupancies at NILS and those trapping sites are considered to be always in equilibrium according to Oriani’s theory [29]; (iii) there are no interactions/communications between those trapping sites; and (iv) they only take the reversible shallow traps into account, and thus can be used to estimate the plastic strain-induced concentration of H trapped at dislocations only. In these models, all the deep traps, such as GBs, are ignored and assumed to be filled with H, which obviously may not be the case in reality;

Type-II: two-way coupling models without considering any underlying grain structures. As an extension of the previous one-way coupling model, in this type of approaches, the effect of H on the material property is added through including the H-induced volume dilatation or flow stress reduction into the framework [30,31,32]. By accounting for the effect of H on the materials’ constitutive response, such models could capture the onset of plastic instability in a H-charged materials under deformation. The predictive capability of them is, however, still limited because, in addition to the same assumptions as that in Type-I models, these models further assumed that H always induces a flow stress reduction and a deformation localization. From the microscopic point of view, this assumption argues that the presence of H will only enhance the mobility of dislocations. According to the recent MD simulations [33], this, however, might not be the case. From the macroscopic point of view, this assumption implies that the presence of H will always cause a materials softening. It is inconsistent with the results from many experimental studies, where both the H-induced hardening [34,35] and H-induced softening [2,36,37,38,39,40] were observed;

Type-III: one-way coupling models for polycrystalline materials with the underlying grain structures being explicitly included. In these models, H transport and the elasto-plasticity are simulated through coupling diffusion with crystal plasticity finite element (CPFE) [41,42,43,44,45,46,47]. Particularly, the effects of hydrostatic stress and plastic strain on the H transport were considered when solving the diffusion equation, while the H-induced material property changes in CPFE were not included. Based on the same assumptions as that in [26,27,28], such approaches share similar limitations with Type-I models but gain in one aspect: H diffusion has been linked to the stress fields calculated from CPFE, rather than that from the standard J2 theory [48,49,50]. With the slip and the elasto-plastic anisotropy being accommodated, these models shed light on assessing the effects of material microstructure on the H diffusion, although only in a qualitative sense due to its ’one-way coupling’ nature and a neglect of the H diffusion along dislocations and GBs;

Type-IV: two-way coupling models with the underlying grain structures being explicitly retained [51,52,53,54,55,56]. If one desires to capture of all the salient features of HIC, a fully two-way coupling model is needed and should at least contain four components: (1) a model that can explicitly characterize the heterogeneous stresses/strains induced by the complex grain structures; (2) a diffusion model that is decorated with local stresses/strain for simulating the H occupancy equilibrium between NILS and trapping sites in materials under deformation; (3) a set of constitutive rules for describing the H-induced property degradation, such as the H-induced material softening/hardening [32,57], GB cohesive strength reduction [51,52,53], and so on; as well as (4) a set of constitutive parameters about the H diffusivity, such as the pipe diffusion along the dislocations and GBs. Nevertheless, to the best of our knowledge, such a fully two-way coupling model equipped with all of those four components does not exist yet. For instance, in a pioneering work by Rimoli and Ortiz [51] and a recent work by Pu et. al. [53], a CPFE model for polycrystals, a Fisher model for GB diffusion, as well as a cohesive zone model for intergranular fracture were all integrated together but with a consideration of neither pipe diffusion nor the H concentration dependence of the CPFE model parameters. As an expansion of [51], in [52], the effects of initial dislocation density and its evolution on HIC were added. Comparing with the traditional CPFE in [51], such a framework indeed gained more physics about how the dislocation slip reacts with the GBs in the presence of H, although with a neglect of the H diffusion along the GBs. By contrast, in a recent diffusion-stress coupled model [54], the GB diffusion under the effect of a normal cohesive traction was explicitly included, but the effects of H on the slip inside the grains were ignored. Obviously, existing two-way coupling models usually include only a few among those four components. Moreover, even with those four components being all equipped, Type-IV models will still require the constitutive rules, such as the H concentration-dependent CPFE model parameters, GB cohesive strength, the local stress-dependent GB diffusivity, and so on, to be carefully calibrated from either experiments or fine-scale simulations. This is non-trivial but is definitely necessary if one desires to achieve a high predictive capability.

Overall, neither continuum nor atomistic models alone can capture all the salient physics underlying the interplay between H diffusion, plastic flow, and the subsequent HIC. Multiscale models are needed to overcome many inherent difficulties in them. To meet this need, a concurrent atomistic-continuum (CAC) methoddeveloped in [58,59,60,61,62,63,64,65] is deployed here. The CAC method is built upon an atomistic field formulation [66]. The governing equations in this formulation are identical in form to that in continuum theories. Thus, continuum modeling techniques, such as FE, can be used to solve them. The utilization of FE in regions where materials deform cooperatively leads to a coarse-grained (CG) model [59,60,62,67]. The combination of CG and atomistic leads to a CAC model. The resulting CAC model was shown to be effective and efficient in predicting: (i) the heterogeneous stress field evolution in polycrystalline SiC with different chemical dopants on the GB [68] and (ii) the formation of a microscale dislocation pileup followed by a atomic-scale dislocation transmission [69], twinning, and phase transformations (PTs) [70,71]. In this paper, taking the plastically deformed bi-crystalline bcc iron sample as a model system, the reaction between a 0.2 $\mu$m long dislocation slip and an atomically structured GB in the presence of H is simulated by CAC. With the interatomic potential being the only constitutive rule, at a fraction of the cost of full MD simulations, the coupled dynamics between the microscale dislocation slip, the long-range internal stress buildup, and the atomic-scale H diffusion at the GBs is characterized. In the rest of this paper, we first briefly introduce the methodology and the computer model set-up in Section 2. CAC simulations results are then presented in Section 3.1. Thereafter, we conclude this paper with a summary of the major findings and a discussion of future research in Section 4.

2. Computer Model Set-up

The CAC model for a bi-crystalline bcc iron sample containing a GB charged with H is shown in Figure 1. One main novelty of such a CAC model is to bridge several length scales relevant to HIC in one single model. In details, it explicitly accommodates the GB at an atomistic resolution while the lagging dislocations away from GB in a coarse-grained (CG) description using FEs. In this way, a majority of the atomistic degrees of freedom (DOF) are eliminated by assuming a collaborative deformation at the lattice-cell level in the CG domains, where the displacements of atoms will be constrained by the FE nodal displacements. Such a computational set-up is multiscale in nature and has two unique features: (i) it requires the interatomic potential as the only constitutive rule, such as the EAM potential for the Fe-Fe and Fe-H interactions [16] used in this work; and (ii) it can scale up to microns and even above but retains an atomistic resolution nearby the critical region (GB in this work) in the materials. A variety of defects, including screw/edge dislocations, twist/tilt GBs, and crystallographic orientation-dependent phase boundaries (PBs) can be constructed as desired. It should be noted that the H segregation from the far region into the GB is not simulated here. As such, a certain number of H atoms at a concentration of 0.5% are directly introduced into a region in the vicinity ( 4nm on both sides of the GB) of the atomically structured GB.

The lattice orientations are chosen to be $x\left[111\right]$, $y\left[10\overline{1}\right]$, and $z\left[1\overline{2}1\right]$ in the left grain, and $x\left[111\right]$, $y\left[1\overline{2}1\right]$, and $z\left[10\overline{1}\right]$ in the right grain, respectively. In this way, one of the slip planes, along which the dislocation pileup will be formed in the left grain, is along y direction and perpendicular to the GB plane. The crystalline grains in the domains 4nm away from the GB are discretized into coarse FEs. Each FE contains 512 atoms. These FEs are carefully configured (see details in [58,67]). They can slide with respect to each other for accommodating the dislocation slip along the FE boundaries. Such a CG description of dislocation was shown to be a first model of its kind that can capture the collective behavior of many µm-long dislocation lines without smearing out the atomic-level kink dynamics along each of them [67]. For the sample with a dimension of 95 nm×300 nm×85 nm, such a CAC model contains 381,000 FEs in the two grains and 5.9 million atoms on the GB, which is equivalent to 201 million atoms in a full MD model. In this way, the number of the computational degree of freedom of such a CAC model will be only 5% of that in full MD simulations. This in turn, provides us with opportunities on quantifying how a pileup containing a few to tens of dislocations may affect the atomic-level H diffusion at the GB.

A queue of dislocations with a uniform separation in between are initially introduced into middle plane of the left grain (see the location of one dislocation nearby the GB in the inset picture of Figure 1). As demonstrated in our previous work [58,59,60,61,62], because the CG domain in a CAC model can accommodate dislocations without smearing out their atomistic natures, the dislocations are built into the CG domain through displacing FE nodes according to the elasticity-based solution fof the displacement field around a screw dislocation. Different from traditional FE models in which the neighboring elements are connected by sharing the FE nodes, the FEs in CAC models across the slip plane are disconnected (see the zoom-in view of an inset picture in Figure 1). Thus, the displacement jump induced by a dislocation along the slip plane will be allowed. After dislocations are introduced into the model, the upper and bottom boundaries of the model is moved along $x\left[111\right]$ but in an opposite direction to impose a shear on the sample. Thereafter the increase of each shear strain, the top and bottom ends of the sample are then fixed for a finite duration to equilibrate the dislocation configuration in the pileup. Since the H diffusion at the GB is concerned here, a temperature of 300K is imposed on the system and maintained as a constant. It should be noted that, although a phonon dynamics-based finite-temperature CG [67]) has been recently developed, the constant temperature in this work is realized through a velocity rescaling scheme (the FE nodal velocities in the CG domain and the velocities of the atoms in the GBs are scaled together towards a desired temperature) for simplicity. A time step of 0.001ps is deployed for the time integration in both the atomistic and CG domains.

3. Simulation Results

3.1. Dislocation pileup-induced internal stress nearby a H-charged GB

After the introduction of dislocations, the model is equilibrated under a constant shear with its top and bottom boundaries being fixed. Due to the large geometric mismatch between the two grains, the random GB under consideration here acts as a strong barrier to dislocation motion. As a consequence, the equilibration of the current model under a shear leads to the formation of a pileup. The separation between dislocations in this pileup is at several nanometers nearby the GB but gradually increases up to tens of nanometers in the region away from the GB. Not surprisingly, a progressive pileup leads to an obvious local stress concentration (Figure 2a,b), where the front and top views of an instantaneous shear stress (${\tau }_{yz}$) distribution are presented, respectively. Although the dislocation pileup is formed only under a very low applied shear, once formed, a stress magnitudes significantly higher than the applied shear can develop at its spearhead. To characterize the dislocation pileup-induced internal stress profile, the material domain ahead of the slip-interface intersection is divided into a series of finite-sized volume elements, each of which is in a dimension of 0.5nm × 0.5nm × 0.5nm. The stresses acting on each volume element can be then calculated and plotted as a function of its distance away from the pileup tip (Figure 2c). Similar to our previous results from 2D simulations, two main observations here are: (a) CAC successfully captures the local stress concentration at a level up to 3GPa when the number of the dislocations in the pileup, noted as n, increases from 5 to 15; (b) this local stress concentration does not decay to zero at a site 80nm away from the pileup tip when n = 15, and is expected to span an even longer range when n > 15.

The stress profiles in Figure 2c can be then fit into an equation resulting from the theory by Eshelby, Frank, and Nabarro [72]:

$${\tau }_{yz}={\tau }_0+\frac{K}{\sqrt{(r+r_0)}}$$

(1)

where K describes the intensity ahead of the slip-GB intersection, r is the distance away from the pileup tip at the GB; ${\tau }_0$ is introduced to accommodate an uncertainty of stress state at the reference point where the stress state was relatively low and uniform; $r_0$ is introduced to indicate the position where the maximum stress concentration occurs. A fitting of the CAC simulation data into Equation (1) leads to a quantification of the stress intensity factor, K, and its dependence on the number of the dislocations (see the fitted values of K in Table 1).

3.2. H atom diffusion nearby the slip-GB intersection

As mentioned before, CAC has a capability of quantifying how the long-range plastic flow affects the atomic-level H diffusion nearby the GB because: (i) the CG model for the material domain far away from the GB enables it to accommodate the long-range plastic flow mediated by a large number of dislocations; and (ii) the fully atomistic resolution at the GB enables us to explicitly track the motion of each H atom, from which the GB diffusivity can be extracted. As a demonstration of such a capability, Figure 3a,b present the trajectories of three H atoms in a domain nearby the slip-GB intersection when zero and 15 dislocations are piled up in the region far away from the GB, respectively. It should be noted that, for comparison purpose, the trajectories in both Figure 3a,b are recorded within a same duration. They are all projected on an initial atomistic configuration containing no dislocations nearby the slip-GB intersection. In this way, the H diffusion trajectories with respect to the same atomistic structure are clearly visible and displayed. Three main observations from Figure 3a,b are: (i) within a same time window, the H atom behind the pileup tip in the left grain diffuse a significantly longer distance when n = 15 that it does when n = 0 (see the left trajectories in Figure 3a,b); (ii) by contrast, the H atom ahead of the pileup tip in the right grain diffuse considerably less when n = 15 than it does when n = 0 (see the right trajectories in Figure 3a,b); and (iii) the H atom within the GB (see the middle trajectories in Figure 3a,b) does not diffuse much when comparing with the H atoms in left and right grains no matter n = 15 or n = 0. This is consistent with common wisdom of considering the GB as a deep trapping sites of H.

To determine the H diffusion coefficient at the slip-GB intersection when different number of dislocations participate in th pileup, we measure the mean square displacement (MSD) of those H atoms, calculate the H diffusion coefficient, and connect it with the pileup tip stress intensity factor at n = 0, 5, 10, and 15, respectively. In details, we calculate the MSD of the H atoms ahead of the pileup tip over a finite duration and then the diffusion coefficient, D, as below:

$$D=\underset{t\to \infty }{lim}\frac{1}{6N}\sum _{i=1}^N{<|\mathbf{x}\left(t\right)-\mathbf{x}\left(0\right)|}^2>=\frac{1}{6}\underset{t\to \infty }{lim}\frac{1}{N}\sum _{i=1}^N{|{\mathbf{x}}^{\left(i\right)}\left(t\right)-{\mathbf{x}}^{\left(i\right)}\left(0\right)|}^2$$

(2)

where N is the number of the H atoms ahead of the pileup tip, $\mathbf{x}\left(t\right)$ is the position of the ith H atoms at time t with respect to a reference position of at $\mathbf{x}\left(0\right)$.

Figure 3c show the MSD curves resulting from such calculations within a duration of 100 ps. It is recalled that, when a particle (the H atoms here) is in a diffusive motion, its MSD will be linearly proportional to the diffusion coefficient. Using the MSD curves in Figure 3c, we measure the local H diffusion coefficient as the slope of each MSD curve. The extracted H diffusion coefficients are summarized as below. Clearly, the more dislocations piled up at the GB, the higher stress intensity and the lower H diffusion ahead of the pileup tip. One consequence of such a reduced diffusion of the H atoms ahead of the dislocation pileup can be a clustering of them. To confirm it, we choose several pairs of H atoms and monitor the relative distances between those two atoms in each pair. Results show that the separation between the two H atoms in each pair decreases upon the increase of the number of the dislocations in a pileup. Figure 3d,e present the typical trajectories of two H atoms in one pair nearby the GB when n = 0 and n = 15, respectively. It is seen that those two H atoms get significantly closer to each other when n = 15 than they do at n = 0.

3.3. H Diffusion Heterogeneity and Its Dependence on Local Stresses

Different from many existing full MD or continuum models/theories [26,27,28] which either only considers very short dislocation segments or formulated in a plane-strain framework by assuming a homogeneity along the dislocation line, the present CAC model is fully three-dimensional and thus can be used to characterize the heterogeneity of the stresses/diffusivity along very long dislocation lines. As a preliminary demonstration, Figure 4a presents an in-plane view of an instantaneous local shear stress distribution resulting from CAC simulations of the pileup of fifteen 95 nm-long dislocations at an atomically structured GB. The H diffusion and stress heterogeneity along the dislocation line direction are then characterized through: (i) dividing the pileup tip material domain along the line direction into 9 bins, each of which is in a width of 10 nm; (ii) calculating the MSD of the H atoms in each bin as shown in Figure 4b and extracting the the slope of each of these MSD curves the local diffusion coefficient; (iii) fitting the stress profiles ahead of each of these 9 bins into Equation (1); and then (iv) correlating the obtained local diffusion coefficients, the local stress intensity factors, and the positions of those 9 bins with each other. The data resulting from the above four steps shows: (a) both the local diffusion coefficients and stress intensity factors can fluctuate more than 20% around the averaged value of them, although such a fluctuation might be suppressed/amplified if a larger/smaller bin width is deployed in the above step-1; (b) for all the 9 positions under consideration here, there exists a good correspondence between local stress intensity and local H diffusion coefficient, i.e., the higher stress intensity, the lower diffusion coefficient. This finding motivates us to develop a mechanics model as below for describing such a local stress-dependent diffusivity ahead of the pileup tip nearby the GB.

Similar to the Oriani theory [29], considering the material domain ahead of the pileup tip as H-trapping sites, the H diffusion coefficient within them can be expressed as

$$D_{\mathrm{pileup}-\mathrm{tip}}=D_{0}f(\sigma ,T),$$

(3)

where $D_0$ is the H diffusion coefficient at NILS under zero stress, $\sigma$ is the local stress, and T is the temperature. The ratio function of $f(\sigma ,T)$ is used here to include the effects of $\sigma$ on $D_{\mathrm{pileup}-\mathrm{tip}}$ by considering the dislocation pileup as an Eshelby inclusion. It can be written as

$$f(\sigma ,T)=\frac{1}{1-K_d+K_d\exp \left(E_{\mathrm{Binding}}\right)/RT}$$

(4)

where R is Boltzmann constant, T is the temperature, $K_d$ is the pileup tip trap density (the ratio between the number of the pileup tip trap sites and the total number of the trapping sites). $E_{\mathrm{Binding}}$ is the H binding energy at the pileup tip. The H-binding energy at such trapping sites is assumed to follow the form of a Eshelby transforming energy dissipation [73,74,75]. If only the first-order interaction between H and trapping site is considered, $E_{\mathrm{Binding}}$ can be expressed as ${\sigma }_{kk}{V}_{incl}$, where ${\sigma }_{kk}$ is the hydrostatic stresses. Here, to include the effect of the shear stress into the consideration, $E_{Binding}$ is expanded as

$$E_{Binding}=({\sigma }_{ij}^{\infty }+\frac{1}{2}{\sigma }_{ij}^T){\epsilon }_{ij}^{*}{V}_{\mathrm{incl}}$$

(5)

with ${\epsilon }_{ij}^{*}$ being the local strain associated with the finite-sized volume $V_{incl}$ in a material domain with a dimension of $5_{nm}\times 5_{nm}\times 5_{nm}$ ahead of the pileup tip, ${\sigma }_{ij}$ and ${\sigma }_{ij}^{\infty }$ are the local stresses and the stress (if there is any) in the far field, respectively. An insertion of Equation (5) into Equation (4), and then Equation (4) into Equation (3) leads to an expression of the pileup tip H diffusion coefficient as below

$$D_{\mathrm{pileup}-\mathrm{tip}}=\frac{D_0}{1-K_d+K_d\left[({\sigma }_{yz}^{\infty }+1/2*{\sigma }_{yz}){\epsilon }_{yz}^{*}{V}_{\mathrm{incl}}\right]/RT}$$

(6)

with ${\sigma }_{yz}$ being the shear components of the stresses ahead of the pileup tip and all other stress components being negligible in the current computational set-up.

Since $D_{\mathrm{pileup}-\mathrm{tip}}$, ${\sigma }_{yz}$, and ${\epsilon }_{yz}$ can be all directly measured from CAC simulations, a fitting of the CAC simulation-predicted $D_{\mathrm{pileup}-\mathrm{tip}}$ at chosen sites (such as Position-5 in Figure 4a) ahead of the tip of a pileup containing different number of dislocations into Equation (6) will enable us to determine all the model parameters in Equation (6). An incorporation of such a CAC simulation data-informed mechanics model as shown in Equation (6) into the higher length scale models, such as the continuum level CPFE model, will be tested and reported in a separate work.

Figure 5. CAC simulation-predicted local H diffusion coefficient and its dependence on the local stresses ahead of the tip of a pileup containing different number of dislocations. Here, only the correlation between the stress intensity and the H-diffusion at Position-5 as indicated in Figure 4a has been considered as an example.

Figure 5. CAC simulation-predicted local H diffusion coefficient and its dependence on the local stresses ahead of the tip of a pileup containing different number of dislocations. Here, only the correlation between the stress intensity and the H-diffusion at Position-5 as indicated in Figure 4a has been considered as an example.

4. Summary and Discussion

To summarize, in this work, a multiscale computational analysis of the effects of a long-range dislocation pileup on the H diffusion nearby an atomically structured GB in a plastically deformed bi-crystalline bcc iron sample is performed. With the interatomic potential being the only constitutive relation, at a fraction of the cost of full MD, the present simulations go beyond not only nanoscale MD model by scaling up in length by accommodating tens of long dislocations using a novel CG description away from the GB , but also continuum models by retaining a fully atomistic resolution near the GB. Our several main findings are: (i) the local shear stress concentration induced by a pileup containing a large number of screw dislocations can be as high as 3GPa although it might be slightly relaxed by the atomic-level H diffusion nearby; (ii) the H diffusivity behind the pileup, at the pileup tip, and ahead of the pileup tip should be differentiated. It diffuses faster behind the pileup tip due to a pipe diffusion along the dislocations arriving at the GB, gets trapped at the pileup tip due to the presence of deep trapping sites within the GB, and diffuses slower ahead of the pileup tip due to a notable stress concentration there. The slow H diffusion ahead of the pileup tip implies a plasticity induced clustering of H atoms (PICH). This may be one important intermediate step during the process of HIC; (iii) both the local stresses and H diffusivity ahead of the pileup tip exhibit a strong heterogeneity along the dislocation line direction. And most importantly, the local stress-dependent GB diffusivity ahead of a pileup tip can fit into a model with its parameters being well calibrated from CAC simulations. This is one main outcome of this research because such a local stress-based model for the long-range pileup-affected GB diffusivity is usually considered in neither MD simulations at the atomic scale nor the continuum-level simulations at the macroscale.

In spite of the above interesting findings, the present CAC simulation results should be also taken with great caution due to its limitations in the following several aspects: (1) although our CG model was recently shown to be able to accommodate not only dislocations but also H diffusion along the dislocation [76], neither the H atoms carried by the dislocations away from the GB nor their segregation from the far region into the GB region is considered here. The reason is that such a H segregation or deposition process is significantly slower than the dislocation pileup. A capture of the dislocation pileup formation together with the H migration from the far region into the GB region necessitates a combination of CAC with the kinetic Monte Carlo (kMC), which is still under development and will be reported in the future; (2) the constant temperature of 300K is simply realized through scaling the velocities of the atoms near the GB together with the FE nodal velocities in the CG domain away from the GB, which may induce a mismatch of the phonon dynamics between the atomistic and the CG domain. The effect of such a phonon dynamics mismatch on the H diffusivity at the GB is unclear. It will be carefully quantified in our future work through implementing either the lattice dynamics-based FE shape function [77] or phonon density of states-based finite temperature algorithm in the CG domain [67]; (3) the observed clustering might result from the artifacts of the EAM potential deployed here. A confirmation of it requires the implementation of high-fidelity potentials, such as the machine learning (ML)-based potential [78], for describing the Fe-Fe, Fe-H, and H-H interactions at an accuracy comparable to that in quantum mechanical calculations. An implementation of the ML-based potential into CAC has been recently implemented in [79], and will be expanded for the problem under consideration in this work.

A further expansion of CAC along those three aspects and a connection of it with the previously reviewed two-way coupling continuum models will be attempted in the future. Such efforts may lead to a predictive multiscale computational platform that can be used to address the full complexity associated with the interaction between plasticity, diffusion, and failure in many engineering materials, such as oil pipeline steel, H-storage materials, and so on.