## Abstract

A gradient-enhanced crystal plasticity model is presented that explicitly accounts for the evolution of the densities of geometrically necessary dislocations (GNDs) on individual slip systems of deforming crystals. The GND densities are fully coupled with the crystal slip rule. Application of the model to two distinct and technologically important crystal types, namely hcp Ti and ccp Ni, is given. For the hcp crystals, slip is permitted with *a*-type slip directions on basal, prismatic and pyramidal planes and *c*+*a*-type slip directions on pyramidal planes. First, a single crystal under four-point bending is simulated as the uniform strain gradient expected in the central span provides a good validation of the code. Then, uniaxial deformation of a model near-*α* Ti polycrystal has been analysed. The resulting distributions of GND densities that develop on the various slip system types have been compared with independent experimental observations. The model predicts that GND density on the *c*+*a* systems is approximately an order of magnitude lower than that for *a*-type systems in agreement with experiment. For the ccp case, slip is considered to take place on the <110>{111} slip systems. Thermal loading of a single-crystal nickel alloy sample containing carbide particles of size approximately 30 μm has been analysed. Detailed comparisons are presented between model predictions and results of high-resolution electron backscatter diffraction (EBSD) measurements of the micro-deformations, lattice rotations, curvatures and GND densities local to the nickel–carbide interface. Qualitatively, good agreement is achieved between the coupled and decoupled model elastic strains with the EBSD measurements, but lattice rotations and GND densities are quantitatively well predicted by the coupled crystal model but are less well captured by the decoupled model. The GND coupling is found to lead to reduced lattice rotations and plastic strains in the region of highest heterogeneity close to the Ni matrix/particle interface, which is in agreement with the experimental measurements. The results presented provide objective evidence of the effectiveness of gradient-enhanced crystal plasticity finite element analysis and demonstrate that GND coupling is required in order to capture strains and lattice rotations in regions of high heterogeneity.

## 1. Introduction

It is becoming increasingly desirable to calculate strain and stress response at progressively smaller length scales. Major drivers for this are in, for example, development of understanding of deformation, fatigue and texture evolution (Rugg *et al*. 2007). Here, we confine ourselves to length scales of order grain size and hence relevant to polycrystalline materials for which there is significant interest in slip transfer, slip localization, grain boundary sliding, twinning, fatigue crack nucleation and micro-texture (Dunne *et al*. 2007*a*; Dunne & Rugg 2008; Gong & Wilkinson 2009; Bache & Dunne 2010; McDowell & Dunne 2010; Wang *et al*. 2010). Also of great interest is the ability to calculate accurately the evolving densities of dislocations, be they geometrically necessary (GND) or statistically stored (SSD), because of the potential roles played in latent hardening, initiation of recrystallization and nucleation of fatigue cracks, for example. There are now many modelling techniques available in crystal plasticity finite element (CPFE) analysis, reviewed in Roters *et al*. (2010), and for the determination of GND density, and their inclusion in modelling studies related to indentation size effects (Britton *et al*. 2010*a*; Wilkinson & Randman 2010), facet fatigue (Dunne *et al*. 2007*a*;) and micro-deformation (Cheong *et al*. 2005; Cheong & Busso 2006), for example, is widespread. However, what is currently less clear is how well current modelling techniques, and, in particular, CPFE methods, properly and accurately capture independent experimentally observed behaviour; that is, at the level of individual grains, the strain and stress distributions, the lattice rotations and the densities of GNDs. In addition, limited work has been reported on the modelling of GND development in hcp polycrystals. Specifically, the slip anisotropy in hcp crystals is significantly different from that in fcc and bcc systems because of the *c*- and *a*-type slip having differing Burger’s vectors and vastly different critical resolved shear stresses (CRSSs). Slip in such systems occurs on basal, prismatic and pyramidal systems (*a*-type), and on first- and second-order pyramidal systems (*c*+*a*-type). Experimental observations (Britton *et al*. 2010*b*; Wilkinson *et al*. 2010) show that *c*+*a*-type slip may occur at a resolved shear stress of 2.5–3.0 times that for *a*-type slip (Gong & Wilkinson 2009). A consequence of this is that the densities of GNDs developing on the different slip systems during deformation have been measured to be orders of magnitude different (Britton *et al*. 2010*b*). This poses challenges for crystal plasticity-coupled GND modelling techniques that have not as yet been addressed. Hence, the primary aim of this paper is to present detailed comparisons between length-scale-enhanced crystal plasticity calculation and high-resolution electron backscatter diffraction (EBSD) measurement of these quantities in an attempt to assess how well the laudable ambitions referred to above are currently achievable.

The whole area of polycrystal modelling and experimental integration is currently, rightly, receiving much attention; an example of this is the recent 140th annual meeting of the Minerals, Metals & Materials Society, USA, which contained this area as a major theme. One of the earlier detailed studies of CPFE capabilities versus experimental observation was carried out by Kalidindi *et al*. (2004), who made use of pseudo-two-dimensional polycrystal experiments with which to carry out detailed strain comparisons with crystal calculations, a technique used later by others (Dunne *et al*. 2007*b*). Wang *et al*. (2011), Yang *et al*. (2011) and Zhao *et al*. (2008) have similarly carried out detailed and interesting studies of polycrystal deformation with crystal plasticity but, for surface-based measurements, the problem of comparing results with crystal calculations in the absence of knowledge of subsurface crystallography remains an issue. There is little doubt, however, that much progress is still to be made in order for CPFE predictions to be considered good (Wang *et al*. 2011; Yang *et al*. 2011), though work reported so far is based on conventional crystal plasticity as opposed to that which is gradient enhanced. In particular, the significant difficulty of capturing the important features of grain boundary slip transfer has been demonstrated (Wang *et al*. 2011).

Our approach in the present paper is to attempt to enable as close as possible a tie-up between the experimental model and that used in gradient-enhanced CPFE analyses. This paper begins by addressing GND development in hcp polycrystals. The physically based crystal plasticity model of Dunne *et al*. (2007*c*) is used, which is fully length-scale-dependent and contains explicitly the coupling between GND density and slip rate. We present a simple methodology for determination of the GND densities on the individual basal, prismatic and pyramidal slip systems and, importantly, present for the first time detailed comparisons of GND activity on the different systems with experimental EBSD observations for a Ti alloy polycrystal.

We then focus on an fcc nickel single crystal which, after heat treatment, contains distinct and reasonably separated carbide particles of size approximately 30 μm. The material when subjected to stress relief by heat treatment at 870^{°}C and subsequent reduction of temperature to 20^{°}C undergoes significant thermal strain at the single crystal–particle boundaries by virtue of the differing thermal expansivities of the nickel single crystal and the carbide particles. The difference is large enough to cause the development of plastic strain with gradients and hence lattice curvatures and densities of GNDs. Our methodology is to adopt high-resolution EBSD (Wilkinson *et al*. 2006, 2010) in order to measure the elastic strains developed at the particle together with determination of the distribution of GND density from lattice curvature measurements (). We also model a nominally identical material and loading using our GND-enhanced crystal plasticity approach, and present detailed comparisons of the results obtained. Finally, we present the results of the same analysis in which the GND development is decoupled from the slip rule in order to assess the importance of the coupling. Firstly, in §2, we describe the crystal plasticity formulation, the coupling between the density of GNDs and the slip rate, and the determination of GND density from the plastic strain gradients.

## 2. Crystal plasticity finite element implementation

Crystal plasticity kinematics rely on the deformation gradient ** F**, which is usually multiplicatively decomposed into elastic and plastic parts (Kroner 1960; Lee 1969),
2.1where the plastic part of the velocity gradient

*L*^{p}incorporates the crystallographic slip from the active slip systems

*i*, with normal vector

*n*^{i}and slip direction vector

*s*^{i}, and is computed according to a slip rule. The slip rule used here was developed by Dunne

*et al*. (2007

*c*), and is used in the following form: 2.2with Δ

*V*

^{i}=

*lb*

^{i2}, where where

*ρ*

^{m}

_{S}and

*ρ*

^{s}

_{S}are the mobile and sessile SSD densities, respectively,

*ρ*

_{G}is the overall GND density,

*ν*is the frequency of attempts (successful or otherwise) by dislocations to jump the energy barrier,

*b*

^{i}is the Burger vector magnitude for the slip system

*i*, Δ

*F*is the Helmholtz free energy,

*k*is the Boltzmann constant,

*T*is the temperature in Kelvin (K),

*τ*

^{i}is the resolved shear stress,

*τ*

^{i}

_{c}is the CRSS (a slip system is considered active when

*τ*

^{i}≥

*τ*

^{i}

_{c}),

*γ*

_{0}is the shear strain that is work conjugate to the resolved shear stress, Δ

*V*is the activation volume,

*l*is the pinning distance and

*ψ*is a coefficient that indicates that not all sessile dislocations (SSDs or GNDs) necessarily act as pinning points. The superscript

*i*indicates the quantity appropriate to the

*i*th active slip system (for example, the Burger vector for hcp

*a*-type slip is different from that for

*c*+

*a*-type slip). The subsequent implicit integration of constitutive equations and the determination of the consistent elastic–plastic tangent stiffness are also detailed in Dunne

*et al*. (2007

*c*).

### (a) Geometrically necessary dislocations in hcp and fcc crystals

The determination of the GND density has been discussed by many authors (e.g. Nye 1953; Arsenlis & Parks 1999; Acharya & Bassani 2000; Cermelli & Gurtin 2001) and here we equate the closure failure (Nye 1953) written in terms of density of GNDs and the cumulative Burger vector (Acharya & Bassani 2000) given with respect to the deformed configuration to give
2.3where *F*^{e} is the elastic deformation gradient, the summation is over all active slip systems and ** G** is introduced for convenience.

In crystals with high degrees of symmetry, the geometrical constraints on dislocation density, imposed by the plastic slip gradient field, can be satisfied with many different dislocation configurations (Arsenlis & Parks 1999). This is the well-acknowledged non-uniqueness problem. In crystals of such symmetry, the number of distinct dislocation types may exceed the nine independent components of the Nye tensor, hence a unique solution for the dislocation density may not be possible. Therefore, the GND density may be obtained only by solving equation (2.3) with an imposed constraint such as minimization of stored energy or dislocation line length, for example.

Here, we address the problem for hcp crystals for which slip occurs on basal, prismatic and pyramidal slip systems, as shown in figure 1, in which the usual convention of a single slip direction associated with a given slip system is assumed. *a*-type slip occurs on the basal, prismatic and *a*-pyramidal slip systems, while *c*+*a*-type slip occurs on the primary (first order) and secondary (second order) *c*+*a*-pyramidal slip systems as shown in figure 1. The Burger vector magnitudes *b*^{i} for *a*-type slip () and *c*+*a*-type slip () differ by the *c*/*a* ratio. Similarly, the CRSS required to initiate slip on the *a*-type (*τ*^{a}_{c}) and *c*+*a*-type slip systems differs. Experimental observations show that *c*+*a*-type slip may occur at a resolved shear stress of 2.5–3.0 times that for *a*-type slip (Gong & Wilkinson 2009). ** G** in equation (2.3) is then written as
2.4where the first and second terms on the left-hand side of equation (2.4) are the contribution to the dislocation content from

*a*-type and

*c*+

*a*-type slip, respectively, and

*N*

_{a}and

*N*

_{c}are the corresponding number of active slip systems, so that the total number of active slip systems

*N*is given as

*N*=

*N*

_{a}+

*N*

_{c}.

*t*^{i}is defined from the slip direction

*s*^{i}and normal

*n*^{i}as

*t*^{i}=

**s**^{i}×

*n*^{i}. Depending on the active slip systems, the available number of screw-type systems varies. The maximum number of screw-type dislocation systems for the hcp crystal type is nine, if dislocations of opposite sign are not differentiated. The number of edge-type dislocation systems is always equal to the number of active slip systems

*N*.

Once the active edge- and screw-type systems have been established, equations (2.3) and (2.4) may be written as a linear set of equations ** Ax**=

**in which**

*b***(**

*A**m*,

*n*) has dimension

*m*=9 and

*n*=

*n*

_{e}+

*n*

_{s}, where

*n*

_{e}and

*n*

_{s}are the number of edge types (

*n*

_{e}=

*N*) and screw types, respectively. This gives a problem of (

*m*=9) linear equations with (

*n*) unknowns.

*n*may be such that

*n*≤

*m*or

*n*>

*m*depending on the number of active slip systems.

For the screw components, a scheme has to be adopted to enable the screw density component for a given slip direction to be shared among the active slip systems having the corresponding slip direction. In this case, the sharing is performed linearly among the systems according to their activity, where the extent of activity *η* for slip system *i* is defined as —a slip system is not active if *η*^{i}<1. For example, if three systems 1, 2 and 3 have the same slip direction and have extents of activity *η*^{1}≥1, *η*^{2}≥1 and *η*^{3}<1, and if the screw density obtained for this slip direction is *ρ*_{s}, then the screw densities on the three slip systems are *ρ*^{1}_{s}=*ρ*_{s}*η*^{1}/(*η*^{1}+*η*^{2}), *ρ*^{2}_{s}=*ρ*_{s}*η*^{2}/(*η*^{1}+*η*^{2}) and *ρ*^{3}_{s}=0.

Arsenlis & Parks (1999) suggested two ways of solving the non-uniqueness problem, one geometrically and the other energetically motivated. These involve minimizing *L*^{1} or *L*^{2}, which are defined as follows by Kysar *et al*. (2010):
2.5where the *w* are possible weightings. Britton *et al*. (2010*b*) used the energy-motivated linear programming scheme to calculate GND densities from measured (i.e. experimental) curvatures using Matlab. The least-squares scheme minimizes the sum of squares of the resulting dislocation densities. It is a special case of the *L*^{2} minimization for which *w*=1. For the minimization, the Lapack libraries may be used or the same result can be achieved via the Moore–Penrose pseudo-inverse. Whatever the choice, one obtains the exact solution to ** Ax**=

**for**

*b**m*≥

*n*. The case for

*m*<

*n*gives the solution to min (∥

**−**

*Ax***∥**

*b*_{2}), and if there are multiple such solutions, then the solution ∥

**∥**

*x*_{2}is returned. Geometrically, this means finding the shortest possible dislocation line lengths. In the present work, the least-squares scheme has been implemented into the CPFE code. Arsenlis & Parks (1999) argue that the differences between this approach and the linear programming scheme are likely to be small. Kysar

*et al*. (2010), however, show that differences can become significant under conditions of plane strain.

Finally, the evaluation of ** G** from the curl terms in equation (2.3) is carried out as any other continuum variable within the finite element method with the added difficulty of accuracy for terms, which are second derivatives of displacements. The theoretical problem of interpreting

**at boundaries (e.g. a grain boundary) crucially depends on length scale, but, within the continuum-level state variable approach, all variables, including**

*G***, are considered continuum, and a grain boundary is simply being represented as a boundary between two regions of differing crystallographic orientation. Interpretation of**

*G***(or indeed any other continuum quantity) at a grain or other boundary must therefore be carried out with full cognizance of these assumptions.**

*G*## 3. Crystal plasticity analysis of hcp single and polycrystal geometrically necessary dislocation development

In this section, we analyse the development of geometrically necessary dislocations within an hcp Ti alloy polycrystal using the GND-coupled crystal plasticity formulation described earlier. We then compare the results obtained with those from independent experimental observations (Britton *et al*. 2010*b*), in particular addressing the GND generation during deformation on the various hcp *a*- and *c*+*a*-type slip systems.

The material considered is taken to be a near-*α* polycrystal Ti alloy. The material properties arising in the slip rule in equation (2.1) are detailed in appendix A, in which the *a*- and *c*+*a*-type systems are differentiated where appropriate (a *c*/*a* ratio of 1.593 is used). The CRSSs for Ti have been obtained from micro-mechanical bend tests (), and the activation energy chosen to ensure that uniaxial stress–strain behaviour resulting for a polycrystal sample under strain-controlled loading is reproduced, including the moderate strain-rate sensitivity observed in these materials even at room temperature (Dunne *et al*. 2007*c*), and *γ*_{0} has been fixed for simplicity to for slip system *i*, where *d* is a characteristic length set to be the average grain size of 30 μm. Burger vector magnitudes, frequency of dislocation jumps and initial density of SSDs are obtained from standard property data and included in appendix A. Briefly, before considering polycrystal behaviour, we apply the crystal model above for the purposes of validation to a simple problem; namely a single-crystal beam subjected to four-point bending in order to produce a linear variation of plastic strain and hence develop a constant GND density field of known slip activity.

### (a) Four-point bending

A four-point single-crystal bending analysis provides a validity test of the model because of the known form of the strain and stress fields developed between the loading points. The linear variation of strain through the beam section is anticipated to lead to a constant plastic strain gradient and consequently GND density within the plastic zone. The (single) crystal orientation is specified such that the *c*-axis is normal to the beam axis, so that prismatic slip only is anticipated. The HCP single-crystal beam, with dimensions of 200×20×20 μm, is supported and loaded (under displacement control) at distances of 0.5 *h* and 2 *h*, respectively, from the end, where *h* is the height of the beam. A deflection of 2 μm is applied, and the beam is constrained in the *z*-direction such that a state of plane strain is imposed.

The GND density was found to be constant through the thickness in the plastic regime in the pure bending region. This is expected because the plastic strains are linear through the thickness so that the GND density, which is derived from their gradients, would be expected to be constant. The magnitudes in the pure bending section of the beam are within the expected range. Fleck *et al*. (1994) gave a one-dimensional single slip approximation for the density of GNDs as follows:
3.1The maximum effective plastic strain within the beam centre is approximately 0.01, the depth of the plastic region is approximately 10 μm and *b*=3.2×10^{−4} μm. Hence, . The corresponding value of *ρ*_{G}, which has been determined at the beam centre (i.e. away from the loading points) to be consistent with this calculation, is about 5 μm^{−2}. Hence, the predictions are physically reasonable. The GND density is calculated and stored according to the contribution from each individual slip system. In this problem, only the prismatic slip systems contribute to the density for the prescribed conditions, as expected.

### (b) hcp titanium polycrystal deformation

The GND-coupled crystal plasticity model is applied to the study of GND development in a model Ti polycrystal. The polycrystal is subjected to compressive, uniaxial straining in the *x*-direction up to a strain of about 12 per cent and constrained not to move in *x* on the *x*-face, in *y* on the *y*-face and in *z* on the *z*-face. The resulting generation of GND densities on the various hcp slip systems are compared with those obtained independently from high-resolution EBSD carried out on a similarly near-*α* Ti alloy previously subjected to deformation by rolling (Britton *et al*. 2010*b*). Clearly, the precise deformation histories applied for the model polycrystal and in the experiment are different, but it remains useful to carry out qualitative comparisons of the results obtained. The model polycrystal comprises 69 grains with random initial orientation, as shown in figure 2*a*. The variable shown is one of the nine components of the orientation matrix ** g** (specifically,

*g*

^{−1}(1,2)). Figure 2

*b*shows the effective plastic strain developed in the polycrystal showing strain magnitudes up to about 50 per cent with considerable variations indicating the existence of strong plastic strain gradients.

Predicted GND density results for the hcp polycrystal are shown in figure 3, in which the densities developing on the specific differing slip system types are shown together with the overall GND density. It is observed that the basal, prismatic and *a*-pyramidal slip systems generate the highest GND densities, and that the *c*+*a* slip systems provide a much smaller contribution than the others. This is in good qualitative agreement with experiments reported by Britton *et al*. (2010*b*) and A. J. Wilkinson, B. Britton & P. Littlewood (2009, personal communication), who carried out a statistical analysis of their experimental results detailing the (normalized) frequencies of GNDs by slip system type, and these results are reproduced in figure 4 together with an equivalent analysis of the predicted results. The model polycrystal analysis also enables the evolutions of the GND statistics to be obtained, and these are shown in figure 5 at four levels of strain through the loading history.

It can be seen that the predicted and experimental results show GND densities for *c*+*a*-type slip to be an order of magnitude lower than those for *a*-type. During deformation, predictions show that GND densities progressively increase, as does the material volume over which GND activity occurs. Both experimental and predicted results show that the distribution of dislocation densities for each slip system type is in effect skewed normal. The experimental data show more significant *c*+*a* GND activity than is observed in the predictions. Both experiment and model show that *a*-pyramidal edge dislocations are the dominant type. The *a*-screw dislocations appear significant too, but it should be noted that they include all screw contributions from all *a*-slip systems: *a*-pyramidal, basal and prismatic systems. Basal and prismatic systems contribute by about the same amount in both the experimental data and predictions. The histograms show that the predicted GND densities fall within the physically expected range reported in the literature and as shown in the experimental results.

### (c) Discussion of titanium model and experiment

In computing the GND density from the experimentally measured curvatures (Britton *et al*. 2010*b*), it was assumed that all slip systems were active, so that, in the density calculation, the number of unknowns (the densities) was always assumed to be higher than in reality. In the crystal plasticity model, only the active slip systems are accounted for. This potentially provides an explanation for the observation that the *c*+*a* coverage is higher for the experimental results than that from the predictions. During deformation, the texture changes in such a way as to become less favourable for *c*+*a* slip. This change is not captured by assuming all systems are always active, as in the experimental analysis. In addition to the above, in the crystal model, there is full knowledge of the strain gradients (i.e. all the nine components of Nye’s tensor are known); in the experiment, only six of the nine components can be obtained (Britton *et al*. 2010*b*) and the assumption of plane stress has to be made. No such assumption or approximation is needed in the modelling.

An unexpected result is the dominance of the *a*-pyramidal slip in contributing to the GND density. In a separate numerical test, it was observed that if *a*-pyramidal systems are assigned the same CRSS as the *c*+*a* pyramidals then they would rarely become activated. Energy considerations make it reasonable that they should have a similar CRSS to that of the basals and prismatics, and this is what has been assumed in the current modelling. The greater influence, therefore, of the *a*-pyramidals over the prismatic or basal systems may be due to the fact that *a*-pyramidal slip can more readily be activated for orientations suitable for both basal and prismatic slip, so that slip is possible over a wider range of crystal orientations than for either of basal or prismatic slip.

A discussion of GND density is incomplete without reference to the length scale over which such gradients are measured. Kysar *et al*. (2010) argue that it is necessary to choose a characteristic length scale for the volume under consideration, so that the calculated value of the geometrically necessary dislocation density converges to a meaningful value. In the experiments, the pixel (point) size is 0.28 μm, whereas the characteristic element size in the simulation is 4 μm. The finite element resolution-related size effect was checked by Mayama *et al*. (2009), who studied the organization of the GND structure in an fcc bicrystal subjected to cyclic loading. Their results showed that a higher mesh resolution led to a higher predicted GND density but qualitatively the spatial distribution was not affected.

In §4, we address GND development in an fcc nickel system and provide more detailed and quantitative assessments of micro-strain, lattice rotation, curvature and GND density by means of detailed comparisons with experiments.

## 4. Crystal plasticity analysis of fcc single-crystal nickel containing a carbide particle

The material studied is a directionally solidified nickel alloy provided by Rolls-Royce plc. Samples of dimensions 13.5×3×3 mm were cut from the stock provided and subjected to an initial heat treatment of 1100^{°}C for 1 h followed by ageing at 870^{°}C for 16 h. This solutionizing and ageing process was applied twice and resulted in significant grain growth, and subsequent EBSD analysis revealed that in fact single-crystal samples had been obtained containing carbide particles of various sizes up to approximately 30 μm. EBSD also enabled determination of the crystallographic orientation and, in the sample studied here, it was found to be *ϕ*_{1}=83.5^{°}, *ϕ*=91.5^{°}, *ϕ*_{2}=0^{°} in Bunge notation. Following cooling from 870^{°}C, high-resolution EBSD was again used in order to quantify the elastic strain distributions that resulted from the mismatch in thermal expansivities of the nickel single-crystal material and the carbide particle. Full details of the EBSD technique and the determination of the full strain field tensor components, lattice rotations and GND density may be found in Karamched & Wilkinson (2011). The particular particle/matrix combination to be studied is shown in figure 6*a*. In the experimental study, the samples were subsequently also subjected to four-point bend tests, but we confine ourselves in this paper to the thermal loading problem for the specific particle/matrix combination shown in figure 6*a*. An advantage of analysing such a system is that, in many respects, it becomes pseudo-two-dimensional in that there is no (or little) variation of crystallographic orientation through the depth. In addition, because it is single crystal, the often-encountered issue of three-dimensional polycrystal surface measurements is largely eliminated, though the potential variation of the particle geometry through the thickness of the sample remains.

The corresponding geometrical finite element representation of the problem is shown in figure 6*b*, in which the carbide particle extends through the thickness of the model. A full three-dimensional model is developed comprising 11 000 20-noded quadratic three-dimensional finite elements, and results obtained from the front, free-surface are representative of those for a state of plane stress. The same assumption has been used in order to determine the full strain tensor from the EBSD analysis. In the finite element model, the boundaries are all assumed to be free, though analysis has shown that changing the conditions (e.g. such that all planar boundaries remain so) has limited impact on results obtained local to the particle, which are those of interest here. The thermal loading is applied in the model as shown in figure 7*a*. The temperature range imposed is identical to that in the experiments, but the time scale used in the modelling is arbitrary because, although the crystal model incorporates rate sensitivity, these effects are likely to be negligibly small because, for most of the cooling, the matrix material behaviour remains elastic. It is only at the end of the cooling process that the thermal strains developed are large enough to cause the onset of plasticity.

The values of the material properties for the nickel alloy studied, including the thermal expansivities of the single-crystal nickel and the carbide, are given in appendix A. The carbide particle is assumed, reasonably, to deform elastically only.

In the slip rule in equation (2.1), the density of GNDs is coupled through the activation volume term Δ*V* such that any increase in GND density impacts directly on the slip rate. In the analyses presented later, results are obtained in using the fully coupled slip rule and also the case in which the GND density is decoupled. That is, for the decoupled analysis, in equation (2.1), while the GND density is calculated, *ρ*_{gnd} is taken to be zero in the slip rule. The GND density for the fcc nickel system is determined for each active and independent slip system, taking full account of the edge or screw nature of the dislocation. The techniques used are described in §2*a*.

### (a) Strain, lattice rotation and geometrically necessary dislocation density analyses of nickel sample

On cooling the particle/Ni matrix material from 870^{°}C to 20^{°}C, the mismatch in thermal expansivity between the single-crystal matrix and carbide particle leads to the establishment of thermal strains and hence stresses. In fact, the difference in expansivity (see appendix A) is large enough, given the temperature change, to cause the development of plastic strains local to the particle. The extent of the predicted plastic straining, obtained from the CPFE modelling, is shown in figure 7*b*. As a result of the heterogeneous nature of the particle geometry, quite considerable variations and indeed gradients in plastic strain develop and it is therefore anticipated that corresponding variations in lattice rotation and GND development will as a result be observed. An important objective of this paper is to establish to what extent the detailed experimental observations can be captured accurately by the gradient-enhanced crystal model presented above. In figure 8, the results of high-resolution EBSD to determine the elastic strains local to the particle, on the material-free surface, are shown (figure 8*a*–*c*) together with the corresponding CPFE results (figure 8*d*–*f*). The colour scale bars in figure 8 relate to the experimental measurements, and there is no one-to-one correspondence with the colour scales in the predicted elastic strains. In the latter, the range of strain has been chosen to show to best effect the calculated strain distributions, therefore providing, firstly, a qualitative comparison with the measurements. Note that EBSD lattice distortion and rotation measurements, and hence strains, are not possible within the carbide particle. Accordingly, predicted strains *within the carbide* are not shown in the field plots. Quantitative comparison is presented shortly. The thermal strains developed are, in the absence of geometrical constraint, dilatational in nature, and are not therefore captured in the EBSD-determined strain fields. Hence, in order to compare like with like, in the calculated results, the purely dilatational strain has been subtracted from the total (for the direct components; this is not necessary for the shear strains) in order to give the differential strains. In addition, the EBSD-measured strains are, of course, the elastic strains so that again, in ensuring we compare like with like, the plastic strains have been removed from the total strain components in order to give just the elastic strains. Comparison of the figures indicates good agreement in a qualitative sense between measured and predicted elastic strains. Figure 9 shows the experimentally measured and predicted in-plane lattice rotations. Again, good qualitative agreement is observed. Note that the high-resolution EBSD technique allows a resolution of lattice rotation of approximately 2×10^{−4} radians (Wilkinson *et al*. 2006), corresponding to about ±0.01^{°}.

In order to carry out a more quantitative assessment, two lines have been introduced in figure 6*a*. All the elastic strain components and lattice rotations along the vertical and horizontal lines have been extracted from both the experimental measurements and CPFE model predictions and are shown in figure 10*a*–*h*. It is apparent that there are some aspects of the strain and lattice rotations which are in good agreement but clearly others for which this is not the case. For example, the direct strains along the vertical line in figure 10*a*,*c* are found to be significantly over-predicted in magnitude, particularly very close to the particle. In contrast, along the horizontal, the predicted xx elastic strains are under-predicted close to the particle. Some of the shear strain distribution characteristics are captured, others are clearly not. Lattice rotations, however, shown in figure 10*g*,*h* are quite well captured by the model in both magnitude and distribution along the two lines shown in figure 6*a*. A particular implication of this is for the CPFE model predictions for the GND densities, which are measured directly from the gradients of the lattice rotations. The predicted GND density distribution developed as a result only of the thermal excursion is shown in figure 11 together with that determined from EBSD measurements taken from the sample-free surface. Again, reasonable qualitative agreement is obtained but a more quantitative comparison is shown in figure 12, in which, again, the variation of the experimentally determined and CPFE-predicted GND densities along the vertical and horizontal lines shown in figure 6*a* are given. Here, reasonable agreement can be seen in terms of both GND density magnitude and the expected nature of the distribution.

### (b) Comparison between the geometrically necessary dislocation coupled and uncoupled models

It is sometimes argued in the literature that the coupling between the evolving GND density and its effect on subsequent plastic slip is crucial in ensuring the accuracy of crystal plasticity predictions. It is therefore of interest to test this assertion particularly in the light of the quality (or otherwise!) of gradient-enhanced CPFE predictions and their comparison with experimental data at a micro-mechanical level. Hence, the CPFE analysis of the particle thermal problem discussed earlier has been repeated but, on this occasion, the GND coupling with the slip rate shown in equation (2.1) is removed. That is, the GND density is calculated from equation (2.3) just as it was in the earlier analysis, but this is not reflected in the slip rule in which, now, the activation volume, Δ*V* , is assumed to be constant and unaffected by the GND density such that its value remains fixed at
4.1where the terms take the values previously used and given in appendix A. Thus, the results obtained are those for a conventional crystal plasticity formulation as opposed to one which is gradient enhanced. The comparisons of the results obtained for the elastic shear strains, lattice rotations and effective plastic strains along the vertical and horizontal lines shown in figure 6*a* for the Ni single crystal/particle material are given in figure 13*a*–*f*.

The elastic shear strains are shown in figure 13*a*,*b*, and, interestingly, the differences observed for the coupled and uncoupled analyses are really quite small. The lattice rotations, however, particularly those close to the particle/matrix interface, do show quite substantial differences (almost a factor of 2) from the coupled model, predicting substantially smaller rotations at the interface. The experimentally measured lattice rotations at the interface, shown in figure 10*g*,*h*, are in fact much closer to the coupled predictions than the uncoupled model results, and this suggests that the GND coupling is important for lattice rotation (and hence GND density) determination at locations of significant heterogeneity, e.g. at grain boundaries, triple points and so forth. It has long been argued that conventional crystal plasticity has been successful at capturing bulk grain deformation but has been shown to be problematic at grain boundaries (e.g. Kalidindi *et al*. 2004). The results presented here provide direct evidence that gradient effects at locations of considerable heterogeneity do at least contribute to some of the differences observed in the past.

Figure 13*e*,*f* shows the distributions of effective plastic strain in which the effect of the coupling can be seen to reduce the size of the highest plastic strains while the region containing plastic deformation is slightly increased in size, particularly close to the particle/matrix interface. This, of course, is the location, from figures 11 and 12, at which the highest density of GNDs is seen to develop and naturally, through the slip rule coupling, this inhibits subsequent plasticity development relative to the uncoupled model. Quite substantial differences in magnitudes of plastic strain can be seen to be obtained from the two models indicating that at regions of heterogeneity, conventional crystal plasticity is likely to overestimate the extent of plasticity.

### (c) Discussion of nickel model and experiments

In this paper, we have attempted to provide detailed comparisons of the results of gradient-enhanced crystal plasticity predictions with EBSD measurements in a model heterogeneous material system under thermal loading. GND effects are now included in a large range of model analyses for crack nucleation and growth, for the initiation of recrystallization and for the nucleation of twinning, for example. It is therefore of considerable importance to demonstrate that the quantities of relevance at the appropriate length scale are in fact being determined correctly and ideally accurately. Hence, we have aimed to test the ability of gradient crystal plasticity at a particular microstructural heterogeneity.

Some of the observed differences in the measured and predicted results inevitably originate from the modelling assumptions made. For example, while free-surface modelling results have been compared appropriately with free-surface EBSD measurements, in the model, the particle is assumed prismatic subsurface; this is unlikely to be strictly the case in reality. In addition, in the experimental Ni alloy material, while it is known that the carbide particles are well dispersed such that inter-particle effects are likely to be small, they cannot be eliminated and the possibility always exists that additional particles, not visible in figure 6*a*, exist below the free surface and unfortunately close to the analysed particle. Further, the matter of boundary conditions is always important in these analyses; the region of the material (figure 6*a*) modelled in figure 6*b* is in fact contained within a larger material sample, and the correct boundary conditions are therefore difficult to reproduce. However, studies with various different boundary conditions show that the effects close to the particle are in fact small.

In summary, we find that the prediction of elastic micro-straining is qualitatively reasonable but not completely accurate. Lattice rotations, however, and GND densities, are found to be better predicted, in a quantitative sense, by gradient-enhanced crystal plasticity and that conventional, non-gradient plasticity is likely to be more in error particularly where gradient effects are important at, for example, phase interfaces, grain boundaries and triple points.

There remains the important need to continue to challenge and to test the abilities of crystal plasticity modelling techniques by means of detailed comparison of predicted results with micro-mechanical testing at the appropriate length scale. Only in this way can confidence in its use for problems of significance be developed.

## 5. Conclusions

Existing numerical techniques that address the well-known non-uniqueness problem have been developed for the particular use of hcp crystals in which GND density is assumed to develop on all the active independent slip systems. A technique has been presented for accounting for the slip systems shared by screw-type dislocations. GND density is fully coupled with the crystal slip rule adopted.

Analyses have been carried out of a model near-*α* Ti polycrystal subject to deformation. The development of GND densities on the hcp *a*- and *c*+*a*-type slip systems (basal, prismatic and pyramidal) has been compared with corresponding experimental measurements. In particular, the model prediction that *a*-type densities are an order of magnitude higher than *c*+*a*-type is entirely in keeping with experimental observations. The predicted distribution of *a*-type basal, prismatic and pyramidal as well as *c*+*a*-type pyramidal GND densities is found to show reasonable agreement with experimental observations.

Crystal plasticity analyses carried out for a single-crystal fcc nickel alloy containing a carbide particle under thermal loading show that plastic straining occurs local to the carbide–nickel interface, leading to the development of plastic strain gradients and densities of geometrically necessary dislocations. Quantitative comparisons of predicted results with those from high-resolution EBSD show that micro-strain and lattice rotation distributions are reasonably well captured, as are the distributions of GNDs. Coupling of the GND density with the crystal slip rule has been shown to be important in locations of high heterogeneity (such as phase interfaces, grain boundaries, triple points) where lattice rotations and strains can be significantly in error if gradient effects are neglected.

## Acknowledgements

We are grateful to the Engineering and Physical Sciences Research Council (grant no. EP/C509870/1), and Rolls-Royce plc for financial support.

## Appendix A. Material properties

(*a*) *Slip rule*

(*b*) *Elastic moduli*

Ti: *E*_{11}=104 GPa; *E*_{33}=143 GPa; *G*_{13}=47 GPa; *v*_{12}=0.48; *v*_{13}=0.2

Nickel single crystal: *E*=207 GPa; *ν*=0.28

Carbide particle: *E*=2070 GPa; *ν*=0.28

(*c*) *Thermal expansivities*

Nickel single crystal: 20≤*θ*^{°}C≤800; 13.0×10^{−6}≤*α*^{°}C^{−1}≤15.0×10^{−6}

Carbide particle: 20≤*θ*^{°}C≤800; *α*=4.5×10^{−6} ^{°}C^{−1}

- Received January 27, 2012.
- Accepted May 5, 2012.

- This journal is © 2012 The Royal Society