Finite element-based model for crack propagation in polycrystalline materials

Finite element-based model for crack propagation in polycrystalline materials^*

N. Sukumar^I; D.J. Srolovitz^II

^IDepartment of Civil and Environmental Engineering University of California, Davis, CA 95616 E-mail: [e-mail]

^IIPrinceton Materials Institute, Bowen Hall Department of Mechanical and Aerospace Engineering Princeton University, Princeton, NJ 08544 E-mail: [e-mail]

ABSTRACT

In this paper, we use an extended form of the finite element method to study failure in polycrystalline microstructures. Quasi-static crack propagation is conducted using the extended finite element method (X-FEM) and microstructures are simulated using a kinetic Monte Carlo Potts algorithm. In the X-FEM, the framework of partition of unity is used to enrich the classical finite element approximation with a discontinuous function and the two-dimensional asymptotic crack-tip fields. This enables the domain to be modeled by finite elements without explicitly meshing the crack surfaces, and hence crack growth simulations can be carried out without the need for remeshing. First, the convergence of the method for crack problems is studied and its rate of convergence is established. Microstructural calculations are carried out on a regular lattice and a constrained Delaunay triangulation algorithm is used to mesh the microstructure. Fracture properties of the grain boundaries are assumed to be distinct from that of the grain interior, and the maximum energy release rate criterion is invoked to study the competition between intergranular and transgranular modes of crack growth.

Mathematical subject classification: 74R10.

Key words: microstructure, grain boundaries, crack discontinuity, partition of unity, extended finite element method, kinetic Monte Carlo, convergence.

1 Introduction

Microstructural features in polycrystalline materials play an important role in determining the failure mechanisms as well as the macroscopic mechanical response in these class of materials. In monolithic, polycrystalline materials, tailored microstructures can induce crack bridging and kinking, leading to improved toughness and failure resistance [1]. The role and influence of microstructural features (e.g., grain size and shape, and grain boundary characteristics) on fracture is well-established [2]. Experimental evidence has demonstrated that the fracture resistance of polycrystals can be enhanced by increasing the fraction of so-called special boundaries and by modifying the grain boundary network and topology [3-5]. However, theoretical and numerical models have not yet emerged that can fully explain these experimental findings. In this paper, we use a finite element-based mesoscale fracture model that accounts for the grain structure and the differences in the critical fracture energy of the grain boundaries vis-à-vis that of the grain interior.

Failure modeling in disordered (heterogeneous) materials has been approached using lattice (spring-network) models [6-10] Voronoi cell-based finite element method [11], and cohesive surface formulations within finite elements [12]. Fracture strength in disordered materials is governed by weakest-link statistics [13]. Spring-network models are intuitive in nature, and can be easily implemented; however, it is difficult to obtain both elastic homogeneity and grid-insensitive crack propagation directions on random lattice networks [14]. Discrete crack growth modeling with finite elements is the prevailing standard. The modeling of arbitrary geometries with varying material properties render finite elements as a powerful computational tool. Even though crack growth modeling with remeshing has reached a mature stage of development [15,16], it is not readily amenable to a microscopic description of failure modeling and hence has not received wide attention.

Recently, there has been significant progress in the direction of discrete crack growth modeling - the development of a partition of unity finite element method for crack modeling (coined as the extended finite element method or X-FEM) [17,18] has provided an accurate and robust numerical method that removes the need to mesh the crack surfaces in static or quasi-static crack growth simulations. In the X-FEM, the framework of partition of unity [19] is used to enrich the classical displacement-based finite element approximation with a discontinuous function (Heaviside step function) and the two-dimensional asymptotic crack-tip fields. This provides a means to model the crack independently of the underlying finite element mesh. As opposed to previous developments in modeling strong (displacement) discontinuities within finite elements [20], some of the noteworthy advantages of the X-FEM are: a single-field variational principle is used with no incompatibilities introduced between elements; the symmetry and sparsity of the stiffness matrix are retained; and the crack discontinuity can be totally arbitrary with respect to the finite element mesh.

It is now well-recognized that failure modeling at the microstructural (mesoscopic) scale in heterogeneous materials is an appropriate and potentially fruitful path [21]. In this study, the incorporation of microstructure within a continuum-based partition of unity finite element method represents a significant point of departure from prior finite element models for fracture. The topology of polycrystalline materials and in particular the grain structure needs to be accounted for in any realistic failure modeling endeavor. The size of the grain (e.g., grain diameter) provides a natural length-scale that is embedded within the mesoscopic fracture model. Such an approach can lead to better integration of experiment and numerical modeling towards the development of a useful tool for fracture simulation in complex heterogeneous materials.

In [22], the computational algorithm for brittle fracture simulations in polycrystalline microstructures was presented; here, the same algorithm is used to elucidate some of the key features of the crack growth model. In Section 2, the extended finite element method is introduced, and in Section 3, computer simulation of polycrystalline microstructures and the microstructural meshing algorithm are summarized. The weak form and discrete equations used in the X-FEM are presented in Section 4. In Section 5.1, a convergence study of the X-FEM for an edge crack problem is carried out. The crack growth procedure is outlined in Section 5.2, and numerical simulation results are presented in Section 5.3. The main conclusions obtained from this study are given in Section 6.

2 Extended finite element method

The partition of unity finite element method [19] is a generalization of the standard Galerkin finite element method. In the literature, numerical techniques such as the extended finite element method (X-FEM) [17,18], generalized finite element method [23], or the element partition method [24] are all particular instances of the partition of unity method. In the X-FEM, the emphasis has been on modeling discontinuities (such as cracks) with minimal enrichment. Even though the idea of adding special functions to the finite element approximation is not new [25], unlike previous attempts, the partition of unity framework satisfies a few important properties which makes it a powerful tool for local enrichment within a finite element setting:

(1) can incorporate application-specific basis functions to better approximate the solution;

(2) automatic enforcement of continuity (conforming trial and test approximations); and

(3) point or line singularities in 2-d, and surface discontinuities in 3-d can be modeled without the need for the discontinuous surfaces to be aligned with the finite element mesh.

In [17], the linear elastic near-tip asymptotic fields were used as enrichment functions, whereas in [26], the asymptotic near-tip fields for a bimaterial interface crack were adopted.

We summarize some of the essential concepts related to 2-d crack modeling in isotropic media [17]. For further details on the X-FEM implementation, the interested reader can refer to References [17] and [27]. Consider a single crack in 2-dimensions. Let G_c be the crack surface (interior) and L_c the crack tip - the closure _c = G_cÈL_c. The enriched displacement (trial and test) approximation for 2-d crack modeling is of the form [17]:

where u_I is the nodal displacement vector associated with the continuous part of the finite element solution, a_I is the nodal enriched degree of freedom vector associated with the Heaviside function H (assumes the value +1 above the crack and -1 below the crack), and is the nodal enriched degree of freedom vector associated with the elastic asymptotic crack-tip functions. In the above equation, is the set of all nodes in the mesh; _G is the set of nodes whose shape function support is cut by the crack interior G_c; and _L is the set of nodes whose shape function support is cut by the crack tip L_c (_GÇ_L = Æ):

The crack-tip enrichment functions in isotropic elasticity are:

where (r,q) are polar coordinates in the local crack-tip coordinate system. In Fig. 1a, the nodes enriched by the Heaviside function (open circles) and the crack-tip functions (filled circles) are shown for an edge crack, whereas in Fig. 1b, the shape function support w_I as well as the local coordinate system for the crack-tip enrichment functions are illustrated.

In Reference [17], the algorithm for the selection of nodes for enrichment is presented. When a crack is nearly coincident with a finite element edge, care must be taken so as to avoid linear dependencies between basis functions; a tolerance criterion is used in the algorithm. If the crack intersects the finite elements, then element partitioning is also carried out; a higher-order (6 ×6) quadrature rule is used in partitioned elements for accurate numerical integration of the weak form. The need for partitioning and its distinction from remeshing is discussed in greater detail in [27]. The above steps ensure that the linear system of equations is well-conditioned.

3 Polycrystalline microstructure

In order to simulate quasi-static crack propagation in a polycrystalline material, a realistic microstructure was first produced using the Potts model for grain growth [28]. In this model, a continuum microstructure is mapped onto a regular two-dimensional square lattice containing N sites. Each lattice site is assigned a number s_i, which corresponds to the orientation of the grain in which it is embedded. The number of distinct grain orientations (spins) is Q. Lattice sites that are adjacent to neighboring sites having different grain orientations are regarded as being adjacent to a grain boundary, whereas a site surrounded by sites with the same grain orientation is in the grain interior. The grain boundary energy is specified by associating a positive energy with grain boundary sites and zero energy for sites in the grain interior, in accordance with the Potts Hamiltonian:

where J is a constant proportional to the grain boundary energy per unit length, and d_ij is the Kronecker delta. In Eq. (5), the summation on i is over all the sites in the lattice, whereas that on j is over the first and second nearest neighbors nn(i) of site i. The kinetics of the boundary motion are simulated via a zero-temperature Monte Carlo technique in which a lattice site is selected at random and its orientation is randomly changed to one of the other grain orientations. The change in energy associated with the change in orientation is then evaluated. If the change in energy is less than or equal to zero, the reorientation is accepted; if the energy is raised, the reorientation is rejected. Microstructures are produced by initially assigning a random value s_i of the grain orientation (1 < s_i < Q) to each site. Time is measured in units of Monte Carlo steps: one MCS corresponds to N attempted changes, with the time increment Dt = 1/N MCS after every reorientation. The Monte Carlo procedure is executed until the desired grain size is produced.

In the crack propagation simulations, the initial finite element mesh is based on the microstructure generated by the Potts model. In the literature, mesh generation algorithms for material microstructures are available [29,30]. In order to perform crack propagation simulations in 2-d, we required a Delaunay triangulation meshing scheme in which the polycrystalline material is represented by distinct grains and the grain boundaries are associated with one-dimensional segments (edges of the finite elements). To meet our needs, we developed a Delaunay algorithm to mesh the polycrystalline microstructure [22]. The input to the meshing algorithm is a polycrystalline microstructure produced by the Potts model, with known spins s_i (1 < s_i < Q, 1 < i < N). A grain boundary conforming finite element mesh is constructed using a Delaunay triangulation algorithm [31], with a provision for local mesh refinement. Mesh refinement is based on a comparison between the actual local length scale l (e.g., element circum-radius) and the desired length scale specified by a scalar variable r called the length density function.

4 Governing equations

4.1 Strong form

Consider a body W Ì

², with boundary G. The boundary G = G_uÈG_t

, where G_u is the essential boundary, G_t the natural boundary, and are the internal cracks. The field equations of elastostatics in the absence of body forces are:

where Ñ_s is the symmetric gradient operator, u is the displacement vector, e is the small strain tensor, s is the Cauchy stress tensor, and C is the tensor of elastic moduli for a homogeneous isotropic material. The essential and natural boundary conditions are:

where n is the unit outward normal to W, and are prescribed displacements and tractions, respectively, and m is the number of cracks. Note that the third equality in Eq. (7) imposes the condition that the crack is traction-free.

4.2 Weak form and discrete equations

The weak form (principle of virtual work) for linear elastostatics is:

where u^h Î U^h and du^h Î are the approximating trial and test functions used in the X-FEM, and d is the first variation operator. The space U^h is the enriched finite element space that satisfies the essential boundary conditions, and which include basis functions that are discontinuous across the crack surfaces. The space is the corresponding space with homogeneous essential boundary conditions. The finite element domain is W^h and the traction boundary conditions are imposed on ¶, which is a subset of ¶W^h.

The trial function u^h and the test function du^h used in the X-FEM are of the form given in Eq. (1). On substituting the trial and test functions in the above equation, and using the arbitrariness of nodal variations, the following discrete system of linear equations is obtained:

where d is the vector of nodal unknowns, and K and f are the global stiffness matrix and external force vector, respectively [22].

5 Numerical results

We first test the X-FEM on a benchmark problem to check the accuracy of the stress intensity factors (SIFs) and to establish the rate of convergence of the method. The convergence of the numerical method is studied by imposing the exact near-tip displacement field on the boundary of an elastic plate that contains a crack that extends from the perimeter to the center of the specimen. Then, we describe the simulation procedure for crack growth through a microstructure and present crack growth simulations for varying grain boundary toughness.

5.1 Convergence study

Consider an elastic plate of dimensions (-1,1) ×(-1,1) with a crack of unit length that extends from (-1,0) to (0,0). The Cartesian components of the near-tip displacement field [2] corresponding to K_I = 1 and K_II = 1 are imposed on the boundary of the specimen:

where m = E/2(1+n) is the shear modulus, E is the Young's modulus, and n is the Poisson's ratio of the material. The constant k = 3 - 4 n under plane strain conditions. The mixed mode SIFs are computed using the domain form [32] of the interaction integral [33]. The radius of the domain r_d = r_kh_e, where h_e is the mesh spacing and r_k is a scalar multiple. All elements that lie within a radius of r_d from the crack-tip are selected to form the 2-d domain that is used in the domain integral computations; for further details on the SIF computations in the X-FEM, see [17].

The stress intensity factors and the relative error in the energy norm are computed as the mesh is refined. The exact energy norm and the error in the energy norm are defined as:

Four different meshes are considered: 10×10, 20×20, 40 ×40, and 160 ×160; a sample mesh (20 ×20 elements) is shown in Fig. 2a. To impose the essential boundary conditions, the crack is explicitly meshed over one element (AB), whereas the remaining part of the crack (BC) is modeled by the X-FEM (Fig. 2a). The scaling factor r_k = 4 is used in the domain integral computations. The results of the convergence study are listed in Table 1 and in Fig. 2b, the relative error in the energy norm is plotted as a function of the mesh spacing h_e (log-log plot). In Fig. 2, the rate of convergence R is also indicated. The numerically computed SIFs are in agreement with the exact solution, and a rate of convergence of one-half is realized, which matches the theoretical convergence rate of the finite element method in the presence of a dominant -singularity [25].

5.2 Simulation procedure

Polycrystalline microstructures are obtained using the Potts grain growth model outlined in Section 3, and a grain boundary conforming finite element mesh of the microstructure is constructed. The problem domain is a square of edge length L. An initial pre-crack of size a (typically 0.02L) is introduced along a grain boundary that emanates from x₁ = 0.5L on the top surface. The top and bottom surfaces are traction-free; uniaxial strain is applied in the x₁-direction by fixing the left edge and imposing displacement boundary conditions on the right edge (Fig. 3).

Let us denote the critical fracture energy of a grain boundary by and that of the grain interior by . The validity of Griffith-Irwin fracture mechanics is not lost at the microstructural level [2], and hence one can use the classical notion of mechanical energy release rate for crack growth. The energy release rate under plane strain conditions is related to the stress intensity factors through Irwin's relation:

The crack will propagate along a grain boundary or through the grain interior depending on which has the larger value of / (k is either gb or i). In the grain interior, the crack is assumed to propagate preferentially in the the maximum hoop (circumferential) stress direction (K_II = 0) [34]. The crack growth procedure is summarized below:

1. Pre-crack the sample and apply load (strain-controlled)

2. If the crack is within a grain, then

• Determine in the maximum hoop stress direction

• If > , propagate crack in that direction

• Move crack-tip a small fraction of grain size or up to a boundary

3. If the crack is along a grain boundary, then

• Determine in the maximum hoop stress direction and in bothdirections along the grain boundary

• Move the crack-tip in the direction of maximum > (k = gb or i) if any are larger than unity

• Move the crack-tip a small fraction of the grain size or grain edge length, without crossing a boundary or junction

4. Increment load and go to step 2

5.3 Crack growth simulations

A 100 ×100 lattice with Q = 100 is used to generate the microstructure. The kinetic Monte Carlo algorithm is executed until a time of 10⁴ MCS. In all the crack growth simulations, plain strain conditions are assumed with Young's modulus E = 10⁵ and Poisson's ratio n = 0.3.

We first mesh a given microstructure using two different values of the length density function r: with increasing r, smaller grains tend to be merged with larger grains. To ensure the presence of all the grains in the mesh, r should be selected to be a fraction of the average grain diameter d. Next, we study the influence of the mesh size on the crack path. For a fixed ratio / = 0.4, we conducted crack growth simulations for different values of r (Fig. 5). The results reveal that for r < 0.03, there is little bias in the crack path. Lastly, we present crack growth simulations using r = 0.03 for varying grain boundary toughness values: / = 0.4, 0.5, 0.6, 0.8. The simulation results are illustrated in Fig. 6 with the percentage of the crack path that is intergranular indicated within braces. As the toughness of the grain boundary is increased, the crack path tends to be transgranular-dominated; for / = 0.6, the crack path is mixed mode.

6 Conclusions

In this paper, we first established that the rate of convergence of the extended finite element method (X-FEM) for crack problems was one-half , which matches the theoretical convergence rate of the FEM in the presence of a dominant -singularity [25]. Subsequently, through numerical simulations, we demonstrated the promise and potential of the X-FEM for crack growth studies in polycrystalline materials.

Acknowledgments

The authors thank Dr. Timothy Baker for developing the microstructural meshing algorithm and Professor Mark Miodownik for providing the code for the Potts grain growth model. Helpful discussions with Professor Jean Prévost are also acknowledged.

Received: 14/IV03. Accepted: 16/IX03.

#571/03.

• [1] A.G. Evans and K.T. Faber, The crack growth resistance of microcracking brittle materials. Journal of the American Ceramic Society, 67 (4) (1984), 255-260.
• [2] B. Lawn, Fracture of Brittle Solids Cambridge University Press, Cambridge, England, second edition, (1993).
• [3] T. Watanabe, The impact of grain boundary character distribution on fracture in polycrystals. Materials Science and Engineering A, 176 (1994), 39-49.
• [4] M. Kumar, W.E. King and A.J. Schwartz, Modifications to the microstructural topology in f.c.c. materials through thermomechanical processing. Acta Materialia, 48 (9) (2000), 2081-2091.
• [5] C.A. Schuh, M. Kumar and W.E. King, Analysis of grain boundary networks and their evolution during grain boundary engineering. Acta Materialia, 51 (3) (2003), 5687-700.
• [6] P.D. Beale and D.J. Srolovitz, Elastic fracture in random materials. Physical Review B, 37 (10) (1988), 5500-5507.
• [7] W.A. Curtin and H. Scher, Brittle fracture of disordered materials. Journal of Materials Research, 5 (3) (1990), 535-553.
• [8] H.J. Herrmann and S. Roux, editors. Statistical Models for the Fracture of Disordered Media North-Holland, Amsterdam, The Netherlands, (1990).
• [9] W.H. Yang, D.J. Srolovitz, G.N. Hassold and M.P. Anderson, Microstructural effects in the fracture of brittle materials. In M.P. Anderson and A.D. Rollett, editors, Simulation and Theory of Evolving Microstructures, pages 277-284, The Metallurgical Society, Warrendale, PA, (1990).
• [10] J.E. Bolander Jr. and S. Saito, Fracture analyses using spring networks with random geometry. Engineering Fracture Mechanics, 61 (1998), 569-591.
• [11] S. Ghosh, K. Lee and P. Raghavan, A multi-level computational model for multiscale damage analysis in composite and porous materials. International Journal of Solids and Structures, 38 (2001), 2335-2385.
• [12] P.D. Zavattieri and H.D. Espinosa, Grain level analysis of crack initiation and propagation in brittle materials. Acta Materialia, 49 (2001), 4291-4311.
• [13] W.A. Curtin, Size scaling of strength in heterogeneous materials. Physical Review Letters, 80 (7) (1998), 1445-1448.
• [14] A. Jagota and S.J. Bennison, Spring-network and finite-element models for elasticity and fracture. In K.K. Bardhan, B.K. Chakrabarti and A. Hansen, editors, Nonlinearity and Breakdown in Soft Condensed Matter. (Springer Lecture Notes in Physics 437), pages 186-201, Springer, Berlin, (1994).
• [15] L.F. Martha, P.A. Wawrzynek and A.R. Ingraffea, Arbitrary crack representation using solid modeling. Engineering with Computers, 9 (1993), 63-82.
• [16] J.B.C. Neto, P.A. Wawrzynek, M.T.M. Carvalho, L.F. Martha and A.R. Ingraffea, An algorithm for three-dimensional mesh generation for arbitrary regions with cracks. Engineering with Computers, 17 (2001), 75-91.
• [17] N. Moës, J. Dolbow and T. Belytschko, A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46 (1) (1999), 131-150.
• [18] C. Daux, N. Moës, J. Dolbow, N. Sukumar and T. Belytschko, Arbitrary cracks and holes with the extended finite element method. International Journal for Numerical Methods in Engineering, 48 (12) (2000), 1741-1760.
• [19] J.M. Melenk and I. Babuska, The partition of unite finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 139 (1996), 289-314.
• [20] J.C. Simo, J. Oliver and F. Armero, An analysis of strong discontinuities induced bystrain softening in rate-independent inelastics solids. Computational Mechanics, 12 (1993), 277-296.
• [21] A.M. Stoneham and J.H. Harding, Not too big, not too small: The appropriate scale. Nature Materials, 2 (2) (2003), 77-83.
• [22] N. Sukumar, D.J. Srolovitz, T.J. Baker and J.-H. Prévost, Brittle fracture in polycrystalline microstructures with the extended finite element method. International Journal for Numerical Methods in Engineering, 56 (14) (2003), 2015-2037.
• [23] T. Strouboulis, K. Copps and I. Babuska, The generalized finite element method. Computer Methods in Applied Mechanics and Engineering, 190 (32-33) (2001), 4081-4193.
• [24] C.A. Duarte, O.N. Hamzeh, T.J. Liszka and W.W. Tworzydlo, The element partition method for the simulation of three-dimensional dynamic crack propagation. Computer Methods in Applied Mechanics and Engineering, 119 (15-17) (2001), 2227-2262.
• [25] G. Strang and G. Fix, An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, N.J., 1973.
• [26] N. Sukumar, Z. Huang, J.-H. Prévost and Z. Suo, Partition of unity enrichment for bimaterial interface cracks. International Journal for Numerical Methods in Engineering, 59 (8) (2004), 1075-1102.
• [27] N. Sukumar and J.-H. Prévost, Modeling quasi-static crack growth with the extended finite element method. Part I: Computer implementation. International Journal of Solids and Structures, 40 (26) (2003), 7513-7537.
• [28] D.J. Srolovitz, M.P. Anderson, G.S. Grest and P.S. Sahni, Grain growth in two dimensions. Scripta Metallurgica, 17 (1983), 241-246.
• [29] W.C. Carter, S.A. Langer and E.R. Fuller Jr., Object-Oriented Finite Element Analysis (OOF) Available at http://www.ctcms.nist.gov/oof, National Institute for Standards and Technology, Gaithersburg, MD, (1998).
• [30] S. Ghosh and S.N. Mikhopadhyay, A material based finite-element analysis of heterogeneous media involving Dirichlet tessellations. Computer Methods in Applied Mechanics and Engineering, 104 (2) (1993), 211-247.
• [31] T.J. Baker, Triangulations, mesh generation and point placement strategies. In D.A. Caughey and M.M. Hafez, editors, Frontiers of Computational Fluid Dynamics, page 101-115, New York, NY, 1994. John Wiley & Sons.
• [32] B. Moran and C.F. Shih, Crack tip and associated domain integrals from momentum and energy balance. Engineering Fracture Mechanics, 27 (6) (1987), 615-641.
• [33] J.F. Yau, S.S. Wang and H.T. Corten, A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. Journal of Applied Mechanics, 47 (1980), 335-341.
• [34] F. Ergodan and G.C. Sih, On the crack extension in plates under plane loading and transverse shear. Journal of Basic Engineering, 85 (1963), 519-527.

*

Revised submission, July 2003.

Publication Dates

• Publication in this collection
  11 July 2005
• Date of issue
  Dec 2004

History

• Received
  14 Apr 2003
• Accepted
  16 Sept 2003