Computational and experimental study of crack initiation in statistical volume elements

. Fatigue crack formation and early growth is significantly influenced by microstructural attributes such as grain size and morphology. Although the crystallographic orientation is a primary indicator for fatigue cracking, the neighbourhood conformed by the first and second neighbour grains strongly affect the fatigue cracking driving force. Hence, two identical grains may result in different fatigue responses due to their interactions with their microstructural ensemble, which determines the fatigue variability. Naturally, macroscopic samples with millions of grains and thousands of competing microstructural neighbourhoods can effectively resemble a representative volume element in which fatigue failure may seem deterministic. However, when considering systems in which fatigue failure is controlled by hundreds or less of grains, fatigue failure is stochastic in nature and the samples are not a representative but a statistical volume. This work studies fatigue crack nucleation in micron-scale Ni beams that contain a few hundred grains. This work presents 3D crystal plasticity finite element models to compute stochastic distribution of fatigue indicator parameters that serve as proxies for crack nucleation in statistical volume elements. The integration of experiments with models provides a method to understand the irreversible deformation at the grain level that leads to fatigue cracking. Our results explain the role of grain morphology of crack nucleation distribution


Introduction
The past decade have seen developments in the manufacturing of miniaturised components. A good example of such components are micro-electro-mechanical systems (MEMS), which have widespread application in industries with stringent reliability requirements such as automotive, aerospace, medical devices etc. [1]. Furthermore, metallic rather than silicon MEMS have emerged as a technological solution for component miniaturization under detrimental environments (e.g., high temperature). As a consequence, there is an increasing interest in the study of the mechanical reliability of microscopic metallic components, among which fatigue failure is a leading concern [2].
Microscopic devices have structural dimensions that are comparable to the characteristic length scale of the material microstructure, and they are thus affected by size effects [3]. In addition, a small number of microstructural units (e.g., grains, twins, grain boundaries) are present on each device, which turns it into statistical rather than representative volume elements (SVE vs RVE) [4]. Hence, two devices cannot be assumed to be identical, which raises particular concerns about their fatigue reliability [5][6][7][8].
Recent years have seen the development of microstructure-sensitive crystal plasticity computational methods [9,10] to represent the fatigue response at the grain scale. These approaches characterise the fatigue damage driving force at the microstructure scale by means of Fatigue Indicator Parameters (FIPs) [10][11][12]. The use of FIPs provide a framework for quantifying the likelihood of crack nucleation at the microstructure scale [13]. Among many FIP formulations, Fatemi and Socie proposed a shear strain based parameter [14], which correlates with early crack nucleation in metals under multiaxial loading [15]. Castelluccio and Mcdowell [10,16] employed a Fatemi-Socie parameter for each octahedral slip system as surrogate driving force for early crack initiation. Their proposed FIP was implemented along a 3D crystal plasticity mesoscale model that was used to investigate microstructure variability on fatigue behaviour at the micro scale. This paper employs the 3D crystal plasticity mesoscale model proposed by Castelluccio and Mcdowell [10,16] to study the influence of grain shape on the fatigue crack initiation in MEMS under bending. We recreate multiple synthetic microstructures to inform finite element models, which effectively correspond to SVE. These synthetic microstructures reproduce elongated grains that are statistical equivalents of the real microstructure. In addition, the fatigue behaviour of the micro beam is investigated for equiaxed microstructure. The results demonstrate the role of grain morphology on fatigue crack initiation and present good agreement with experiments.

Micro-beam fatigue damage
We consider the micro-resonator technique developed by Pierron and co-workers [17] , which applies cyclic loading to a cantilever micro-beam under fully reversed loading (Fig. 1). The micro-resonators consists of a micro-beam connected to a large mass with two comb drives on each side. Stress controlled fatigue tests are performed on the micro-beams by driving the specimen at an in plane resonance frequency of, fo 8-10 kHz. The result is a fully reversed cyclic bending of the micro-beams under large enough amplitude of rotation. The micro-beam is made of electroplated nickel fabricated using the LIGA technique. The dimensions of the micro-beam are 11.4μm at its thinnest section, a thickness of 20 μm, and a length of about 60 μm. The cross section of the micro-beam at the gauge section consists of 5-10 grains across its width, and 2-4 grains across its thickness [18]. The small size of the micro-beam contain a few tens to hundreds of grains on the surface, which effectively makes the beam a SVE for fatigue failure. 3 Computational fatigue damage prognosis

Synthetic microstructures
We constructed a 3D beam model with synthetic microstructures by informing Dream 3D [19] with multiple layers of Electron Back Scatter Diffraction (EBSD) images. The micro-beam has a strong (001) crystallographic texture and consists of large columnar grains with smaller equiaxed grains filling the gaps between the columnar grains. Studies have shown that fatigue cracks tend to initiate from the relatively larger grains [20]. Hence, we need to consider the role of the morphology of the columnar grains on fatigue cracking. The EBSD data are processed to extract grain morphological and crystallographic statistics using Dream 3D with a misorientation threshold of 2 deg. Grains that are less than 230 voxels are excluded and added to neighbouring grains that meet the characterisation criteria. The corresponding histograms showing grain size distributions, misorientation, and the texture pole figure are shown in Fig. 2. An in-house Matlab script was integrated with Dream 3D and Abaqus python scripts to generate and mesh the microbeam finite element model. Dream 3D generated a spatial distribution of grains for a voxelised synthetic microstructure with geometric dimensions equal to that of the gauge section of the micro-beam and same texture as the real microstructure. The resulting grain diameter was approximated by a log normal distribution with a mean and standard deviation of 1.19μm and 0.26μm, respectively. We recreated ten equivalent microstructural realizations containing the micro-beam gauge section that represents the microstructure and texture from experiments. In addition, we recreated ten artificial realizations with similar grain size distribution and texture, but equiaxed grain morphology. Fig 3 shows the statistical distributions of both grain shape shapes sizes as modelled in Dream 3D.
To mitigate mesh sensitivity [21], we consider the average fatigue indicator parameters over volumes corresponding to the microstructure characteristic length scale. On each grain and for each slip orientation we choose to average the FIP over the elements that correspond to the volume occupied by slip bands. These sets are one-element-thick slabs that span the entire grain cross-section and are oriented parallel to each slip systems. All four slip systems normal directions are considered, which means that the band sets can fully reconstruct each grain four times. These bands serve as non-local volume to quantify the fatigue driving force [21].   The matlab script creates Abaqus input files that includes the micro-beam geometry, mesh, loading conditions, and grain morphology. The script also identifies the elements that make up the slip bands volume and write these into the input file as element sets. The micro-beam is meshed with cuboidal reduced integration 8-node linear brick elements (C3D8R). The micro-beam gauge section is discretised with 17500 elements for all instantiations. The loading sequence consist of applying an equivalent displacement controlled simulation for a fully reversed bending rotation angle of 25 mrad to the nodes at the micro-beam free end with the support end fixed. A rotation angle of 25 mrad is equivalent to an experimental microbeam fatigue test ran at ̴ 470 MPa with a fatigue crack nucleation life, that lies between 1.42 x10 7 and 2.34 x10 7 cycles in a vacuum environment [22].

Crystal plasticity model
We consider a multiscale crystal plasticity model proposed by [23], which was validated for nickel single and polycrystal. The rate-dependent flow rule follows parametrization proposed by [24] i.e., (1) in which is the density of mobile dislocations, is the mean advance distance of a dislocation between obstacle bypass events, B is the Burgers vector ~2.5 x 10 -10 m, and is an estimate of the attempt frequency ~10 12 s -1 . Also, is the energy required to bypass dislocation barriers, is the effective stress, is the directional intragranular contribution to the slip system back stress, is Boltzmann constant, and T is the absolute temperature in Kelvin. The activation enthalpy, , pertaining to the mean effective barrier, has been parameterised [24,25] in the form

(2)
Here F0 and p, q are the activation energy and profile parameters, respectively, for dislocation barrier bypass at =0, is thermal slip resistance at 0 K, and μ and μ0 represent the shear modulus (C44) at temperatures T and 0K, respectively. In addition, the effective stress that drives the barrier bypass at the slip system level is given by, where is the local slip system resolved shear stress, is the threshold stress and are the Macaulay brackets. The threshold stress required to bow a dislocation against a dislocation pile-up between dislocation walls is given by (4) in which corresponds to an intrinsic lattice friction stress, the second term corresponds to the stress required to bow out a dislocation, which depends on the dislocation line energy, , and dislocation wall fraction, . The third term corresponds to the strength interaction between collinear dislocations in pile-ups, which depends on the strength of the self-interaction coefficient, .
The wall fraction follows a parametrization based on the maximum plastic shear strain range per cycle, where, , , and are material constants that can be estimated by computing from TEM images.
The back stress evolution is based on an Eshelby mean field approach that incorporates dislocation substructure effects [23], (6) is the Eshelby tensor for a prolate spheroid, η corresponds to the normalized dislocation mean free path. The back stress evolution relation represents the homogenized response of a quasiperiodic dislocation structure that induces residual stress in the material due to the coexistence of a plastically deformable matrix surrounded by an essentially linear elastic structure of dense dislocation walls.
The mean free path for dislocation glide, , depends on the local strain range and it is computed following the scheme proposed by [23]. Thus, the mean advance distance is determined by identifying the appropriate dislocation mesoscale substructures.

(7)
Here, is the characteristic dislocation substructure [26,27], i.e., The mobile dislocation density on slip system α results from a balance between dislocation multiplication, annihilation and cross slip (9) Here kmulti = 1 for cell structures and kmulti = 2 for parallel walls such as PSBs and labyrinth, γ=1 if or ; otherwise and . The rate of cross slip of mobile screw dislocations from plane α to ζ is estimated using the probabilistic formulation proposed by several researchers [28][29][30], but considering the local shear stress to account for screening effects from dislocation structures, i.e., (10) Here, and are the cross slip characteristic frequency and dislocation length, respectively, is the critical cross slip stress, and is the magnitude of the local stress on the cross slip plane of a screw dislocation in system α. The cross slip activation volume is inversely proportional to the slip system shear stress, i.e.
The critical cross slip stress follows, (12) where corresponds to the annihilation distance of screw dislocations. Parameters are summarized in Table 1 to   Table 2. Parameters related to atomistic scale unit processes for pure Ni.   The model was implemented as a user-material sub routine (UMAT) with the use of the Newton Raphson incremental line search and backward -Euler implicit integration method in Abaqus [31].

Fatigue life damage prediction model
To quantify the driving force for trans-granular fatigue crack nucleation we compute a crystallographic version of the Fatemi-Socie fatigue indicator parameter (FIP). The FIP is evaluated for each octahedral slip system. (13) where is the cyclic shear strain range at each slip system, is the peak stress normal to the slip system, the cyclic yield strength of the polycrystal, and k=0.5 as proposed by Fatemi-Socie [14]. To regularise the finite element discretisation and also represent the non-local influence of the fatigue damage zone, the FIPs are computed on each integration point, and are averaged along mesoscale bands parallel to slip planes across an entire grain as shown in Figure  4. We consider a crack nucleation formulation that predicts the minimum number of cycles it takes for the first grain to crack ( ) as a function of the dimension of the mesoscale bands , where is an irreversible coefficient and is a reference scale considered for extending the crack and is given by, (15) where is the sum of the length of the band in the ith grain, and are the length(s) of bands in neighbouring grains having low disorientation ith grain, with as the disorientation factor which is computed as (16) where is the disorientation between the crystallographic orientations of the crystallographic orientations of the bands, and the Macaulay brackets are defined according to (a) = a if a>0, (a) =0 if a ≤ 0. Castelluccio and Mcdowell [10,16] provides a detailed explanation of the mesoscale fatigue damage prediction model as used in this work.

Results and discussions
A total of 10 microstructure realisations were considered for each grain shape to compare the fatigue crack nucleation life, . For each realization we compute , and Fig. 5 represents the boxplot showing the fatigue crack nucleation life for both grain shapes. The nucleation life for the micro-beam instantiations of the grain shapes represent the distribution of extreme values as each crystal plasticity finite element (CPFE) simulation only contributes a single value to each grain shape statistics, that is the grain with the lowest nucleation life. The boxplot shows that nucleation life in the lower, medium and upper quartiles of the equiaxed grain is about 10% larger than the elongated grain, and with a larger dispersion in its spread. The prediction of a higher nucleation life for the equiaxed grains compared to the elongated grains is in agreement with an experimental study by [32] that was carried out for nickel based super alloy.

Conclusions
The fatigue behaviour of nickel micro-beams is investigated for the influence of grain shape on its crack nucleation life using CPFE modelling. The CPFE simulation results show that that the fatigue life at the mesoscale level is influenced by grain shape. The model correlates well with experimental data on the number of cracked grains during cyclic loading. The results agree with fatigue behaviour trends from even experiments reported in literature. The mesoscale fatigue model showed its capability to model the influence of grain shape on the crack nucleation life.