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
. The presence of
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
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
,
, and
in the left grain, and
,
, and
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
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.
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.
Author Contributions
Conceptualization, LX; methodology, YP, RJ, TP, and LX; software, YP, RJ, TP, XC, and SX; validation, AB and LX; formal analysis, NZ, SX, AB, and LX; investigation, YP, RJ, and TP; resources, AB and LX; data curation, YP, RJ, and LX; writing—original draft preparation, YP and LX; writing—review and editing, YP, RJ, TP, XC, NZ, SX, AB, and LX; visualization, YP, RJ and LX; supervision, LX; project administration, AB and LX; funding acquisition, NZ, AB and LX. All authors have read and agreed to the published version of the manuscript.