Realistic microstructure-based modelling of cyclic deformation and crack growth using crystal plasticity

Using crystal plasticity, ﬁnite element analyses were carried out to model cyclic deformation for a low solvus high refractory (LSHR) nickel superalloy at elevated temperature. The analyses were implemented using a representative volume element (RVE), consisting of realistic microstructure obtained from SEM images of the material. Monotonic, stress-relaxation and cyclic test data at 725 (cid:1) C were used to determine the model parameters from a ﬁtting process and their sensitivity to RVE size and random grain orienta- tion. In combination with extended ﬁnite element method (XFEM), the crystal plasticity model was further applied to predict surface crack growth, for which accumulated plastic strain was used as a fracture criterion. Again, realistic microstructure, taken from the cracking site on the surface of a plain fatigue specimen, was used to create the ﬁnite element model for crack growth analyses. The prediction was con- ducted for a pseudo-3D geometrical model, resembling the plane stress condition at specimen surface. The loading level at the cracking site was determined from a viscoplasticity ﬁnite element analysis of the fatigue specimen. The proposed model is capable of predicting the variation in growth rate in grains with different orientations. (cid:3) The Authors.


Introduction
Polycrystalline metals and alloys possess a heterogeneous nature of grain microstructure and exhibit anisotropic response during mechanical deformation. These crystallographic orientation effects cannot be captured by traditional isotropic homogeneous models. Therefore, physically-based crystal plasticity models were developed and have been successfully applied in literature to predict the anisotropic mechanical response of individual grains in polycrystalline materials. Combination of crystal plasticity and finite element method has the ability to model the global and local stress-strain response of crystalline materials under various types of loading regimes including fatigue [1][2][3].
In the efforts to consider the effects of the material's microstructure on the mechanical behaviour of nickel based superalloy, a large number of studies were conducted using crystal plasticity [4][5][6]. These studies aimed to predict the global and local stress-strain response as well as grain texture evolution and micro-structural crack nucleation under various loading conditions. For example, Fedelich [7] employed a crystal plasticity model with implicit dependency on precipitate size and volume fraction to predict the influence of microstructural parameters such as lattice misfit on deformation behaviour of a single crystal nickel based superalloy. The model was also used to calculate the back stress as a function of local deformation state, i.e., distribution of plastic strain in the microstructure. Kumar et al. [8] adopted crystal plasticity-based finite element approach to study the effects of microstructural variability on cyclic deformation of a single crystal nickel base superalloy. This study considered the effect of size and volume fraction of the c 0 precipitates on fatigue resistance of the material. Shoney et al. [9] presented a microstructure sensitive crystal plasticity model which has the capability to capture the variation in stress-strain response with changes in grain size, distribution of precipitate size and precipitate volume fraction as well as dislocation density for each slip system. In this regard, Lin et al. [10] used 2D crystal plasticity model to study the effect of grain microstructure on localised stress-strain distribution of a polycrystalline nickel based superalloy. The study was subsequently extended to 3D to simulate the global stress-strain response of the material under fatigue loading [11]. In addition, effects of grain microstructure on crack tip deformation and crack growth path were also explored using the crystal plasticity model. However, these models were limited to the use of artificial grain http microstructure generated by employing the Voronoi tessellation algorithm, and cannot capture the local stress and strain behaviour precisely, an important indicator for nucleation of fatigue cracks. In order to capture the actual localised stress-strain response, 2D model based on realistic microstructure were developed by Choi et al. [12] to predict the localised stress in polycrystalline Mg alloys using crystal plasticity finite element method. Similarly, Wang et al. [13] studied grain level heterogeneous deformation in polycrystalline a-Ti material using crystal plasticity model based on 2D realistic microstructure. However, the work is limited to simulation of local deformation, without consideration of global deformation. Also to our best knowledge, there is hardly any work devoted to prediction of fatigue crack growth in polycrystalline materials using the realistic grain microstructure and crystal plasticity model.
Fatigue failure in engineering materials has been an important subject for several decades. A number of computational techniques have been developed and used in the literature to predict fatigue crack initiation and propagation, especially the finite element method (FEM) [14]. However, the major shortcoming of FEM is that the crack surfaces need to be in line with the edges of the elements, which presents a challenge for FEM in modelling the growth of discontinuities and cracks. In order to overcome this difficulty, a novel approach, called extended finite element (XFEM), has been developed recently [15]. In this technique, the crack geometry is mesh-independent and there is no need for the re-meshing to accommodate the crack as it grows. This is achieved by incorporating the enrichment functions into standard finite element approximation space [16,17]. With an appropriate fracture criterion, crack propagation can be modelled without the introduction of a predefined crack growth path, a unique feature of the XFEM technique. In fact, the method has already been widely applied to model a variety of crack problems such as crack propagation in thinwalled structures, crack branching and dynamic crack growth [18][19][20].
A variety of fracture criteria have been proposed in the last decades to predict crack propagation, which can be classified as either stress or strain based approach. For example, Erdogan and Sih [21] proposed the maximum principal stress criterion as the crack growth driving force. They suggested that the crack would propagate in the direction perpendicular to that of maximum principal stress when the stress value at the crack tip reached a critical value. Similarly, Moes et al. [22] used the circumferential stress criterion for crack growth. The direction of crack propagation was determined by setting the shear stress equal to zero as the circumferential stress was a principal stress in the direction of crack propagation. So, this was essentially a principal stress-based criterion. Newman [23] used a strain criterion for crack advancement and the crack started to grow when the strain immediately ahead of the crack tip reached a critical strain level. McDowell and Dunne [5] used the concept of accumulated plastic strain to reveal the fatigue crack formation in Ni-based superalloys. Using a crystal plasticity model, it was shown that the site of crack nucleation observed in experiments matched the location of maximum accumulated plastic strain obtained from simulations. Recently, accumulated plastic strain has been explored as a criterion to predict fatigue crack growth. This was based on the observation of strain accumulation, from both modelling [24][25][26][27] and experiments [28,29], near a crack tip under fatigue loading conditions. Also, such a criterion has been utilised to predict fatigue crack propagation in a Ni-based superalloy at elevated temperature using the XFEM technique [30]. With the assistance of a cyclic viscoplasticity model, the predicted crack propagation path and rate were in very good agreement with experimental data [30]. However, there is very limited work in predicting the full path and rate of crack propagation by considering the explicit grain microstructure, and this can be attempted by combining the XFEM technique, a crystal plasticity model and a suitable fracture criterion such as strain accumulation.
In this paper, crystal plasticity model, in combination with XFEM, has been applied to study cyclic deformation and fatigue crack growth in a nickel-based superalloy LSHR (Low Solvus High Refractory) at high temperature. The first objective of this research was to develop and evaluate a RVE-based finite element model with the incorporation of a realistic material microstructure. The second objective of this work was to determine the parameters of a crystal plasticity constitutive model to describe the cyclic deformation behaviour of the material by using a user-defined material subroutine (UMAT) interfaced with the finite element package ABAQUS. The model parameters were calibrated from extensive finite element analyses to fit the monotonic, stress relaxation and cyclic test data. The third objective was to predict crack growth by combining the XFEM technique and the calibrated crystal plasticity UMAT, for which accumulated plastic strain was used as the fracture criterion.

Crystal plasticity model
The theoretical framework used here is based on large deformation and rate-dependent crystal plasticity theory presented in [8,31,32]. In the theory, shear flow along slip systems is solely responsible for the plastic deformation in the crystals. The shear deformation along slip planes is caused dislocation generation and motion, with the associated inelastic velocity gradient (L P ) given as: where F P is plastic deformation gradient, _ c is shear strain rate of slip system a, and unit vectors s a and m a refer to the direction of shear slip and the normal to the slip plane, respectively.
The shear strain rate ( _ c a ) for each slip system can be expressed in terms of the resolved shear stress (s a ) using an exponential function: where S a is the slip resistance, B a is the back stress, l and l 0 are shear moduli at absolute temperature and 0 K, respectively, and j is the Boltzmann constant. In Eq. (2), F 0 ,ŝ 0 , p, q and _ c 0 are material constants which need to be calibrated for proper simulation of the material behaviour. The slip resistance and the back stress are two internal state variables which are associated with the current dislocation network introduced at the slip level. Specifically, the slip system resistance (S a ) is a measure of the impedance of dislocation motion on a slip system by short-range interactions between all dislocations [33]. The slip resistance for a slip system can be described by the following equation: where S a 0 is considered as initial value of the slip resistance for a slip system, and h s and d D are the material constants linked to static and dynamic recovery terms, respectively.
The back stress is associated with dislocations bowing between obstacles, such as precipitates or pinning points, and can be expressed in terms of a hardening constant (h B ) and a dynamic recovery constant (r D ), as: The above constitutive relationship was implemented as a UMAT subroutine, based on the fully implicit (Euler backward) integration algorithm, to interface with a finite element package (ABAQUS) for modelling of cyclic deformation and crack growth in this work.

Model development
The material used in this study was powder metallurgy LSHR superalloy, provided by NASA and used for turbine blades and rotors in the hot section of gas turbine engines. The material possesses excellent high temperature tensile strength and creep performance as well as good characteristics due to low c 0 solvus temperature [34]. The material was supersolvus heat treated to yield a coarse grain microstructure, yielding a range of grain sizes, i.e., 10-140 lm, as shown in Fig. 1a. An average grain size of 36.05 ± 18.07 lm was obtained by line intercept measurements.
In order to simulate the mechanical behaviour of a polycrystalline material at grain level, a representative volume element (RVE) consisting of sufficient number of grains is required in order to represent the global material behaviour. In this paper, RVE was developed to have realistic microstructure of the material rather than the artificially constructed using Voronoi tessellation technique. Starting from the material's microstructure obtained from SEM, coordinates of grain boundaries in Fig. 1b were determined using an in-house image-processing code developed within MATLAB. In order to obtain the actual dimensions of all the grains, SEM image was processed according to the actual scale to obtain the coordinates of all grain boundaries. A python code was developed to generate the geometry using the grain boundaries coordinates obtained from MATLAB. Using the python code, those coordinates were used as input in ABAQUS CAE, which were then connected to give the FE model with realistic grain microstructure. A flowchart showing steps to generate RVE based on realistic microstructure is shown in Fig. 2. The model developed here using the aforementioned technique is based on surface images of the material and thus considers the realistic microstructure in only 2D while neglecting the heterogeneity of the material in the third dimension.
Periodic boundary conditions were applied on the model to realise the repetitive deformation of the RVE. In order to apply the periodic boundary conditions, RVE was meshed in such a way that the nodes on the opposite edges of RVE are equal in both number and space, which was performed manually. Then, multiple constraint conditions in ABAQUS were used to enforce parallel deformation of the opposite edges of RVE as shown in Fig. 1c. Rigid body motion was removed by fixing the vertex A in both x and y directions as well as constraining vertices B and C in x and y directions, respectively. In order to simulate the uniaxial straincontrolled loading conditions, displacement was applied to vertex D in x-direction. Additionally, in order to study the effects of multiple point constraints (MPCs) and periodic boundary condition (PBCs) on the modelling results, simulations by applying simple boundary conditions (i.e. fixing the RVE at the left and bottom sides in the x and y direction, respectively, and applying displacement to the right side in the x direction) was performed. The results obtained were almost identical to those with MPCs/PBCs, in terms of local as well as global stress/strain response. This confirmed that MPCs/PBCs neither reduce the effectiveness of the RVE model nor affect the results.

RVE size and mesh sensitivity
In order to obtain the true representative properties of heterogeneous material, it is necessary to determine the appropriate RVE size [11]. In this paper, the size of RVE was evaluated in terms of global stress-strain response by carrying out a series of FE analyses. Homogenisation based on averaging theorem over the whole RVE was used to obtain the stress-strain curves of the material. This homogenisation technique was implemented using a postprocessing numerical subroutine written in computer language FORTRAN and interfaced with ABAQUS.
RVEs with realistic material microstructure and containing different number of grains of various sizes were developed using the aforementioned technique. In this work, five RVEs were created with a sizes of 300 Â 300, 450 Â 450, 600 Â 600, 750 Â 750 and 1000 Â 1000 lm 2 respectively. The corresponding number of grains for each RVE was 42, 90, 169, 250 and 445, respectively. FE analysis of each RVE was performed by applying a monotonic tensile loading condition up to 1% strain at a strain rate of 0.0083%/s. The results obtained for each RVE is compared in Fig. 3, which shows that the convergence of stress-strain curve was obtained for 250-grain RVE with a dimension of 750 lm Â 750 lm. Further FE analysis of the RVE containing 250 grains exhibited the negligible perturbation in the stress-strain curves for different sets of random grain orientations (see Fig. 4). These results confirmed that the RVE with 250 grains is large enough to give the effective mechanical properties of LSHR.
The mesh convergence study, in terms of element number, was also investigated for a 250-grain RVE. Three different FE meshes were considered for this study to examine the convergence of stress-strain response of the RVE. The number of elements for the three meshes are 7871, 79,284 and 149,502, corresponding to an average element size of 0.71 lm, 0.071 lm and 0.037 lm, respectively. The results showed a good local convergence achieved for 79,284-element mesh (Fig. 5), which was subsequently used in our simulations of the mechanical behaviour of LSHR.

Model development
This part of work was inspired by a fatigue test on a plain specimen, with dimensions of 50 mm Â 4 mm Â 4 mm, in vacuum (1 Â 10 À5 mbar) and at high temperature (725°C). The test was load-controlled with a trapezoidal loading waveform (1s-1s-1s-1s). The maximum load was chosen to be 1.5 kN, with a load ratio of R = 0.1. Under fatigue, a short crack was observed to emerge from the centre location of the top surface and grew through a number of grains, as shown by the EBSD mapping (Fig. 6).
To simulate the growth of such a crack, an FE model with the same grain microstructure and orientation as obtained from EBSD mapping of the region containing the fatigue crack was developed   using the aforementioned technique (see Section 3.1.1). In this model, a transgranular pre-crack (about 1 lm) was introduced to the grain in which crack initiation was observed experimentally. In order to determine the load to be applied in the microstructural model, a viscoplasticity model was used to determine the deformation level at the centre of the plain specimen (see Fig. 7) subjected to the same loading condition as the test. The details of viscoplasticity model and model parameters are given elsewhere [30]. The viscoplasticity FE analysis showed that the centre of the top surface was found to experience a maximum and a minimum strain of 0.14% and 0.81%, respectively, when subjected to the fatigue loading condition of the experiment. In order to simulate the crack growth, the maximum and the minimum strain values, calculated by the viscoplasticity model, were consequently applied to the microstructural FE model with the same load ratio and waveform as the fatigue test. This technique is equivalent to the sub-modelling technique. In order to prevent rigid body movement of the model, the corner node on edge K (see Fig. 6) was fixed in the x and z directions.

XFEM procedure
In order to simulate the growth of the surface crack, the XFEM approach in combination with crystal plasticity was adopted. In the XFEM, a crack is represented by enriching the classical displacement-based finite element approximation through the framework of partition of unity [27]. A crack is modelled by enriching the nodes whose nodal shape function support intersects the interior of the crack by a discontinuous function, and enriching the nodes whose nodal shape function supports intersect the crack-tip by the two-dimensional linear elastic asymptotic neartip fields. XFEM enrichment was applied to the whole model to allow the crack to propagate through a solution dependent path determined by local material response. Thus it allows the crack propagation to be predicted without pre-defining the crack growth path (a big advantage over cohesive element approach). Another advantage of XFEM is its less dependency on mesh. In conventional approach of FE modelling of crack growth, mesh should be conformed to geometric discontinuity (i.e. crack path), which requires mesh refinement. However, with special enriched function combined with additional degrees of freedom, the presence of discontinuity is possible for XFEM with less dependency on mesh. Also, the average element size used in our XFEM is around 2.43 lm, which is sufficiently small for both stationary and growing crack analyses based on our experience with crack deformation modelling for nickel alloys (convergence is normally achieved for mesh size less than $3 lm) [11,25,26].
A strain accumulation criterion was adopted in this work. It was assumed that the crack starts to grow when the accumulated plastic strain reaches a critical value at the crack tip. The accumulated inelastic strain is referred to the accumulation of equivalent inelastic strain, whose rate is defined as _ _ F P is the rate of inelastic deformation gradient and I is the identity matrix. It is assumed that crack growth direction is orthogonal to that of the maximum principal strain. During the simulation, an energy-based criterion was used to define the evolution of damage until eventual failure. In XFEM, the damage energy refers to critical energy release rate for fracture, which was calculated as 56 N/mm from the fracture toughness of LSHR at 725°C (in the range of 100 MPa p m). In the present work, three different values of accumulated plastic strain (0.01, 0.1 and 0.23) were chosen for crack growth prediction. The experimental observation of crack growth in this research was performed on the surface of the specimen, for which plane stress analysis is suitable due to negligible influenced from the material beneath. Nevertheless, the psudo-3D model (see Fig. 16) used in this work has stress-free surfaces and is equivalent to plane stress model. Also the crystal plasticity UMAT in this study shows better numerical convergence for 3D solid element.

Calibration of model parameters
In this work, the 12 octahedral slip systems are considered to be potentially active considering LSHR has a FCC crystal structure. The model presented in Section 2 has a number of materials constants that need to be determined from experimental data. The components of stiffness, C 11 , C 12 and C 44 , were first evaluated according to the elastic response of the material from the FE analysis of RVE under monotonic loading. The value of these stiffness components were found to be C 11 = 491.12 GPa, C 12 = 251.47 GPa and   After an initial estimate of the unknown material parameters, subsequent iterations were performed until a consistent description of the experimental data set was obtained. The material constants associated with the internal slip system variables, namely, h S , d D , h B and r D , were fitted in a similar manner to that outlined above. Overall, extensive FE analyses were carried out in order to obtain a good description of monotonic, relaxation and cyclic deformation. The fitted parameters are given in Table 1.

Stress-strain response and distributions
Comparisons of the experimental data and the model simulations are given in Fig. 8 for monotonic tensile and stress relaxation behaviour (strain rate = 0.0083%/s, maximum strain 1%) and in Fig. 9 for stabilised loop of a strain controlled cyclic test (strain rate = 1%/s, strain range De = 1% and strain ratio = 0). Clearly, the model predictions are in close agreement with the test data. Furthermore, the shape of the hysteresis loop is also well reproduced by the crystal plasticity model. Figs. 10 and 11 show contour plots of the normal stress and strain in x-direction (direction of loading) at saturation stage for the RVE with realistic microstructure. A heterogeneous stress and strain distribution can be clearly seen from these images. The maximum value of the local stress is about 5 times the global stress and the maximum local strain is about 1.5 times the global strain. Stress and strain concentrations are located in different grains, generally at grain boundary regions. A closer inspection shows that most of the stress concentrations are at the triple or multiple points where smaller grains encounter larger grains. The potential reason for high stress concentration lies in the fact that the mis-   match of the mechanical properties between the neighbouring grains is confined to a small area, which generates pronounced stress concentrations. In order to study the local convergence, stress-strain responses of individual grains were compared across RVEs with varying size. It was verified that local stress-strain response is basically dominated by grain orientation. As long as the orientations of all surrounding grains are unchanged, stressstrain response of individual grain is very consistent and insensitive to RVE size. Local convergence was also studied by Przybyla and McDowell for Statistical Volume Element (statistical representation of underlying microstructure) approach [35][36][37]. It was reported that convergence can be easily achieved for an individual grain as long as it is enclosed by two layers of neighbouring grains within the SVE [35]. In our study, we also confirmed that RVE size is not particularly an issue as long as the orientations of surrounding grains remain consistent across the RVEs of different size. A few grains with high stress and strain concentrations were selected, and the stress-strain response is plotted in Fig. 12 for each individual grain. Some grains (e.g. grains 1 and 4) with high strain concentration have low strain hardening in the stress-strain response and are normally referred to soft grains. On the other hand, other grains (e.g. grains 2 and 3) with high stress concentrations have high strain hardening in the stress-strain response and are normally referred to hard grains (see Fig. 12). As demonstrated in Fig. 12, the stress-strain loops are much wider for the soft grains than for the hard grains, indicating the development of more inelastic deformation in the soft grains RVE models with realistic microstructure are often developed from surface-based images obtained from EBSD or SEM. One simple method is to overlay the microstructural map on the regular finite element mesh and assign the lattice orientation of each grain to the corresponding element integration point. As shown in Cheong et al. [38], the generated model has the capability to predict the regions prone to crack initiation and growth, but with some false predictions due to its inability to incorporate precise microstructural information such as the presence of grain boundaries. In this paper, the series of RVE models were developed based on the real grain microstructure with delineated grain boundaries, which is more accurate but requires strenuous effort. This is also the method predominantly adopted in literature (e.g. Vaxelaire et al. [39], Delaire et al. [40] and Sistaninia and Niffenegger [41]), where efforts were made to study intragranular strain/stress heterogeneities, deformation localisation and fatigue life of the relevant materials, but only with a small number of grains. In this paper, we particularly evaluated the effect of RVE size on the simulated global stress-strain response by exploring a wide range of the number of grains. The results highlighted the importance of considering the sufficient number of grains if a RVE with realistic grain microstructure is used to capture the representative global response of polycrystalline alloys. This is further explored in the following Section through comparison of RVE models with artificially generated grain microstructures.

Comparison with an artificial microstructure model
In order to help understand further the results obtained, a RVE with artificial grain microstructure was developed using the Voronoi technique. There are several parameters used to construct the artificial microstructure similar to the realistic one [42,43]. These parameters include the number of contiguous neighbour grains, grain volume and orientation/misorientation distribution functions. The number of contiguous grains and grain volume are more relevant for 3D grain microstructure and unable to be quantified in this research. In this paper, for simplicity, the artificial RVE was constructed only based on the average grain size, without considering the orientation/misorientation distribution for which statistical analysis will be required and is a limitation of this research work. Again, random orientations were assigned to each Voronoi cell representing each individual grain of the polycrystal. The distribution of stress and strain obtained from the artificial RVE is shown in Fig. 13. This also confirms the heterogeneous nature of the deformation. An interesting observation found by comparing Figs. 8 and 13 is that considerably higher values of localised stress are obtained for the RVE with a realistic microstructure as compared to those of the artificial one. A possible reason for this phenomenon is the irregular edges of some grain boundaries forming very sharp curves which are common in realistic microstructure in contrast to the straight grain boundaries in the artificial microstructure. Also the local variation of grains sizes, i.e. larger and smaller grains, magnified the property mismatch between individual grains, especially the neighbouring ones. However, the converged stress-strain response of the artificial RVE (150 grains), which is shown in Fig. 14, is very close to that given by the RVE with the realistic grain microstructure.
These results suggest that in order to have a reliable interpretation of microstructural deformation, it is important to consider the realistic grain microstructure in the RVE-based simulation of mechanical deformation of polycrystalline materials. RVEs with realistic microstructure possess considerably higher localised stress concentrations than artificial ones, signalling the chance of earlier crack initiation and material failure. In order to study the convergence behaviour of RVEs with an artificial microstructure, four FE models consisting of 50, 100, 150 and 200 grains corresponding to RVE sizes of 255 Â 255 lm 2 , 360 Â 360 lm 2 , 441 Â 441 lm 2 and 509 Â 509 lm 2 , respectively, were created using the Voronoi tessellation. For all RVEs, the grains have an average size of 36.05 lm. To realise repetitive deformation of the RVEs, periodic boundary conditions were applied to the RVE in the same way as mentioned in Section 3.1.1. An interesting observation is that the RVEs with artificial microstructure converged nicely at 150 grains (see Fig. 14) as compared to 250 grains for RVEs with realistic microstructure (Fig. 3). The results suggest that RVEs with artificial microstructure lead to faster convergence during simulation, but it may give misleading prediction of localised stress-strain response and behaviour as well as associated fatigue failure even if the global stress-strain behaviour is predicted reasonably well. Therefore, it is important to always consider the realistic microstructure for proper simulation of both global and local Fig. 11. Strain contour at the maximum strain of the stabilised cycle, (1% strain range, 0.0083%/s strain rate, T = 725°C). mechanical response of the material. It should be noted that random crystallographic orientation is used for both realistic and artificial RVEs, due to the lack of EBSD data on grain orientations. As shown earlier (Fig. 4), the effect of grain orientation was basically negligible for the converged RVE size (750 Â 750 lm 2 ). Based on this, it is supposed that the stress-strain response of the RVE would remain about the same when using the actual crystallographic orientations, which are generally random for a sufficiently large RVE, measured from EBSD analysis. So, the general conclusion is still valid. However, this will be verified in our further study when the full EBSD data are available.

Prediction of surface crack growth
Using the XFEM and crystal plasticity, the predicted crack length is plotted in Fig. 15 as a function of the number of fatigue cycles. Obviously, the increase in crack growth rate was observed for the lower values of accumulated plastic strain. For different values of accumulated plastic strain, fluctuation in crack growth rate was always observed. This fluctuation in crack growth is linked with the presence of grain microstructure. When the propagating crack reached a grain boundary and grew into the next grain, it changed, depending on the grain orientation. The crack propagation path predicted by the simulations is shown in Fig. 16 in comparison with the experimentally observed crack path shown in Fig. 6(b). Experimentally observed cracks showed both an intergranular and transgranular path whereas only a transgranular crack path was observed in the simulations. A potential reason for this discrepancy lies in the fact that the vacuum was not perfect during the experiments, which leads to intergranular crack growth.
An interesting observation is that the crack propagation rate is not constant across different grains (see the crack length versus the number of cycles for Grains A and B in Fig. 15) although it generally increases with the number of loading cycles. For instance, hard grains may act as obstacles to crack growth, and the growth of a crack that is controlled by the accumulated plastic strain at the   crack tip decreases when it encounters a hard grain. This mechanism yields an oscillating crack growth rate that is the typical characteristic of microstructurally short fatigue crack growth. The scattering is clearly evidenced from Fig. 17, which shows the predicted crack growth rates da/dN against Da. Typically, crack growth rate is presented in terms of DK for fatigue loading; however, in our case, it is not convenient to calculate an accurate DK at the grain level, and therefore Da is used in Fig. 17. The source of the scatter in crack growth rates is the mismatch in material properties of the neighbouring grains. A straight line was also fitted for crack growth rate for each critical accumulated plastic strain, and the crack growth rate tends to increase with the decrease of critical accumulated plastic strain. The fitted straight lines are presented merely to show the trend of crack growth rate as a function of crack length which is difficult to follow due to the high scatter.
As seen from Fig. 17, for different values of critical accumulated plastic strain, fluctuation or scattering in crack growth rate was always observed, due to the nature of random grain orientations. When the propagating crack reached a grain boundary and grew into the next grain, it experienced either retardation or acceleration, depending on the grain orientation. In order to study the effect of orientation on crack growth, simulations considering only single property (i.e. single orientation) were carried out. Three single orientations were chosen, namely, [0 0 1], [1 1 1] and a random one (i.e., a random choice of the three Euler angles). The results are given in Fig. 18, in comparison with that obtained for a polycrys-talline model (taken from Fig. 15). The critical value of accumulated plastic strain was chosen to be 0.1 for the simulations. It is observed that crack growth is highly dependent on grain orientation, which can be captured by crystal plasticity model. The predicted crack growth rates da/dN against Da are shown in Fig. 19 for single orientation, in comparison with that for polycrystal model (taken from Fig. 17). As expected, there is no scattering in crack growth rate for single orientation. Despite the importance of characteristic fracture strain in the instantaneous crack growth, it is not well explored in literature [25]. The plastic strain accumulation criterion is linked with the cyclic deformation, thus making it a more appropriate approach to quantify the crack tip mechanics and physical process of material damage under fatigue [30]. The present work also confirmed the suitability of this criterion for crack growth prediction. To the authors' knowledge, it is the first time that the accumulated plastic strain, in combination with XFEM and crystal plasticity, is used to predict the crack growth at microstructural level in polycrystalline alloys under high temperature fatigue loading.
Although the model is capable of simulating the crack growth path and rate in a polycrystalline, still there are limitations of this model. One of the major limitations is the lack of crystallography which would allow consideration of slip planes for growing crack analysis. Slip planes play a significant role in the crack growth as it is likely to grow along the preferred slip planes, resulting in a meandering crack path [44,45]. Consequently, a reason for the sim-   ulated crack path along nearly straight line (Fig. 16) as opposed to the experimental observation can be attributed to the fact that the slip planes were not considered. Another limitation of this work lies in the crack growth criterion. Since alloy LSHR is a facecentred-cubic (f.c.c.) polycrystalline metallic alloy with random grain orientations, with each grain orientation exhibiting a different mechanical behaviour under same loading conditions, a single value of critical accumulated plastic strain cannot capture the intrinsic fluctuations in crack growth caused by orientation mismatch and grain boundaries. This is also the limitation imposed by ABAQUS as it only contemplates single value for each crack growth criterion used. Further work is required in this direction, including a direct comparison with experimentally measured crack growth rates and digital mapping of strain accumulation along the crack path (including the RVE model) -this can be accomplished by using high resolution EBSD analysis and digital image correlation.

Conclusion
A RVE, with realistic grain microstructure and random grain orientations, was developed for crystal plasticity modelling of cyclic deformation of a nickel based superalloy at elevated temperature. Model parameters were calibrated using monotonic, stressrelaxation and cyclic test data at 725°C. Deformed RVE exhibited severe heterogeneous distributions of stress and strain due to property mismatch of individual grains, with magnitudes considerably higher than those in an artificially constructed RVE making it essential to always consider the realistic microstructure for polycrystalline materials. In combination with XFEM, the calibrated crystal plasticity model was further applied to predict crack growth at the grain level using an accumulated plastic strain fracture criterion. Results showed that there was a significant variation of crack growth rate in different grains, indicating crack growth acceleration and retardation imposed by the intrinsic material microstructure.