Modeling Crack Propagation in Polycrystalline Microstructure Using Variational Multiscale Method

Crack propagation in a polycrystalline microstructure is analyzed using a novel multiscale model. The model includes an explicit microstructural representation at critical regions (stress concentrators such as notches and cracks) and a reduced order model that statistically captures the microstructure at regions far away from stress concentrations. Crack propagation is modeled in these critical regions using the variational multiscale method. In this approach, a discontinuous displacement field is added to elements that exceed the critical values of normal or tangential tractions during loading. Compared to traditional cohesive zone modeling approaches, themethoddoes not require the use of any special interface elements in themicrostructure and thus canmodel arbitrary crack paths. The capability of the method in predicting both intergranular and transgranular failure modes in an elastoplastic polycrystal is demonstrated under tensile and three-point bending loads.


Introduction
Efficient microscale modeling tools are needed to compute microstructure-dependent properties of advanced structural alloys used in aerospace, naval, and automotive applications [1].Microstructural effects become an important consideration in regions of stress concentrations such as notches, cracks, and contact surfaces.While crystal plasticity finite element (CPFE) method [2,3] has emerged as an effective tool for simulating the mechanical response of the microscale containing aggregates of few hundred metallic crystals, simulation of "macroscale" components that contain millions of grains is a challenging task even when using current stateof-the-art supercomputers.Multiscale methods developed to address this problem can be broadly classified into two categories: concurrent and hierarchical methods.In hierarchical methods, microscale models that contain features significantly smaller than the coarse scale mesh resolution are simulated to recover the material properties at the macroscale.An example is the method of computational homogenization in which stresses in a polycrystal are volume averaged to compute the macroscopic stress state [4].Homogenization approaches are inadequate in problems where there is no clear separation of length scales, for example, when simulating localization or failure phenomena [5].
In concurrent methods, both micro-and macroscale problems are solved together in a single domain using multiresolution meshes.Most popular among these methods are quasicontinuum methods, where resolution of meshes from atomistic to continuum scales has been considered in the past to study crack propagation [6].In concurrent methods, treatment of interfaces between fine and coarse scale regions becomes problematic when multiphysical domains (e.g., crystal plasticity at fine scale and continuum plasticity at coarse scale) are considered.For example, Dunne et al. [7] employed a combined continuum crystal plasticity finite element approach for modeling a three-point bend test.One of the potential difficulties identified in the model was the inconsistency of yield stress at the interface between the domains.To minimize the effect, the choice of the critical resolved shear stress in the crystal plasticity model was made such that polycrystal yielding occurred at a value close to the continuum yield stress.In this paper, the interface problem is solved through the use of a single physical 2 Mathematical Problems in Engineering model at both the fine and coarse scales.This is achieved (i) by explicitly resolving single crystal deformation at the microscale (crack tip) using a crystal plasticity model and (ii) by employing a statistical model at the macroscale (far field), where polycrystal aggregates are represented using an orientation distribution function (ODF).Under an applied deformation, microstructure evolution in the far field is simulated by numerically evolving the ODF (rather than the microstructure itself) using conservation laws [8].Such statistical approaches are significantly (several orders of magnitude) faster than full order (CPFE) methods [9,10] and balance the trade-off between accuracy and computational efficiency in far-field regions.
The paper presents the first use of this efficient multiscaling framework for prediction of crack trajectories in a polycrystalline microstructure.Towards this end, a cohesive zone model (CZM) [11,12] is introduced at the microscale for modeling the crack kinematics.In this approach, the crack is represented using special interface elements that obey a relationship between the interface traction and the crack separation.At the microscale, these relationships have been motivated from the results of atomistic simulations of grain boundaries [13][14][15].Cohesive element-based approaches have been previously used for modeling intergranular crack paths [16][17][18].However, in order to simulate both transgranular and intergranular crack paths in an elastoplastic material, these interface elements are needed practically at every element-to-element boundary.This both makes the problem computationally expensive and also results in unsatisfactory mesh convergence properties [19].Other computational methods have emerged in the past two decades that incorporate displacement discontinuities using enrichment functions that are classified into element enrichment methods (e.g., variational multiscale method (VMM) [20,21]) and nodal enrichment methods (e.g., Extended Finite Element Method (XFEM) [22,23]).Sukumar et al. [24] employed XFEM to simulate intergranular cracks in elastic polycrystals.Rudraraju [25] and Rudraraju et al. [26] recently explored the use of VMM for prediction of crack propagation in elastic composites and show good comparison with experiments.However, the use of these methods has not yet been demonstrated in elastoplastic regimes, in particular, for cases involving intergranular to transgranular failure mode transitions at the microstructural level.In this work, a numerical framework based on VMM is proposed for addressing such a case.The paper is organized as follows.In Section 2, VMM scheme is briefly introduced.In Section 3, the polycrystalline constitutive model will be briefly introduced followed by the statistical (ODF) modeling approach in the far field.In Section 4, examples illustrating crack propagation at the microstructural level are presented, while concluding remarks are presented in Section 5.

Introduction of Variational Multiscale Method
The microstructural simulations (in 2D) described here employ an updated Lagrangian framework wherein the configuration at the previous time step ( =   ) serves as a reference configuration for the next time step ( =  +1 ).Let F  refer to the total deformation gradient in the reference configuration (Ω  ) at time   with respect to the initial undeformed configuration (Ω 0 ) at time  = 0. Similarly, the total deformation gradient in the current configuration (Ω +1 ) at time  =  +1 is written as F +1 .Then F  denotes the relative deformation gradient between the two configurations; that is, A crack in the reference configuration is represented as a surface Γ   .The displacement jump in the crack is denoted by ⟦u⟧.The crack surface in this configuration has a normal direction n and tangential direction m as shown in Figure 1.
Standard finite element procedure for deformation problems involves solving for the nodal displacements d in the finite element mesh.In the variational multiscale approach, an additional degree of freedom ⟦u⟧ is introduced within cracked elements such that the overall corrected displacement in the cracked element is u = d−⟦u⟧.The finite element equations for the crack element are described next.

Model for the Cohesive Zone. Traction continuity on the crack surface is given by
where t  is the nominal traction vector on the crack surface and (   ,    ) are the normal and tangential components of the nominal traction in the reference configuration.Here, the Piola Kirchhoff I stress P in the reference configuration at time   is expressed in terms of the Cauchy stress  in the current configuration at time  +1 as P = (det F  )F −  .The traction components are assumed to be linearly related to the components of the crack displacement ⟦u⟧ in the reference configuration via the softening moduli H: where H  , H  are the softening moduli for the mode I and mode II crack and   0 ,   0 are the critical values of normal and tangential nominal tractions that lead to the formation of a new crack.The linear traction separation law is shown in Figure 2.
Substituting (3) into (2), the residual r for the crack element is obtained as Linearizing the residual with respect to the displacement unknown gives Here, P = (det F  F −  ) is the increment in the PK-I stress which can be simplified as The increment in the Cauchy stress can be computed in an implicit form  = C[F  ] using the crystal plasticity constitutive model as shown in Appendix B. Note that the increment in deformation gradient is of the form F  = (∇  d − ∇  ⟦u⟧) which can be substituted in (6) such that (5) is obtained in terms of the unknown displacement increments.Note that d is unknown and thus the crack element equation (see (5)) is to be solved in a coupled manner with the standard finite element equations (explained next).

Nonlinear Finite Element Model.
For any kinematically admissible test function ũ, the weak form of the virtual work equation for displacement-controlled problems in the absence of body forces is of the following form: The Newton-Raphson iterative scheme is employed to solve the above equation.The increment G due to infinitesimal changes in the displacements is computed as In the variational multiscale method, the test function is written in terms of virtual displacements d; that is, In addition, the increment in the PK-I stress can be written in an implicit form 6). Thus, the linearized equation reduces to 2.3.Coupled Solution Procedure.Let u  +1 be the displacement field at the th iteration; then the linearized equations are solved for the increments u = u +1 +1 − u  +1 .The analysis procedure involves the solution of two coupled equations for the increments in displacements d and ⟦u⟧. Let and ∇  ⟦u⟧ = G{⟦u⟧} denote the gradients in terms of finite element matrices (B d and G) and the nodal unknowns.Each NR step involves solving for increments d, ⟦u⟧ using the system of equations given below: where where H is a matrix computed by rearranging the component of the crack normal n.In our implementation, for example, where   and   are the components of the crack normal n.

Finite Element Representation of Cracks.
Our present implementation employs three-node triangular elements.The displacement d is interpolated using the usual linear shape functions [ 1 ,  2 ,  3 ] for the triangle element.
The displacement discontinuity ⟦u⟧ is interpolated using a specially constructed shape function that ensures that the displacement jump is localized to the element.The fine scale shape function can be considered to be the difference between a linear shape function and a Heaviside jump function centered at the crack.Two possible cases of such a shape function are shown in Figure 3.One of these two shape functions is chosen by comparing the outward normal direction of the element edge not intersected by the crack and the normal direction of the discontinuity within the element.In the first case shown in Figure 3(a), n and n 1 point in the same direction and the shape function is the difference between  1 and a Heaviside function.In the second case in Figure 3(b), n and n 1 point in the opposite directions and the shape function is the difference between (1 −  1 ) and a Heaviside function.
When using this description, the displacement gradient (written in a vector format) is computed using the standard finite element shape functions as where {⟦u⟧} is a 2 × 1 unknown vector for each cracked element.In this equation, index  corresponds to the node common to the element edges intersected by the crack.The correct sign of G is determined by the choice of one of the two shape functions shown in Figure 3.When using the shape function shown in Figure 3(a), a positive sign is used for G.If the shape function defined in Figure 3(b) is used, a negative sign is employed.In addition, the crack normal (n) in all the elements forming a contiguous crack is aligned to the same side of the crack path so that the discontinuity is consistent across these elements.

Constitutive Modeling of Single Crystal Response
A rate independent crystal plasticity theory [27] will be used to model deformation response of particles within each grain.
The theory is based on the notion that plastic flow takes place through slip on prescribed slip systems.For a material with  = 1, . . .,  slip systems defined by orthonormal vector pairs (m  0 , n  0 ) denoting the slip direction and slip plane normal, respectively, at time  = 0, the constitutive equations relate the following basic fields (all quantities expressed in crystal lattice coordinate frame): the deformation gradient F defined with respect to the initial undeformed crystal which can be decomposed into elastic and plastic parts as F = F  F  (with det(F  ) = 1), the Cauchy stress , and the slip resistances   > 0. In the constitutive equations to be defined below, the Green elastic strain measure E  = (1/2)(F  F  − I) defined on the relaxed configuration (plastically deformed, unstressed configuration B) is utilized.The conjugate stress measure is then defined as T = det F  (F  ) −1 (F  ) − .The constitutive relation, for stress, is given by , where L  is the fourth-order anisotropic elasticity tensor.It is assumed that deformation takes place through dislocation glide and the evolution of the plastic velocity gradient is given by where S  0 = m  0 ⊗ n  0 is the Schmid tensor and γ  is the plastic shearing rate on the th slip system.The resolved stress on the th slip system is given by   = T ⋅ S  0 .A rate independent algorithm is employed to solve the single crystal model.The resolved shear stress   is taken to attain a critical value   (the slip system resistance) on the systems where slip occurs.These active systems have a plastic shearing rate γ  > 0. There is no plastic shearing rate ( γ  = 0) on inactive slip systems where the resolved shear stress does not exceed   .The evolution of slip system resistance is given by the following expression: where ℎ  is the slip system hardening term, γ  is the plastic shearing rate on the th slip system, and   0 is the initial slip system resistance on the th slip system.
The algorithm for computing the plastic shear increment Δ  from this model is found in Appendix A. Subsequently, the plastic part of the deformation gradient is updated using (14); the elastic part is computed from F = F  F  .The conjugate stress measure T is then computed from T = L  [E  ] and converted to the PK-I stress for use in the weak form.The slip resistances are also updated at the end of the time step using (15).Finally, the tangent modulus /F for use in the weak form is computed using a fully implicit algorithm described in Appendix B.

Constitutive Model in Macroscopic
Regions.Microstructural effects become important at regions of stress concentrators such as notches, cracks, and contact surfaces.A multiscale model is employed that efficiently captures microstructural details at such critical regions.The approach is based on a multiresolution mesh that includes an explicit microstructure representation at critical regions where stresses are localized.At regions farther away from the stress concentration, a reduced order model that statistically captures the effect of microstructure is employed.The configuration of the multiscale model is shown in Figure 4.
At the macroscale, each integration point models the response of aggregates of several hundred grains probabilistically using an orientation distribution function (ODF).The ODF, A(g), gives the probability density of finding an orientation g in the microstructure.The ODF satisfies the following conservation equation at all times during deformation: where g is a differential volume element (the invariant measure) of the orientation space.In this work, a finite element (FE) mesh is used to model the fundamental region and the ODF is defined at the nodal points of this mesh [8,28].The probability values between nodal points are obtained as a result of interpolation using finite element shape functions.The evolution of ODF is given by the conservation equation ( 16) as where the ODF at time , A(g  ), represents the probability density of crystals with orientation g  and g  represents the volume element in the undeformed (initial) ODF mesh which becomes volume element g  at time .A Jacobian (g  , ) = det(G) gives the ratio of elemental volumes such that g  = (g  , )g  .Here, G is the reorientation gradient given as G(g  , ) = g  /g  which can be computed from the single crystal constitutive model.Using the Jacobian, a map of the current mesh (at time ) to the reference mesh (at  = 0) can be made: The quantity written as Â(g  , ) is the volume density A(g  ) plotted over the corresponding orientation (g  ) in the initial mesh.Thus, Â(g  , ) gives the Lagrangian representation of the current ODF in the initial mesh.If the integrand is continuous, a localized relationship of the following form can be used to update the ODF at any time : During deformation, it is assumed that every orientation is subject to the same macroscopic deformation gradient F. Once the stress  at each orientation is computed using the single crystal model in Section 3, the stress ⟨⟩ at the macroscale is obtained by averaging over the ODF as follows: where σ(g  , ) is the stress ((g  ) = (g  + Δg)) plotted over the corresponding orientation (g  ) in the initial ODF mesh.The average tangent modulus for the macrostructure is obtained in a similar way: ⟨/F⟩ = ∫(σ/F)A(g  ,  = 0)g  .

Numerical Examples
We consider planar polycrystals characterized by a twodimensional rotation R that relates the local crystal lattice frame to the reference sample frame.A parametrization of the associated rotation group is where  is the angle between the crystal and sample axes, E is the two-dimensional alternator ( 11 =  22 = 0;  12 = − 21 = 1), and I is the identity tensor.A general planar crystal with symmetry under rotations through  is considered here.Under the symmetry, crystal orientations can be described uniquely by parameters drawn from a simply connected fundamental region (0, ).Due to symmetry, the orientation 0 is exactly the same as orientation .The crystal reorientation velocity (v = g/) is obtained by taking a derivative of relation (21): where Ω is the spin tensor defined as Ω = Ṙ R  .Here, R  is evaluated through the polar decomposition of the elastic deformation gradient F  as F  = R  U  .The texturing of the polycrystal is tracked by computing the change in orientation Δg = vΔ of crystals at each time step.Using Δg, the Jacobian  can be computed for evolving the ODF in (19).
For modeling the single crystal response, two slip systems at orientations −/6 and +/6 were considered.The values in the elastic stiffness matrix are taken as  11 = 2 GPa and  12 =  44 = 1 GPa.The particular hardening law in (15) was chosen as follows: )  (no sum on ) .The slip system hardening term (ℎ  ) includes latent hardening through parameter .In this term, ℎ   is the hardening coefficient and    is the saturation resistance of slip system .The values for the slip hardening parameters are taken to be identical for both slip systems and are listed below: In order to find the direction of crack growth, we assume that the crack grows perpendicular to the direction of maximum tensile stress motivated by recent microscale experiments in Dunne et al. [7].However, instead of searching for the maximum tensile stress direction in the immediate vicinity of the crack tip, nonlocal criteria are employed [25].This involves finding the angle along which the "averaged" tensile stress over a finite distance ahead of the crack tip is the largest.Since VMM is an element-based method, we search across multiple elements ahead of the crack tip.The search is performed along the same angle but from start points at the center of each element edge encountered along the path as shown in Figure 5. Once the direction of maximum tensile stress is found, all elements along this path whose tensile stress exceeds the critical normal tractions are chosen and the fine scale variable ⟦u⟧ is augmented.Subsequently, both the normal and shear tractions in these elements are controlled using the traction separation laws in (3).The loading displacement is maintained at the same value and the system is reequilibrated.After reequilibration, the direction of crack propagation is again searched from the location of the new crack tip, without increasing the loading displacement.More crack elements are assigned on elements whose tensile stress exceeds the critical normal traction.The next loading step is applied only after no more cracks can form in the previous loading stage.In the case of failure of elastoplastic elements, we assume that once the crack has formed, we assume that the element begins to unload following an elastic unloading curve with the initial elastic modulus as shown in Figure 5(b).

Isotropic Model Tests.
In the numerical examples, we consider the case of a homogeneous material (without explicit representation of grains) to study mesh convergence and to compare with results published in literature.To model the isotropic material, uniform ODF with all orientations equally weighted, corresponding to a random microstructure, is employed.
A tensile linear elastic specimen is first simulated to test the convergence of the approach (in mode I) as shown in Figure 6(a).Plasticity is switched off by assigning a very high critical resolved shear stress for activation of slip systems.A loading displacement in the -direction is imposed on the right edge.Isotropic stiffness matrix constants corresponding to the random ODF are  11 = 2.25 GPa and  12 =  44 = 0.75 GPa.The cohesive strength   0 is taken as 10 MPa and H  is 0.25 GPa/mm.Figure 6 The crack path search method employed involves searching over a finite distance ahead of the crack tip for the angle along which the average tensile stress is the largest.In the implementation, the search is done across several elements ahead of the crack tip in the same direction but from start points at the center of the element edge.(b) After the crack is formed in an elastoplastic element, the unloading is assumed to be elastic following the initial elastic modulus.conditions employed in a specimen of dimension 2 mm by 1 mm (loaded in plane stress).The crack is initiated by specifying a starting element at the bottom edge from where the crack path search is initiated.A loading rate of 0.003 mm per step was employed once the crack starts to propagate.The resulting elastic load-displacement curves for different mesh densities, for unstructured grids with 442, 700, 930, and 2500 elements, respectively, are shown in Figure 6(b).

P, Δ
In order to compare with Rudraraju [25],  and Δ have been normalized with   = 20 MPa and Δ  = 0.2 mm.The result in Figure 6(b) indicates that the load-displacement response does not show any significant mesh dependency for the two cases considered and the result is indeed identical to that published in literature [25].Figures 6(c) and 6(d) show the crack paths and displacement contours for a coarse and a fine mesh.It can be seen that the right half carries all the displacements at this point while the displacement varies along the crack path.
In order to verify that the approach indeed gives the correct crack trajectory for an elastoplastic case, a three-point bending test was employed as shown in Figure 7(a).Plasticity was switched on using an initial critical resolved shear stress of  0 = 20MPa for both slip systems.Two cases were simulated: the first one involves a centrally located loading (Figure 7(a)) and the other involves an eccentric loading at the top surface (Figure 7(b)).A 2 mm by 1 mm specimen was employed and an initial notch 0.4 mm high and 0.02 mm wide was specified at the center bottom surface in order to nucleate the crack.The two constraints at the bottom surface were placed at 0.33 mm from the left and right edges, respectively.
A 2226-element mesh was used in the simulation.In the case of three-point bending, the crack is expected to propagate towards the point of application of force at the top surface of the specimen.As seen from Figures 7(c) and 7(d), the overall crack direction is indeed predicted as moving towards the point of application of the loading in both cases.Note that, for a better visualization of the crack path, the crack path elements are removed from the simulation images as shown in Figures 7(c) and 7(d).

Polycrystals: Tensile Test with Intergranular Failure Mode.
Here, the entire computational region is divided into two levels: macroscale (ODF-based) and microscale.At the microscale, aggregates of grains are explicitly modeled to track propagation of cracks.A Voronoi tessellation approach was used to build a polycrystalline microstructure and random orientations were assigned to each crystal.The mesh used in this work is shown in Figure 8(a) depicting a 0.32 mm by 0.16 mm domain with the color of each grain in the microscale region mapping to a particular crystal orientation from 0 to  radians.There is a gradual transition of mesh density from microscale to the macroscale wherein the farfield elements are much bigger than the microscale elements.A uniform ODF with all orientations equally weighted, corresponding to a random microstructure, is employed in the far-field region.The loading and boundary conditions used in these examples are the same as those shown in Figures 6 and 7 corresponding to a tensile test and a threepoint bending test, respectively.A loading rate of 0.001 mm per step was employed in these simulations.In these simulations, the crack is initiated by specifying a starting element along a grain boundary.The crack is advanced in the microstructural region and the simulation is stopped when the crack reaches the macroscale element.Since VMM is an element-based method and there are no special interface elements involved, the grain boundaries are 0.25 GPa/mm.The critical traction components for the grain interior elements were assigned to be much higher such that the grain boundary elements fail exclusively.The simulation was stopped when the crack reached the macroscale element.The critical normal stress employed is the same for both the elastic and the elastoplastic case.-Δ response for the tensile test case in Figure 8(b) shows that the elastic case failed earlier than the plastic case and the crack grew in two discrete increments.In the elastic case shown in Figure 8(c), the crack propagation continued at the same loading step (displacement corresponding to the first large drop in -Δ response) until the triple junction is reached.At the triple junction, the crack paused until the next loading step wherein the stresses exceed the critical normal traction in one of the grain boundaries leading to the next large drop in -Δ response.The elastoplastic case (Figure 8(d)) follows the identical crack trajectory until the triple junction, beyond which another grain boundary branch is chosen.In the elastoplastic case, the specimen can be loaded to a much larger displacement for the same crack length.The progression of cracking is slower and only a few grain boundaries fail when the critical stress is first reached.Subsequently, the specimen takes up a large displacement loading following an elastoplastic loading curve until all the grain boundaries fail in the final loading step leading to a sharp drop in -Δ response.
VMM approach gives us the ability to simulate arbitrary crack paths within the polycrystal.This ranges from intergranular cracks following the grain boundaries to transgranular cracks following the interior of the grains.In addition, features like crack branching can be simulated.In the above example, when the crack reaches a triple junction along one of the grain boundaries, tensile stresses of elements along the other two grain boundaries are compared.If normal traction along one direction is much higher than those along the other direction, the former boundary is chosen for the crack to advance.Figure 9(a) shows a case where both grain boundaries at a triple junction exceed the critical stress level in which case the method simulates crack branching.In addition, by giving comparable critical stresses to grain boundaries as well as the elements in the grain interior, one can expect a transition to transgranular failure mode as shown in Figure 9

Polycrystals: Three-Point Bending Test with Intergranular
to Transgranular Transition.In this example, a three-point bending test was chosen with the same mesh as in the previous example.Two simulations were performed.In the first case, the critical stresses for the grain interior elements are much higher such that only the grain boundary elements fail.In the second case, elements in the grain interior were given critical traction values comparable with those for the grain boundary elements to test conditions when the crack transitions from intergranular to transgranular mode.A crack in the three-point bending case was initiated by assigning a start element at the bottom edge of the bar.-Δ response for the three-point bending test in the intergranular failure mode is shown in Figure 10(a) and shows more softening and progressive failure response compared to the tensile test case in Figure 8(b).Figures 10(b) and 10(c) show the intergranular failure mode in the elastic and elastoplastic case, respectively.Compared to the elastic case, the elastoplastic failure tends to show a larger fracture toughness and a larger displacement to failure as expected.Crack branching was not observed.
The second case allows transition from intergranular to transgranular crack as illustrated in Figures 11(b) and 11(c).To allow for intragranular cracks, elements in the grain interior were assigned lower values of critical traction components   0 = 40 MPa and   0 = 10 MPa.The softening moduli for normal and tangential directions for these elements were assigned as 2 GPa/mm and 0.5 GPa/mm, respectively.These values are higher than those assigned to elements adjoining grain boundaries, which implies that the crack always starts in the intergranular failure mode.However, as the crack reaches a triple junction, it is seen that the transgranular mode can become favorable.In the transgranular mode of crack propagation, the crack propagated perpendicular to the direction of maximum normal traction.When the transgranular crack reached another grain boundary with an easier propensity to form a crack, the crack transitioned back to intergranular mode.The crack paths in the elastic and elastoplastic cases were different.In the elastic case, the crack showed an increasing tendency to follow the grain boundary.-Δ response in Figure 11(a) again shows that the crack is formed at a slower rate in the elastoplastic case.

Mathematical Problems in Engineering
A comparison of -Δ response of a purely intergranular mode in Figure 8(a) and a combination of inter-and transgranular modes in Figure 11(a) in the elastoplastic case reveals several interesting features.The failure is faster when transgranular paths are available for the crack even though the critical stresses for these paths were higher.Further, the crack tends to form in fewer steps and more abruptly when transgranular modes are available.This aspect is especially seen when a crack reaches a triple junction, where a transgranular mode becomes more favorable for the crack to propagate to the next grain boundary rather than wait for the stresses to exceed the critical tractions in the grain boundaries of the triple junction.

Conclusions and Future Directions
Modeling failure at microstructural scales is valuable for predictive modeling of component life and ensuring structural integrity in aerospace structures.In this paper, variational multiscale method (VMM) is used to model crack propagation in polycrystalline materials.In this approach, a discontinuous displacement field is added to elements that exceed the critical values of normal or tangential traction during loading.This additional degree of freedom is represented within the cracked element using a special discontinuous shape function that ensures that the displacement jump is localized to that particular element.The finite element formulation and code implementation details were presented.Compared to traditional cohesive zone modeling approaches, the method does not require the use of any special interface elements in the microstructure.
Since failure is usually localized at regions of stress concentrators such as notches, cracks, and contact surfaces, a multiscale model is developed in this work which is based on a multiresolution mesh.The model includes an explicit microstructure representation at critical regions, while at regions farther away from the stress concentration, a reduced order model that statistically captures the effect of microstructure is employed.To demonstrate applicability of the methodology, a polycrystalline grain structure was monotonically loaded in tension and three-point bending modes.The primary contribution in this work is to demonstrate the use of the methodology in predicting crack growth paths that are arbitrary (transitions from intergranular to transgranular paths and vice versa).In VMM method, the crack path is computed on the fly during the simulation and thus arbitrary crack paths can be simulated.Both elastic and elastoplastic cases were compared, and the inclusion of crystal plasticity leads to increased fracture toughness and more progressive failure over a larger number of loading steps compared to the elastic case.By controlling the critical normal stress for the grain boundary elements and the elements in the grain interior, transitions from intergranular to transgranular failure modes can be simulated.One aspect to be pursued in a future work is to extend the method towards 3D problems.In addition, experimental identification of the parameters in the cohesive model is of interest.This includes identification of critical tractions, the fracture energy, and criteria for identification of crack growth direction in elastoplastic crystals.In situ fracture experiments within a scanning electron microscope are being targeted towards achieving these goals.

Figure 1 :Figure 2 :
Figure 1: Representation of crack opening ⟦u⟧, domain Ω  , and crack surface Γ   .Crack surface has normal direction n and tangent direction m.Nominal traction on the crack surface is denoted as t  .

Figure 3 :
Figure 3: Shape functions for representing the displacement jump.Two cases are possible by comparing the outward normal direction of the element edge not intersected by the crack and the normal direction of the discontinuity within the element.(a) n and n 1 point in the same direction: the shape function is the difference between  1 and a Heaviside function centered at the crack.(b) n and n 1 point in the opposite directions: the shape function is the difference between (1 −  1 ) and a Heaviside function centered at the crack.

( 23 )Figure 4 :
Figure 4: The configuration of multiscale model.The model is divided into two scales.Microscale contains explicit representation of crystals with assigned orientations in a region where crack is expected to propagate.The macroscale contains a statistical representation of the polycrystal (ODF) at each of its integration points.

Figure 5 :
Figure 5: Aspects of the numerical implementation.(a)The crack path search method employed involves searching over a finite distance ahead of the crack tip for the angle along which the average tensile stress is the largest.In the implementation, the search is done across several elements ahead of the crack tip in the same direction but from start points at the center of the element edge.(b) After the crack is formed in an elastoplastic element, the unloading is assumed to be elastic following the initial elastic modulus.

Figure 8 :
Figure 8: Tensile test with the multiscale model.(a) Boundary conditions, (b) load-displacement response, elastic versus elastoplastic material behavior, (c)   contours for the elastic case, and (d)   contours for the elastoplastic case showing the crack trajectory.

Figure 9 :
Figure 9: (a) Branching of a crack at a triple junction.(b) Intergranular to transgranular crack transition.

Figure 10 :
Figure 10: Results from a three-point bending test using the multiscale model with failure constrained to the intergranular mode.(a) Force displacement curve showing elastic versus elastoplastic material behavior.(b, c)   contours for elastic and elastoplastic response. (b).

Figure 11 :
Figure 11: Results from a three-point bending test using the multiscale model allowing for transgranular mode in addition to the intergranular mode.(a) Force displacement curve showing elastic versus elastoplastic material behavior.(b, c)   contours for elastic and elastoplastic response.