**1. Introduction**

Crystal plasticity is a class of constitutive equations that can be used to describe deformation mechanisms of polycrystalline materials at the grain scale [1–3]. It has been used in many studies to understand the evolution of heterogeneous deformation fields within individual grains in "real" time scales [4–7]. A significant number of crystal plasticity studies have focused on modeling deformation via slip, which is controlled by the movement of dislocations on a particular plane known as the slip plane and in a particular direction known as the slip direction. The strength of the obstacles interacting with dislocations determines the resistance of a slip system to the movement of the dislocations and the consequent material hardening. These obstacles can be categorized into short-range and long-range barriers [8]. The primary short-range obstacles are generally assumed to be the other dislocations that intersect the slip plane. The long-range obstacles may include the elastic stress field induced by far field dislocations or grain boundaries. It is generally assumed that the long-range obstacles are affected by the density of the geometrically necessary dislocations (GNDs), while the evolution of both statistically stored dislocations (SSDs) and GNDs during crystallographic slip equally increase the short-range interactions [9,10]. When conventional crystal plasticity models are used, the difference between these two obstacles is generally ignored and the critical resolved shear stress (CRSS) for the movement of the dislocations does not depend on the deformation state in the neighboring points. However, the plastic response of a polycrystal depends on

**Citation:** Sedaghat, O.; Abdolvand, H. Strain-Gradient Crystal Plasticity Finite Element Modeling of Slip Band Formation in α-Zirconium. *Crystals* **2021**, *11*, 1382. https://doi.org/ 10.3390/cryst11111382

Academic Editor: Wojciech Polkowski

Received: 28 October 2021 Accepted: 10 November 2021 Published: 12 November 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/).

both the local state of deformation and the gradients of the plastic strain. Incorporation of the plastic strain gradient in crystal plasticity constitutive equations makes the model response dependent on the neighboring elements—this is generally investigated in the non-local crystal plasticity approach [11–15]. The strain-gradient plasticity theories are formulated in two different forms, i.e., the lower-order and higher-order approach [16]. In the lower-order approach [17–19], the plastic strain gradient terms are only incorporated in the hardening laws. However, in the higher-order approach, the equilibrium equations are also adjusted [20–23]. This paper focuses on examining the methods available in the literature for calculating GND densities in the lower-order crystal plasticity framework by providing a direct comparison between the results of each method to those measured experimentally.

GND densities are inherently linked to the lattice curvature and can be determined using the gradients of plastic strain [24]. The derivation of the dislocation densities from the "curl" of the deformation gradient in the Nye equation has been accompanied with inconsistencies in the literature, some of which have been reported by Das et al. [25]. Generally, two methods are used in the literature for extracting GNDs using Nye equations. In the first method, the contribution of each slip system, e.g., α, to GND density is assumed to be proportional to the plastic shear accommodated on the same slip system α. Hence, the number of linear equations that should be solved is equal to the number of unknowns, i.e., GNDs on each slip system α. This is called the direct method, where the density of GNDs on each slip system is determined unambiguously. The direct method was initially proposed by Dai [26] and was subsequently implemented in crystal plasticity finite element (CPFE) models [10,18,27–29]. In the second method, the cumulative contribution of plastic shears on all slip systems to the total GND densities on all slip systems is assumed to be proportional to the Nye tensor. This usually results in an under-determined system of equations where the number of unknowns is more than the number of equations. Two main approaches are suggested for determining the unknowns. In the first approach, known as L1, a solution is found by minimizing the dislocation line energy. In the second approach, known as L2, a solution is found by minimizing the sum of the squares of the dislocation densities [30]. Das et al. showed that the use of the L2 method leads to a solution where dislocation densities are evenly distributed on all slip systems, while the use of the L1 method leads to capturing some variations of GNDs between the slip systems [25]. Negligible difference was found between the magnitudes of the total dislocation densities obtained from the two minimization approaches. The L<sup>2</sup> minimization approach has been widely used in the CPFE models [31–34]. In contrast to the direct method, which provides a unique solution for the GND density on each slip system, both minimization approaches lead to non-unique solutions. Therefore, an investigation is required to compare the magnitudes and distributions of the obtained GND densities from these two methods.

The measurement of internal stress and dislocation densities helps validate dislocationbased crystal plasticity models. Diffraction-based experimental methods are mainly used for measuring internal lattice strains [35–38] and localized stress fields [38,39] as well as the density of dislocations [40]. For example, the dislocation densities in the vicinity of slip bands in HCP polycrystals were measured using micro X-ray diffraction [41]. Although a three-dimensional view of the dislocation densities is provided, only a few grains can be studied in this technique. A high angular resolution electron back scatter diffraction (HR-EBSD) technique can be used to measure "relative" elastic strains and rotations and "absolute" GND densities for many grains, yet close to the sample surface. Diffraction patterns are collected and then cross correlated using a reference pattern. The deformation gradient tensor is then calculated and used to determine the relative elastic strain tensor [42–45] and GND densities using the Nye tensor [46]. Many studies have used HR-EBSD to compare the numerically obtained GND densities from CPFE models to the measured ones [33,47–49].

The objectives of this paper are to compare the GND densities calculated using the direct method to those obtained from the minimization-based method and investigate the capability of the GND-based CPFE models in capturing the formation of slip bands in an α-zirconium polycrystal. Hence, the direct and the minimization-based methods for calculating GND densities are formulated and implemented in a lower-order strain-gradient CPFE model. After calibrating single crystal parameters, the evolution of GND and SSD densities predicted using each method is studied in detail using single crystal FE models. The capability of the developed non-local models in capturing the experimentally observed slip bands of a deformed polycrystalline α-zirconium is subsequently investigated. This is followed by comparing the results of the models for GND densities, stresses, and elastic rotations to those measured via HR-EBSD.

#### **2. Sample Preparation and Experimental Set Up**

Blocks of α-zirconium material were firstly annealed at 700 ◦C and then air cooled to relieve residual stresses as much as possible. A dog bone sample was subsequently made with a gauge length of 20 mm. The sample was mechanically polished down to 4000 grit, followed by polishing in a 50 nm colloidal silica suspension and electropolishing in a −30 ◦C solution of 90% methanol and 10% perchloric acid for 60 s at 25 V. The cross section of the sample after polishing was measured to be 0.5 × 0.5 mm2.

The sample was uniaxially deformed at room temperature under strain control to a tensile strain of 1.2% using a strain rate of 2.6 × 10−<sup>5</sup> s<sup>−</sup>1. The EBSD measurement was performed in a Zeiss MERLIN field emission gun scanning electron microscope (FEG-SEM) with 20 keV beam energy, 15 nA probe current, and at the working distance of 18 mm. Diffraction patterns were collected every 0.5 μm using a Bruker high resolution EBSD detector. Following the method developed by Wilkinson et al. [42], elastic strain and lattice rotations were calculated by cross correlating the collected 800 × 800 pixels Kikuchi patterns [50]. An EBSD map of the deformed sample with the coordinate system used is shown in Figure 1a. As shown, the X-axis coincides with the loading direction, the Z-axis is along the electron beam direction and sample thickness, and the Y-axis is defined by the cross product of the other two axes. In addition, as shown in Figure 1b, the c-axis of the HCP crystals is oriented towards the z-axis with a little spread towards the x-axis. Since the specimen was only deformed to 1.2% and the c-axis of most of the crystals is perpendicular to the loading direction, deformation twinning was not active in this experiment. No twin was observed in the EBSD and the ARGUS images taken from sample after deformation. As such, the effects of twinning are ignored in all mathematical formulation and CPFE modeling.

**Figure 1.** (**a**) An inverse pole figure-Z of the sample with the corresponding pole figure shown in (**b**).

#### **3. Crystal Plasticity Formulation and Input Model**

In this section the constitutive equations implemented for calculating GND densities using both the direct and minimization-based methods are described. Here, the direct and minimization-based methods are called "Method I" and "Method II", respectively. These equations are implemented in a User MATerial (UMAT) subroutine originally developed by Abdolvand et al. [51] and recently updated by Sedaghat and Abdolvand [47]. After describing the implemented constitutive equations, the steps followed to prepare the FE input model are presented.

#### *3.1. Crystal Plasticity Constitutive Equations*

The total deformation gradient (*F*) can be decomposed to the elastic (*Fe*) and plastic (*Fp*) parts:

$$F = F^\sharp F^p \tag{1}$$

The total velocity gradient tensor (*L*) in the current configuration can be divided into the elastic (*Le*) and plastic (*Lp*) parts as:

$$L = \dot{F}\mathcal{F}^{-1} = \dot{F}^{\mathcal{E}}\mathcal{F}^{\epsilon-1} + F\dot{F}^{\mathcal{P}}F^{p-1}F^{\epsilon-1} = L^{\epsilon} + L^{p} \tag{2}$$

The total velocity gradient tensor can be divided into its symmetric part, i.e., deformation rate tensors (*D<sup>e</sup>* , *Dp*), and asymmetric part, i.e., spin tensors (**Ω***<sup>e</sup>* , **Ω***p*). The plastic part of the velocity gradient tensor is calculated using the following equation [52]:

$$L^p = D^p + \Omega^p = \sum\_{a=1}^N \dot{\gamma}^a \overset{\rightarrow}{m}^a \otimes \overset{\rightarrow}{n}^a \tag{3}$$

where <sup>→</sup> *mα* , → *n α* , and . *γ<sup>α</sup>* respectively represent the slip direction, the slip plane normal, and the shear rate on the *<sup>α</sup>*th slip system for the *<sup>N</sup>* number of active slip systems. <sup>→</sup> *mα* <sup>⊗</sup> <sup>→</sup> *n α* is known as the Schmid tensor of the slip system *α*. The shear rate on the slip system *α* is calculated based on the resolved shear stress (*τα*) acting on the same slip system [52]:

$$\dot{\gamma}^{a} = \dot{\gamma}\_{0}^{a} \text{sign}\left(\frac{\tau^{a}}{\mathcal{S}^{a}}\right) \left|\frac{\tau^{a}}{\mathcal{S}^{a}}\right|^{n} \tag{4}$$

where . *γ*<sup>0</sup> is a reference shear strain rate, *n* represents the sensitivity of the material to a strain rate, and *g<sup>α</sup>* is the CRSS of the slip system *α*. The resolved shear stress acting on the slip system *α* is proportional to Kirchhoff stress (*ψ*) through the following equation [53]:

$$
\tau^a = P^a \; : \; \psi \; \tag{5}
$$

where *P<sup>α</sup>* is the symmetric part of the Schmid tensor. The Jaumann rate of the Kirchhoff stress tensor (*ψ*ˇ) is related to the elastic part of the deformation rate tensor (*De* ) and the elastic stiffness tensor (C) of the HCP crystal after rotation to the deformed configuration as:

$$\Psi = \mathbb{C} \, : \, \mathcal{D}^{\varepsilon} \, \tag{6}$$

The elastic modulus of zirconium HCP crystals used in this study is the one determined by Fisher and Renken [54]: C11 = 143.5 GPa, C33 = 164.9 GPa, C12 = 72.5 GPa, C13 = 65.4 GPa, and C44 = 32.1 GPa. The objective rate of Kirchhoff stress, in Equation (6), is defined as:

$$
\psi = \dot{\psi} - \Omega^{\varepsilon} \psi + \psi \Omega^{\varepsilon} \tag{7}
$$

The strength of each slip system (*gα*) in Equation (4) follows a dislocation-based hardening law [10–12,27]:

$$\mathbf{g}^{a} = \mathbf{g}^{a}\_{\*,0} + \frac{H^{a}}{\sqrt{D}} + \tilde{\xi}Gb^{a} \sqrt{\sum\_{\beta=1}^{N} q^{a\beta} \left(\rho^{\beta}\_{GND} + \rho^{\beta}\_{SSD}\right)} = \mathbf{g}^{a}\_{\*,0} + \frac{H^{a}}{\sqrt{D}} + \tilde{\xi}Gb^{a} \sqrt{\rho\_{T}} \tag{8}$$

where *g<sup>α</sup>* is the updated CRSS, *g<sup>α</sup>* <sup>∗</sup>,0 is the initial CRSS, *D* is the equivalent grain diameter, *H<sup>α</sup>* is the Hall-Petch parameter for slip system *α*, *ξ* is a material constant, *G* is the shear modulus, *b<sup>α</sup>* is the size of the Burgers vector for the slip system *α*, and *ρ β GND* and *ρ β SSD* are

the GND and SSD densities on the slip system *β*. A hardening matrix *qαβ* is introduced to consider the effects of self- and latent-hardening. For the sake of simplicity in deriving a solution, the term under the radical is called *ρT*. It is assumed that the GND and SSD densities have the same hardening matrix *qαβ*. The shear modulus, *G* in Equation (8), is assumed to be the average value of C44, C55, and C66 [55]. The material constant *ξ* is assumed to be 0.5.

It is assumed that the evolution of SSD densities follows the equation below [10,56]:

$$\dot{\rho}\_{SSD}^{a} = \frac{\left| \dot{\gamma}^{a} \right|}{b^{a}} \left( K^{a} \sqrt{\rho\_{SSD}^{a} + \rho\_{GND}^{a}} - 2 y\_{c}^{a} \rho\_{SSD}^{a} \right) \tag{9}$$

where *K<sup>α</sup>* is the dislocation accumulation constant and *y<sup>α</sup> <sup>c</sup>* is the dislocation annihilation length. The first term in the right-hand side of the Equation (9) is the dominant term at the early stages of plasticity, while with further loading, the effects of the second term will be non-negligible. In the following sub-sections, the two methods used for calculating GND densities are presented.

#### 3.1.1. Method I

In the first method, the GND density of the slip system *α* is uniquely correlated with the resolved shear strain accumulated on the same system [26]:

$$\dot{\rho}\_{\rm GND}^{a} = \frac{\dot{\gamma}^{a}}{b^{a}} \left| \text{Curl} \left[ \stackrel{\rightarrow a}{n\_{j}} \stackrel{p}{\mathbf{F}\_{j\bar{i}}^{p}}{\mathbf{F}\_{j\bar{i}}^{a}} \right] \right| = \frac{\dot{\gamma}^{a}}{b^{a}} \left| \epsilon\_{ik\eta} \frac{\partial (n\_{j}^{\rightarrow \mathcal{A}} \stackrel{p}{F}\_{j\eta}^{p})}{\partial \mathbf{x}\_{k}} \right| \tag{10}$$

where  *ikq* is the permutation tensor. Hence, *ρ<sup>α</sup> GND* has a unique solution. For implementing Equations (8)–(10) into the UMAT, the time derivative of Equation (8) is firstly calculated and coupled with Equations (9) and (10) to solve for the increments of shear strains using the Newton-Raphson iterative algorithm. The evolution of CRSS can be reformulated by taking the time derivative of Equation (8):

$$\dot{\mathbf{g}}^a = \mathbf{\tilde{g}} Gb^a \, \frac{\dot{\rho}\_T}{2\sqrt{\rho\_T}}, \text{ where } \dot{\rho}\_T = \sum\_{\beta=1}^N q^{a\beta} \left( \dot{\rho}\_{GND}^\beta + \dot{\rho}\_{SSD}^\beta \right) \tag{11}$$

By substituting Equations (9) and (10) into Equation (11) we have:

$$\dot{\mathbf{g}}^{a} = \sum\_{\beta=1}^{N} \left[ q^{a\beta} \frac{\mathbf{g}^{\tau} \mathbf{G} b^{a}}{2b^{\delta} \sqrt{\rho\_{T}}} \left( \left| \epsilon\_{ikq} \frac{\partial \left( \stackrel{\rightarrow a}{n\_{j}} \mathbf{F}\_{jq}^{p} \right)}{\partial \mathbf{x}\_{k}} \right| + K^{\delta} \sqrt{\rho\_{SSD}^{\theta} + \rho\_{GND}^{\theta}} - 2 \mathbf{g}\_{c}^{\theta} \, \rho\_{SSD}^{\theta} \right) \right] \dot{\gamma}^{\theta} \tag{12}$$

The complexity of implementation of this method in the UMAT subroutine is calculation of the gradient term, which was originally developed by Abdolvand [57].

## 3.1.2. Method II

In the second method, the density of GNDs on the slip system *α* is defined as [30]:

$$\left(\operatorname{curl}\left(\mathbf{F}^{p^{\mathrm{T}}}\right)\right)^{\mathrm{T}} = \sum\_{a=1}^{N} \left(\rho\_{\mathrm{GND},s}^{a} \stackrel{\rightarrow}{b}^{a} \otimes \stackrel{\rightarrow}{m}^{a} + \rho\_{\mathrm{GND},ct}^{a} \stackrel{\rightarrow}{b}^{a} \otimes \stackrel{\rightarrow}{t}^{a} + \rho\_{\mathrm{GND},cn}^{a} \stackrel{\rightarrow}{b}^{a} \otimes \stackrel{\rightarrow}{n}^{a}\right) \tag{13}$$

*ρα GND* for each slip system can be decomposed into three components, one for screw type dislocation *ρ<sup>α</sup> GND*, *<sup>s</sup>* and two edge type dislocations, *<sup>ρ</sup><sup>α</sup> GND*,*en* and *<sup>ρ</sup><sup>α</sup> GND*,*et* with *t <sup>α</sup>* = *m<sup>α</sup>* × *nα*, respectively. In contrast to Method I, Method II has infinite solutions where a minimization scheme is usually used to find a solution. Equation (13) can be solved by minimizing the sum of the squares of GND densities [19]:

$$\{\rho\_{\rm GND}^a\} = A^T (AA^T)^{-1} B \tag{14}$$

where *ρα GND* is 3N <sup>×</sup> 1 column vector including the components of GND for slip system *α*, *B* is a 9 × 1 vector, containing the curl function components in Equation (13), and *A* is a 9 × 3N matrix with the basis tensors of Equation (13). The time derivative of Equation (13) is determined as:

$$\left(\operatorname{curl}(\dot{\mathbf{F}}^{\mathcal{P}})\right)^{T} = \sum\_{a=1}^{N} (\dot{\rho}\_{\text{GND},s}^{a} \stackrel{\rightarrow}{b}^{a} \otimes \stackrel{\rightarrow}{m}^{a} + \dot{\rho}\_{\text{GND},\text{ct}}^{a} \stackrel{\rightarrow}{b}^{a} \otimes \stackrel{\rightarrow}{t}^{a} + \dot{\rho}\_{\text{GND},\text{cn}}^{a} \stackrel{\rightarrow}{b}^{a} \otimes \stackrel{\rightarrow}{n}^{a}) \tag{15}$$

By substituting Equations (15) and (9) into Equation (11), we have:

$$\dot{\mathbf{g}}^{a} = \sum\_{\beta=1}^{N} \frac{q^{a\beta} \mathfrak{J} G b^{a} (K^{\beta} \sqrt{\rho\_{\rm SSD}^{\theta} + \rho\_{\rm GND}^{\theta}} - 2 \dot{y}\_{\varepsilon}^{\theta} \, \rho\_{\rm SSD}^{\theta})}{2 b^{\beta} \sqrt{\rho\_{\rm T}}} \cdot \left| \dot{\gamma}^{\beta} \right| + \sum\_{\beta=1}^{N} \frac{q^{a\beta} \mathfrak{J} G b^{a} (\dot{\rho}\_{\rm GND}^{\theta})}{2 b^{\beta} \sqrt{\rho\_{\rm T}}} \tag{16}$$

Similar to Method I, the first term in the right-hand side of Equation (16) shows an explicit correlation between the shear strain rate . *<sup>γ</sup><sup>β</sup>* and . *g α* , but the second term represents the effects of GND evolution rate on . *g α* . Therefore, both GND density and the rate of GND density evolution should be determined in Method II. Further, a minimization approach must be followed to determine the two tensors. Accordingly, the implementation of Method II is relatively more complicated and computationally costly compared to Method I. Here, the method originally developed by Abdolvand [57] was adopted to calculate the curl of *Fp*, which was subsequently implemented in the UMAT as described in [47].

The single crystal parameters used for HCP zirconium for each method are provided in Table 1. The single crystal parameters for the Method II model were calibrated using a comprehensive data set measured in a series of in-situ neutron diffraction experiments conducted on HCP zirconium polycrystals [58]. Measured macroscopic stress-strain curves, lattice strains, and texture development were used to calibrate the model parameters as explained in [47]. The model parameters for Method I are calibrated such that they reproduce identical stress-strain curves to those from the model with Method II. The initial SSD density is set at 10<sup>10</sup> m−<sup>2</sup> for both methods. The effect of twinning is ignored in the simulations since no twin was observed in the experiment. Here, only the single crystal parameters for Method I require calibration. The process of calibration is discussed in Section 4.1.


**Table 1.** Single crystal parameters used in the hardening law for HCP zirconium from Sedaghat and Abdolvand [47].

*3.2. Input Models*

3.2.1. Single Crystal

The evolution of GND and SSD densities as well as the CRSS of each slip system is studied using the single crystal model. This is to provide a comparison between the two

methods used. A cube-shaped single crystal with a side length of 20 μm is used for the FE input model. The cube is deformed uniaxially to 5% applied strain at the strain rate of 10−<sup>5</sup> s<sup>−</sup>1. Uniaxial tensile strain is applied on the surface DCGH along the X-direction (Figure 2a). Fixed displacement boundary conditions are applied on the ABFE, AEHD, and EFGH surfaces along X, Y, and Z directions, respectively. The model is discretized with 1000 quadratic brick elements with reduced integration points (C3D20R), with 10 elements in each direction. The c-axis of the HCP crystal was rotated with respect to the loading direction from 0◦ to 90◦ at the step size of 22.5◦. As shown in Figure 2a, the misorientation between the loading direction and crystal c-axis is presented by θ.

**Figure 2.** (**a**) Applied boundary conditions for the single crystal model. (**b**) The EBSD-measured microstructure imported into the finite element solver with the meshed map shown in (**c**).

## 3.2.2. Polycrystal

The measured EBSD map shown in Figure 1a is imported into the ABAQUS FE solver. The imported polycrystal model is shown in Figure 2. The as-measured strain and strain rates are applied to the FE model. The boundary conditions used are *ux* = 0 and *uy* = 0 on the ABCD and CGHD surfaces, respectively, and *uz* = 0 on the AD and DH edges. The model is extruded in the z direction to the thickness of 50 μm to be consistent with the measured average grain size of the specimen. After conducting a mesh sensitivity study to ensure that FE results are converged, the model was discretized to 65,157 C3D20R brick elements. Three elements along the model thickness are used. The EBSD-measured grain orientations are assigned to each modeled grain. The model was deformed in two steps; in step 1, 1.2% uniaxial strain was applied to the model in the global x-direction, and in step 2, the model is unloaded and allowed to relax. In all simulations presented, the maximum time increment allowed is 5 s to ensure that the calculated GND densities are independent from the size of the selected time increment and are converged.

#### **4. Results**

A single crystal model is firstly used to calibrate the single crystal parameters for Method I of the CPFE model. This is followed by comparing the results of the CPFE models to those measured using HR-EBSD for stress, lattice rotations, and GND densities within the polycrystal. The capability of each method in determining slip activity and capturing slip bands is also discussed.

#### *4.1. Single Crystal*

The single crystal model shown in Figure 2a is used to calibrate the single crystal parameters for Method I of the CPFE model. The calculated stress of this method is compared to the one from Method II of the CPFE model. In addition, the calculated slip activities and GND densities are compared between the two methods. This is to understand the differences between each method at the single crystal level and to establish a foundation for the differences observed in the polycrystalline simulations.

In Figure 3a, the stress-strain curves calculated with the CPFE model using both Method I and II at five different crystal orientations are shown. In addition, the calculated relative activities of each slip system from both methods are provided in Figure 3b. Relative activities are the resolved shear strain (*γα*) for each slip set divided by the total shear strain at the applied strain of 5%. It is shown that identical stress-strain curves and slip activities are obtained for different crystal orientations.

**Figure 3.** (**a**) The average stress-strain curve, (**b**) relative slip activity, and (**c**) relative GND density calculated for a single crystal with different c-axis misorientation with the loading direction (θ). Comparisons are made for Method I (M I) as well as Method II (M II) of the CPFE model.

Figure 3c compares the calculated relative GND density of each slip system from both methods of the CPFE model. Relative GND densities are determined by dividing the total GND density calculated for each slip set by the total GND density calculated for all slip systems. The calculated relative GND densities from Method I follow the same trends captured for the slip activities. However, the relative GND densities of Method II follow a nearly horizontal line, which is independent from the misorientation angle θ and the calculated slip system activity. Moreover, the relative GND density of the pyramidal slip system is almost four times higher than those from prism and basal slip systems. The relative GND densities calculated for prism and basal slip systems from Method II are identical. These two curves overlapp in Figure 3c. The relative GND densities calculated using Method II for each slip set is proportional to the number of variants for each slip set. That is, GND densities are evenly distributed on all variants of a slip system when Method II is used.

The evolution of the average GND and SSD densities of each slip set from the two methods implemented in the CPFE framework is compared in Figure 4. For each slip set, the total calculated dislocation density, either GND or SSD, is divided by the number of variants in the same slip set, i.e., *<sup>M</sup>* in <sup>∑</sup>*<sup>M</sup> <sup>α</sup>*=<sup>1</sup> *<sup>ρ</sup><sup>α</sup> GND <sup>M</sup>* is 3, 3, and 12 for prism, basal, and pyramidal slip systems, respectively. Results are for three different misorientation angles (θ). In addition, the calculated total dislocation densities, GNDs plus SSDs, from the two methods are compared. Although both methods predict identical average stress-strain curves (Figure 3a), a significant difference in the predicted GND densities is observed. The SSD and GND densities from Method II have the same order of magnitude. However, in Method I, SSD densities are at least two orders of magnitude higher than GND densities. This means that the hardening of slip systems from Method I is mainly driven by SSD densities as the contribution of GND densities is almost negligible. The cumulative dislocation densities calculated from both SSDs and GNDs are almost equal to those calculated from Method II. This results in calculating the same evolution for CRSS using both methods.

**Figure 4.** The evolution of the average GND density, average SSD density, and total dislocation density of prism, basal and pyramidal slip systems at different crystal c-axis misorientations with the loading direction (θ). Results are for CPFE model using Method I and II. M is 3, 3, and 12 for prism, basal, and pyramidal slip systems.

#### *4.2. Polycrystal*

Results presented in this section are from the CPFE model using both Method I and II. The CPFE results for GND densities are presented and compared to those measured via HR-EBSD. In addition, the capability of both methods in simulating formation of slip bands is evaluated by comparing the numerical results to those from SEM images of the deformed specimen. Finally, the distribution of calculated stresses and elastic lattice rotations are compared to those measured using HR-EBSD.

#### 4.2.1. GND Density and Slip Activity

The distributions of total GND densities after the unload step for Method I and II are shown in Figure 5a,b, respectively. The total GND density from Method I is about an order of magnitude smaller than that calculated using Method II. For example, in grains G2, G3, and G5, marked in Figure 5a, a significant portion of the grains have GND densities lower than 10<sup>9</sup> m−<sup>2</sup> when Method I is used. Moreover, the maximum GND density from Method I is about 4 × 1012 m<sup>−</sup>2, but it is ~7 × 1014 m−<sup>2</sup> when Method II is used. Although there is a noticeable difference in the magnitudes of calculated GNDs, the predicted trends are quite similar between the two methods. For example, most of the grain boundaries have high GND densities. In addition, both methods predict highly localized GND concentration sites at the top-right sections of grains G1 and G2, upper side of grain G3, and right side of grain G6. These areas are shown by the white circles in Figure 5a. HR-EBSD measurement also

shows high GND density at these regions, as shown in Figure 5c. However, the magnitudes of GND densities from HR-EBSD are higher than those from both CPFE methods, but much closer to those from Method II. In comparison to Method I, the magnitude of GND densities from Method II is closer to those measured with HR-EBSD. This is not surprising because in the calculation of GND densities from HR-EBSD measurements, Method II is used along with the L1 optimization scheme. As a result, it is expected to have a better agreement between HR-EBSD measurements for GNDs and those from Method II of the CPFE model. The effects of using L<sup>1</sup> and L2 optimization schemes in determining GND densities from HR-EBSD are discussed elsewhere [59,60]. In addition, Method I cannot be used for HR-EBSD measurement as the plastic resolved shear strain on each slip system (*γα*) cannot be directly extracted with this technique (see Equation (10)).

**Figure 5.** The distribution of total GND density from (**a**) CPFE model using Method I, (**b**) CPFE model using Method II. and (**c**) HR-EBSD. (**d**) SEM-ARGUS image of the deformed sample. In (**d**) slip bands are shown with the black dashed lines and crystal orientations are shown using HCP crystals.

Figure 5d shows an ARGUS image of the deformed polycrystal, in which slip bands can be observed in most of the grains. For better visualization, the orientations of slip bands are shown with the black dotted lines. Interestingly, CPFE results from both Method I and II also show GND density localization in the form of slip bands. For example, in grains G1, G2, G3, and G4, a distinct slip band is formed that spans from one side of the grain to another, as shown with the white dotted lines in Figure 5b. The orientations of these bands are in agreement with those observed in the SEM image in Figure 5d. In addition, both CPFE methods have captured the crossing of the two separate slip bands observed in grain G7. The two slip bands can also be observed in the ARGUS image. The CPFE results reveal that these slip bands are associated with the localization of plastic shear on the prism planes. This is also consistent with the measured grain orientations as shown in Figure 5d. For example, in grains G2, G3, and G4, the observed slip bands of the ARGUS image are parallel to ones of the variants of the prism planes of the illustrated HCP crystals. Therefore, both Method I and II are equally effective when predicting the direction of slip bands.

In order to determine the relationship between the observed slip bands and active slip systems, the average resolved shear strain and GND density of each slip set are extracted from the CPFE results and are shown in Table 2. These values are calculated by taking the volume average of the resolved shear strains, or GNDs, calculated at all integration points assigned to each grain. The calculated average value for each slip set is the summation over all variants, i.e., 3 variants for basal or prism, and 12 variants for pyramidal <c + a>. Table 2 summarizes the amount of slip and GND density accumulated on the prism, basal, and pyramidal slip systems of grains G1 to G5. The predominant slip system with the highest value is presented in bold. Both methods confirm that the predominant slip system for all five selected grains is slip on the prism planes. For example, for grain G4 almost all plastic deformation is accommodated by prism slip systems whereas for grain G3 both basal and prism slip systems are active with prism being the most active one. Although these trends are also observed for the GND densities obtained from Method I, the pyramidal slip system accommodates most of the GNDs calculated in Method II. In fact, the GND density from method II is almost evenly distributed on all eighteen variants, but since the pyramidal slip system has 12 variants, it always has the highest GND density once the summation is done over all 12 variants. These results are consistent with the trends observed for the single crystal model. The heterogeneous distribution of the GND densities on the different slip variants in Method I is more reasonable compared to Method II, where the GND densities on individual slip variants are always equal, regardless of the applied loading conditions.


**Table 2.** Average resolved shear strains and GND densities from the CPFE models. Results are for grains G1 to G5 of the polycrystalline model at 1.2% applied strain. The predominant slip system is shown in bold.

> Figure 6 shows the distribution of resolved shear strains on the first and second most active slip variants. Results are for grains G2, G3, and G4. All variants are related to the slip on prism, except for grain G3 where its second dominant variant belongs to the basal system (see Table 2). Shear strain localization in the form of slip bands is observed in all three grains. These slip bands are shown with a red dotted lines in Figure 6. These bands are parallel to the ones observed in the GND maps in Figure 5a,b as well as those observed in the ARGUS image of Figure 5d.

#### 4.2.2. Stress

The normal stress component along the loading direction (*σ*11) and the in-plane shear stress component (*σ*12) of grains G3, G4, and G5 are shown in Figures 7 and 8, respectively. The distribution of *σ*<sup>11</sup> is almost identical between Method I and II. After unloading, the normal stress drops significantly. The distributions of *σ*<sup>11</sup> from HR-EBSD are also shown in the last row of Figure 7. In an HR-EBSD measurement, a reference point is selected within each grain to measure "relative" stresses and elastic lattice rotations with respect to this

point. These reference points are shown with the red dots in Figure 7. Since stresses from HR-EBSD are relative, the calculated stresses at all integration points assigned to a grain are reduced from those calculated at the integration point that coincides with the HR-EBSD reference point. This is to provide a like-to-like comparison between CPFE results and the HR-EBSD measurement.

**Figure 6.** The distribution of resolved shear strain on the most active (**first row**) and the second most active (**second row**) slip systems. Results are for the CPFE model with Method I shown in the right column and Method II in left column. Grain IDs are shown on the top of each grain.

**Figure 7.** The normal stress component along the loading direction for grain G2, G3, and G4. A comparison between Method I (**right column**) and Method II (**left column**) of the CPFE model and HR-EBSD (**last row**).

**Figure 8.** The in-plane shear stress component of grains G2, G3, and G4. A comparison between Method I (**right column**) and Method II (**left column**) of the CPFE model and HR-EBSD (**last row**).

For grain G2 of Figure 7, a tensile stress field is predicted in the lower half of the grain where the stress magnitude increases toward the grain boundary. This is in agreement with HR-EBSD measurement. In addition, two compressive stress fields are observed in the upper half of the grain. In comparison to the measurement, the size of these two fields is overestimated in CPFE results, although the trends are the same. Grain G3 can be divided into three regions, i.e., top, bottom, and middle regions. Calculated stress *σ*<sup>11</sup> is tensile in the top and bottom regions, while it is compressive in the middle region. This is in agreement with the HR-EBSD measurement. For grain G4, the left side of the grain has a compressive stress field, and the bottom side of the grain has a tensile stress field in both the model and experiment.

In Figure 8, it is shown that the distribution of shear stress *σ*<sup>12</sup> is nearly identical between Method I and II. For grain G2, a negative shear stress band is observed in the CPFE results on the top-left side of the grain. This band is also observed in the HR-EBSD measurement. In addition, a positive *σ*<sup>12</sup> region, at the top-right side of the grain, is calculated in both CPFE models, which is also observed in the HR-EBSD measurement. However, there are regions in grains G3 and G4 where the calculated shear stress *σ*<sup>12</sup> does not match with those from the HR-EBSD measurement. The magnitude and distribution of calculated normal (Figure 7) and shear stresses (Figure 8) are quite similar between the Method I and II.

#### 4.2.3. Lattice Rotations

The in-plane elastic rotation component of grains G2, G3, and G4 are shown in Figure 9 and the results of Method I and II are compared with those from HR-EBSD. Since measurements were done at the surface, only the in-plane component of the lattice rotation (Ω*<sup>e</sup>* 12) is shown. The relative values with respect to the reference point are presented in the last row of the simulation results. No difference is observed between the two methods in all the three grains after the loading and unloading steps. This means that although different GND densities are calculated in these two methods, the parameters are calibrated well, such that the resulting deformation tensors are almost identical from the two methods. The relative values with respect to the reference point of the EBSD follow the same trend as the HR-EBSD measurement. For example, both CPFE and HR-EBSD results show positive elastic rotation Ω<sup>12</sup> in the left half of grain G2. In addition, in both the model and experiment, the top side of grain G3 has negative elastic rotation Ω12. Furthermore, the CPFE results show negative Ω<sup>12</sup> on the top-right side and positive Ω<sup>12</sup> on the lower left side of grain G4, which agrees with the HR-EBSD measurement.

**Figure 9.** The elastic rotation component Ω<sup>12</sup> of grains G2, G3, and G4. A comparison between Method II (**left column**) and Method 1 (**right column**) of the CPFE model and HR-EBSD (**last row**).

#### **5. Discussion**

A lower-order strain-gradient crystal plasticity finite element model was developed in which two different methods were used to describe the evolution of GND densities. In Method I, a definite value was calculated for the magnitude of GNDs on each slip system, whereas in Method II, a minimization approach was used to solve for an under-determined system of equations. The performance of CPFE models was examined both for single crystal models and for polycrystalline models. In addition, the as-measured microstructure of a deformed α-zirconium sample was imported to the FE solver to compare the calculated grain-scale stresses, elastic lattice rotations, and GNDs to those measured using HR-EBSD.

Results shown in all previous sections were for single or polycrystalline models with an average grain size of more than 20 μm. The results of the strain-gradient CPFE model, however, are length-scale-dependent and are highly sensitive to the size of the studied grains [47]. Such a size dependency originates from the accumulated GND density, which is a function of plastic strain gradients. Since the magnitudes of GND densities from Method I and II are very different, it is important to check the effects of grain size on the results obtained. To do so, the single crystal model used in Section 3.1 is also used here, but the model is scaled so that different grain sizes can be studied. The model is scaled to ten different sizes from 0.125 μm to 64 μm and the calculated stresses, GND, and SSD densities from Method I and II are compared. The misorientation between the c-axis of the grain and the loading direction is set at 22.5◦, 45◦, and 67.5◦. The macroscopic strain applied to

the polycrystal is 5%. The Hall-Petch effects are switched off to only study the geometrical effects on the material hardening.

Figure 10a–c show that the average stress for the single crystal with a misorientation angle of 22.5◦, 45◦, 67.5◦, respectively. It is shown that the stress magnitude is highly affected by the model size when Method II is used, particularly when the grain size is smaller than 10 μm. Although Method I also shows some dependency on the grain size, the magnitude of the stress is only significantly affected by the grain size for the grains smaller than 1 μm. This is the reason why the stresses obtained from the two methods for the polycrystal, presented in Section 4.2.2 with an average grain size of 50 μm, are not that different. Figure 10d–f show that the average GND and SSD density for the single crystal with a misorientation angle of 22.5◦, 45◦, 67.5◦, respectively. It is shown that the calculated GND density has a linear relation with the grain size in the log-log scale. Interestingly, the slopes of these two lines are the same for the two methods used. In contrast to GNDs, SSDs are almost unaffected by the grain size. Therefore, the hardening of bigger grains is SSD driven, while that of smaller grains is GND driven. This means that SSDs are as important in determining stress development. Hence, when the formulation provided here is used to study fracture of polycrystals at grain scale [31,32], it is very important to fully characterize the contribution of both GNDs and SSDs to material plastic deformation.

**Figure 10.** A comparison between the results of the CPFE model, using Method I and II, for the single crystal model in Figure 2a at ten different grain sizes. Stress along the loading direction for the misorientation angles of (**a**) 22.5◦, (**b**) 45◦, and (**c**) 67.5◦. GND and SSD densities for the misorientation angles of (**d**) 22.5◦, (**e**) 45◦, and (**f**) 67.5◦. The values on the x-axis represent the size of the single crystal. Hall-Petch effects are switched off.

A lower-order strain-gradient crystal plasticity model is used here, where the material hardening law is modified by including the strain gradient terms. One of the disadvantages of this approach is the formation of unrealistic, localized deformation fields as described in [16]. Further, when lower-order formulation is used, CPFE results are mesh-size dependent; on the other hand, when higher-order formulation is used, CPFE results are almost mesh independent [61].

The simulation results obtained from CPFE models showed several discrepancies with respect to those from HR-EBSD. The observed discrepancies could be due to ignoring the sub-surface grain structures in the CPFE results. As shown in Figure 2, a columnar grain structure was assumed in the CPFE model by extruding the measured EBSD map along the z-axis. This is because the EBSD measurement provides the information about the sub-surface grains. Although ignoring the sub-surface grain structures may affect the stress values at the free surface in the CPFE results [62–64], Zhang et al. showed that the calculated trends in the GND densities and strains are not much affected by the sub-surface microstructure [48].

Table 3 summarizes the results of the comparisons made between Method I and II. The effectiveness of the two methods can be assessed using two different perspectives: implementation and performance. From the implementation point of view, Method I is more straightforward and can be used for developing higher-order non-local CPFE models. This method is also computationally efficient since no optimization is required to determine the density of GNDs. In terms of performance, the use of Method II leads to closer numbers to those measured with HR-EBSD for GND density but results in uniform distribution of GNDs on all slip variants. Finally, when larger grain sizes are used, typically more than 10 μm, there is no difference between the two methods.

**Table 3.** A comparison between Method I and II for determining GND densities.

