A microstructure-sensitive driving force for crack growth

A stored energy based methodology for calculating the driving force for crack growth is introduced which can capture the highly local microstructural sensitivity. This has been implemented in the context of crystal plasticity finite element simulations with explicit representation of the crack with the eXtended Finite Element Method (XFEM), with nonlocal approaches for both stored energy and J-integral calculation. The model is shown to have good agreement with discrete dislocation plasticity (DDP) models in terms of the crack-tip dislocation configurational energy, and with experimental observations of long and very short (microstructurally-sensitive) cracks for both fracture toughness and crack growth rate data. The method is shown to capture the microstructural sensitivity, in contrast with the widely used J-Integral method. By modelling different crack lengths, the diminution of the microstructural sensitivity with increasing crack length is quantified and a critical length defined above which the microstructural sensitivity is insignificant. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The ability to predict engineering component lifetimes relies on the accurate prediction of crack growth from the smallest size to the point of failure. This requires an intimate knowledge of the processes working at all stages of the life of a crack. In ductile metals, the focus of this work, large cracks -many times the microstructural scale -do not interact strongly with the microstructure and can be well characterised experimentally using elastic-plastic fracture mechanics. However, in many applications, cracks spend a great deal of their lifetime at the microstructural scale where they are strongly influenced by the local microstructure, taking tortuous paths and growing at highly variable rates as they pass through different grains. These cracks with length on the order of a few grains are referred to in this work as small cracks. It is therefore desirable to develop a methodology through which the driving forces for crack growth can be characterised and which captures this microstructural sensitivity. In order to be predictive of a wide variety of materials without the need for empirical fitting, it is essential that this model is physically based and ideally applicable to the entire lifetime of the crack.
The focus in this work is therefore on building a general methodology based on the mechanisms and concepts that are equally applicable to the full range of ductile metals in which microstructural sensitivity is observed. The main observations regarding microstructurally small crack growth, occurring both under fracture and fatigue loading conditions, can be summarised as follows: 1. Crack propagation can occur with loading under the stress intensity factor threshold for long crack growth. Or, conversely, for a given applied stress intensity factor, the growth rate is higher for a small crack. 2. There is significant scatter in observed growth rates due to microstructural sensitivity. 3. This microstructural sensitivity diminishes with increasing crack length with a threshold on the order of about 3-10 grains.
Since these effects are common to both fatigue and fracture loading, it is informative to discuss both when considering the microstructural effects on the driving force for crack growth. Both are therefore compared here, but it should be noted that the models in this work are restricted to monotonic loading. Future work will address fatigue loading conditions. HCP metals are of particular interest due to their strong anisotropy which not only magnifies these effects but also strongly effects long crack growth and fracture. The key factor for the long crack growth effect is the alignment of the crystal c-axis (the longer axis of the unit cell, along which HCP metals typically exhibit a higher stiffness) with the loading direction and so the texture can be used together with the applied load to characterise the driving force for long crack growth. A new method of calculating the crack driving force should reproduce these findings for a long crack in addition to the general microstructurally small crack observations noted above.
To make sure that the parameters are physically realistic we use pure zirconium as a representative HCP alloy. Due to the safety critical nature of their main application as fuel rod cladding in nuclear reactor cores, knowledge of the resistance of these materials to crack growth and failure is essential. There has therefore been extensive characterisation of the resistance, primarily as fracture toughness measurements for long cracks, over the past 40 years.
Although microstructurally sensitive crack growth has been shown to occur in reactor grade Zr (see for example the tortuous paths observed by Snowden and Stathers (1977) , to the authors' knowledge there is no publicly available experimental data for microstructurally small crack growth rates in either zirconium or zirconium alloys. Microstructural sensitivity in the growth rates of similar cracks has however been observed in a wide range of ductile metals including magnesium ( Tokaji et al., 2004 ), aluminium, nickel ( Gao et al., 2007 ) and titanium ( Tokaji et al., 1994 ), indicating that they are fundamentally due to the geometry and morphology rather than due to chemistry or material specific slip characteristics. Microstructural sensitivity of the growth rate is therefore also expected in zirconium and so must be captured by the quantity used to measure the driving force for growth in this model.
Energy is the fundamental quantity governing failure at the smallest scales and so this is the quantity of focus in this work. In the context of fatigue, it is clear that the energy build-up which leads to failure is not all provided in one loading cycle, but rather builds over the course of many cycles. This is in line with previous suggestions that, for small crack propagation, more of the energy comes from internal stresses and dislocation structures than from the energy directly applied within a single cycle ( Sadananda, 1997 ). However, more recent dislocation dynamics work has shown that, dependent on there being suitable nucleation sites and obstacles, fracture can also proceed through plastic dissipation ( Cleveringa et al.,20 0 0 ), indicating that long crack growth may also be governed by the same processes. The methodology presented here thus aims to calculate the total energy available for crack propagation.

Microstructural sensitivity for small cracks
Given the paucity of published experimental work at the microstructural level in Zirconium, it is necessary to look to studies in other metals. A significant difference that occurs in nickel superalloys when going from physically large cracks to small cracks is that large cracks tend to grow along grain boundaries whereas small cracks tend to grow through grains. This behaviour is potentially modulated by temperature and environment. These small cracks have been noted to grow preferentially along slip planes when growing through the grains ( Zhai et al.,20 0 0 ) and often have highly tortuous paths. Importantly from a fatigue perspective, small cracks propagate at far lower stress intensity factors than large cracks. There are many factors which have been identified to affect the rate and direction of crack growth, including grain-boundary interactions, bridging of ductile regions, crack deflection due to hard particles and phase transformations ( Kumar and Curtin, 2007 ).
There are two related effects of grain boundaries on small crack growth: deflection of the crack path and retardation of the growth rate. The deflection of the crack path is illustrated by the tortuous path taken by the crack in Fig. 1 . There has been significant discussion about the reasons for the direction of crack deflection, but we restrict ourselves here to discussion of the driving force for crack growth. The retardation of the growth rate at grain boundaries has been the subject of further studies ( Gao et al., 2007;Schäf et al., 2013 ), but the interaction is not simple as it is a function of both the grain boundary misorientation and the inclination of the grain boundary. High angle grain boundaries have been found to slow growth more than low angle grain boundaries. This deceleration gives an improved small crack growth resistance in finer microstructures ( Gao et al., 2007 ), the opposite finding to that of large crack growth studies. The Zhai model has had success in characterising these grain-boundary effects ( Zhai et al.,20 0 0 ) by relating the crack growth retardation to the tilt and twist angles between the crack planes on each grain. Work by Holzapfel et al with 3D microstructures has highlighted that some deceleration at low-angle grain boundaries, an apparent discrepancy with the Zhai model, can be explained due to the presence of sub-surface precipitates impeding crack growth ( Holzapfel et al., 2007 ).
Work in mild steel indicates that cracks are less likely to cross grain-boundaries ( Schäf et al., 2013 ) and failure tends to occur by the agglomeration of smaller cracks rather than a single crack growing to a critical size. Given the differences noted  Gao et al. (2007) . between FCC and BCC materials, it seems likely that further differences might be expected for the highly anisotropic HCP materials. Fundamentally, the modulation of the growth rates at grain boundaries tells us that the driving force is affected near the grain boundary. This will be discussed in more detail later.
A major barrier to studies of grain boundary effects has been the lack of sub-surface grain-boundary information. This has led to previous studies simply extrapolating from the surface information ( Gao et al., 2007 ). The development of focused ion beam tomography has opened the possibility of investigations with 3D characterisation of the grain boundaries ( Kubis et al., 2004 ;Uchic et al., 2007 ), but this still has the drawback that it is a destructive technique which can only be done after the crack growth analysis. More recently, detailed 3D datasets are becoming available from X-Ray synchrotron experiments. The power of these large datasets is demonstrated in work such as that of Rovinelli et al. (2018) , who showed reasonable agreement in terms of crack growth rates using a probabilistic model.
It has been previously noted in theoretical studies that dislocations at persistent slip bands produce stress concentrations under the same assumptions as LEFM ( Brown, 2013 ). This suggests that there may be some similarities between the crack growth problem and that of crack nucleation. Given that energy-based methods have proven able to accurately predict crack nucleation sites it is reasonable to assess the applicability of these methods to crack growth. Taking this approach gives the same advantage of a more physically based model.
Recent works using crystal plasticity, both studying nickel alloys, have taken this approach, including that of Wan et al. (2016) and that of Farukh et al. (2016) . Farukh et al have used a critical value for the plastic strain, which showed good growth rate agreement in previous work using a viscoplastic material model ( Farukh et al., 2015 ) albeit without explicit consideration of microstructure. The crack growth direction was orthogonal to the maximum principal stress direction. The growth rate was clearly dependent on the grain crystallographic orientation, however the path was purely transgranular, whereas the experimental path followed a mixed transgranular and intergranular path. Wan et al used a maximum principal stress criterion together with an energy-release rate based traction-separation response to study crack growth in real ferritic steel microstructures ( Wan et al., 2016 ). This work highlighted the importance of including elastic anisotropy to capture the highly-localised stress states at the grain boundaries. A further finding was that the explicit inclusion of grain boundary effects was found necessary to correctly capture the crack path.
The crack growth rate appears to depend upon the energy release rate in the material ahead of the crack tip. Thus, by modelling the energy around the crack it may be possible to predict the growth of the crack. An advantage of using an energy-based approach to this problem is that the relative energies associated with the generation of crack surfaces at grain boundaries and within the grain are quantifiable using atomistic methods, giving a firm physical basis to the predictions. This has not yet been attempted but may be able to provide a unified approach to the crack growth problem. A further issue with both methods described above is that neither one considers the crystallography directly with regards to the path, even though cracks at this scale are often observed to grow along slip planes.
Concerning long cracks, there is some data for Zirconium alloys. However, very few fatigue studies have been conducted and thus fracture toughness data are essentially the only available measurements of the resistance to crack growth. Stress intensity factor studies have indicated a critical value of 43.9 MPa √ m for β treated Zircaloy 4 ( Kreyns et al . , 1996 ) which, assuming isotropic elasticity with a Young's modulus of 88,0 0 0 MPa, corresponds to a Griffith energy of 19.4 kJm −2 . Grigoriev et al. (1996) found a fracture toughness of 100 kJm −2 in cold-worked Zircaloy 2 with similar values after irradiation and slightly higher values for annealed samples. Hsu and Tsay (2011) found a similar value of 97 kJm −2 for Zircaloy 4. The most widely used standard for fracture testing specifies plane strain conditions; however, as pointed out by Raynaud et al. (2008) , the geometry of the tubes gives specimens resulting in plane stress conditions except in the case of through-thickness cracks. This means that the stress state at the crack tip may be significantly different from that found from experimental testing of axial or radial crack growth.

Driving forces for crack growth
There are many approaches that have been taken for modelling crack growth, which can broadly be categorised according to the parameter used for the driving force. This can be stress, as is common in fracture and long crack growth models, plastic strain, energy or dislocations. As will become clear, these parameters are tightly coupled, and hence they can give comparable results.
The first quantitative analyses of crack growth driving forces were based upon isotropic elastic models and were initially focussed on fracture. The Griffith failure criterion gives the point of maximum free energy, after which failure is energetically favourable. This was later modified by Irwin to include plastic dissipation giving where σ f is the stress at fracture, a is the length of the crack, E is the Young's modulus, γ is the surface energy and G p is the plastic dissipation energy. γ and G p have units Jm −2 . A related development was that of stress intensity factors.
These factors characterise the stress state at the crack tip under three different loading modes. Mode I is opening, with loading perpendicular to the plane of the crack, mode I I is in-plane shear, with loading parallel to the plane of the crack and perpendicular to the crack front and mode I I I is tearing with a shear stress applied parallel to both the crack plane and the crack front. For a crack oriented with the crack plane perpendicular to the y axis the stress intensity factors for the three modes are defined as where K x is the stress intensity factor for mode x , r is the distance perpendicular to the crack front and σ ii is a stress component. For plane stress problems, this is simply related to the strain energy release rate by where G is the strain energy release rate with units Jm −2 . The stress intensity factors have been used extensively to characterise crack growth, with a critical value of the mode I stress intensity factor, K IC , often being given for fracture in brittle materials. The K IC is dependent on the morphology, texture, chemistry as well as other factors and so must be experimentally determined along with any orientation dependence for each new set of these parameters. For example, fracture toughness testing of an extruded magnesium alloy by Somekama et al showed significant deviation in the measured K IC when loaded either parallel to or normal to the extrusion direction ( Somekawa and Mukai, 2005 ). This variability of K IC represents a major limitation. K IC is not a fundamental material parameter and so it's applicability in predictive modelling is highly limited without extensive experimental characterisation. These concepts have also been utilised to characterise fatigue crack growth. Paris' growth law ( Paris et al., 1961 ) gives the rate of growth of the crack as a function of these stress intensity factors where a is the crack length, N is the number of fatigue cycles, K max and K min are the maximum and minimum stress intensity factors and C and m are material constants. Presentation of fatigue crack growth rate data is now predominantly done in terms of K. Whilst this gives reasonably reliable results for long cracks, for small cracks it is not refined enough to capture the local effects that lead to crack propagation. Although this approach is not appropriate at the microstructural scale, any new approach would still need to agree with K approaches in the limit of longer cracks.
Whilst LEFM gives a good description of very brittle materials such as glass and other materials at large length scales, in ductile materials the assumptions of elasticity break down as the plastic region becomes significant -i.e. when the crack size is microscopically small. When the plastic region is small, the J-integral gives a value equivalent to, and for linear elastic materials identical to, G . It has been shown that the integral is path independent but must travel from one face of the crack to the other. The J-integral, as introduced by Rice (1968) is given by where is the path around the crack tip, W is the strain energy density of the material, T is the traction vector, u is the displacement vector, s is the arc length along the contour and x and y are the distances parallel to and normal to the crack respectively. The units of J are conventionally given as kJm −2 .
The J integral is usually used to quantify the fracture toughness, but provides an important measure of the energy required for crack propagation. For application to fatigue problems, the J-Integral has been further extended to the 'cyclic J-Integral' ( Lamba, 1975;Chow and Lu, 1991 ). This has been applied in the case of large cracks, with Banks-Sils et al finding that the cyclic J integral provided results which were comparable with K approaches under elastic or small-scale yielding conditions ( Banks-Sills and Volpert, 1991 ), but hasn't been applied to microstructurally small cracks due to concerns about the applicability under large scale yielding.
Another measure of the driving force is the crack-tip displacement, δ. This is related to the values of the orthogonal crack-tip opening displacement ( CT OD ) and the crack-tip sliding displacement ( CT SD ) via Although inherently a phenomenological measure, δ is closely related to the J integral, and for an elastic-perfectly plastic material the relation is simply where m is a constant ( ∼ 1 ), σ ys is the yield stress and δ is the crack tip displacement. Use of the CT D has been justified by linking to dislocation pile-ups in atomistic simulations ( Farkas et al., 2005 ). This has shown that the CT D is intrinsically linked to the build up of dislocations at the crack tip which has also been demonstrated as a mechanism of crack growth.
In agreement with the experimental observations noted earlier covering driving forces at grain boundaries, the CT OD has been demonstrated to be strongly affected around grain boundaries in crystal plasticity models by Potirniche and Daniewicz (2003) . These observations thus support the idea that the dislocation pile-ups are an essential part of crack propagation. This was also found to be true in the context of fatigue, using the CT OD -the difference between maximum and minimum CT OD across a loading cycle. Köster et al. (2010) have also utilised the CT D , with a Paris type law to model crack growth in Ti6Al4V and showed a remarkable agreement with experimental results. However, this implementation relied on explicit consideration of dislocation motion which is infeasible over larger length scales in 3D.
Since, as noted above, dislocations have been demonstrated to provide a viable atomistic mechanism for crack growth it is logical to attempt to use dislocations directly in the crack growth model. We now move our focus to models at the dislocation scale. These models are prevalent in both crack growth and fatigue crack growth studies.
The double-slip mechanism proposed by Neumann (1969) for stage II fatigue crack growth shows that the combination of two slip systems can give rise to crack extension in a direction which does not correspond to a single crystallographic slip system. Work by Wilkinson et al. (1998) demonstrated a link between the dislocation activity at the crack tip, in particular the slip irreversibility, and the driving force for stage I fatigue crack growth. This has also been developed by Deshpande et al. (2002) for fatigue cracks growing along metal ceramic interfaces. Importantly, this approach was one of the first to be applicable to both monotonic and cyclic loading. It was noted in this paper that the material properties can strongly influence the relevance of the dislocations to crack propagation under monotonic loading -materials with large numbers of dislocation sources tend to undergo crack-tip blunting, whereas low numbers lead to brittle fracture. As with that work, the focus here is on the intermediate regime where the dislocations contribute to crack growth. Work by Olarnrithinun et al. (2013) has shown the importance of the equilibrium dislocation dipole spacing as a key factor in determining the fracture toughness.
In the face of the clear evidence that a major driver of crack advance is the build-up of dislocations, it may seem necessary for the model to explicitly include dislocations. However, a key limitation of discrete dislocation dynamics simulations is that only small length scales are accessible due to computational expense. Thus, although the explicit inclusion of dislocations into a model would be a more rigorous method for addressing crack advance, the modelling of polycrystals which are more representative of experiments and components is intractable with current computing resources. It is therefore desirable to develop a continuum scale model that is able to capture the broad effects of the processes occurring at the dislocation scale, at a reasonable computational expense, without necessarily capturing the fine detail.
In this vain, many more global approaches have been introduced. The problem encountered by all of the global approaches is that they are not capable of detecting the highly localised effects of plasticity which are essential to accurately model cracks at the microscopic level. In the context of fatigue, these limitations have driven the development of Fatigue Indicator Parameters (FIPs). These methods use a variety of parameters to predict crack growth, but the focus is very much on achieving good agreement with experiments, or mesoscale measures such as the CTD or plastic deformation fields, rather than a rigorous interpretation of the underlying physics, thus they are likely to be far more material dependent. The Fatemi-Socie FIP ( Fatemi and Kurath, 1988 ) is based on the maximum cyclic plastic shear strain γ p max and the peak normal stress on the maximum plastic shear range plane σ max n to model the failure of the material under fatigue loading: Castelluccio and co-workers utilised the Fatemi-Socie FIP to model cracks of varying size in copper ( Castelluccio and McDowell, 2012 ). This was evaluated against the CT D in order to assess the effectiveness of the FIP. They found that the two parameters are correlated and that the implementation of localised slip bands had a significant effect on both. However, the correlation of the FIP and the CT D was found to be dependent on the loading mode. More practically, the mesh refinement was found to have a strong effect. Since the crack was explicitly included in the model and the mesh it is difficult for this model to include the effect of propagation.
Another implementation based on the FS FIP is that of Musinski and McDowell (2012) , in which the effect of notches on the driving force for crack growth through the microstructurally-sensitive crack (MSC) growth period was assessed. The MSC growth regime is treated statistically, without an explicit representation of the crack, to give a probabilistic prediction for crack growth to the point at which LEFM techniques are applicable.
Although more refined than the global approaches, the links between FIPs and the underlying physics of crack propagation are often not well developed and the reliability and precision has also been questioned ( Rovinelli et al., 2017 ). Thus they may be sufficient for generating a statistical description of crack growth in a characterised material, but developing a fully predictive framework would require explicit consideration of the driving forces. An approach which combines both a knowledge of the dislocations and the energies for propagation would thus provide a model with a stronger mechanistic foundation.
In the context of crack nucleation, and building on the concepts of dislocation pile-ups leading to crack propagation, a stored energy approach has been proposed by Wan et al. (2014) and Chen et al. (2017) . This links the energy available for crack growth with energy stored in the dislocation structures formed during plastic deformation. These structures are energetically unfavourable and thus a crack propagating through these structures, and annihilating them, releases this energy. For this crack propagation to occur, the energy of these structures at the crack tip must be raised to some threshold level at which propagation is energetically favourable. The stored energy is thus a highly local quantity analogous to the Griffith energy -far lower than the global energy dissipation given by the J-Integral.
The stored energy is calculated from the density of the plastic strain energy that is not dissipated. This energy, E ρ , is stored in dislocation structures. This quantity is calculated by taking a proportion ( ξ = 5%) of the total plastic work This calculates the energy density at each material point, but to predict crack nucleation or growth we need to calculate the amount of energy that can be provided at the crack tip. The critical length scale over which this energy is transferable to the crack tip is argued to be that of the dislocation mean free path, which can be calculated as λ This is shown diagrammatically in Fig. 2 . The stored energy that is available for growth at the crack tip per unit of new crack area can be determined by the integration over the loading history: In terms of crack nucleation, this approach is clearly related to the magnitude of plastic strain, which has previously been shown to be indicative of crack nucleation. However, the stored energy also gives consideration to the stress under which the deformation occurred and the dislocation densities. The stored energy has been shown to be predictive of crack nucleation sites and cycles to nucleation in Nickel alloys ( Wan et al., 2014;Chen et al., 2017 ).
The stored energy therefore seems to have strong potential as a measure of the driving force for crack growth. In this paper we aim to address the following: 1. Is the stored energy valid and useful as a measure of the driving force for crack growth? If so, it would need to capture crack advance for averaged macro-scale behaviour (ie replicate the J integral), as well as behaviour dependent on microstructural features local to the crack tip. 2. How does the driving force change with increasing crack length, and when does the effect of local microstructure on crack growth become negligible?
A non-local form of this stored energy as a measure of the driving force for crack growth is presented. The paper is divided into two sections. The first section includes an introduction to the methods and models used in the work, along with the essential theoretical and mechanistic basis. In the second section representative results which show the validity of the stored energy as a measure of the crack growth driving force are shown, followed by an investigation of the microstructural sensitivity; how the microstructural sensitivity changes with crack length and when it becomes negligible.

Methods
This section outlines how the key methods have been implemented within the framework used for this investigation: the material constitutive model, which uses crystal plasticity; the stored energy implementation with associated non-local averaging scheme; a brief overview of XFEM, the method by which cracks are included, the crystal-level J-Integral calculation scheme and finally the discrete dislocation plasticity (DDP) framework with the dislocation configurational energy density explicitly calculated.

Crystal plasticity
The deformation gradient F can be decomposed into the elastic and plastic parts; F e and F p respectively.
In the crystal plasticity method, the plastic velocity gradient is calculated by summing the contributions from all slip systems which are activated according to Schmid's rule: where N s is the total number of slip systems, ˙ γ i is the slip rate on slip system i and n i and s i are the corresponding slip plane normal and slip direction respectively. The slip systems present in zirconium are shown in Fig. 3 . The slip rule originally proposed by Dunne et al has been shown to accurately represent the deformation of various materials including Zirconium:  Table 2 Critical resolved shear stresses for the slip systems in pure zirconium.  where ρ m is the density of mobile dislocations, ν the frequency of attempts of dislocations to jump obstacle energy barriers, b the Burgers vector, F the thermal activation energy, k the Boltzman constant, T the temperature (295 K), τ i and τ i c the resolved shear stress and critical resolved shear stress on slip system i respectively, and V CP is the activation volume. This model is derived from the notion of dislocation movements being pinned. A thermodynamic driving force causes the dislocation to escape these pinning points which leads to rate sensitive slip. The parameters used here for pure zirconium are listed in Table 1 and the critical resolved shear stresses for each slip system in Table 2 .
The critical resolved shear stresses used in this work are taken from micro-cantilever tests by Gong et al. (2015) , except the < a > type pyramidal, which was assigned the same CRSS as the prismatic system, and < c + a > 2nd order pyramidal, which was given a value 5x that of the 1st order since it is not observed experimentally. Hardening on these systems due to the evolving GND density is introduced according to: where ρ SSD and ρ GND are the sessile densities of statistically stored dislocations (SSDs) and geometrically necessary dislocations (GNDs) respectively. The GND density is determined from the relationship between the Nye dislocation tensor and the plastic strain gradient ( Dunne et al., 2012 ), and the sessile SSD density is taken here to be constant. To validate the parameters above, a cuboidal representative volume element was created with 9 × 9 × 9 grains, each made up of 5 × 5 × 5 elements as shown in Fig. 4 . Grain orientations were randomly assigned from the experimental texture. The resulting stress-strain response matches bulk uniaxial compression results well.

Stored energy, non-local methodology, XFEM and J-integral implementation
The local stored energy density introduced above in Eq. (15) is determined within the crystal plasticity formulation from This has been implemented as a discrete sum so that at each integration point and for each step the increment of stored energy G is calculated in the following manner: These increments are summed through time at each integration point: Since the GNDs are only calculated at the end of each time increment, there is a small error introduced since the value of ρ GND used in the calculation of the stored energy is that of the previous time increment. However, this error can be minimised by using very small time increments. Since the increments with a crack present are necessarily small, the change in ρ GND over a single increment is less than 1% and hence the error in the calculated stored energy is much less than 1%.
G is a local finite element quantity and thereby suffers potentially from problems of mesh refinement and dependence. In order to remove such effects, a non-local calculation method is introduced utilising the approach of Nguyen et al. (2015) and full details are provided in Appendix A . The eXtended Finite Element Method (XFEM) has been utilised to introduce cracks into the model. This has been done using the XFEM implementation in Abaqus software with the propagating crack formulation which does not include crack tip enrichment. It is therefore effectively identical to a cohesive zones approach when the crack is stationary, which is the case in all simulations presented in this work. However, these models can be used in future for investigations of cracks propagating in arbitrary directions without any changes to the model. Full details of the XFEM formulation are given in Appendix B .
J-integral calculations are needed in order to relate microstructure-sensitive stored energy with macro-level and experimental fracture toughness test data. In the context here, the J-integrals are calculated at crack tips within both single crystals and polycrystals where the latter introduces grain boundaries and discrete regions of elastic anisotropy, and yet it needs to be demonstrated that the J-integral calculation remains path independent. Full details and demonstration of path independence and agreement with analytical calculations for the given material system and model are given in Appendix C . It should be noted that, although the J-integral measured here is sufficiently reproducible in this material system and microstructure, this is not necessarily generally applicable and so care should be taken to ensure this before application in other systems.

Discrete dislocation plasticity formulation
The two dimensional discrete dislocation plasticity formulation is used to simulate the dislocation activities with the presence of a crack. The classic superposition method ( Giessen and Needleman, 1995 ) is used to resolve the boundary conditions where the strain and stress field are written as The dislocation thermally activated escape from obstacles ( Zheng et al., 2016 ) is implemented in the DDP model such the room temperature strain rate sensitivity of Zirconium alloys can be captured. Each obstacle is assigned a thermal activationrelated, stress-associated time parameter t obs , which is equal to the inverse of the successful jump frequency, i.e. t obs = 1 / obs . The frequency of successful jumps is governed by a similar rule to that used in the crystal plasticity model as] where ν D is the frequency of attempts of dislocations to jump the energy barrier, l obs the average obstacle spacing, τ the local shear stress act on the pinned dislocation and V DD is the activation volume in DDP model. More details about the DDP formulation are well described elsewhere ( Zheng et al., 2016 ).
The free energy in DDP is introduced as by Benzerga et al. (2005) . with which consists of the elastic stored energy W e , dislocation self-energy and core energy E l and the dislocation configurational energy E con f which maintains the dislocation structure under the present stress state. From this starting point, we propose that the dislocation configurational energy can thus be calculated from where D is the stiffness matrix. Note that the dislocation configurational energy is defined here as the energy difference between the considered dislocation structure and the energy associated with the dislocations in the structure were they isolated (hence non-interacting) in the structure, hence the configurational energy can be either positive or negative depending on whether the dislocation structure is more or less stable than isolated dislocations. The physical motivation here is that the propagation of a crack requires energy sufficient to generate new crack surface area, so what is the origin of this energy? Plastic deformation results from slip activation through dislocation glide. The progressive (monotonic or cyclic) development of slip leads to dislocation pinning, to the establishment of lattice curvature carried by geometrically necessary dislocations, and the interaction of glissile dislocations with those which are sessile. Dislocation pile-ups develop, often forming structure in the form of cells, ladders, or more complex forms. The dislocation structures so formed store energy during plasticity, and this energy is what we term the dislocation configurational energy. It is the energy associated with the interaction of dislocations driven by externally applied loading. It is argued that the majority of the energy associated with dislocation structure is generated by geometrically necessary dislocations. If so, considering the dislocation structures formed at a crack tip, energy is stored locally by the GNDs supporting the lattice curvature. It is conceivable that sufficient energy could be released by reducing or eliminating the lattice curvature supported in this way to provide a sufficient driver for new surface generation and hence crack extension. This in turn would lead to a substantial GND density reduction because of the release of the lattice curvature. Hence in summary, our hypothesis is that crack extension occurs because it is energetically favourable for new surfaces to be generated by the diminution of local lattice curvature supported by GNDs. The dislocation configurational energy may be calculated using a discrete dislocation analysis, and this is carried out here in order to demonstrate that its magnitude and spatial distribution can be reasonably captured by a higher length scale crystal plasticity quantity; that of the local stored energy introduced earlier in Eqs (15) and (20) . Thus an objective is to demonstrate that the stored energy density calculated using CP approximately replicates the lower length scale quantity dislocation configurational energy determined from DD analysis. In this way, a mechanistic basis will have been established for the CP stored energy density as a driver for crack extension.
The stored dislocation configurational energy is a calculable, physical quantity that can provide the energy required for crack propagation. For this, the local density of this energy, which is available to provide energy to the crack tip, must be considered.
The local configurational energy area density which is transferable through dislocation structure and is available for crack growth is defined as where is the average configurational energy volumetric density over the volume V L and λ d = 1 / ρ dis ( V L ) is the corresponding dislocation mean free path in V L . The dislocation mean free path depends on the volume size over which the dislocation density ρ dis ( V L ) is averaged, hence an appropriate grid size (here 0.3 μm ) needs to be determined such that the characteristics of the dislocation distribution are well captured. The dislocation configurational energy density is that quantity which the crystal plasticity stored energy density, G ( Eq (20) ), attempts to capture, but from a higher length scale approach, hence this is assessed later.

Microstructure, texture and crack length studies
This section addresses the investigation of the crack-tip stored energy density as a measure of the driving force for crack growth. This is realised through multiscale modelling, first comparing the dislocation and stored energy fields against those predicted using discrete dislocation dynamics. The relationship with the J-integral is then investigated and model predictions compared with experiments in the case of both single crystal and textured polycrystal models. The section also addresses the relationship between crack length and microstructure, and the bounds on the microstructural sensitivity of crack growth. Throughout this section key quantities are investigated with respect to the applied stress intensity factor ( K app ) as opposed to the applied far-field stress ( σ f f ). The former has been calculated using the standard analytical solution for an edge crack under applied loading σ as follows where the lengths a and b are as shown in Fig. 5 . In this way, the effects of specimen size on the elastic field are eliminated so that crack length effects may be independently addressed.
One of the key problems with using stress-intensity factor or the J-Integral as a measure of the toughness of a materialthe resistance to crack growth -is that these quantities are not microstructurally sensitive in a predictive sense. For example, polycrystal texture is known to influence fracture toughness which therefore has to be determined empirically. Similarly, single crystals, fine or large grained polycrystals, or multi-phase materials all give differing fracture toughnesses which need to be obtained from experiment. Furthermore, small cracks exhibit different behaviour to longer cracks. Hence critical values of the stress intensity factor and J-integral energy are not intrinsic material properties. The stored energy density is argued to be a more fundamental property of the material which provides the driving force for small crack and large crack behaviour, and one which enables microstructural sensitivity to follow directly. Since for long cracks LEFM results are characteristic, we can assure that the stored energy is predictive by comparing the predictions to experimental results in the long crack regime. In this regime, the grain size is negligible compared with the crack length and cracks do not exhibit significant interactions with the microstructure.

Comparison with discrete dislocation analysis
A key question to answer before using the CPFE implementation of the stored energy is whether the dislocation activity, and in particular the dislocation stored configurational energy, is being adequately captured to be predictive of the propensity for crack growth. Since CPFE modelling captures a continuum scale description of the dislocation activity within the model, it will not match the dislocation-level detail, but a qualitative match would be sufficient to make an informed prediction of the toughness.
A DDP model representing zirconium was developed using the same experimental parameters derived from the microcantilever testing. As can be seen in Fig. 6 , the agreement in the stress-strain response of a single crystal under uniaxial loading between the DDP model and the CP model is very good. This DDP model and the CPFE model were therefore applied to identical single crystal models which include an edge crack under uniaxial loading. Two crystal orientations were chosen: one hard, with the c-axis parallel to the loading direction (loaded to 10 0 0 MPa), and one soft, with the c-axis perpendicular to the loading direction (loaded to 300 MPa). Contour plots showing the dislocation configurational energy area density and the corresponding stored energy density for the DDP and CP models respectively are shown in Fig. 7 a- There are a number of important similarities to highlight here. Most generally, it can be seen that the shape of the stored energy distribution is similar across both models. The stored energy field in the soft orientation is diffuse, extending far away from the crack tip, whilst in the hard orientation the energy field is focused near the crack tip, including a significant portion behind the tip. Looking at the intensity of the field, there is significantly more dislocation activity and more stored energy in both models for the soft orientation, despite the much lower load.
In terms of the slip system activity, it can be seen that the identity and distribution of slip systems activated in both models are also similar: exclusively prismatic slip occurring in the soft orientation, with the same locations of dominance, and < c + a > pyramidal ahead of the crack with < a > basal behind in the hard orientation.
Fundamentally, the stored energy is predicted to be greater in the soft orientation rather than the hard orientation in both the DD and the CPFE analyses. Since the CPFE model is able to capture this difference, as well as the distributions of the key quantities that are thought to drive crack growth, this method could form the basis of a continuum scale model for  addressing crack growth. To validate this model further, we must also compare with experimental results. We therefore now move to compare the CPFE predictions with experimental results at higher length scales.

Edge-cracked single crystal behaviour
The loading and boundary conditions for the analysis of cracked single crystals are shown in Fig. 8 . The model consists of a 30 μm thick plate subjected to a monotonically increasing, uniaxial, stress-controlled loading in the X direction. For each model an edge crack extending through the entire thickness in the z-y plane is introduced from the top (positive y direction) at the centre point in the x direction. The crystallographic orientations chosen for the single HCP Zr crystals are indicated in Fig. 8 . Stored energy Eqs (20) -(22) and the J-integral Eqs (8) -( (10) ) are evaluated at the crack tip as the remote loading is increased, and both quantities are determined numerically non-locally utilising the methodology described above and in Appendix A .
The results are shown in Fig. 9 from which there are three key points to draw. First, both the J-Integral and stored energy approaches predict a necessarily higher crack driving force for crack growth in the S crystal orientation (c-axis perpendicular to loading direction) compared with the H orientation (c-axis parallel to the loading direction). Second, the effect of crystallographic orientation only becomes significant for the J-Integral at higher loads greater than about 10 MPam 0.5 ( Sadananda, 1997 ) -when the plastic zone becomes large -whereas the orientation effect is much larger for the stored energy approach from much lower loads of about 4 MPam 0.5 ( Sadananda, 1997 ). The J-integral gives energies per unit area which are three orders of magnitude higher than those for the stored energy density ( Eq (20) ). This last point highlights the key difference between the two parameters; the J-Integral captures energy associated with both the applied loading and the plastic work that has been done, whereas the stored energy is more focused on the portion of the plastic work done locally at the crack tip which is available for crack growth. Since much of the plastic energy measured by the J-integral is well removed from the crack tip, it cannot be directly relevant to the crack growth mechanism until the crack length is  considerable and the plastic zone size much bigger than the local microstructural feature size. Fig. 9 (a) and (b) enable the relationship between local stored energy and the J-integral energy to be quantified (since the value of each is known for given applied stress intensity) and these results are shown in Fig. 10 . This figure then facilitates the link between experimental macro-level J-integral measurements and the corresponding local stored energy.
Once the stored energy reaches a critical value, it is argued that it is energetically favourable for crack growth to occur. This point has been extensively demonstrated for long cracks within the framework of fracture mechanics in terms of the plane strain K IC and plane stress K C values, which are related to each other by in which ν is Poisson's ratio. The critical values of K IC for varying texture have been determined for long cracks in strongly textured samples of zircaloy4 by Walker and Kass (1974) and in zircaloy2 and zircaloy4 by Cockeram and Chan (2009) , Cockeram and Chan (2013) . Cockeram and Chan (2013) analysed the mechanism of fracture and showed that in alphaannealed zircaloy-2 at the K C point (and below the stress intensity factor for steady-state growth) fracture was controlled by slip and slipband cracking, upon which this model is based, rather than through other mechanisms involving voids or particles.
Texture is often described in zirconium alloys by means of Kearns factors, which give an estimation of the alignment of grain c-axes in a polycrystal to a set of orthogonal sample axes. The Kearns factor can be calculated as by Cockeram and Chan (2009) using f ≡ π 0 I ( ϑ ) cos 2 ϑ sin ϑdϑ π 0 I ( ϑ ) sin ϑdϑ (30) where f is the Kearns factor in a given direction and I(ϑ ) is the average pole figure intensity at an angle ϑ from the given direction. The factor takes a value between 0 and 1, with 1 indicating perfect alignment of the basal pole parallel to the given direction, and a value of 0 indicating all basal poles are perpendicular to the given direction. In the experiments of Cockeram et al Cockeram and Chan (2013) , the critical stress intensity for the particular zircaloy2 polycrystal tested with Kearns factor of 0.29 in the loading direction was determined to be 40.7 MPam 0. 5 ( Sadananda, 1997 ). The corresponding critical J-integral value may then be determined from to give 17.76 kJm −2 (using E = 93.3 GPa). Since zirconium is elastically anisotropic, an effective value of E can be taken from the elastic portion of the J-Load curve. Because the experimentally tested samples were argued to be strongly textured, we make the simplification that the textured polycrystal with Kearns factor 0.29 in the loading direction may be considered approximately equivalent to a single crystal with the corresponding c-axis direction for the purposes of determining its fracture toughness (the demonstration that this is reasonable is given in a later section). In particular, if the critical experimental J-integral value for a given crystal orientation is known, then Fig. 10 allows the intrinsic critical stored energy to be determined. Hence, with the critical J-integral value from experiment of 17.76 kJm −2 , the intrinsic critical stored energy is determined from Fig. 10 to be 7.25 Jm −2 . With knowledge of the intrinsic critical stored energy, G c , it is argued that the critical fracture toughnesses K IC for all textures then become predictable from stored energy. The predicted results for zircaloy2 for a critical stored energy of G c = 7.25 Jm −2 are shown in Fig. 11 along with the experimental measurement for zircaloy2 and zircaloy4. Agreement between predicted K IC and experiment is good for zircaloy2, and reasonably good for zircaloy4 (taking the same critical stored energy for zircaloy2) with the exception of one anomalous looking experimental point for zircaloy4 with the Kearns factor of 0.675. Contour plots showing plastic slip, stored energy, XX stress and GND density for a sub-set of the crystal orientations are shown in Fig. 12 a, loaded to the point at which the J-integral is 5 kJm −2 . The slip systems that are active are shown in Fig.  12 b and are different across the set of orientations, broadly in line with the Schmid factor in the loading direction. For the S orientation the slip is almost exclusively prismatic, whereas the H orientation has a significant amount of < c + a > pyramidal slip activated.
There are multiple effects at work here. One is the shape of the slip field -a field that is more focused (with a higher intensity of slip per unit area) ahead of the crack tip will give a lower toughness. Second is the amount of slip, which is determined by the critical resolved shear stress, the (local) Schmid factor and the slip rule -more slip leads to higher stored energy and thus a lower toughness. Third, the GND evolution -a lower GND density gives rise to a higher rate of accumulation of slip, but to a lower local stress (since there's reduced lattice curvature). The slip field for the S orientation leads to the stored energy distribution being more extensive, whereas for the 60 orientation it is more focused ahead of the crack tip (see Fig. 12 a). The differing slip system critical resolved shear stresses and Schmid factors lead to higher slip in the lower angle orientations. The final effect arises from the GND density; in the S orientation, the GND density at the crack tip is much lower than in the other orientations, which leads to uninhibited slip and hence high stored energy and low toughness.
Overall, the single crystal results thus show good agreement with long crack experiments, but of more interest in this study is the microstructural sensitivity of crack length resulting from crack interactions with grain boundaries and polycrystal grain orientations, which are assessed next.

Microstructural sensitivity and texture
For the purposes of this study, uniformly sized square grains were introduced into the central region of an edge-cracked polycrystal model. This ensures that investigation of crack length is not obscured by other effects such as grain morphology, or proximity of crack tip to grain boundary effects, and provides consistency and repeatability of conditions. Grain crystallographic orientations are assigned from an experimentally measured texture in order to be representative of the original orientation and misorientation distributions. An outer region taken to be plastically isotropic is also included in order to eliminate far field boundary effects. The model, orientation distributions and boundary conditions are shown in Fig. 13 . The particular grain in which the crack tip is located was assigned to be either in the S (c-axis normal to loading) or H (caxis parallel to loading) orientation as described above, in order to examine in some detail the role of the local crack-tip microstructural sensitivity. This was done for three (macro) textures: a texture representative of that found in service, the  Edge-cracked polycrystal model loading and boundary conditions, and the three orientation distributions used to investigate textured behaviour. Texture A gives a higher proportion of crystal c-axes parallel with loading direction (a 'hard' texture) as does C but to a lesser extent, even though it is more random; texture B has a higher proportion of crystal c-axes normal to the loading (a 'soft' texture). same texture but rotated with respect to loading direction, and a random texture -labelled A, B and C respectively. For each texture, three realisations were generated in order to investigate scatter.
The three textures differ significantly in terms of the average c-axis alignment with the applied loading direction. The c-axes in texture B are most likely to be aligned perpendicular to the loading axis (average angle 73.4 degrees, Kearns factor 0.13), favouring < a > type slip, whereas texture A exhibits slightly more alignment with the loading axis (average angle 54.4 degrees, Kearns factor 0.34), favouring < c + a > type slip. Texture C is more random and has a slightly lower proportion of c-axis alignment with the loading axis (average angle 56.6 degrees, Kearns factor 0.35). These differences mean that B is a 'softer' texture while A and C are 'harder' textures. Due to the dependence of the J-Integral on both the elastic and the plastic material response it is expected that its value increases in cases of increased plasticity. Thus 'harder' textures -in this case textures with increased c-axis alignment with the loading direction -would be expected to show a lower J-Integral value for a given applied load.
The calculated J-Integral results for the set of three representative textures are presented in Fig. 14 a. The anticipated relationship between the J-integral and the c-axis alignment is observed in the model predictions; that is, as can be seen in Fig. 14 a, the J-Integral values for 'soft' texture B are higher than that for 'hard' textures A and C. If the J-Integral represents the energy driving crack propagation at the crack tip, this suggests that increased c-axis alignment with the loading direction imparts increased crack growth resistance to the material.
This shows that varying the texture can have significant effect on the energy concentration at the crack tip. This agrees with experimental work, which has shown that the critical value for crack propagation is dependent upon the alignment of the c-axis with the loading direction. Critically however, although the crack length is just 270 μm and some microstructural sensitivity would be expected, this is not observed in the J-integral behaviour. As seen in Fig. 14 a the J-Integral does not change significantly when the crystal orientation of the grain at the crack tip is changed from hard to soft orientations. Thus, while the J-integral might be sufficient to investigate macroscale effects due to the texture, it does not predict significant microstructural sensitivity and so is not appropriate as a criterion for microstructurally-sensitive crack growth.
The stored energy is also affected by the change in texture, as shown in Fig. 14 b. Since it is a much more local measure, the orientation of the grain at the crack tip has a much greater effect on the stored energy, which can be seen to differ considerably depending on the crystal orientation of the grain in which the crack terminates. Since a growing crack potentially needs to propagate through hard and soft orientated grains as it is progressing, with the nature of the grains strongly determined by the texture, the microstructural sensitivity is anticipated to magnify the differences due to the surrounding grains when comparing two textures, for the case when the crack growth is microstructurally sensitive. In order to investigate the predicted sensitivity of the critical J-integral value (for fracture) to texture, an arbitrary critical stored energy of 1 Jm −2 is considered, and the corresponding predicted critical J integral values for varying loading direction Kearns factor are shown in Fig. 15 . Considering the average value obtained for each polycrystal texture, it can be seen that the single crystal results are reasonably representative of the average behaviour even with a small sample size. This analysis therefore supports the earlier case made that a single crystal with c-axis orientation chosen to match that direction with the highest number of c-axes in a strongly textured polycrystal is reasonable, providing confidence in the model predictions shown in Fig. 11 .
To assess the sensitivity of the driving force to the distribution of crystal orientations in the grains in the neighbourhood of the crack tip it is useful to compare the variability resulting from the different texture realisations. As can be seen in the spread of the individual realisations of each texture in Fig. 15 , the observed J integral at the critical stored energy shows significant variability with regards to the distribution of orientations in the neighbouring grains even when the grain orientation at the crack tip remains unchanged. This shows that not only is the grain immediately at the crack tip important for small crack growth, but so too are the grains immediately surrounding that grain. The implication of this is that the misorientation distribution is also key in determining the crack driving force.

Crack length and microstructural sensitivity
The crack length has been shown to be a very significant factor in the microstructural sensitivity of small crack growth rates. This microstructural sensitivity is greater for smaller cracks and diminishes with increasing crack length. Another experimental observation is that small crack growth occurs below the thresholds established for long cracks. Here, the dependence of the stored energy at the crack tip on the length of the crack is investigated. The boundary conditions utilised are identical to those in the above texture study, but a random crystallographic texture is modelled, and the model is summarised in Fig. 16 a. An edge crack of varying length is introduced (indicated by the colour bar) with respect to fixed microstructure and texture. For the purposes of investigation of microstructural sensitivity, again the crystallographic orientation of the grain containing the crack tip is assigned either to be hard (H) or soft (S) with respect to the loading direction. Fig. 16 b shows the calculated J-Integral versus applied loading K app for differing crack lengths. For a given applied stress intensity, the J-Integral transpires to be higher for shorter crack lengths, with this effect diminishing with increasing crack length. This agrees with experimental observations that cracks grow at applied loading below the long crack K threshold when cracks are small. However, when comparing the results with a hard grain at the crack tip with those for a soft grain, it is clear that the J-integral again predicts very little microstructural sensitivity.
However, moving on to the stored energy calculations for the same crack and microstructural configurations, it can be seen in Fig. 16 c that the driving force at a given K app is again predicted to be higher for shorter cracks, but the difference between the two crack tip stored energies for differing local crystal orientation is now large. At lower loads this microstructural sensitivity is even larger than the difference due to increasing crack length. This is due to the differences being driven by the plasticity, which is less significant at lower loads. The crack length, L μ , above which the crack's microstructural sensitivity diminishes to a sufficiently small level is practically very important since it determines, for a given microstructure, the crack length for which classical fracture mechanics, including Paris law for fatigue crack growth, becomes appropriate. For crack lengths less than L μ , the behaviour is microstructurally sensitive such that critical stress intensity and J-integral methods cease to be appropriate and more fundamental, microstructurally 'aware' methods are needed. Hence the present microstructurally-sensitive stored energy density methodology is compared against existing experimental observations for crack growth rates.
The ratio between the stored energies calculated with either the hard or soft grain orientations at the crack tip provides a dimensionless measure of the degree of microstructural sensitivity. This ratio does have a dependence on K app , and has been calculated for K app = 5 MPa √ m over a range of characteristic crack lengths, and the resulting variation in predicted microstructural sensitivity is shown in Fig. 17 by the circular and square symbols (for large and small sized samples to show negligible size effect). The characteristic crack length accounts for the remaining ligament length and so is equivalent to the length of an edge crack in a semi-infinite plate that would give the same stress intensity factor for a given load. It is well established that much of the scatter in growth rates stems from crack interaction with grain boundariesparticularly in the case of arrest. However, the purpose of the present study is to identify the effect of the crack length on the microstructural sensitivity. As such, the calculations were done with the crack at a fixed distance from the next grain boundary and only the relative change in the driving force for each orientation is of interest.
Experimental studies in pure titanium by Tokaji et al. (1994) show a decrease in the coefficient of variation of the growth rate with increasing crack length. These measurements were done with two different grain sizes (fine and coarse grain diameters of 73 and 115 μm respectively) with both giving similar results for a given crack depth (normalised by the grain size). These results are overlaid on Fig. 17 . The results were normalised against the sensitivity at a crack length a = 1 . 5 d, where d is the grain size, to allow comparison between experimental and model results. The relative microstructural sensitivity follows the same trend in both the titanium experiments and in the current zirconium predictive model. Again, as noted before, the absolute magnitude of the sensitivity is likely to be greater than that calculated here due to other features not considered such as grain boundary interactions and crack path deflection. However, the results presented provide evidence that the relative reduction in microstructural sensitivity is almost entirely a function of the characteristic length and the grain size.
It can be seen that there is no clear length L μ at which the crack becomes completely insensitive to the microstructure. However, in the absence of a clear line, we can define L μ to be the relative crack length at which the microstructural sensitivity is half as strong compared to its original length. Defining the critical length in this way gives a critical relative length for the present HCP material of ∼20 grains below which the crack may be deemed to be microstructurally sensitive, and above which the microstructure becomes less important and for which conventional fracture mechanics approaches become relevant.
Experimental studies of microstructurally small fatigue crack growth tend to be done at constant stress amplitude, meaning that the stress intensity factor is a function of the crack length. If an edge crack is imagined to be growing through a specimen of large dimensions, the stress intensity factors can be calculated using K app = 1 . 12 σ √ π a . For a given applied loading the stress intensity factor at each crack length can then be calculated. The predicted stored energy for each of these stress intensity factors can then be used to predict the growth rate of a crack under that loading. By overlaying the curves  Tokaji et al. (1994) for microstructural sensitivity in pure titanium as a function of crack size (relative to the grain size). The experimental values have been uniformly scaled to match the plot axes. Fig. 18. Stored Energy at 200 MPa and at 300 MPa plotted against stress intensity factor for increasing crack lengths. For each given load, the shaded area shows the difference between the stored energy when the grain at the crack tip is in the S orientation vs the H orientation -the microstructural variability of the driving force. generated when the grain containing the crack tip is in both orientations, the scatter can be estimated. Given that there are many different factors which contribute to the crack growth rate, such as grain boundaries, the estimates of the scatter presented here would be expected to underestimate the amount of scatter present in real data. However, Fig. 18 shows the stored energy range (scatter) expected for cracks growing under constant loading at two different stresses. Since the stress intensity factor is dependent on both the far-field stress and the crack length the variability seen for a given K app for the 300 MPa far-field stress is representative of a shorter crack than that for the equivalent for the 200 MPa case. This again demonstrates that for a given applied loading the microstructural variability is lower for a longer crack. Similar experiments done in titanium ( Tokaji et al., 1994 ) show that for a given value of K app there is an increased average crack growth rate at higher levels of applied loading. This observation is also predicted by this model.

Conclusions
A microstructurally-sensitive stored energy density has been introduced as a measure of the energy available for crack growth. This has been shown to be able to capture the distribution and magnitude of the dislocation configurational energy as calculated using discrete dislocation dynamics, and is also predictive of the experimentally determined critical stress intensity K C as a function of polycrystal texture (Kearns factor) when using homogenised single crystal models.
These results have been further confirmed using textured models, for which results averaged over a number of realisations of the texture show agreement with the single crystal results, whilst also exhibiting strong microstructural sensitivity.
The experimentally observed increase in driving force coupled with a reduction in the microstructural sensitivity with increasing crack length is predicted using the stored energy. The relative reduction in the microstructural sensitivity with increasing crack length has been shown to be quantitatively comparable with that observed experimentally.
The stored energy density provides the determination of the critical crack length L μ above which microstructural sensitivity diminishes such that classical fracture mechanics approaches (eg critical stress intensity and J-integral measures) become appropriate. Conversely, for crack lengths less than the critical length, it provides a quantitative determination of the stresses needed to drive crack growth for a given microstructure.

Appendix B. XFEM
XFEM allows the introduction of arbitrarily aligned discontinuities, for example cracks, in the displacement field of an element. This is simply implemented through the addition of an enrichment term to the displacement field and the standard FE basis functions are still used ( Belytschko and Black, 1999;Belytschko et al., 2001;Moës et al., 1999 ). In the original implementation of XFEM the enrichment functions for a given element are chosen according to three possible conditions: 1. a crack passes completely through the element, 2. a crack tip is within the element or 3. the element is not cracked. For case 3 clearly no enrichment is required. For case 2 the enrichment functions contain crack tip functions however these are not used in the Abaqus implementation of XFEM and so are not considered here. The enrichment functions are thus required for elements that are completely bisected by the crack.
To illustrate the enrichment functions for bisected elements consider the displacement field for a single 1D element with nodes positioned at u 1 = 1 and u 2 = 2 . The original FE displacement field, shown in Fig. B1 , is given by where u ( X, t ) is the displacement field, n is the number of nodes and N I and u I are the shape functions of and displacement at node I. The XFEM method adds enrichment functions q to equation ( 34 ) in order to generate a discontinuity. These enrichment terms use the same shape functions as the standard XFEM terms but contain Heaviside functions ( H ) which give rise to the discontinuity and set its position at X = a . An example set of enrichment functions with a = 0 . 5 are shown in Fig.  B2 . The displacement field is then given by Equation ( 35 ); the original XFEM formulation as proposed by Moës et al. (1999) .
The Abaqus implementation of XFEM utilises the phantom-node formulation proposed by Song et al. (2006) . To reach this formulation, equation ( 35 ) is rearranged to give equation ( 36 ). To understand this rearrangement, notice that the   Heaviside functions cause the terms to be active exclusively on one side of the crack or the other. These two distinct sets of terms can thus be thought of as two separate, superposed elements each of which is then only integrated up to the crack from one side. Element 1 has nodes at u 1 and u 2 − q 2 , whilst element 2 has nodes at u 1 + q 1 and u 2 . The new nodes are referred to as 'phantom nodes'. Fig. B3 shows the resulting displacement field from using the displacement field of element 1 on X < 0 . 5 and of element 2 on X > 0 . 5 .
In higher dimensions, use is made of the level-set method to track the position of the crack. This uses a function f( X) in the place of ( X − a ) in the 1D example.
XFEM has been shown to accurately reproduce analytical solutions for stress states around crack-tips in 2D, with a reasonably low mesh-dependence ( Moës et al., 1999 ). The main benefit of this method over other methods is that the crack is not constrained to grow via mesh gridpoints nor does the mesh need to be regenerated after each growth step. Further testing has shown good agreement of the stress intensity factor ( Rannou et al., 2010 ) and crack growth rate with experimental measurements ( Ferrié et al., 2006 ) when coupled with a Paris-type growth law.
This method has been applied to the modelling of crack nucleation and growth. The enrichment functions are used to introduce a crack into the system and to simulate the development of stresses around the crack. The two crack growth studies introduced earlier made use of this methodology ( Wan et al., 2016;Farukh et al., 2016 ).

Appendix C. J-Integral calculation
The J-integral calculation has been implemented as a post-processing script in python. In order to make the code as general as possible, it is evaluated at a circular series of points around the crack tip and ignores the mesh. Since these points do not coincide with node or integration point locations, all parameters are calculated by linearly interpolating the nodal values. Due to limitations with the Abaqus implementation of XFEM, only linear elements can be used with a crack present. This presents a problem for the displacement gradient calculation since the errors at each point can be significant. The gradients are therefore calculated over a few elements rather than using the gradient of the element on the path.
The calculation steps are shown in Fig. C1 . Typically around 100 points were used for the integration, which was found to be sufficient for good agreement with the theoretical calculations. To check the validity of the implementation, it was applied to an isotropic elastic model for which the theoretically correct value can be calculated. As can be seen in Fig. C2 , the agreement between this implementation and the theoretical value is extremely close.

Path and Mesh Independence
The J-Integral is a global variable and as such, for the value obtained to be useful, the J-integral should be path independent. If the value is path-dependent, then it is very difficult to determine the applicability of any measurement made. This path-independence has been proven in the case of a non-linear elastic material under proportional loading and it is valid with small-scale plasticity, but not for the case of anisotropic elastic-plastic behaviour with grain boundaries. Thus it is important to demonstrate that the J-integral gives a consistent value.
To determine whether the integral was path independent in an elastically anisotropic material with crystal plasticity, the J-integral was measured with different contour radii. This was first done for a single crystal model and a polycrystal model with 50 μm grains and a 10 μm element size. The J-Integral calculated at different radii for a single crystal model is shown in Fig. C3 a. It can be seen that the J-integral is essentially insensitive to the path over which the contour integral is calculated.
There is a small difference in the value at 26 μm, which is mainly due to the fact that the path only goes through ∼20   C2. Comparison of the J-Integral calculated using the method presented here and that calculated using the implementation built in to Abaqus with the theoretical value for an isotropic elastic material. elements. Since the element basis functions are required to be linear, the gradients calculated at this small radius are less accurate.
The same plot for a polycrystal model is shown in Fig. C3 b. Again, it can be seen that the path has little effect on the calculated J-Integral. Finally, the J-Integral was calculated using a model of the same polycrystal with a refined mesh around the crack tip. Refined and original meshes are shown in Fig. C3 c. The result is shown in Fig. C3 d, with the J-Integral calculated at 50 μm on the unrefined model shown for comparison. These results show that the J-Integrals being calculated are converged with respect to mesh refinement at the crack tip as well as the radius of the contour over which they are calculated. Fig. C3. J-Integral calculated along contours at different radii around crack tip for a) single crystal model with c-axis perpendicular to loading direction, b) polycrystal model. c) Fine and coarse mesh around crack tip, with some example contours plotted in the coarse mesh with radius 26, 51 and 101 μm d) J-Integral calculated along contours at different radii around crack tip for polycrystal model, with the fine (fi) and coarse (co) meshes.