A mechanistic modelling methodology for microstructure-sensitive fatigue crack growth

A mechanistic methodology for simulating microstructurally-sensitive (tortuosity and propagation rate) fatigue crack growth in ductile metals is introduced which utilises the recently introduced dislocation configurational stored energy as the measure of the driving force. The model implements crystal plasticity finite element simulations using the eXtended Finite Element Method (XFEM) to represent the crack. Two methods of predicting the direction of growth (based on the crystallographic slip or the maximum principal stress) are compared. The crystallographic slip based direction model is shown to predict microstructurally-sensitive fatigue crack growth in single crystals which displays many features of path tortuosity that have been observed experimentally. By introducing a grain boundary, the crystallographic model is shown to capture behaviour similar to that observed experimentally including crack deflection and retardation at the grain boundaries. Finally, two experimental examples of fatigue cracks growing across three grains are analysed, and the model is shown to capture the correct crystallographic growth paths in both cases.


INTRODUCTION
Accurately quantifying the resistance of ductile metals to fatigue crack growth is an essential part of predicting component lifetimes. Current models can capture this resistance when the cracks are many times longer than the microstructural scale, however, a sizeable portion of the crack lifetime for many components is spent with the crack length of the same order as the grain size. At this scale, the crack interacts strongly with the surrounding microstructure, taking a tortuous path and growing at highly variable rates. Three important features that significantly affect the fatigue lifetime are the behaviour at grain boundaries, the effect of crystal orientation, and the presence of other phases.
Traditional fracture mechanics techniques are not able to capture these effects. Furthermore, the focus for fatigue crack growth modelling at this length scale has been on phenomenological models which ignore the path and aim to capture the growth rate as a function of some set of indicators. This is partly due to the lack of a consistent mechanistic model for crack growth.
To be able to accurately predict the lifetime of a crack in this regime, as a function of microstructure including texture, a physically based model must be developed which is able to explicitly capture the growth of the crack. Numerous experiments, such as those of Herbig et al. [1], have shown that the rate and path of crack growth are strongly coupled.
There are therefore two essential components of such a model: 1. A measure of the driving force for crack growth, to determine the rate of growth.
2. An understanding of the growth mechanism, to determine the crack path.
Despite the requirement for both of these, there is currently no established measure for either. Recent work has shown that the dislocation structure stored energy [2] gives a physically-based, microstructurally-sensitive measure of the driving force for crack growth. In terms of the direction of propagation, the mechanistic basis for using the stored energy as a measure of the driving force for crack growth suggests a crystallographic path. Indeed, in titanium and its alloys, there is considerable evidence of crack propagation along slip planes which are favourably oriented for slip [3]. However, at longer crack lengths, crack growth normal to the maximum principal stress direction gives a reasonable prediction of the crack path. Even in the microstructural regime, the normal to the maximum principal stress direction can also give reasonable agreement, with some specific exceptions [4].
In this work, the mechanistic basis of determining crack growth direction in which growth rate is controlled by local stored energy density is studied. Single crystal crack growth is considered first, in order to address the important issue of crack propagation direction with respect to crystallographic orientation. Bicrystal models are then considered to allow the behaviour of cracks growing at grain boundaries to be studied and explained. Finally, we assess the capability of the mechanistic understanding obtained by studying direct comparisons between simulated crack paths and experimentally observed cracks at the microstructural length scale.
To begin, we briefly review some of the important existing mechanistic models for fatigue crack growth, and modelling methods, providing the context for the work presented in this paper. The present work does not address chemical effects such as oxidation, nor are materials containing significant numbers of voids considered. Clearly these will be important in some cases, but the focus here is on crack growth arising from internal processes under mechanical loading.

PROPAGATION MECHANISMS FOR FATIGUE CRACK GROWTH
Cracks are strongly microstructurally-sensitive when the length of the crack is below ~20 grains. In this regime, the plastic zone is large and encompasses the entire area within which the processes leading to crack propagation occur. It seems reasonable, therefore, that processes such as slip are considered in the model. One result of the large plastic zone is that Paris' law based models are not appropriate for microstructurally short cracks, since the stress-intensity factor is not a reliable measure of the conditions at the crack-tip (which is clear from the considerable scatter present in da/dN vs ΔK plots). A satisfactory mechanistic model is therefore one which is able to bridge between the nucleation of the crack through to the Paris-regime (i.e. Stage I and stage II crack growth), ideally being able to capture both. Some of the most important aspects of current mechanistic understanding are assessed here, but a more in-depth review of the literature has been given by Chowdhury and Sehitoglu [5].
Dislocation-based models, such as those of Neumann [6] have long been used to explain the facets associated with stage II cracks. From this basis of cracks being formed from the alignment of dislocations, several dislocation-level models for fatigue crack growth (FCG) have emerged. Two key ideas have driven progress in this field. The first is that fatigue crack propagation is driven by the nucleation of dislocations at the crack tip, which then contribute to crack growth. The second is that FCG is not governed by the emission, but by slip irreversibility associated with these dislocations.
An early model for crack propagation-based dislocation nucleation was that of Gerberich [7]. Here, the key criterion was the stress required to nucleate a dislocation, which then extends the crack by one burgers vector. Slip irreversibility models are based on the idea that fatigue crack growth occurs through strain-localisation followed by separation. Wilkinson et al [8] showed crack advancement through the accumulation of vacancy dipoles at the crack tip. Related to the slip irreversibility, the energy stored in dislocation structures has recently been shown to influence the nucleation of fatigue cracks [9]. Since nucleation involves the growth up to a measurable length, the same buildup of energy under slip irreversibility may also be responsible for the propagation of the crack at the microstructural length scale. Recent modelling work has confirmed that this gives a microstructurally-sensitive measure of the driving force for crack growth [2]. It is thus proposed that it is the energy associated with these dislocation structures, built up over the loading history, which drives fatigue crack growth. As the crack propagates, the dislocation structures are annihilated, releasing the energy. This is supported by the fact that the energies typically developed ahead of the crack are of the same order as the surface energy, and that the majority of the energy is stored in geometrically necessary dislocations supporting local lattice curvature.
With regards to the direction of growth, dislocation based models implicitly suggest crystallographic propagation direction. There is also considerable experimental evidence showing that this is the case. Work by Herbig et al. [1] has shown than even when at a higher length scale the path does not appear to follow a single slip plane and some cracks grow along alternating slip planes. This allows an overall crack growth direction which does not correspond to a single crystallographic plane when observed at a high length scale but is very clearly created through crystallographic growth. As pointed out by Herbig et al., this could be explained by the double-slip mechanism proposed by Neumann [6]. A further behaviour observed was a coupling between the rate and the apparent growth direction of the crack -cracks growing parallel to slips planes had higher growth rates compared with those which changed planes. Capturing the behaviour observed by Herbig et al. may be the key in developing a model covering both 'crystallographic' and 'principal stress' directed fatigue crack growth.
Another aspect that a microstructurally-sensitive crack growth model must address is that of the grain boundary behaviour. Fatigue crack growth can be retarded at grain boundaries or experience little change. Navarro et al [10] were able to capture this in an analytical model. By considering slipblockage at grain boundaries, which leads to retardation of growth until the stress rises to drive slip in the following grain, they defined the rate of growth as a function of the distance from the cracktip to the next grain boundary. This model was able to capture the decreasing variation of the growth rate as the crack transitions into the long crack regime. This is an important behaviour to capture, but to be more predictive it is desirable to use a more fundamental quantity governing the rate rather than one which is phenomenological such as the distance to the following grain boundary.
The key addition to literature in terms of the directional behaviour at grain boundaries is the Zhai model [11], which gives a somewhat heuristic approach to explaining the crack path and retardation at grain boundaries based on the tilt and twist angles between the crack paths in the two grains. There are cases where this relationship between the tilt and twist angles with the growth rate holds, for example the finding that increased twist character inhibited crack transmission by Ludwig et al [12]. However, it has been noted in recent literature that there are a number of cases where this model seems to break down. For example Rovinelli et al assessed several slip transmission criteria in a BCC material that experienced pencil glide and found that all those based on the relative slip plane orientations were poor predictors of the crack transmission [13]. Another example has been given by Zhang et al [3], who found that transmission with a poor alignment between the incoming and outgoing planes was possible provided the slip system in the outgoing grain is favourably oriented for slip. Given that the most favourable slip across a grain boundary will often be aligned, and the favourability of the outgoing plane for slip is a better predictor of the direction, it seems likely that the most active slip system may be the best determinant of the growth direction of the crack. We also note that there are some cases where crack deflection occurs within a grain, for example in the ferritic steel results of Wan et al [4]. It is therefore desirable to develop a model that is able to capture the wider variety of behaviours both at grain boundaries and away from them.
Underlying all of the above observations, there are some significant complications to the understanding and modelling of experimental data. One is that some of the changes in direction experimentally observed with surface measurements can be caused by sub-surface features such as grain boundaries [14], [15]. Another, fundamental to the Zhai model, is that the 3D nature of the cracks is important. Hence the interpretation of experiments must therefore be considered in the context of there being a considerable number of unknowns which could affect the observations. An example of more fundamental experimental observation of fatigue crack growth is that of Mine et al [16], who looked at fatigue crack growth in single crystal commercially pure titanium. In this work based on CT specimens, crack growth was found to be heavily dependent on the crystal orientation. For specimens with the c-axis oriented out of the plane of the sample the crack was found to grow alternating between two prismatic planes, although this was complicated by crossslipping on the <a> pyramidal planes. The oscillations occurred at a scale that made them largely imperceptible when viewing the crack at full scale. The crack growing in crystals with the c-axis in plane tended to grow closer to the mode-1 growth direction, but again this was complicated by twinning. An issue in applying the results from this work to microstructurally short crack growth is that the initial notch was 25% of the thickness of the CT specimen, which may make the stress state at the crack tip different from that experienced in microstructurally short cracks due to the large influence of the initial notch.
The mechanisms discussed so far have focused on the effect of dislocations. A dislocation scale model would therefore likely be best able to capture the full behaviour, but the high computational overhead makes analysing larger models over many fatigue cycles intractable. Hence the focus here is on developing grain-scale models, for which recent efforts are assessed next.

GRAIN-SCALE MODELLING
There is a huge number of models which predict the rate of crack growth as a function of a multitude of parameters. We explicitly focus here on models that are intended to capture the fatigue crack below the Paris' regime. Two approaches which have been used are stress-based and (plastic) strain based respectively. In addition to these, we introduce the dislocation stored energy criterion used in this work.
The use of stress as a key criterion for crack growth in fracture mechanics has led to its use in a considerable number of studies for fatigue crack growth. An example of this is work by Potirniche et al. [17], which showed that the driving force for crack growth (measured using both CTOD and the ratio between the opening stress and the peak stress at the crack-tip) could both be increased or decreased as the crack approaches a grain boundary, depending on the crystal orientations. They noted that the slip resistance of the second grain affected both the stresses and the CTOD. This was also found by Ferrie et al [18].
Forest et al [19], [20] defined a damage indicator based upon the resolved shear stress, normal stress and the slip on each slip system, with the direction leading to the greatest damage point being used for propagation. Single crystal analyses using this model gave a crack path similar to the crystallographic direction in the case of single slip, and an alternating path in dual slip systems, reminiscent of that seen experimentally. Extension into full 3D microstructures has shown 'short crack behaviour', but has yet to be definitively linked with experimental results [21]. It is pointed out here that the microstructure-scale modelling is able to compute many of the quantities relevant to FCG, such as resolved shear stresses and accumulated slip -which could aid the understanding of what is driving the growth.
The strength of microstructure-scale modelling has been demonstrated in a study by Rovinelli et al [22], in which a Bayesian network was applied to investigate the correlations between a large number of variables and the observed rates and directions of growth. This demonstrated that, in the context of a BCC material which exhibits pencil-glide mechanics, the rate of growth was most strongly correlated with the crack length, as might be expected, followed by the maximum accumulated plastic shear along a given slip direction mis-alignment of the primary and secondary principal stresses with the slip direction. Meanwhile, the direction of growth was shown to be correlated most strongly with the mis-alignment of the maximum principal stress (showing a maximum around 60°) and the resolved shear stress. Analytical models for crack growth based on these insights, which were generally stress based, showed significant improvement over other fatigue metrics.
Generally, the stress is essential for determining the conditions at the crack tip for a given sample/model, but the fundamental problem for its use as a measure of the driving force for microstructurally-sensitive fatigue crack growth is that it is not able to capture the level of microstructural sensitivity that is observed [2]. Plastic strain accumulation was utilised by Farukh et al [23] who were able to show fluctuating rates as the crack grew through successive grains. However, the fluctuations were not as significant as those seen experimentally. Crack arrest at grain boundaries, a commonly observed behaviour in short crack experiments, was not predicted by the model.
The fundamental quantity that leads to the creation of new surfaces is energy. Under fatigue conditions, a portion of the plastic strain energy gradually builds up ahead of the crack tip and is stored in the form of dislocation structures. Of this stored energy, only that energy which is within a critical length scale of the new crack surface is released when the crack extends. The stored energy, G, is therefore calculated from: where the proportion of plastic strain energy that is stored in dislocation structures is given by (taken here to be ~5%) and the critical length scale over which the energy contributes to crack extension is given by the mean free path of dislocations which can be calculated from the dislocation density using = 1 √ ⁄ . This stored energy provides a physically-based measure of the driving force for crack growth through the annihilation of dislocation structures leading to energy release. This has been shown to be sensitive to the microstructure. A more detailed discussion of the theory and implementation is given in [2].
In stage II, the crack growth direction tends to be normal to the maximum principal stress. This approach has been assessed in previous work studying stage I FCG, and shows decidedly mixed agreement with experiments at the microstructural scale [4]. Although this gave some agreement in some cases, it is clear that the stress state in an anisotropic elastic model is not sufficient to predict the crack path with good accuracy. It is therefore apparent that the mechanism determining crack path is much more than that resulting simply from the elastically anisotropic stress state; while the stress is important, there are other contributors which derive from the meso-scale dislocation activity, which may be accessible from crystal plasticity. The mechanistic criteria therefore warrant more attention.
It is well established that fatigue cracks often have at least facets which are aligned with crystallographic planes, and many materials exhibit a significant amount of growth along slip planes. This also occurs during quasi-static fracture loading [24]. This has informed the methods in some previous work such as that of Lin et al. [8]. In this work the crack was sequentially grown along the slip trace corresponding to the highest accumulated slip through each grain, to the next grain  [3] boundary. By updating the crack and remeshing at each grain boundary, they were able to get a prediction of the crack path, which was found to be generally perpendicular to the external load. This has the limitation that the transgranular path is linear, and so does not capture more complex changes in rate and path around grain boundaries. For example, from the experimental results in Figure 1 [3], it is apparent that crack deflection within grains can be significant.
There remains the absence of linkage between the mechanistic understanding, which is confined to the dislocation scale, and FCG models, which are generally at the continuum scale, and which are generally focused on matching experimental growth rate results, forcing lifing models to rely on statistical methods for incubation up until the crack enters the Paris regime. Previous work has shown selective agreement with experimental results using either crystallographic or maximum principal stress based criteria. There is therefore a need to develop a model which is informed by the mechanistic understanding at the dislocation scale, but which is capable of explicitly modelling crack growth from nucleation through to the long crack stage where conventional growth models are applicable. In this work we utilise the stored dislocation structure energy criterion as the driver for crack growth to investigate FCG rates and growth directions. In previous work, a dislocation-based crystal plasticity formulation for an HCP zirconium alloy was utilised together with the dislocation structure stored energy in order to address microstructurallysensitive drivers for crack growth. The link through to the macroscale J-integral was demonstrated, and the key contribution was the establishment of a microstructure-sensitive driver for crack growth which tended to the J-integral for cracks of a length beyond the microstructurally sensitive. This work did not address the growth of microstructurally-sensitive cracks, nor the path tortuosity which are considered in the present paper. Hence, a brief summary is given of the key methods utilised, but readers are referred to the earlier paper [2] for full details.  We again consider zirconium as a representative HCP material in the following two sections, but argue that the mechanistic approach is general in nature and applicable across other crystalline metals. The elastic constants used here (all in GPa) are E1 = 98.3, E3 = 123.3, G12 = 35.1 and G13 = 32.0. The slip rule used which has been shown to accurately represent the deformation of various materials including Zirconium is: where is the density of mobile dislocations (which is assumed to be constant), the frequency of attempts of dislocations to jump obstacles, the Burgers vector, the thermal activation energy, the Boltzman constant, and the resolved shear stress and critical resolved shear stress on slip system respectively, and is the activation volume. is the temperature (295K). The properties/parameters used here for pure zirconium are listed in Table 1 and the critical resolved shear stresses for each slip system in Table 2 are taken from micro-cantilever tests by Gong et al [25]. Dislocation hardening on these systems due to the evolving SSD and GND density is given by: where SSD and GND are the sessile densities of statistically stored (SSDs) and geometrically necessary (GNDs) dislocations respectively. Determination of the GND density (reported in detail elsewhere, e.g. [26]) relies on the relationship with the curl of deformation gradient : where is the number of slip systems (denoted ), and and are the cumulative Burgers vector and GND density respectively on slip system . The curl of the deformation gradient is calculated in this model geometrically over each element. The sessile SSD density is taken here to be constant at 1 × 10 −2 µ −2 . The choice of a constant SSD density is done on the basis of previous experimental work [25], which demonstrated an absence of significant hardening in the single crystal

{0001}
P a g e | 9 response. In order to demonstrate the capability of the model to capture macroscale polycrystal behaviour, a cuboidal representative volume element was created with 9x9x9 grains, each made up of 5x5x5 elements as shown in Figure 2. Grain orientations were assigned in order to represent the experimentally measured texture. Figure 2 shows that good representation of the experimental data is achieved.
The microstructurally-sensitive driver for crack growth introduced in [2] is calculated at the crystal plasticity length scale but its mechanistic basis is that the energy for crack growth originates from that released by dislocation structures in order to generate new free surfaces. In particular, the GND structures established at the crack tip, and which carry the local lattice curvature, dominate (over SSDs) and the majority of the energy released is argued to derive from the release of the local curvature. Discrete dislocation analyses for edge cracks undergoing cyclic loading support this argument and demonstrate that the stored energy density in equation (1) calculated at the dislocation-based crystal plasticity level (provided geometrically necessary dislocations are incorporated) is a good quantitative measure of the true dislocation configurational energy [2]. The local stored energy density introduced above in equation (1) is determined within the present crystal plasticity formulation through an integral calculated pointwise throughout the entire deformation history from Hence its value at any particular loading cycle, or within any individual cycle, is known. G is a local finite element quantity and thereby suffers potentially from problems of mesh refinement and dependence. These effects are minimised using a non-local methodology, which calculates a weighted average of the stored energy with a radius of 30 µm (which is around 3 element widths). The implementation here is done utilising the approach of Korsunsky et al [27], and full details of the non-local formulation utilised are provided in [2].
The eXtended Finite Element Method (XFEM) has been utilised in the present work to introduce initial cracks into the single crystal, bicrystal, and polycrystal models to allow crack propagation according to the mechanistic rate and direction criteria introduced. This has been done using the XFEM implementation in Abaqus software with the propagating crack formulation which does not include crack tip enrichment. The latter effect has been shown to be very small in previous work [2] by comparison with analytical calculations. The propagation criterion is checked in the element ahead of the crack tip at every increment in the calculation (typically on the order of 10 -6 s during cracking). Since the loading cycles in this work are of duration 20 s each, there is essentially no applied restriction in terms of crack growth in a given cycle. Further algorithmic details are provided in Appendix 1.

CRACKS GROWING IN SINGLE CRYSTALS
Initially we consider a crack growing in a single crystal under uniaxial, cyclic, R=0, stress controlled loading (far field stress 200 MPa -which corresponds to a stress intensity factor of 7.73 MPa√m) as shown schematically in Figure 3. The crystal orientation was assigned to be one of 7 orientations (A-G, shown in Figure 4). These were chosen to be representative of the conditions under which stable FCG occurs. To speed the calculation, allowing higher numbers of fatigue cycles, only the central 750µm was modelled with the full CPFE. The outer region was modelled using an isotropic elastic plastic model (elastic properties identical to the crystal in orientation A, first yield at 425MPa, with linear hardening to 500MPa at 1.6% strain) chosen to replicate the bulk yield response shown in Figure 2. The CPFE region is sufficiently large such that the isotropic elastic plastic region is essentially fully elastic throughout the growth. The elastic discontinuity is both small enough and far enough away that it has little effect on the conditions at the crack tip. The crack was introduced using XFEM as a 270 µm edge crack as shown in Figure 3.
The cracks were grown within the crystal plasticity finite element approach by one element size when the stored energy ahead of the crack tip reached a critical value. The direction of growth was chosen to be either along the slip plane with the most slip ahead of the crack tip (the crystallographic model) or along the direction normal to maximum principal stress (the MPS model).
For the slip plane calculation, the shear strain on each slip system was integrated over the course of the simulation, and an average value calculated across all integration points within each element. When the stored energy reached the critical value in a given element, the crack was extended through the element, parallel to the slip system with the highest averaged total shear strain. In the MPS model, the maximum principal stress direction is calculated at each integration point within the element, and again averaged across the element. The crack is then grown in the plane normal to this. In order to get sufficient resolution and mesh refinement in the growth region, the model was a single element in thickness (30 µm). Conditions of plane stress applied, as on a free sample surface. Due to the absence of resolution in the Z direction (out of plane) the crack front was constrained to remain perpendicular to the x-y plane for both crystallographic and principal stress models.
The critical stored energy was chosen to be 0.75 Jm -2 , which is lower than the critical value for fracture identified in previous work of 7.2 Jm -2 [2]. This was chosen such that the number of cycles to failure (i.e. crack growth) was low enough to be computationally tractable, but long enough to give sufficient differentiation of rates across the different crystal orientations and, for later sections, at grain boundaries. This may affect the stored energy field ahead of the crack tip compared to the case where many more cycles are needed when a higher critical stored energy is employed for crack extension. However, it is noted in previous studies that the field changes most significantly in the first loading cycle [28], and then only gradually in subsequent cycles. Hence the mechanistic drivers for crack propagation are unlikely to be significantly changed by virtue of the low critical stored energy utilised since the differences in the key fields (slip, stress and GND density) are not expected to be great.
Field plots of the stored energy, plastic strain and GND density at the point of the start of crack growth are shown in Figure 4 for each of the crystal orientations modelled. For the purposes of this study, a range of crystal orientations has been chosen that represent extremes in terms of the c-axis orientation, and also to cover orientations for single and double slip, but it is not, of course, exhaustive. We consider first the cycling before propagation occurs, focusing on the build-up of the stored energy, and this is followed by assessment of the growth paths and rates predicted by the two growth criteria.

PRE-PROPAGATION CYCLING
Field plots showing the distributions of stored energy, plastic strain and GND density at the onset of crack growth are shown in Figure 4b-d. There are several important points to draw from this. Firstly, the stored energy distributions are very different, being much more concentrated for orientation C, and much more diffuse for orientations A and D. This has important implication for the crack growth rates, in that the distribution will cause an increase in the rate if the crack grows along the directions of high stored energy. This will be discussed in more detail later.
The build-up of stored energy in the element at the crack tip is shown for each crystal orientation in Figure 5, along with the number of cycles until propagation starts. The radius of the non-local averaging is 30 m, vs a crack length of 270 m and so, although not a fully local approach, the controlling quantity is at the crack tip, and much of the distributed stored energy build-up that can be seen in Figure 4 does not contribute to this. Since the critical value of the stored energy has been chosen here to give stable crack growth within the small number of cycles that are possible computationally, of interest is not the specific number of cycles, but rather the comparative differences between the different crystals. The stored energy develops past the growth threshold within the first cycle in crystal E, the third cycle for F, and only after 12 cycles for crystal G. The high rate of stored energy development in crystal E is due to the c-axis alignment with the crack direction -note that this is also observed to be the weakest orientation in experimental fracture toughness measurements of highly textured titanium, such as those of Tchorzewski et al [29], which show that this orientation gives the lowest fracture toughness. Another interesting point is that the higher GND densities generated in the softer orientations A-C cause the stored energy to be reduced relative to that seen in orientation E, despite having a similar level of plastic slip.
Finally, it can be seen that the distribution for orientation D is highly non-symmetrical. This is due to the c-axis orientation being non-perpendicular with the loading direction (orientation F does not exhibit strong asymmetry since the c-axis is twisted out of the plane of the model).

GROWTH PATHS
For the growth paths we consider first orientations A-C, which are shown in Figure 6. These orientations share a common c-axis direction, which is out of plane of the model, and so would give an identical elastic response in the absence of slip. Looking first at the MPS model paths, it can be seen that the cracks in orientations A and B were generally straight and perpendicular to the remote loading, whereas the crack in orientation C grew at an angle. The differences here are due to orientations A and B giving symmetric slip under mode I loading whereas Orientation C is orientated for single slip. This difference means that the maximum principal stress direction at the crack tip is rotated in crystal C.
For cracks growing in single crystal orientations A-C under the crystallographic model, the cracks all grew along prismatic systems. In crystal C the crack grew along the slip system that has the highest Schmid factor, as might be expected. More striking is the difference between orientations A and B.
In these orientations, in addition to the identical elastic response, the plastic response is also similar, showing the distributions of plastic strain and stored energy in Figure 4. Despite these similarities, the crystallographically-driven growth of the crack showed significant differences in path. The crack remained on the same plane in the A orientation but repeatedly switched planes in the B orientation. Note that since the model is symmetric, the initial growth along either prismatic direction (i.e. to the right or to the left) is arbitrary, but the subsequent growth paths, and their differences, are entirely deterministic and occur for the following reasons.
The contour plots in Figure 7 show different regions around the crack coloured according to the slip system on which there is the greatest amount of slip (shear strain). These are the slip systems along which the crack would propagate if it passed through these regions. Before growth begins (see the cycle 7 plots in Figure 7), the lobes were similar in shape and were dominated by two prismatic planes in both crystal orientations. However, comparing the two orientations, the planes corresponding to each region are oriented at 90° to each other. This important difference gave rise to the striking differences observed in the crack paths seen in Figure 6.
In the A orientation the crack grew along the slip system indicated by the blue colour. This plane is orientated such that, after extending the crack parallel to the slip plane, the crack tip remained in a region dominated by the same slip system. This can be seen in Figure 7, where slip in the region around the crack tip is dominated by the same slip system (indicated by the blue colour) at cycles 7, 25 and 31. Since this was always the case through the entire simulation, the crack did not change direction.
In the B orientation, the most active slip plane in each region is oriented at 90° to that in crystal A. Initially (from cycles 7-25), the crack grew along the slip plane indicated by the yellow colour.
Extending the crack along this direction brought the crack tip away from the region in which this slip system was most active (which is located to the right of the crack tip in Figure 7). Eventually, at cycle 25, the crack tip grew into the region where the slip system indicated by the blue colour was dominant. These deflections occured twice in the simulation, at cycles 25 and 31. This gave rise to the highly tortuous (zig-zag) path, shown in Figure 6, that is reminiscent of those seen in similar experimental tests in single crystal Nickel samples [30]. Note that there is likely to be some mesh dependence in terms of the precise location at which these deflections occur -the key point, however, is the qualitative difference of a tortuous path versus a straight path resulting simply from crystallographic orientation.
The tortuosity therefore seems to be fundamentally driven by two features: first is the spatial variation of the slip fields for each individual slip system, and second is the directionality of the crack growth along the slip system, which takes the crack tip into an area which will have different relative activities of the different slip systems. The prediction of these effects is non-trivial however, since spatial variation of the slip fields will be dependent on a number of factors including the loading history, crack loading mode, and of course the crystallography.  The lobes of high activity of distinct slip systems have been noted before by Potirniche [31] and are relatively straightforward to rationalise using analytical stress functions in the case of the A, B and C orientations (since the c-axis is out of plane). Using analytical functions for the stress state around a crack-tip (Equations 6-8), in an isotropic elastic plate, the activity of the three prismatic slip systems can be estimated.
Using this analytical model, the plastic zone size for a given slip system can be estimated by assessing the area ahead of the crack tip for which the resolved shear stress on a particular slip system (calculated using the full stress state given by equations 6-8) exceeds the critical resolved shear stress. The plastic zone sizes predicted for each of the prismatic slip systems using the analytical model are shown in Figure 8, normalised against the largest size across orientations considered. The slip systems predicted to be most active are the same as those predicted in the CPFE model. Furthermore, it is clear that in the C orientation there is a single dominating slip system. The analytical equations assume pure mode I loading and so the fact that the correct slip systems are predicted means that the stress state is not being affected sufficiently by the differing orientations to completely change their identity, but clearly the stress state does change sufficiently to favour one or other of the similarly active systems in orientations A and B.
Moving on to look at the results for crystal D orientation shown in Figure 9, this was the orientation for which the maximum principal stress directed growth exhibited the most significant deviation from linearity and also normality to the remote principal stress direction. This occurred due to the caxis being oriented neither parallel nor perpendicular to the loading direction, which caused the maximum principal stress direction at the crack tip to be rotated compared with the far-field load direction by 29.8°. This change in the stress state also affected the crystallographic growth model, for which the crack did not grow along the basal plane (Labelled D1 in Figure 9a) which had the highest Schmid factor of 0.5 but instead grew along one of two <a> pyramidal systems (Labelled D2 and D3 in Figure 9a) which had identical Schmid factors of 0.31 (the planes can be seen in Figure 9a). The high level of slip on the pyramidal system compared with the basal system can be seen in Figure  9b, along with the predicted crystallographic crack growth path in Figure 9c. These slip planes are collinear in the plane of the model and so the crack progresses along these systems exclusively, before switching onto a <c+a> plane. This occurs because as the crack gets longer the stress state at the crack tip rotates to become more aligned with the far-field loading direction (Figure 9d), as also occurred in the principal stress model, and so the resolved shear stress on the <c+a> plane (shown in Figure 9e) increases at longer crack lengths. The E, F and G crystal orientations give crack paths that were very similar under both crystallographic and principal stress crack growth (although the later crack growth in E was tortuous), being on average perpendicular to the remote loading direction ( Figure 10). The most active slip systems, shown in Figure 11 for each crystal orientation, were different for each orientation and are twisted out of plane in 3D. The growth in crystal orientation E was along prismatic planes, F was dominated by basal slip and G by <c+a> pyramidal slip whose plane is a simple rotation around the axis perpendicular to the loading direction in the plane of the model. Growth along the <c+a> plane in orientation G has been observed experimentally in titanium by Bantounas et al [32].
Comparing the modelling results obtained here with the microstructurally sensitive crack growth seen experimentally, e.g. that shown in Figure 1, it is clear that the crystallographic growth is better able to capture the range of behaviours seen experimentally, with some cracks growing straight along slip planes and others growing along multiple planes. It is also clear that, under the conditions  prescribed in the crystallographic growth model, a small change in the crystal orientation can have a large effect on the direction and type of fatigue crack growth. A good example of this was discussed earlier in which crack growth along prism planes can give rise to completely different crack paths, despite there being no change in the in-plane elasticity properties, resulting from crystallographic growth with respect to the local spatial stress state. The effects of changing the crystallographic orientation are not limited only to the path however and are also potentially important for rate of crack growth, which is addressed next.
Generally, the crack paths in the MPS model were in all cases essentially linear and perpendicular to the overall loading. This is expected since, with the exception of orientation D, the crystal orientations give a symmetric stress response across the crack plane. Orientation C also showed deviation from the medial plane due to the single slip. Comparing with crystallographic growth, the two methods sometimes gave similar paths (E-G), whilst in others the paths were found to be very different (A-D). Furthermore, it can be seen that the prismatic plane crack growing parallel to the caxis (E) was initially straight and later had low amplitude oscillation, whereas the prismatic plane cracks growing perpendicular to the c-axis (A-C) either did not change planes or had higher amplitude oscillations. This increased deviation from the medial plane in those prismatic cracks growing perpendicular to the c-axis, and the associated broader slip band behind the crack front, is something that may be observed in experiments.
In the crystallographic model, growth was seen along a range of slip systems, similarly to experimental observations, and the paths only sometimes aligned with the maximum principal stress direction either through a single slip system or through the combination of multiple slip systems giving an overall path that was perpendicular to the remote loading direction. The crystallographic model is therefore predictive of the kinds of growth path behaviour observed experimentally by Herbig et al [1] with microstructurally short cracks growing parallel to slip systems in some crystal orientations and growing on alternating planes which give an 'average' direction which does not correspond to a given slip plane.
Another general observation in the crystallographic model is that, at longer lengths, the cracks started growing on a wider variety of slip planes. This suggests that as the crack gets longer, growth on alternating slip planes becomes more likely. The 'average' path is therefore less likely to align to a given crystallographic plane as the path is made up of many individual facets on different slip systems. This is again in agreement with experiments, which have shown that cracks which are long relative to the microstructure give paths that do not follow a specific crystallographic plane, but which do have individual facets on the crack surface aligned with specific planes.

GROWTH RATES
Crack extension -measured along the path of the crack -as a function of fatigue cycles for all of the single crystal analyses are shown in Figure 12a-b, for crystallographic and principal stress driven crack paths respectively. Note that the local rate of growth is controlled by the value of the critical stored energy, and recall that a relatively low value (0.75 Jm -2 ) has been selected here to ensure a small number of cycles is required in the computations to lead to significant crack growth ~500 m.
The key insight sought here is not the absolute rate of growth, but rather the comparative rates across the different crystal orientations. It is clear that the crystal orientation has a significant effect on the rate of growth, with uncontrolled growth in orientation E and very slow growth in orientation G under both crack drivers. The variation in rate observed is in qualitative agreement with experimental observations in that the predicted growth rate fluctuates significantly in the microstructurally-sensitive region [33], as well as previous work which has demonstrated the microstructural sensitivity of the stored energy [2]. Furthermore the slowest growth being in orientation G has been observed in titanium [16], [32]. The results also show a strong coupling between the path and the rate of growth. The number of cycles for the crack to grow from 50 to 200 µm in each model is shown in Table 3.
In the cases where the crystallographic and maximum principal stress directed paths aligned (orientations E and F), the number of cycles to extend the crack were identical, whereas when the paths diverged from each other, the rates were significantly affected (differing by up to 1500%).
There are two key factors that affect the rate of growth here: the rate of development of the stored energy at the crack tip and the shape of the slip field. The shapes of the stored energy fields were shown previously in Figure 4. Considering the directions of growth, in the A and B orientations, the crack grows along the lobe of high stored energy in the crystallographic model, but into a region of low stored energy in the MPS model. As such, the crack propagation for both orientations was faster in the crystallographic model. This is shown diagrammatically in Figure 13. The crack grows into an area of high stored energy density, which leads to the increased growth rate relative to the MPS model. The highly directional stored energy field in orientation D also gives a higher rate of growth in the crystallographic model.  The higher rate in C occurs since the maximum principal stress directed crack grows doesn't grow perpendicular to the loading direction, but instead grows at an angle, which takes it along a region of higher stored energy than if it grew perpendicular to the loading direction. Also interesting here is that the rate of growth in both models is similar for the first 100 µm -being just 4% different. Once at 100 µm however, the higher stress intensity at the crack tip in the MPS model caused a significant increase in the stored energy ahead of the crack tip and increased the rate of crack extension (the effect of which can be seen as a large area of high stored energy around the crack in Figure 6). This demonstrates the complexity captured by this crack growth model, since the stress state at the crack tip depends not just on the load and the length of the crack, but also on the path taken.
At a higher length scale, the repeated changing of direction in the crack would give an effective path that would match the principal stress path. The slower growth along the principal stress path therefore agrees qualitatively with the observation of Herbig et al [1] that the cracks that were growing along non-crystallographic directions grew more slowly than those growing along crystallographic directions.

CRACKS GROWING IN BICRYSTALS
One of the biggest factors affecting the rate of crack growth is grain boundaries. Grain boundaries affect crack paths but are also experimentally observed to affect the rate of growth, even for crack tips away from the boundary itself. Since, as noted before, the crack path and growth rate are highly coupled, in this section we consider these effects of the changing crystallography together.
The Zhai model [11] is often used to rationalise the change in plane across grain boundaries, but this has not led to the development of a predictive model for crack growth and in some cases the behaviour at the grain scale is not what would be expected [34]. It is therefore interesting to investigate whether crack plane compatibility is predicted, and why it might not always occur.
The model used for the bicrystal crack investigations is shown in Figure 14. This is the same as in the single crystal model, but with a grain boundary introduced at a distance of 80µm ahead of the initial

µm Crack
Crack tip crack tip position, with differing crystal orientations applied either side of the grain boundary. Since there are countless possibilities here, we focus on some specific examples which in turn give differing and interesting behaviour. In addition, we confine the bicrystal studies to utilise always a sub-set of two of the crystal orientations (A-H) described already.
The specific combinations of bicrystals studied, together with some key results are shown in Figure  15. For each combination of orientations the stored energy distribution is shown at the onset of growth, and then at the instant when the crack reaches the grain boundary. The final crack paths are also shown.

Bicrystals A-B / B-A
This pair of crystal orientations (the first two results columns in Figure 15) have parallel c-axes, normal to the initial edge crack, and are rotated by 30° relative to each other about the c-axis. As the elastic response of both crystal orientations is identical, the initial propagation rate and path, shown in Figure 16, was found to be very similar to that in the single crystal cases. Upon reaching the grain boundary, the cracks grew faster/slower in accordance with the single crystal results for the crystal orientations of the second crystal. In the B-A case, the crack moved on to the prismatic plane which was best aligned with the incoming plane, indicating that if there are two active planes that might be observed in a single crystal, compatibility may be the deciding factor determining the behaviour at grain boundaries. In the A-B case the crack followed the behaviour seen in single crystal B, staying closer to the centre of the sample and switching direction.

Bicrystals A-D / D-A
The results for bicrystals A-D and D-A are shown in third and fourth columns in Figure 15. In the A-D bicrystal the crack path in crystal 'A' followed the same prismatic plane seen in the single crystal model. The dominant slip system overall in crystal 'D' was again the <a> pyramidal, but as the crack approached the grain boundary the activity on the basal plane increased, and the crack initially moved on to a basal plane in crystal 'D' before switching back to an <a> pyramidal plane. The high activity of the basal plane close to the grain boundary is likely to be due to the high compatibility between this and the prismatic plane in crystal 'A'.
In bicrystal D-A, the crack grew in crystal 'D' along the same <a> pyramidal plane seen in the single crystal 'D' model and then in crystal 'A' along the same prismatic seen in the single crystal 'A' model. These two samples differed significantly in the crack growth rate. In bicrystal A-D, the crack grew across the grain boundary without slowing significantly (as can be seen from the lack of any change in gradient of the black line in Figure 17a) whereas in bicrystal D-A, the crack was significantly retarded at the grain boundary, taking a number of cycles to start growing again (grey line in Figure  17a can be seen to abruptly change gradient at the point the crack reaches the grain boundary indicating that the crack growth has slowed). Consideration of the stored energy distribution at the point the crack reaches the grain boundary (shown in Figure 17 b and c) provides the explanation for these observations. In the bicrystal A-D case, significant stored energy was developed in the second crystal, ahead of the crack tip, and so driving the crack such that there was no retardation of the growth. In bicrystal D-A, very little stored energy was developed in the second crystal before the crack reached the grain boundary. This led to the retardation in the crack growth rate at the boundary until sufficient stored energy developed with cycling in the second crystal.

Bicrystal A-G
In this bicrystal (results shown in final column of Figure 15), interestingly the crack grew through the 'A' crystal to the grain boundary much slower than in the single crystal. This is caused by the elastic constraint effect of the second crystal (G) inhibiting slip in crystal A -as can be seen in the reduced slip field observed in the bicrystal A-G vs the single crystal A in Figure 18

Figure 17 -Plot of crack extension against cycles for bicrystals A-D and D-A, with single crystal results shown for comparison. b-c) contour plot of stored energy distribution around crack at point crack tip reaches grain boundary in b) A-D and c) D-A bicrystals. Note that there is significantly more stored energy developed in the second grain in the A-D bicrystal compared with the D-A bicrystal, which causes the retardation in growth rate at the grain boundary in the latter.
loading history. Upon reaching the grain boundary at the second crystal (G), the crack propagation rate was significantly reduced. The contour plot of the stored energy at cycle 61 (the point at which the crack reached the A-G grain boundary) in Figure 18d shows that there was very little stored energy being developed ahead of the crack tip in the second grain, which lead to the significant arrest in growth.
In order to see further growth in the second crystal, the simulation here was extended to 100 fatigue cycles. Interestingly, the eventual growth in crystal G was along an <a> pyramidal plane rather than the <c + a> pyramidal plane seen in the single crystal growth. This plane is better aligned with the incoming prismatic plane in crystal A.

Crack passes grain boundary
The results above indicate that the crack growth model based on a rate determined by the dislocation stored energy density and a path determined crystallographically, captures much of the experimentally observed behaviour associated with microstructurally-sensitive fatigue crack growth. However, no direct comparisons with experimental data have so far been made.
Experimental results showing fatigue crack growth paths in Ti-6Al-4V have recently been published by Zhang et al [3]. In these experiments EBSD imaging was used to track the path taken by the crack through the surface microstructure. Two differing sections of the crack growth were highlighted by Zhang et al (shown in Figure 19) -one where the crack maintains a consistent direction through three grains, and another where the crack deflects strongly at the first grain boundary to grow along a prismatic system and then again at the following grain boundary, growing in the third grain initially along a direction which is not aligned with any one slip system, before growing onto another prismatic plane. The experimentally measured crystal orientations from [3] were adopted for two corresponding CPFE models with three successive grains, as shown in Figure 20, for each of the two sections examined experimentally. The experimental conditions of R=0 cyclic loading were also imposed in the model representation, and the remote peak stress was set at 200MPa. For these simulations CRSS values were modified such that the ratios matched those of titanium [35], but the CPFE parameters were otherwise identical to those used in the preceding sections. The critical value of the stored energy was again chosen to be 0.75 Jm -2 since the primary interest for this study is in crack path dependence on microstructure.
The crack paths predicted using these two microstructural models (A and B) are shown in Figure 20, together with the slip planes along which the simulated cracks grew. The predicted slip planes and crack growth directions are also overlaid on the EBSD images from [3]. As can be seen, there was remarkable agreement between the paths predicted using the model and those observed experimentally. In sample A, the simulated crack grew along precisely the same directions as observed in the experiment, with very little deflection at the grain boundaries. In sample B the simulated crack grew along the experimentally observed direction in the first grain, correctly deflected at the first grain boundary onto another prismatic plane (and not on to the basal plane with the higher Schmid factor), but in the third grain only grew along the slip plane which the crack moved on to after initially growing along a non-slip system direction.
The intial crack growth in grain B3 is aligned with what appears in the EBSD image to be a scratch on the surface, and this potentially explains the differences observed, which only arise over the length of the apparent 'scratch'. An entirely new microstructurally-sensitive and mechanistic fatigue crack growth methodology has been established in which the rate of propagation is determined by stored dislocation structure energy, and the crack path crystallographically such that the direction is along that slip plane with the highest slip. The latter has also been compared with another formulation in which the crack path is that given by the normal to the local maximum principal stress direction. The stored energy formulation for growth rate has been developed at a crystal plasticity length scale in order to replicate the dislocation configurational energy determined using discrete dislocation analysis [2].

µm
The crystallographic growth model has been shown to be able to capture a wide range of behaviours in single crystals and bicrystals which are representative of the range of behaviours observed in experimental studies. Crucially, the mechanistic model was also able to differentiate and correctly predict the crack paths observed in Ti polycrystal studies which examined two very differing crack growth scenarios. A common comparison made in previous studies has been the Schmid factor calculated from the overall loading direction which has often been used to compare the crystallographic crack planes. However, this does not always give the correct crack plane because the stress state is not uniaxial near the crack tip and so other planes may in fact be more active. Furthermore, there are often multiple planes with similar Schmid factors which cannot be differentiated using the global Schmid factor alone. The new mechanistic model presented captures all of these features appropriately.
It has been suggested in previous work that, at grain boundaries, good compatibility between slip systems (i.e. when there is a low angle between the planes) may lead to one system being selected over another for crack growth [3]. However, in this work the key factor appears to be the ease of activation. In the bicrystal studies, the A-B, A-D and D-A models all showed cracks growing across grain boundaries with very poor compatibility between the incoming and outgoing planes. The fact that even without the presence of the grain boundary the crack plane was observed to switch between planes with poor compatibility (e.g. single crystal B) also supports this idea. These observations lead to the conclusion that compatible systems may be often observed simply because two well-aligned slip systems are likely to be activated together, rather than being due to any interaction between the slip systems. In addition to this, the change to mixed mode I and mode II loading when the crack propagates at an angle may also make compatible planes more likely to be observed.
To be useful in developing the understanding of material resistance to fatigue, capturing the rate of growth is perhaps the most important factor. Using the stored energy to predict the propagation gave a strong dependence on the crystal orientation and the path. The retardation of cracks at grain boundaries was seen in the model as well as the observation that some cracks pass through grain boundaries with little to no retardation. The key factor driving this was the ability to drive slip in the outgoing grain, which was strongly affected by the relative orientations of the two crystals. In the bicrystal examples considered, instances of direct crack propagation across a grain boundary occurred, where grain orientations were well-aligned for slip. In addition, when grain orientations were badly aligned for slip, significant retardation of crack growth rate was predicted, in keeping with some experimental observations. Significant changes to crack paths were observed at grain boundaries, depending on crystallographic orientation.
Despite a significant amount of research recently indicating that the subsurface microstructure will have a significant effect on the crack path, the experimental results were accurately reproduced without any consideration of the subsurface. It therefore seems that, in this particular case, the path of the crack at the surface is largely independent and the crack growth is led by surface features. Despite this, the subsurface features would likely have a more significant effect on the rate of growth since it would have the effect of holding the surface together even after the surface has cracked.
The unpredicted growth along a defect observed in the experiment highlights the importance of this in building predictive capability. The models depend on the surface being smooth, which means that a small scratch can cause a significant deviation between the model and the experimental results. When comparing with component lifetimes, the presence of a significant number scratches along the surface could mean that the majority of cracks grow along these defects at the surface rather than along the crystallographic planes. This could potentially cause a significant change in the fatigue lifetime.
In terms of accurately capturing the rate of growth, behaviours observed in experiments including growth promotion, retardation and arrest were all predicted by the model, but have not been quantitatively linked with experimental data here. More precise experimental data showing crack growth through the microstructure is still needed in order to be able to give a direct, quantitative comparison.

CONCLUSIONS
1. Modelling microstructurally-sensitive fatigue cracks by extending along the direction of the most active slip systems produces crack paths which are representative of the paths seen in the experimental literature, capturing behaviours such as cross-slip. By comparison, the paths predicted by extending the crack normal to the maximum principal stress direction have much less dependence on the microstructure than is observed experimentally. 2. Using the stored energy as the criterion for crack growth, the model can capture variations in growth rate due to crystal orientation. The rate and path are highly coupled and so both must be accurately predicted for a model to be able to assess the rate of microstructurally sensitive crack growth -models which do not consider both will fail to capture the full effect of the microstructure. 3. Crack deflection and growth promotion, retardation and arrest, behaviours which have been observed experimentally in many metals, were all predicted at different types of grain boundary. These effects were primarily governed by the ability to develop slip and thus stored energy in the outgoing grain. 4. The complete model was able to accurately capture the crack planes and the deflections at grain boundaries observed in the two experimental cases studied.

ACKNOWLEDGEMENTS
FPED and DW gratefully acknowledge support from the HexMat EPSRC programme grant EP/K034332 , and FPED his Royal Academy of Engineering/Rolls-Royce research chair support. Both DW and FPED wish to thank Mike Martin, Rolls-Royce, for financial support and for many useful discussions and insight.