Journal Club Theme of August 2012: Mesh-Dependence in Cohesive Element Modeling

The classical cohesive zone theory of fracture finds its origins in the pioneering works by Dugdale, Barenblatt and Rice [1–3]. In their work, fracture is regarded as a progressive phenomenon in which separation takes place across a cohesive zone ahead of the crack tip and is resisted by cohesive tractions. Cohesive zone models are widely adopted by scientists and engineers perhaps due to their straightforward implementation within the traditional finite element framework. Some of the mainstream technologies proposed to introduce the cohesive theory of fracture into finite element analysis are the eXtended Finite Element Method (X-FEM) and cohesive elements.

Sukumar et al. [4] first utilized the X-FEM for modeling 3D crack growth by adding a discontinuous function and the asymptotic crack tip field to the finite elements. Subsequently, the method was extended to account for cohesive cracks [5]. While the X-FEM approach can deal with arbitrary crack paths, it becomes increasingly complicated for problems involving pervasive fracture and fragmentation.

On the other hand, the cohesive element approach consists on the insertion of cohesive finite elements along the edges or faces of the 2D or 3D mesh correspondingly [6-9]. Even though this approach is well suited for problems involving pre-defined crack directions, a number of known issues affect its accuracy when dealing with simulations including arbitrary crack paths, e.g., problems with the propagation of elastic stress waves (artificial compliance), spurious crack tip speed effects (lift-off), and mesh dependent effects (c.f. [10] for a complete review). However, the robustness of the method makes it one of the most common approaches for pervasive fracture and fragmentation analysis.

Artificial compliance and lift-off effects can be avoided by using an initially rigid cohesive law [8] or, more elegantly, a discontinuous Galerkin formulation with an activation criterion for cohesive elements [11,12]. However, the problem of mesh dependency, more precisely mesh-induced anisotropy and mesh-induced toughness, is an active area of research.

When a mesh is introduced to represent the continuum fracture problem within the cohesive element formulation, a constraint is introduced into the problem due to the inability of the mesh to represent the shape of an arbitrary crack. If we think of an arbitrary crack as a rectifiable path (2D) or surface (3D), the ability of a mesh to represent a straight line (2D) or plane (3D) as the mesh is refined is a necessary condition for the convergence of the cohesive element approach [14]. In 2-dimensional problems, the ability of a mesh to represent a straight segment is characterized by the path deviation ratio, defined as the ratio between the shortest path on the mesh edges connecting two nodes, and the Euclidian distance between them (see Fig. 1a). It is desirable, then, for the path deviation ratio to be independent on the segment direction (to avoid mesh-induced anisotropy) and to tend to one as the mesh size tends to zero (to avoid, in the limit, mesh-induced toughness). A mesh that satisfies these requirements is said to be isoperimetric.

[img_assist | nid=12896 | title=Figure 1 | desc=(a) 4k mesh, (b) 4k mesh with nodal perturbation, (c) K-Means mesh, and (d) Conjugate-Directions mesh. In all meshes, the path deviation ratio is defined as the ratio between the length of the red solid line and the blue dashed line. | link=none | align=center | width=640 | height=480] 

The work of Radin and Sadun [13] shows that the pinwheel tiling of the plane has the isoperimetric property. Based on this result, Papoulia et al. [14] observed that crack paths obtained from pinwheel meshes are more stable as the meshes are refined compared to other types of meshes.  It is worth noting, however, that pinwheel meshes are hard to generate and there is no known extension to the 3-dimensional case. In a recent work, Paulino et al. [15] analyzed the behavior of 4k meshes modified by nodal perturbation and edge swap operators (see Fig. 1b). Their results show that the expected value of the path deviation ratio (over all possible directions) is decreased for a given mesh size. It is worth noting, however, that even though the mesh induced-toughness is reduced in this way, the meshes under consideration still exhibit a considerable anisotropy [16].

Recently, K-Means meshes generated by the application of Delaunay’s triangulation to nodes resulting from clustering random points (see Fig. 1c), have been proposed as a way of alleviating mesh-induced anisotropy while keeping acceptable triangle quality [16]. For reasonable mesh sizes, K-Means meshes exhibit the same mean value of the path deviation ratio as 4k meshes with nodal perturbation while being perfectly isotropic. In the same article, another type of mesh termed Conjugate-Directions mesh is introduced. Conjugate-Directions meshes are generated by the application of barycentric subdivision to K-Means meshes as depicted in Fig. 1d. In this way, the barycentric subdivision adds new directions to the existing K-Means mesh which tend to be orthogonal to the original directions as the K-Means tends to be smoother (i.e., as the K-Means algorithm is applied over a larger number of random points). This can be interpreted as enriching the set of directions provided by the original mesh with new conjugate directions. Consequently, Conjugate-Directions meshes exhibit the same isotropy observed in K-Means meshes while producing a drastically reduced mean value of the path deviation ratio for identical mesh sizes. Figure 2 shows the polar plot of the path deviation ratio vs. mesh direction for 4k, K-Means and Conjugate-Directions meshes.

[img_assist | nid=12897 | title=Figure 2 | desc= Polar plot of the path deviation ratio (minus 1) as a function of the mesh direction for the meshes under consideration. |
link=none | align=center | width=352 | height=264]  

Moreover, preliminary results show that the convergence of K-Means meshes in the sense of the mean value of the path deviation ratio is similar to that of 4k meshes as the mesh size is reduced [17], see Fig 3a. At the same time, the standard deviation decreases at a fastest rate when compared to 4k meshes as shown in Fig. 3b. However, numerical evidence shows that the path deviation ratio tends to saturate around 1.04 for both meshes. On the other hand, the same numerical experiment shows no indication of saturation for Conjugate-Directions meshes in the studied range of mesh-sizes. In summary, K-Means meshes are isotropic thus not providing preferred crack propagation directions. In addition to being isotropic, Conjugate-Directions meshes exhibit a better convergence behavior in the sense of the path deviation ratio, making them good candidates for cohesive element analysis of crack propagation problems where the crack path is not known a priori.

[img_assist | nid=12898 | title=Figure 3 | desc=convergence of the path deviation ratio; (a) its mean value, and (b) its standard deviation. |
link=none | align=center | width=640 | height=480] 

To conclude, we should emphasize that even though convergence in the sense of the path deviation ratio is a necessary condition for convergence of the cohesive element formulation, when dealing with finite meshes an occasional misalignment of an edge in the cohesive crack might be enough to cause the simulated crack path to diverge from the physical one. Needless to say, the issue of crack path convergence in the cohesive element formulation is still an open problem. The hope of many of the previously mentioned research efforts is that arbitrary crack propagation can be achieved through mesh design (pre-processing).

References

[1] Dugdale, D. S., “Yielding of steel sheets containing slits,” Journal of the Mechanics and Physics of Solids, Vol. 8, No. 2, 1960, pp. 100–104.
[2] Barenblatt, G. I., “The Mathematical Theory of Equilibrium Cracks in Brittle Fracture,” Advances in Applied Mechanics, Vol. 7, 1962, pp. 55–129.
[3] Rice, J. R., “Mathematical analysis in the mechanics of fracture,” Fracture: an advanced treatise, Vol. 2, 1968, pp. 191– 311.
[4] Sukumar, N., Moes, N., Moran, B. and Belytschko, T., “Extended finite element method for three-dimensional crack modeling,” International Journal for Numerical Methods in Engineering, Vol. 48, 2000, pp. 1549-1570.
[5] Moes, N. and Belytschko, T., “Extended finite element method for cohesive crack growth,” Engineering Fracture Mechanics, Vol. 69, 2002, pp. 813-833.
[6] Xu, X. P. and Needleman, A., “Numerical simulations of fast crack growth in brittle solids,” Journal of the Mechanics and Physics of Solids, Vol. 42, 1994, pp. 1397–1397.
[7] Xu, X. P. and Needleman, A., “Numerical simulations of dynamic crack growth along an interface,” International Journal of Fracture, Vol. 74, 1995, pp. 289–324.
[8] Camacho, G. T. and Ortiz, M., “Computational modeling of impact damage in brittle materials,” International Journal of Solids and Structures, Vol. 33, 1996, pp. 2899–2938.
[9] Ortiz, M. and Pandolfi, A., “Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis,” International Journal for Numerical Methods in Engineering, Vol. 44, No. 9, 1999, pp. 1267–1282.
[10] Seagraves, A. and Radovitzky, R., Advances in cohesive zone modeling of dynamic fracture, Springer Verlag, 2009.
[11] Noels, L. and Radovitzky, R., “An explicit discontinuous Galerkin method for non-linear solid dynamics: Formulation, parallel implementation and scalability properties,” International Journal for Numerical Methods in Engineering, Vol. 74, 2008, pp. 1393-1420.
[12] Radovitzky, R., Seagraves, A., Tupek, M., and Noels, L., “A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method,” Computer Methods in Applied Mechanics and Engineering, Vol. 200, No. 1-4, 2011, pp. 326-344.
[13] Radin, C. and Sadun, L., “The isoperimetric problem for pinwheel tilings,” Communications in Mathematical Physics, Vol. 177, No. 1, 1996, pp. 255–263.
[14] Papoulia, K.D., Vavasis, S.A., and Ganguly, P., “Spatial convergence of crack nucleation using a cohesive finite-element model on a pinwheel-based mesh,” International Journal for Numerical Methods in Engineering, Vol. 67, No. 1, 2006, pp. 1-16.
[15] Paulino, G. H., Park, K., Celes, W., and Espinha, R., “Adaptive dynamic cohesive fracture simulation using nodal perturbation and edge-swap operators,” International Journal for Numerical Methods in Engineering, Vol. 84, No. 11, 2010, pp. 1303–1343.
[16] Rimoli, J. J.,  Rojas, J. J., and Khemani, F. N., “On the mesh dependency of cohesive zone models for crack propagation analysis,” in 53rd AIAA Structures, Structural Dynamics, and Materials and Conference, Honolulu, HI, 2012.
[17] Rimoli, J. J. and  Rojas, J. J., “Meshing strategies for the alleviation of mesh-induced effects in cohesive element models,” submitted. Preprint: http://arxiv.org/abs/1302.1161