**1. Introduction**

Grain boundaries (GBs) have significant influence on the mechanical properties of polycrystalline materials by interacting with dislocations and other defects [1–3]. Therefore, the interactions between GBs and dislocations have attracted much attention for many decades [4–9]. Due to the complexity, the detailed interaction between GB and dislocation has not been understood thoroughly to date [3,10,11]. However, it is well accepted that the interactions between dislocations and GBs have a close relationship with the GB misorientation angle [2,3]. Dislocations trend to pile up at high-angle GBs (HAGBs) and pass through low-angle GBs (LAGBs) in popular opinion [12–14]. Therefore, the angle of GBs is employed to estimate the effects of GBs on deformation behavior in many studies.

In recent years, bicrystalline micro-pillar experiments have been widely used to investigate the effects of GBs on the deformation behavior of polycrystals [15–17]. Kheradmand et al. [18,19] have observed slip transfer LAGBs in a bicrystalline micro-pillar, which is consistent with the previous understanding. However, there are some phenomena that contradict the common view of dislocation transmission in the bicrystalline micro-pillar with HAGBs. For example, although the work of [19] has shown no slip transfer in the HAGB in a bicrystalline pillar, as anticipated, the works of [10,11,20] have successfully revealed slip transfer HAGBs in some bicrystalline pillars. Although all the GBs are HAGBs, the orientations of component crystals and the directions of load are randomly distributed, which may induce the different activation behavior of slip systems and further influence the slip transfer [17,21]. This may be the reason for the contradictory effects of HAGBs on slip transfer. It implies that the consideration of GB misorientation angle alone would not lead to fully understanding the slip transfer mode of bicrystalline specimens.

Due to the complicated development of stress and strain formed in the deformation process of bicrystals, simulations have been widely used to investigate the stress and strain evolution [22]. Raabe and Lu et al. [14,23] have employed phenomenological crystal plasticity models to qualitatively analyze the effects of GB misorientation angle on the local plastic deformation of bicrystals and found the blocking effect of GBs on plastic

**Citation:** Zhou, H.; Wang, P.; Lu, S. Investigation on the Effects of Grain Boundary on Deformation Behavior of Bicrystalline Pillar by Crystal Plasticity Finite Element Method. *Crystals* **2021**, *11*, 923. https://doi.org/10.3390/ cryst11080923

Academic Editor: Wojciech Polkowski

Received: 11 July 2021 Accepted: 5 August 2021 Published: 9 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

strain. However, the interactions between dislocations and GBs were not considered in the models. Therefore, those investigations cannot reveal the microscopic mechanism behind the effect of GBs on the plastic deformation behavior [22]. Generally, plastic deformation is the result of dislocation motion [9]. Therefore, it is necessary to establish a model that based on the evolution of dislocations to reveal the underlying mechanism. Zikry et al. [24,25] have established a dislocation density-based crystal plasticity finite element method, and coupled it with dislocation density–grain boundary interaction (DDGBI) schemes [5,26,27]. Using this model, Zikry et al. [24,28] have revealed that the different plastic slips between the component crystals influence the pile-up of dislocations at GBs in bicrystalline specimens [24]. This implies that to understand the effects of GBs on the deformation behavior of specimens, it is necessary to gain an insight into the interactions between the GB and activated slip systems. The activity of slip systems is directly influenced by crystal orientations. However, the effects of crystal orientations on the interactions between slip and GBs in bicrystalline specimens have been overlooked. Therefore, the effects of GBs on the deformation behavior of bicrystalline specimens with HAGBs are still unclear. In fact, this unclear point is also found in polycrystal material.

For a detailed understanding of the effects of GBs on the deformation of bicrystalline specimens, a DDGBI scheme coupled with the dislocation density-based crystalline plasticity finite element method has been established and used to investigate the deformation behavior of three high-purity nickel (99.99%) bicrystalline pillars with different crystal orientations but an identical GB misorientation angle (38.7◦). The global coordinate systems and schematic of the bicrystalline pillar are shown in Figure 1a. The detailed orientations of the component crystals of each bicrystalline pillar are listed in Table 1. The loading direction is along [001] of global coordinate systems. The diameter and height of samples are 1.4 μm and 2.1 μm, respectively.

**Figure 1.** (**a**) Schematic of bicrystalline pillar, (**b**) engineering stress–strain curves of specimens that are used to determine the material parameters. The contours of the accumulative plastic slip on the surface of specimen BC1, (**c**) on the front side and (**d**) on the reverse side.


**Table 1.** Crystal orientations of the three bicrystalline specimens.

#### **2. Crystal Plasticity Simulation Model**

The constitutive model used in this study is based on the theories proposed by Asaro [29] and Zikry [30]. The equations of the shear strain rate, reference stress of slip systems and dislocation density have been introduced in detail in [24,31–33]. Here, only the DDGBI scheme is described briefly. To obtain a detailed understanding of DDGBI, a transmission scalar for each slip system, *λαβ AB*, was proposed [9] as

$$
\lambda\_{AB}^{\;a\;\beta} = \left( n\_A^a \cdot n\_B^\beta \right) \left( \mathbf{s}\_A^a \cdot \mathbf{s}\_B^\beta \right) \tag{1}
$$

where *s* is the current slip vector and *n* is the normal vector. The superscript *α* and *β* correspond to the slip system *α* and *β* in grain *A* and *B*, respectively. Grain A and B are separated by GBs. Due to the orthogonality of the unit vector, *λαβ AB* varies from zero to one, where complete dislocation blockage and complete dislocation transmission correspond to a value of zero and one, respectively. Therefore, *λαβ AB* can be employed as an index to indicate whether the slip on the *α* slip system in grain *A* transmits across the GB to the *β* slip system in grain *B* or not [24].

Additionally, slip transmission only occurs between the activated slip systems [21]. Therefore, the second criterion of slip transmission between slip systems in the adjacent grains can be expressed as

$$\left| \mathbf{r}\_A^{(a)} \right| \ge \left| \mathbf{r}\_c^{(a)} \right| \left| \mathbf{r}\_B^{(\beta)} \right| \ge \left| \mathbf{r}\_c^{(\beta)} \right| \tag{2}$$

where *τ* is the resolved shear stress, and the subscript *c* corresponds to the critical stress to drive the motion of dislocations. Moreover, definition of the slip systems used in the crystal plasticity model and the Schmid factors of the slip system for the three specimens are listed in Table 2.


**Table 2.** Definition of the slip systems used in the crystal plasticity model and the Schmid factors of the slip system for the three specimens.

Figure 1a shows the finite element mesh of the pillar for uniaxial compress deformation and the boundary condition in commercial ABAQUS software. A 3D solid element, C3D8, with eight nodes and one integration point, was employed. The mesh density in the regions adjacent to GBs was slightly higher than other regions. In order to verify the mesh convergence, mesh with different sizes in the regions adjacent to GBs was considered. The mesh size of the model was about 50 nm. The displacement boundary conditions were applied to specimens as shown in Figure 1a. For uniaxial compression, a prescribed displacement was imposed on the top surface of specimens, and the fixed condition (ENCASTRE) was applied on the bottom surface of specimens. Additionally, the strain rate was 0.001 s−<sup>1</sup> during the loading.

To obtain the material parameters to be used in simulation, deformation behavior of the single crystalline pillars, crystal *A*1 and *B*1, was simulated first. Based on the relationships of the dislocation density evolution coefficients [24] and engineering stress– strain curves of the specimens [19], the material parameters were obtained. Thereafter, the accuracy of material parameters and DDGBI model were verified by simulating the deformation of specimen BC1 first. The engineering stress–strain curves of these samples are shown in Figure 1b. It can be observed that the simulations are in good agreement with the experiments. Figure 1c shows the accumulative plastic slip on the surface of specimen BC1 from the front side of the specimen. It is clear that the plastically deformed area in grain *B*1 is larger than that in grain *A*1. Moreover, the area close to the GB in grain *B*1 has higher plastic slip than other areas. The similar phenomenon is also observed from the reverse side (Figure 1d). Those phenomena are consistent with the distribution of slip bands in experiments reported in [19]. Those results suggest that the DDGBI scheme and the material coefficients can be used to study the heterogeneous deformation of bicrystals. All the parameters that were used in the simulation are listed in Table 3.


**Table 3.** Material parameters used in the crystal plasticity finite element simulations.

#### **3. Results**

Figure 2 shows the contours of elastic stress on the longitudinal section of the three specimens after 0.05% engineering strain. For specimen BC1, elastic stress on crystal *B*1 is higher than that on crystal *A*1 (Figure 2a), while the difference in the elastic stress between crystal *A*2 and *B*2 is slightly smaller in specimen BC2 (Figure 2b). For specimen BC3, the contours of elastic stress on crystal *A*3 and *B*3 are almost same (Figure 2c). Those phenomena reveal that crystal orientations have great influence on the distribution of elastic stress in bicrystalline pillars.

**Figure 2.** Contours of stress of the three specimens when the engineering strain is 0.05%. (**a**) BC1, (**b**) BC2 and (**c**) BC3.

In the elastic regime, the stresses distributed on the component crystals are influenced by their elastic moduli in the loading direction, which can be given [34] by Equation (3):

$$
\sigma\_A = \frac{\sigma\_T A\_T}{A\_A + \frac{E\_B A\_B}{E\_A}}, \; \sigma\_B = \frac{\sigma\_T A\_T}{A\_B + \frac{E\_A A\_A}{E\_B}} \tag{3}
$$

where *σ<sup>T</sup>* is the total stress applied to the bicrystalline pillar, *AT*, *AA* and *AB* are the crosssection areas of the bicrystalline pillar and the component crystal *A* and *B*, respectively. *EA* and *EB* are the elastic moduli for crystal *A* and *B* in the loading direction, respectively. In this simulation, *AA* = *AB*. Therefore, Equation (3) can be expressed as

$$
\sigma\_A = \frac{2E\_A \sigma\_T}{E\_B + E\_A \prime}, \sigma\_B = \frac{2E\_B \sigma\_T}{E\_B + E\_A} \tag{4}
$$

The orientation dependence of the elastic moduli E in a uniaxial stress state <sup>→</sup> *χ* direction is given [34] by

$$\frac{1}{E\_{\frac{\cdot}{x}}} = \frac{\mathbb{C}\_{11} + \mathbb{C}\_{12}}{(\mathbb{C}\_{11} - \mathbb{C}\_{12})(\mathbb{C}\_{11} + 2\mathbb{C}\_{12})} + \left(\frac{1}{\mathbb{C}\_{44}} - 2\frac{1}{\mathbb{C}\_{11} - \mathbb{C}\_{12}}\right) \frac{\left(\mathbf{x}\_{1}^{2}\mathbf{x}\_{2}^{2} + \mathbf{x}\_{2}^{2}\mathbf{x}\_{3}^{2} + \mathbf{x}\_{3}^{2}\mathbf{x}\_{1}^{2}\right)}{\left\|\chi\right\|^{4}} \tag{5}$$

Based on Equations (4) and (5), the elastic moduli of the component crystals in a uniaxial loading direction, and the ratios of the stresses on the component crystals to the total stresses of specimens, were calculated and are listed in Table 4. The difference in the distributed stresses between crystal *A* and *B* in specimen BC1, BC2 and BC3 is 0.346*σT*, 0.272*σ<sup>T</sup>* and 0.022*σT*, respectively. Therefore, the difference in the stress distribution between the component crystals in these three specimens decreases gradually (Figure 2). The above analysis illustrates that crystal orientation determines the elastic moduli and further influences stress proportions of the components in the elastic deformation stage.

**Table 4.** Elastic moduli and the ratio of the applied stresses on the component crystals to the total stress.


The different stress proportions of the components in the elastic deformation stage will induce different levels of initiation and development of plastic deformation. To reveal the initiation of plastic deformation, the evolutions of plastic slip at two special nodes close to the GB were analyzed. The locations of the nodes on the three specimens are shown in Figure 3a–c. For specimen BC1, node L1 starts plastic deformation when the engineering strain is 0.14% (Figure 3d). For node R1, plastic deformation begins when the engineering strain is 0.167% (Figure 3e), which is later than node L1. This is because the largest Schmid factor of crystal *A*1 is 0.486 on γ5, and the largest Schmid factor of grain *B*1 is 0.313 on γ<sup>8</sup> and γ11, as listed in Table 2. Combining the different stress proportions with the Schmid factors of the component crystals, the resolved shear stress on γ<sup>5</sup> of crystal *A*1 is higher than those on γ<sup>8</sup> and γ<sup>11</sup> of crystal *B*1. After the initiation of plastic deformation, plastic slips interact with the GB and influence the evolution of the plastic slip. For specimen BC1, the transmission scalars *λ*5,8 *<sup>A</sup>*1,*B*<sup>1</sup> and *<sup>λ</sup>*5,11 *<sup>A</sup>*1,*B*<sup>1</sup> are 0.126, and the other transmission scalars are equal to zero. This means the plastic slip on γ<sup>5</sup> of crystal *A*1 almost cannot transmit to γ<sup>8</sup> and γ<sup>11</sup> of crystal *B*1. Therefore, the distribution of plastic slip on γ<sup>5</sup> of crystal *A*1 is very different from the plastic slip on γ<sup>8</sup> and γ<sup>11</sup> of crystal *B*1 on the two sides of the GB, i.e., the plastic deformation is not continuous in the vicinity of the GB. This is confirmed by the contours of the plastic slip on the GB of specimen BC1, as shown in Figure 4a–c. It implies

that the GB is an obstacle to the development of plastic slip when the transmission scalars of the activated slip systems are low. Therefore, the accumulative plastic slip on specimen BC1 is heterogeneous and there is a sharp change in the distribution of plastic slip between *B*1 and *A*1 (Figure 3a).

**Figure 3.** (**a**–**c**) Field of accumulative plastic slip on the longitudinal section of BC1, BC2 and BC3, respectively, after 0.4% strain. L1 (**d**), R1 (**e**), L2 (**f**), R2 (**g**), L3 (**h**) and R3 (**i**).

Similarly, the combined effects of Schmid factors and the stress proportion lead to the plastic deformation at node R2 beginning earlier than that at L2, as shown in Figure 3f,g. Additionally, the plastic deformation at R2 is mainly contributed to by the plastic slip on γ<sup>12</sup> and γ10, while it is γ<sup>12</sup> and γ<sup>5</sup> at L2. For specimen BC2, the transmission scalars *λ*12,10 *A*2,*B*2 and *λ*12,12 *<sup>A</sup>*2,*B*<sup>2</sup> are 0.887 and 0.826, respectively. Due to the higher transmission scalar, the plastic slip on γ<sup>12</sup> of crystal *A*2 can transmit to γ<sup>10</sup> and γ<sup>12</sup> of crystal *B*2, as marked by arrows in Figure 4d–f, which promotes continuous plastic deformation in the vicinity of the GB. Therefore, the accumulative plastic deformation on the component crystals is nearly homogeneous (Figure 3b).

**Figure 4.** Contours of plastic slip on grain boundaries of the three specimens when the engineering strain is 0.18%. (**a**–**c**) correspond to the plastic slip on γ<sup>5</sup> of crystal *A*1, γ<sup>8</sup> and γ<sup>11</sup> of crystal *B*1, respectively. (**d**–**f**) correspond to the plastic slip on γ<sup>12</sup> of crystal *A*2, γ<sup>10</sup> and γ<sup>12</sup> of crystal *B*2, respectively. (**g**,**h**) correspond to the plastic slip on γ<sup>11</sup> of crystal *A*3 and γ<sup>12</sup> of crystal *B*3, respectively.

For specimen BC3, the plastic deformation begins almost simultaneously on the two component crystals, as shown in Figure 3h,i. The plastic deformation of crystal *A*3 and *B*3 is mainly contributed to by the plastic slip on γ<sup>11</sup> and on γ12, respectively. As the transmission scalar *λ*11,12 *<sup>A</sup>*3,*B*<sup>3</sup> is 0.891, the plastic slip on γ<sup>11</sup> of crystal *A*3 and on γ<sup>12</sup> of crystal *B*3 can transmit the GB easily. Therefore, the accumulative plastic slips on γ<sup>11</sup> of crystal *A*3 and γ<sup>12</sup> of crystal *B*3 are nearly the same, as marked by the arrow in Figure 4g,h, which induces continuous accumulative plastic slip on the GB in the BC3 specimen (Figure 3c). The differences in the accumulative plastic slip among the three bicrystalline pillars reveal that the crystal orientations determine the activation of slip systems and further influence the interaction between plastic slip and GBs, and the angle between the activated slip systems of component crystals are the critical factor influencing the homogeneity of deformation.

Besides these above three specimens, we performed the simulations on other bicrystalline micro-pillar specimens with different degrees of HAGBs, such as 45◦, 55◦, 75◦, etc. Additionally, the same conclusion has been reached that the angle between the activated slip systems of component crystals are the critical factor influencing the homogeneity of deformation. Motivated by the discussion, the underlying mechanism of the effects of GBs on the heterogeneous plastic deformation of specimens can be understood from the evolution of the activated slip systems. For the specimens with LAGBs, the orientations of crystals are close to each other, thus the component crystals are stressed almost identically [19] in the elastic stage. Moreover, the Schmid factors of the same slip systems in the component crystals are nearly the same. The combined effects of the Schmid factors and applied stress induce same activated slip systems of the component crystals. Meanwhile, the angles between the slip systems of the component crystals are low. This ensures that the slip on the activated slip systems of the component crystals can transmit the GBs easily [35]. In contrast, for the specimens with HAGBs, the stress on the component crystals is always unequal and the Schmid factors of the same slip systems of the component crystals are different. The combined effects of applied stresses and Schmid factors may induce different activated slip systems in the two component crystals, as shown in Figure 3. The angles between these activated slip systems of the two components are varied, thus the corresponding transmission scalars are random, just as for the transmission scalars of specimen BC1–BC3. That is why the plastic deformation of the specimens with the same HAGBs are different (Figure 3a–c). This analysis suggests that besides the angle of GBs, the angle between the activated slip systems should also be considered to understand the plastic deformation of specimens.

#### **4. Conclusions**

In summary, the effects of GBs on the heterogeneous deformation of bicrystalline pillars have been investigated by the crystalline plasticity finite element method coupled with a dislocation density grain boundary interaction scheme. It is found that the orientations of the crystals determine the elastic moduli and stress proportions of the component crystals at the elastic deformation stage. The combined effects of stress proportions and Schmid factors further influence the initiation of plastic deformation. Additionally, then, the angle between the activated slip systems of the component crystals determines the interaction of slip with GBs, which finally influences the heterogeneous deformation of specimens.

**Author Contributions:** Conceptualization, H.Z., P.W. and S.L.; methodology, H.Z. and P.W.; software, P.W. and H.Z.; validation, H.Z., P.W. and S.L.; formal analysis, H.Z., P.W. and S.L.; investigation, P.W. and H.Z; resources, H.Z., P.W. and S.L.; data curation, H.Z., P.W. and S.L.; writing—original draft preparation, H.Z. and P.W.; writing—review and editing, H.Z. and P.W.; visualization, H.Z., P.W. and S.L.; supervision, P.W. and S.L.; project administration, P.W. and S.L.; funding acquisition, S.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China, Grant No. 51474203 and the Youth Innovation Promotion Association, Chinese Academy of Sciences, Grant No. 2013126.

**Acknowledgments:** The authors thank the financial support received from the National Natural Science Foundation of China (grant No. 51474203) and the Youth Innovation Promotion Association, Chinese Academy of Sciences (grant No. 2013126). The authors are grateful to Sheng Gang for helping with the subroutine.

**Conflicts of Interest:** The authors declare no conflict of interest.
