Is stored energy density the primary meso-scale mechanistic driver for fatigue crack nucleation?

Fatigue crack nucleation in a powder metallurgy produced nickel alloy containing a non-metallic inclusion has been investigated through integrated small-scale bend testing, quantitative characterisation (HR-DIC and HR-EBSD) and computational crystal plasticity which replicated the polycrystal morphology, texture and loading. Multiple crack nucleations occurred at the nickel matrix-inclusion interface and both nucleation and growth were found to be crystallographic with highest slip system activation driving crack direction. Local slip accumulation was found to be a necessary condition for crack nucleation, and that in addition, local stress and density of geometrically necessary dislocations are involved. Fatemi-Socie and dissipated energy were also assessed against the experimental data, showing generally good, but not complete agreement. However, the local stored energy density (of a Griffith-Stroh kind) identified all the crack nucleation sites as those giving the highest magnitudes of stored energy.


Introduction
The microstructural and mechanistic origins of the drivers of fatigue crack nucleation remain elusive, and yet the process of crack nucleation may often comprise a significant fraction of overall component life in a range of key industries including those in aerospace, power generation and land transport. In modern aero engines, nickel-based superalloys have been widely used for applications in turbines due to their excellent properties, with high strength, corrosion and temperature resistance. Increasingly, some of these components are produced by powder metallurgy (PM) routes due to their near-net-shape forming, high cost effectiveness, efficiency, and more homogeneous microstructure compared to other conventional forming methods (Reed, 2008). However, it has been found that non-metallic inclusions and agglomerates are unavoidable in these alloy systems resulting from the manufacturing processes, and they significantly degrade mechanical properties and service life of these components. Fatigue failure which is driven by the presence of inclusions shows shorter lifetimes to fracture (Texier et al., 2016a), and as a consequence, their role has been investigated by many researchers (Caton et al., 2004;Pineau and Antolovich, 2016;Naragani et al., 2017), as well as for the effects of grain size distribution (Alexandre et al., 2004), twin and twin boundaries (Texier et al., 2016b), the size and shape of agglomerate (Gabb et al., 2008), which further illustrates the importance of fatigue in the presence of inclusions. However, it still remains a key scientific and technological challenge to establish full mechanistic and physical understanding of microstructure-sensitive crack nucleation at inclusions (Sweeney et al., 2013;Findley and Saxena, 2006;Reed et al., 2008;Sangid, 2013;Miao et al., 2012;McDowell and Dunne, 2010).
Over the past few years, the crystal plasticity finite element (CPFE) method, integrated with experimentally characterized https://doi.org/10.1016/j.ijplas.2017.11.005 Received 6 September 2017; Received in revised form 21 November 2017; Accepted 23 November 2017 manufactured via a PM route. Prior to mechanical test, this sample was annealed at 750°C for 7 h to release residual stress and energy stored within the sample during manufacturing processes. The free surface of this sample was polished to obtain a high-quality metallographic finish using a series of grit silicon carbide papers, finished with colloidal silica suspension for 40 min to remove the residual deformation induced during sample preparation processes and obtain the required surface quality for EBSD characterisation.
Interrupted cyclic three-point bending testing was carried out on this prepared nickel sample with a force ratio of R = 0. For each cycle, the peak force was increased, as described in Fig. 1(b). The sample was loaded to 3650N in the first cycle, followed by subsequent cycles with a load increment of 200N until the peak load reached 5900N. At the final stage (unloaded state from 5900N), several micro-cracks were observed to be apparent local to the non-metallic inclusion. At the unloaded state of each cycle, the sample was removed from the testing rig and loaded into the scanning electron microscope (SEM) chamber for SEM imaging, HR-DIC or HR-EBSD measurements to identify and quantify microcracks, as well as to quantify the slip and GND densities developed.
Prior to the mechanical test, the crystal orientations around the inclusion region were characterized by EBSD, as shown in Fig. 1(c). The measured initial microstructure local to the ceramic inclusion is shown with respect to the loading direction and revealed by inverse pole figure (IPF) map. The inclusion is located at the centre of the EBSD map and highlighted by the white line, as shown in Fig. 1(d), which shows a zoomed in view of the inclusion and surrounding nickel matrix that are modelled explicitly using a three dimensional CPFE beam sub-model. It can be seen from the EBSD map that a great number of annealing twins exist in this nickel polycrystal.
The establishment of GND densities and in-plane strains, under three-point bending, are analyzed utilizing HR-EBSD and HR-DIC techniques, which have been reported in detail elsewhere (Jiang et al., 2016a).

Crystal plasticity finite element modelling
In FCC crystals, there are 12 slip systems that can be activated potentially during plastic deformation (listed in appendix A). In general, the sample's crystallography and stress state are the primary factors to determine whether a slip system is active or not. The plasticity that the specimen undergoes during the experiment comes from the contributions of activated slip systems. Critical resolved shear stress is the criterion used to determine the onset of crystallographic slip.
The crystal plasticity slip rule is rate-dependent and implemented in the user material subroutine (UMAT) using ABAQUS standard/explicit analysis. The total deformation gradient F is split into elastic F e and plastic F p parts, where F e is responsible for rigid body rotation and lattice stretch while F p contributes to the movements of dislocations on slip planes The total deformation rate D is used to calculate the total deformation gradient F, which can be expressed using equation (2), where D e the elastic deformation rate is determined using Hooke's law, and D p the plastic deformation rate is obtained from slip rule.
The plastic deformation rate D p is the symmetric part of plastic velocity gradient L P , expressed as where s α and n α are the line vectors along slip direction and slip plane normal of slip system α respectively, and the slip rate on a given slip system is expressed by Here, ρ ssdm is the mobile dislocation density, b the magnitude of Burger's vector, ν the frequency of dislocation attempts to jump energy barriers, successful or otherwise, ΔF the activation energy for dislocation escape from pinning, with corresponding activation volume ΔV, k the Boltzman constant, T the operating temperature, τ α the resolved shear stress for slip system α and τ c α is the corresponding critical resolved shear stress. The frequency of successful attempts to jump the obstacles is given by the term , considering both forward and backward activation events (Dunne et al., 2007b).
The dislocation-based hardening rule is shown in equation (6), where τ c0 is the intrinsic critical resolved shear stress and G 12 is the shear modulus.
The development of sessile SSD and GND densities provides hardening through shear resistance on each slip system. Here, the evolution of ρ SSD is computed from the accumulated plastic strain rate ṗ thus giving rise to isotropic hardening where λ provides the rate of accumulation and is determined from experiments to ensure appropriate polycrystal hardening, and the effective plastic strain rate ṗ is determined from in which ε p is the plastic strain tensor (see Appendix B for the treatment of this for the experimental DIC measurements). The establishment of density distributions of GNDs is the result of strain gradients accommodating lattice curvature and Nye's dislocation tensor Λ is used to compute the densities of GNDs obtained from where ρ Gs α are screw components on slip system α while ρ Get α and ρ Gen α are edge components, b α is Burger's vector on slip system α, and m α , n α and t α are an orthogonal unit vector set. The 36 independent GND components (12 screw components and 24 edge dislocation components) need to be computed from 9 equations, hence equation (9) is expressed in matrix form: where Λ is the 9 × 1 column vector form of Nye's dislocation tensor, A is a tensor of size of 9 × 36 which contains the basis tensors ⊗ b m α α , ⊗ b t α α and ⊗ b n α α , such that the 9 rows in A incorporate the 9 components for each calculated tensor product, for each of the 36 GND components. However, equation (9) cannot necessarily determine the 36 GND components uniquely. Therefore, the L 2 -norm minimization method is utilized to make the sum of the weighted squares of the resulting GND densities minimal (Arsenlis and Parks, 1999;Cheng and Ghosh, 2015). The L 2 -norm is written as the sum of the squares of the ρ Get α , ρ Gen α and ρ Gs α edge and screw dislocation density components which yields a scalar density of GNDs on each individual slip system (Cheng and Ghosh, 2015), given by

Material properties
Uniaxial tensile testing was carried out on polycrystal alloy FGH96, and the experimentally measured macroscopic true stressstrain curve is shown in Fig. 2. The critical resolved shear stress and SSD hardening parameter λ, were identified and obtained using the uniaxial tensile test data. A 3-dimensional finite element model with representative volume element (RVE), comprising 1000 randomly orientated grains, was developed for calibration. The properties so determined are listed in Table 1. According to Frost and Ashby (1982), the activation energy ΔF can be estimated by = Δ F ωGb 3 , where ω is a constant and depends on the strength of the obstacles. In this study, taking = ω 0.024 in order to give undetectable strain rate sensitivity at room temperature, with shear modulus = G 90GPa, burger's vector = b 2.54e-4μm, then ΔF is computed as 3.456e-20 J/atom. The initial mobile dislocation density ρ ssdm is estimated as − 0.01μm 2 (Hull and Bacon, 2001). The activation volume ΔV is associated with the motion of the climbing dislocation segment and has a magnitude of λ bd p , where d is the length of the activation event and is approximated to be equal to burgers vector b, λ p is the obstacle spacing and is taken to be 30 nm (=11.6b), hence the activation volume ΔV is given by = Δ V 11.6b 3 (Manonukul et al., 2002). The frequency of dislocation escape attempts from obstacles ν is taken to be × − 1 10 s 11 1 (Sweeney et al., 2013) which is considerably smaller than the Debye frequency, determined from Debye temperature and the Boltzman and Planck constants as ∼ × − 1 10 s 13 1 (Cottrell, 1967), in order to minimise strain rate sensitivity in this study. The computed stress-strain response obtained from the model is also shown in Fig. 2, giving reasonably good agreement with experiment data.

Crystal plasticity sub-model
The EBSD image of the region of interest of the tested beam, containing the non-metallic inclusion, is shown in Fig. 1(c), with a magnified view of the non-metallic inclusion shown in Fig. 1(d). A microstructural sensitive crystal plasticity finite element submodel, comprising the inclusion and surrounding nickel matrix, was then developed using 30,292 3-dimensional 20-noded userdefined elements with reduced integration (C3D20R). The microstructure representation of the sub-model is shown in Fig. 3(a), consisting of a homogeneous non-metallic inclusion surrounded by fine-grained nickel matrix. Fig. 3(a) demonstrates the appropriate crystallographic orientations implemented in the CP sub-model. In addition, the microstructures are assumed to be prismatic through this small depth due to the absence of knowledge of the 3D grain morphologies and orientations. This is a necessary simplification but nonetheless enables a quantitative comparison of the key metrics with respect to fatigue crack nucleation sites for fixed, given microstructural representation. Comparative results are therefore useful but clearly, any quantitative, absolute results for the experimental microstructure would ideally need to take full account of the sub-surface grain structure also.
In the modelling, the experimentally characterized ceramic inclusion (coloured in green in Fig. 3(b)) is assumed to be elastic only, the fine nickel grains around the inclusion (given by different colours in Fig. 3(b)) are modelled using elastically anisotropic crystal plasticity with the crystallographic orientations identical to the experiment and obtained from EBSD, and the outer nickel matrix region (coloured in red in Fig. 3(b)) is modelled by crystal plasticity but assigned a reference crystallographic orientation (i.e., [100] parallel to the x-direction etc.) in order to make the computations tractable. This is reasonable since the primary zone of interest, where fatigue cracks nucleate, is at the Ni-agglomerate interface at which the full surface grain morphologies and orientations are fully represented.
The determination of the appropriate boundary and loading conditions for the CP sub-model is necessary in order to obtain the correct mechanical and physical behaviour developed in the experiments. Hence, a macroscopic finite element three-point beam model with Mises isotropic hardening was developed to represent the beam test shown in Fig. 1(a). The boundary and loading conditions were applied to reproduce the experimental loading history shown in Fig. 1(b). Under three-point bending conditions, the inclusion region is deformed effectively under uniaxial loading and the shear effects can be neglected as they are anticipated to be very small near the beam free surface. In addition, the agglomerate region is sufficiently small with respect to the beam and applied loading that it becomes very reasonable to argue that locally, the loading is uniform and uniaxial. In this way, the macroscopic finite element beam model in Fig. 1(a) has been used to establish the average peak loading to impose at the sub-model agglomerate scale shown in Fig. 3(b). However, local plasticity at the agglomerate region means that the local mean stress developed by the three-point bending is not known, and so this has been determined using the CP sub-model to ensure that the experimentally (DIC) measured average strains are captured correctly by the model. The corresponding sub-model applied loading is shown in Fig. 4, and the comparisons between the CP computed and DIC-measured average strains at the agglomerate region are shown in Fig. 3(d). In addition, the spread of strains along the edge regions indicated in Fig. 3(a) and (c) for both DIC measurement and CP computation are shown by the error bars in Fig. 3(d).

Cyclic crystal plasticity modelling and experimental measurement
The effective plastic strains at the unloaded states of cycle 1 (3650N), 3 (4100N), 7 (4900N) and 11 (5700N) predicted by the The plastic strain heterogeneity is revealed at the beginning of the cyclic deformation, as shown in Fig. 5(a) and (e). The strain is accumulated progressively and extends from the inclusion into the Ni matrix with subsequent cyclic loads. The surface effective plastic strains indicate a tendency to form bands at ∼45°to the loading axis arising from directions of maximum macroscale shear. However, increasing strain with cycles shows a strong dependence on microstructure. From observation of the CPFE and HR-DIC maps, plastic strain initiated at the inclusion and the strain pattern is established at the beginning of the cyclic deformation and tends to increase with cycling. Two 'butterfly' deformation bands develop in the HR-DIC measured strains showing significant localized within these two bands, adjacent to the inclusion. However, the qualitative (contour plot) strain pattern is slightly different in the CPFE computation compared to those in HR-DIC measurement.
Only one obvious 'butterfly-shaped' strain band is observed at the top of the inclusion, and the strain bands at the top left and bottom right (indicated using red ellipses in Fig. 5(d)) are not as clear nor high in magnitude as those from DIC measurement in Fig. 5(h). However, quantitative assessment of the plastic strains in the region of interest is facilitated by extracting the strain values from  CPFE calculation and HR-DIC measurement. The effective plastic strain evolution at nine arbitrarily selected points, shown in Fig. 6(a), are plotted in Fig. 6(b)-(j). As expected, both CPFE predicted and HR-DIC measured strains tend to increase with increasing cyclic load, with the exception at location 5. However, it can be seen from the qualitative assessment in Fig. 5 and the quantitative result in Fig. 6 that the strain at location 5 is very small (∼0.1%) such that the experimentally measured strain at this location is likely noise. In addition, at location 4, the HR-DIC measured plastic strain develops much more quickly than that indicated by the CPFE model. The qualitative assessment in Fig. 5 indicates that initially higher strained regions (e.g. 'butterfly' bands) continue to develop more quickly than low strain regions with progressive loading. Hence, with the exception of location 4, the crystal model gives strain evolutions with cycling in good agreement with the DIC measured strains. CPFE predicted full-field Mises stress ( Fig. 7(a)-(d)) and hydrostatic stress (Fig. 7(e)-(h)) are shown at the unloaded state (that is, they are local residual stresses) for cycle 1 (3650N), 3 (4100N), 7 (4900N) and 11 (5700N) respectively, which show considerable heterogeneity around the agglomerate. In addition, careful examination shows that the stress hotspots show progressively increasing residual stress with cycling, reflecting the local accumulation of plastic strain shown in Fig. 5. Comparison of stress distributions with the microstructure shown in Fig. 6(a) reinforce the effects of grain boundaries, local crystallographic orientations and the geometry of the non-metallic inclusion on stress distributions. Mises stress discontinuities arise at grain boundaries on account of the differing crystallographic orientations, and of course the stress distributions near the inclusion are strongly affected by its complex geometry. Consideration of SEM images (see later) shows that the fatigue microcracks tend to nucleate at high stress locations, but they do not always match with the highest stress locations. The Mises stress variations around the inclusion can reach∼1500 MPa.
The change in GND density between initial unloaded state and after 12 loading cycles has been determined from both the CPFE computation ( Fig. 8(a)) and HR-EBSD measurement ( Fig. 8(b)), and the overall GND density is observed to increase significantly from both methods, showing considerable spatial variation. It is noted that the Ni regions above and below the inclusion show high GND density, whereas a low density is observed at the left and right of the inclusion, as shown in Fig. 8(a) and (b). The changes in average GND density (calculated over the region of interest) determined by HR-EBSD measurement and CPFE modelling are 1.05 × 10 14 m −2 and 7.8 × 10 13 m −2 respectively. Experimentally measured GND density increases are found generally to be higher than those predicted by the CPFE model, but this has been observed and addressed by other authors. In addition, the GND density is predicted by  the CPFE model to increase rapidly in the first cycle of loading, and at a much lesser rate with subsequent cycles, which is in a good agreement with previous interrupted progressive GND studies near an inclusion (Jiang et al., 2016a).
The GND density is further quantitatively analyzed by comparing model with experimental results at the unloaded state of cycle 12, along the two paths chosen in Fig. 9(a), and the line graphs showing spatial variation of GND density are given in Fig. 9(b) and (c). The localization of GND density at the inclusion interface is observed in both paths AA′ and BB'. In addition, the predicted heterogeneous distributions of GND density along the two paths are also captured rather well. This, together with the good agreement for slip accumulation, provides the confidence to utilize the CP model to investigate the key crack nucleation drivers at the locations where cracks are observed to nucleate in the experiment.
Note that the color scale bars for strain, stress and GND are chosen to maximize contrast, and where comparisons are presented, color scales for model and experiment results are not the same.

Experimental observations of crack nucleation and CP analysis of local slip
Six microcracks were experimentally observed (Jiang et al., 2016a) to occur around the non-metallic inclusion. These tiny cracks nucleate at or in some cases possibly within the inclusion and then extend to the Ni matrix within one grain. It is noted that the observed cracks tend to extend along slip traces within the Ni grains, indicating a strong dependence on crystallographic orientation and slip system activation. The locations of crack nucleation are identified and labelled in Fig. 10. From knowledge of the SEM images shown in previous work (Jiang et al., 2016a), it is noted that crack 1 is observed to nucleate at the end of cycle 4100N (cycle 3), whereas cracks 2-6 have been found to form at cycle 4900N (cycle 7) and cycle 5700N cycles (cycle 11) respectively.
The inclusion in this study is rather porous, containing an apparent large void at its central region. The presence of these pores may be important for the effective stiffness of the inclusion, such that Niinclusion decohesion is not observed and instead, crystallographic slip-driven cracks in the surrounding Ni alloy are found to develop, unlike that for the agglomerate studied in Zhang's et al. work (Zhang et al., 2015(Zhang et al., , 2016, in which the decohesion of Ni matrix/oxide interfaces is argued to be the main mechanism of crack nucleation. The slip lines on the sample surface shown in Fig. 10 provide confirmation of accumulated plastic deformation local to the agglomerate. The six local crack regions have been identified and assessed experimentally and computationally (CP), as shown in Fig. 11. Crack 1. Is observed to extend along active slip planes within grain 1. From the observations of the activated slip systems, both slip systems 7 and 8 are activated, but it is noted from the model prediction that slip system 7 is more heavily activated and hence dominates. The intersections of these two slip systems with the sample surface illustrate that crack 1 is slip dominated.
Crack 2. Is found to extend from the inclusion to the high accumulated slip region within grain 2, whereas crack 3 is nucleated from the inclusion and propagates along a grain boundary (grain 3). The directions of these two cracks are argued to be crystallographic slip system dominated within grain 2 and grain 3.
Crack 4. Is found to extend away from the inclusion along the boundary of grain 4. Slip system 9 is dominant within grain 4, so crack nucleation and growth is crystallographically controlled and growth occurs along this slip system direction.   11. The six micro-cracks which nucleate at the interface between nickel matrix and inclusion, corresponding SEM images and EBSD maps captured at the end of cycle 12 (5900N). CP predicted slip on relevant activated slip systems is shown (colour legend shows slip magnitude), followed by identification of CP predicted (red) and experimentally observed (blue, broken) slip traces on the sample surface around six micro-crack regions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Crack 5. There are multiple slip systems activated in Grain 5, but slip systems 10 and 11 show greater slip than for systems 1 and 3. The propagation direction of crack 5 is observed to be parallel with the observed slip traces for systems 10 and 11. Crack 6. Is observed to nucleate from the boundary of grain 6 and propagates within the same grain. Slip systems 3 and 12 are predicted to be activated within this large grain. From the observations of the local SEM images and the slip traces on the surface, it seems that crack 6 does not propagate along slip system 3 or 12. In addition, there is some slip found to occur within the twin region, but no micro-cracking at the twin is observed.
It is noted that all bar one of the microcracks are crystallographic in nature; that is, they nucleate and grow along activated slip systems. However, the magnitude of slip is not in itself a strong indicator, in that there are examples of higher slip accumulation but the absence of crack nucleation.

Crack nucleation mechanisms
To further assess the relationship between stress, dislocation density, strain and experimentally observed crack locations (Fig. 10), these key quantities from the CPFE model are extracted along the complete path DD' (highlighted by the red line in Fig. 12(a)) forming the perimeter of the inclusion, at the final stage (the unloaded state of 5900N), as all six microcracks are observed to become very apparent by the end of cycle 12 (5900N). The GND density and Mises stress along path DD' have been plotted in Fig. 12(b) and (c) respectively, and the perimeter locations at which the cracks are observed to nucleate have been identified. From the observations of the line graphs, it is apparent that there is no over-riding correlation between crack nucleation sites and the GND or stress hotspots. Although most of these cracks are generated at locations of high stress or high GND density, cracks are also found not to nucleate at locations of higher or highest stress or GND density, also indicated by the field plots shown in Figs. 7 and 8. Hence, these quantities seem to be relevant to crack nucleation, but not to be sufficient in their own rights to be diagnostic. CPFE predicted effective plastic strain along path DD' (Fig. 12(a)) is plotted in Fig. 13 at the unloaded states of cycle 1, cycle 3, cycle 7 and cycle 11. The strains at the experimentally observed crack nucleation sites have been enlarged in the inset figures for clarity. The sites of crack nucleation show much stronger correlation with the high effective strain locations. It is again apparent from the line graphs that the accumulated plastic strain pattern of spatial distribution around the inclusion remains largely unchanged and that the subsequent cycling has the effect of increasing the strain magnitudes at the locations where slip localization has already occurred. This confirms the same result from the strain maps shown in Fig. 5. However, there are some strain hotspots where no crack nucleation is observed to occur. For instance, in Fig. 13, the peaks at the left side of crack 5 and right side of crack 1 are higher than the strain at crack 4, but no crack is found to nucleate at these sites. In addition, it is also found from the strain distributions shown in Fig. 5 that there are some sites far from the inclusion, which have very high strain magnitudes, but do not generate crack nucleation. The relationship between slip accumulation, persistent slip band formation, and crack nucleation is ubiquitous across the literature in fatigue and hence we argue that is a necessary requirement for crack nucleation, but not sufficient.
From the study of Findley et al. (Findley, 2005), fatigue crack nucleation in Ni-base superalloys has been found to be sheardominated microstructurally and also to correlate well with the Fatemi-Socie (FS) fatigue indicator parameter (FIP) P FS , which is defined by where Δγ max p is the maximum range of cyclic plastic shear strain, σ n,max is the maximum tensile stress normal to the plane associated with the maximum shear range, and σ y is the cyclic yield strength. Hence in this work, we assess this criterion, together with that of dissipated energy, against the experimental observations. In equation (12), k is related to material properties and chosen to be 0.5 for Nickel alloy systems (Findley, 2005;Fatemi and Socie, 1988;Shenoy et al., 2007). We calculate the FS measure at the grain/subgrain scale, and its value along path DD' (Fig. 12(a)) at cycle 12 (5900N), is shown in Fig. 14. The results indicate some correlation between the FS measure and the occurrence of fatigue crack nucleation, but it is noticeable that there are multiple locations around path DD′ for which a high value of the FS parameter is determined but for which no cracks are observed in the experiment.
Following the study of Korsunsky et al. (Korsunsky et al., 2007;Cruzado et al., 2017), the dissipated energy (Jm −2 ) per cycle can be calculated over a complete loading cycle from where σ and ε are stress and plastic strain tensors respectively. The energy dissipated at the end of cycle 12 (5900N) has been extracted along path DD' (Fig. 12(a)) and shown in Fig. 15 from which it is seen that the observed fatigue crack nucleation sites correlate reasonably well with the high dissipated energy locations, but it is also apparent that some high dissipated energy sites don't generate crack nucleation. For example, the two energy peaks between the sites of cracks 5 and 6, which have high dissipated energy, do not occur at sites of crack nucleation. A local stored energy criterion has been proposed by Wan et al. (2014Wan et al. ( , 2016 and by Chen et al. (2017) recently for crack Fig. 13. CPFE computed effective plastic strain and along path DD' (Fig. 12(a) for reference), with magnified views shown at the locations of six micro-cracks. B. Chen et al. International Journal of Plasticity 101 (2018) 213-229 nucleation. The criterion is closely associated with slip localization, stress state, SSD and GND densities, and follows from the idea that it is the local elastic sored energy density, driven by dislocation pile-up and structure formation, over a length scale determined by the local dislocation density, which balances with the surface energy needed to establish the new crack. It is therefore Griffith related in the stored energy density sense, but considered at the microstructural scale, and Stroh-related since the local dislocation density, or number per unit area, contributes to the local stress as the driving force.
In the crystal plasticity implementation, the dissipated plastic energy ΔU (Jm −2 ) per cycle can be calculated as in equation (13). A fraction (arguably ∼5%) of this dissipated energy is stored as dislocation pile-ups and structures, and the storage volume may be expressed in terms of a critical length and an area by where ′ λ is the mean free distance of resultant SSD and GND dislocations. Hence, the stored energy density per cycle (ΔG) can be written obtained as in which ζ represents the fraction of dissipated plastic energy stored in dislocation structure. It is reasonable that this quantity achieves a steady state when evolution rates of accumulated slip and dislocation reach a steady state. The precise magnitude of ζ remains to be determined but is often taken to be approximately 5% of total energy associated with plastic deformation, as the empirically accepted fraction of dissipated energy emerging as heat is well documented to be about 95% (Kamlah and Haupt, 1997;Hodowany et al., 2000;Rusinek and Klepaczko, 2009). In this study, the energy density which is stored in the first and subsequent cycles local to the inclusion are assessed and shown in Fig. 16, from which the pattern of spatial distribution is observed to remain throughout the cyclic loading, as was the case for accumulated plastic strain. The highest stored energy points are located around the inclusion, and the experimentally observed crack nucleation sites are found to coincide with local peak value of stored energy. This is confirmed by the line graph shown in Fig. 17, in which the stored energy density is extracted from the CPFE calculation along path DD' (Fig. 12(a)) around the inclusion at the end of cycle 1, cycle 3, cycle 7 and cycle 11. Again, with cycles, the spatial distribution remains largely unchanged but the magnitudes of the stored energy density increase with increasing load. It is found that the six observed experimental cracks are generated at the six locations at which the highest stored energy density magnitudes are found to occur along the inclusion perimeter. The stored energy density is the only quantity of all of those considered for which this has been found to hold. Comparing Figs. 15 and 17 for the dissipated energy and the stored energy density respectively shows strong agreement, but it is notable that for the stored energy density, the high magnitudes obtained between cracks 1 and 2 are lower than that for crack 4. In contrast, the dissipated energy in Fig. 15 indicates that cracks would preferentially nucleate at multiple other sites in preference to that for crack 4 (since they give  higher dissipated energies). The stored energy criterion therefore, which explicitly accounts for dislocation density, provides a more accurate prediction over all the previous indicator quantities, including Fatemi-Socie and the dissipated energy, for crack nucleation.

Discussion
In this study, microstructure-sensitive CPFE modelling with integrated experimental studies are utilized to provide new insights in understanding the physical and mechanical basis of crack/defect nucleation in PM Ni polycrystalline alloy FGH96, containing a nonmetallic inclusion. It is ubiquitous from many experimental studies that crack nucleation requires the localization of strain, in turn resulting from slip, for fatigue micro-cracks to nucleate. The present experimental and computational studies reinforce this conclusion but with the proviso that the slip accumulation is a necessary condition, but by no means sufficient. Crack nucleation in Ni polycrystalline FGH96 containing a non-metallic inclusion was observed to appear at the Ni matrix/inclusion interface, and extend into the Ni matrix with increasing loading cycles. Crystallographic slip was found to be a major factor at the location of crack nucleation and in influencing the subsequent crack propagation direction. The CP modelling performed well in capturing the experimentally observed slip activation where the cracks were nucleated, as recently reported in Guan et al., 2017), and more generally at arbitrarily chosen points in the microstructure in successfully reproducing experimental DIC strain measurements. In addition, the quantitative comparison of HR-EBSD measured GND densities with the Nye-based formulation for its determination from CP computations also showed very good agreement, in keeping with previous studies (Zhang et al., 2016).
In this paper, the initial GND density near the inclusion quantified by HR-EBSD is consistent with other observations in similar alloy systems. For instance, the average GND density at the interface of carbide/nickel matrix measured using HR-EBSD in a directionally solidified Ni alloy subjected to thermal and mechanical loading was found to be ∼3 × 10 14 m −2 (Karamched and Wilkinson, 2011). An average GND density of ∼10 14 m −2 in Ni polycrystalline with a non-metallic inclusion was also observed and reported in the study of Zhang et al. (2015;. This is found to be consistent with the related work by Jiang et al. (2015bJiang et al. ( , 2013, in which the GND density is measured in deformed copper under cyclic loading. It can be seen from the field plots of GND density shown in Fig. 8 that the geometry of the inclusion and grain boundaries play important roles in the dislocation distribution. The GND density near the inclusion in the current study determined by HR-EBSD is also ∼10 14 m −2 . Following deformation, high GND densities were found at the locations generating heterogeneity such as grain boundaries, triple junctions, sharp corners of the inclusion/matrix interface and twin boundaries where large strain gradients are expected. From the observations of GND distributions (Fig. 8), while high GND density is strongly associated with crack nucleation site, there are many locations with high GND densities where fatigue cracks do not nucleate. This observation is also confirmed by other studies (Zhang et al., 2015;Zhang et al., 2014Zhang et al., , 2016Jiang et al., 2016a;Jiang et al., 2015a). This illustrates that crack nucleation is not uniquely implied from high GND density alone.
Hence consideration of the local stored energy density as a mechanistic driver for fatigue crack nucleation was argued to be reasonable given that all the quantities upon which it depends (local slip, stress, and GND density) have all been demonstrated to be well captured by the CP model. A significant range of experimental anaysis has brought to light the potential roles of microstructurelevel stress, slip accumulation and density of GNDs in fatigue crack nucleation. Some studies have in isolation implicated, for example, the role of GND density in capturing crack nucleation site and subsequent crack growth path (Jiang et al., 2016a). However, (many) other studies reinforce the role of accumulated slip, and yet others that of stress. Hence there is considerable collective evidence that all of these quantities are potentially important in fatigue crack nucleation. The results reported in this paper may provide a unification of the experimental observations through the concept of the local stored energy density. The stored energy density has been shown to be highest, above all other locations, where the experimentally observed cracks are found to occur. The other fatigue indicator parameters (accumulated plastic strain, Fatemi-Socie, dissipated energy) all give one or more predicted locations of high magnitude at which cracks are not observed in the experiments (where the magnitude of the FIP is higher than that for locations where experimental cracks are observed). However, these FIPs do show some success in capturing some of the observed

Conclusion
Fatigue crack nucleation in a powder metallurgy nickel alloy containing a non-metallic inclusion has been investigated through integrated small-scale bend testing, quantitative characterisation using HR-EBSD and HR-DIC, and computational crystal plasticity which replicates the experimental microstructures.
Multiple sites of fatigue crack nucleation were observed to occur at the nickel matrix-inclusion interface and the cracks were found to be driven by local crystallography and slip. Crack growth direction was found to be crystallographic in all cases but one, and where multiple slip systems were active, that with the highest slip magnitude controlled the crack direction.
Crystal plasticity modelling was used to replicate the microstructure in terms of the nickel polycrystal morphology and crystallographic orientations, and the non-metallic inclusion. The CP model was demonstrated to reproduce the experimentally DICmeasured strains and their distributions, and the HR-EBSD measured GND density distributions. At the nickel matrix-inclusion interface where crack nucleation occurred, the CP model provided quantitative detail of slip activation, GND density, stress, dissipated energy, the Fatemi-Socie FIP, and stored energy density so that a quantitative assessment of the drivers of the observed crack nucleation could be obtained. Slip activation, stress, and GND density were all found to be involved, but no one of them proved to be uniquely identifiable as the primary driving force, although slip was found to be a necessary, but not sufficient condition for cracking. Highest magnitudes of local stored energy density, however, were found to correlate with the experimentally observed crack locations, and comparisons were presented with other fatigue indicators (eg Fatemi-Socie and dissipated energy).