Abstract
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, the method does not require the use of any special interface elements in the microstructure and thus can model 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.
1. 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 state-of-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 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–15]. Cohesive element-based approaches have been previously used for modeling intergranular crack paths [16–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.
2. 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 (). Let refer to the total deformation gradient in the reference configuration () at time with respect to the initial undeformed configuration () at time . Similarly, the total deformation gradient in the current configuration () at time is written as . Then 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 . 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 in the finite element mesh. In the variational multiscale approach, an additional degree of freedom is introduced within cracked elements such that the overall corrected displacement in the cracked element is . The finite element equations for the crack element are described next.
2.1. Model for the Cohesive Zone
Traction continuity on the crack surface is given bywhere 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 in the reference configuration at time is expressed in terms of the Cauchy stress in the current configuration at time as .
The traction components are assumed to be linearly related to the components of the crack displacement in the reference configuration via the softening moduli :where , are the softening moduli for the mode I and mode II crack and , 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 for the crack element is obtained as
Linearizing the residual with respect to the displacement unknown gives
Here, 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 using the crystal plasticity constitutive model as shown in Appendix B. Note that the increment in deformation gradient is of the form which can be substituted in (6) such that (5) is obtained in terms of the unknown displacement increments. Note that 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).
2.2. 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 due to infinitesimal changes in the displacements is computed asIn the variational multiscale method, the test function is written in terms of virtual displacements ; that is, . In addition, the increment in the PK-I stress can be written in an implicit form using (6). Thus, the linearized equation reduces to
2.3. Coupled Solution Procedure
Let be the displacement field at the th iteration; then the linearized equations are solved for the increments . The analysis procedure involves the solution of two coupled equations for the increments in displacements and .
Let , , and denote the gradients in terms of finite element matrices ( and ) and the nodal unknowns. Each NR step involves solving for increments , using the system of equations given below:wherewhere is a matrix computed by rearranging the component of the crack normal . In our implementation, for example, where and are the components of the crack normal .
2.4. Finite Element Representation of Cracks
Our present implementation employs three-node triangular elements. The displacement is interpolated using the usual linear shape functions for the triangle element. The displacement discontinuity 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), and point in the same direction and the shape function is the difference between and a Heaviside function. In the second case in Figure 3(b), and point in the opposite directions and the shape function is the difference between and a Heaviside function.
(a)
(b)
When using this description, the displacement gradient (written in a vector format) is computed using the standard finite element shape functions aswhere is a 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 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 . If the shape function defined in Figure 3(b) is used, a negative sign is employed. In addition, the crack normal () 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.
3. 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 slip systems defined by orthonormal vector pairs () denoting the slip direction and slip plane normal, respectively, at time , the constitutive equations relate the following basic fields (all quantities expressed in crystal lattice coordinate frame): the deformation gradient defined with respect to the initial undeformed crystal which can be decomposed into elastic and plastic parts as (with ), the Cauchy stress , and the slip resistances . In the constitutive equations to be defined below, the Green elastic strain measure defined on the relaxed configuration (plastically deformed, unstressed configuration ) is utilized. The conjugate stress measure is then defined as .
The constitutive relation, for stress, is given by , where 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 bywhere 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 .
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 . There is no plastic shearing rate () 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 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 . The conjugate stress measure is then computed from 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 for use in the weak form is computed using a fully implicit algorithm described in Appendix B.
3.1. 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, , gives the probability density of finding an orientation in the microstructure. The ODF satisfies the following conservation equation at all times during deformation:where 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) aswhere the ODF at time , , represents the probability density of crystals with orientation and represents the volume element in the undeformed (initial) ODF mesh which becomes volume element at time . A Jacobian gives the ratio of elemental volumes such that . Here, is the reorientation gradient given as 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 ) can be made:The quantity written as is the volume density plotted over the corresponding orientation () in the initial mesh. Thus, 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 . 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 is the stress () plotted over the corresponding orientation () in the initial ODF mesh. The average tangent modulus for the macrostructure is obtained in a similar way: .
4. Numerical Examples
We consider planar polycrystals characterized by a two-dimensional rotation that relates the local crystal lattice frame to the reference sample frame. A parametrization of the associated rotation group iswhere is the angle between the crystal and sample axes, is the two-dimensional alternator (), and 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 . Due to symmetry, the orientation is exactly the same as orientation . The crystal reorientation velocity () is obtained by taking a derivative of relation (21):where is the spin tensor defined as . Here, is evaluated through the polar decomposition of the elastic deformation gradient as . The texturing of the polycrystal is tracked by computing the change in orientation of crystals at each time step. Using , the Jacobian can be computed for evolving the ODF in (19).
For modeling the single crystal response, two slip systems at orientations and were considered. The values in the elastic stiffness matrix are taken as GPa and GPa. The particular hardening law in (15) was chosen as follows: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 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).
(a)
(b)
4.1. 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 GPa and GPa. The cohesive strength is taken as 10 MPa and is 0.25 GPa/mm. Figure 6(a) shows the boundary 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). In order to compare with Rudraraju [25], and have been normalized with MPa and 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.
(a)
(b)
(c)
(d)
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 MPa 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).
(a)
(b)
(c)
(d)
4.2. 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 to radians. There is a gradual transition of mesh density from microscale to the macroscale wherein the far-field 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 three-point bending test, respectively. A loading rate of 0.001 mm per step was employed in these simulations.
(a)
(b)
(c)
(d)
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 represented by choosing elements to one side of the grain boundary with the crack normals assigned to be parallel to the grain boundary normal. The values of critical traction components for the grain boundary elements and were picked as 30 MPa and 5 MPa and softening moduli for normal and tangential directions are 1.5 GPa/mm and 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(b).
(a)
(b)
4.3. 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.
(a)
(b)
(c)
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 MPa and 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.
(a)
(b)
(c)
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.
5. 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.
Appendix
A. Constitutive Update Scheme
The quantities at the current time step are denoted by subscript . The deformation gradient is known at the current time step and the update procedure given here is used to compute PK-I stress , where is the constitutive model described in Section 3.
An Euler-backward time integration procedure for (14) leads to the following approximation:Substituting (A.1) into the multiplicative decomposition results in the following:where is the trial elastic deformation gradient and is given as .
The Green elastic strain measure is computed using (A.2) aswhere and are defined as
Using (A.3) in the constitutive relation for stress leads to the following:where .
A trial resolved shear stress is then computed. A potentially active set of slip systems can be identified based on the trial resolved stress as the systems with .
During plastic flow, the active systems are assumed to follow the consistency condition: . Increment in shearing rates at each time step is obtained by solving the following equation obtained by resolving (A.5) along slip directions:where .
A system of equations is obtained of the following form:whereIf for any system , then this system is removed from the set of potentially active systems. The system is repeatedly solved until for all systems .
B. Implicit Computation of Tangent Moduli
The linearization process of Cauchy stress is given by
The above expression requires the evaluation of and using the constitutive model in terms of . In order to obtain , we consider the linearization of (A.5) to obtainThis computation of requires the evaluation of , obtained by linearization given byUsing the definition of , the above set of equations yield an implicit form for use in (B.2). Next, is obtained from
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
The computations have been carried out as part of research supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under Award no. DE-SC0008637 that funds the PRedictive Integrated Structural Materials Science (PRISMS) Center at University of Michigan.