User login

You are here

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

Julian J. Rimoli's picture

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).



[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:


Arash_Yavari's picture

Dear Julian:

Thank you for the excellent post. So, for a problem with a complicated geometry and no known solution (not knowing where the crack will go) what would be a practical solution? Trying a few different meshes and see if there are any agreements between the predicted paths?


Julian J. Rimoli's picture

Dear Arash:

Thank you very much for your comments. Let me give you my understanding on this:

In general, for problems with no known crack path, I would first make sure that the mesh has no structure. If the mesh is structured it will certainly provide preferential (low energy) directions for crack propagation, thus generating a bias in the solution. It is very important to have this aspect in mind when using commercial software packages, since they generally generate meshes with some sort of structure as a result of their meshing technique (advancing front, etc.). For this purpose, K-means, blue-noise, or maximal Poisson sampling meshes would do the job.

The second aspect to take care of is the convergence issue. Meshes have an implicit “roughness” which in general decreases as the mesh is refined, as seen in Figure 3a. Thus, a first step would be to use adaptive refinement near the crack tip (or tips) following some criterion. The effect of this refinement would be to decrease the path deviation ratio locally. The work by Parks et al. (Parks, Paulino, Celes, and Espinha, IJNME, 2011) is a good example of this approach. However, since this “roughness” saturates about 4% error for most meshes, there is only so much you can do through mesh refinement. That is the point where I think Conjugate-Direction meshes have great potential, but that remains to be seen.

That being said, the necessity for lack of structure (randomness) makes the problem a statistical one, so I would not be happy by just looking at one specific mesh for a given problem. I think there is a necessity to look at multiple realizations of the same problem and look at the convergence in a statistical way. This is, in fact, what you are suggesting. You can read, for example, Joe Bishop’s work in this area (Bishop and Strack, IJNME, 2011).

I hope this clarifies your question.



Hi Julian,

Interesting topic and a clear presentation!

Also, I enjoyed your answer to Arash's excellent query above; especially relevant here are these excerpts: "... I would first make sure that the mesh has no structure. If the mesh is structured it will certainly provide preferential (low
energy) directions for crack propagation, thus generating a bias in the
solution." and "the necessity for lack of structure (randomness) makes the problem a statistical one."



Do you know if anyone has defined any statistical parameters with the aim of measuring the structure of a given mesh? And, to a bit deeper, still: what does the term "structure of a mesh," really mean? What kind of quantitative properties or attributes?


Here, I have a feeling that researches involving stereology might have at least some very nice and relevant clues if not already well worked out or sophisticated answers.

Now, speaking about myself, by now, it's been almost two decades since I stopped working in that field. As a matter of fact, I seem to have forgotten all of my stereology. However, back in the early 1990s, I was a PhD student at the University of Alabama at Birmingham, and I do remember my guide Prof. B. R. Patterson showing me an examiner's copy of a PhD thesis written by a student of Prof. Arun Gokhale's; this thesis (at GATech) was concerned, if I recall it right, with stereological characterization of the digitally captured microscopic features of metals and ceramics, perhaps also including the grain boundary network.

Coming back to the present issue of quantitatively characterizing crack-paths, I think stereologists might have worked on it, even if, to them, "mesh" must have meant in the first instance something like: "a network of grain boundaries."     

Thus, I very definitively think that there must be valuable inputs to this question coming from stereology; and so, it would be wonderful if someone could look up the stereological literature, or discuss it with the stereology researchers, and provide the computational engineering community with some kind of a guidance on defining quantitative parameters to measure the structure of a mesh. Also, other things like characterizing the geometry of the crack path surface, and the geometry of its propogation in simulations.


Another point. Apart from stereology, can there be any other ways to measure the structure of the geometric entities like the computational mechanics mesh?

Perhaps, based on the Fourier analysis? Here, I (very vaguely) recall of the Fourier analysis being mentioned during discussions on the last year's Chemistry Nobel (Prof. Schechtman's work on quasicrystals). The Fourier analysis might not yield as detailed---i.e. acute---parameters to capture the various fine aspects concerning the geometrical structure of a mesh. However, given that FFT is computationally relatively fast, with excellent libraries (FFTW) now already available, who knows, characterizations based on the Fourier analysis might be of help to some extent. 

In contrast, the measures based on stereology seem to be all: definitely relevant, perhaps already worked out to a fair degree of completeness/maturity, and practically relevant even for issues and concerns from computational mechanics/engineering.


My 2 cents.


- - - - -


Julian J. Rimoli's picture

Dear Ajit,
Thank you very much for your comments. I am glad you enjoyed the post! Let me try to express my thoughts on the very interesting points you raised:
1-    People working on this area usually look only at the mean value of the path deviation ratio for a given mesh. While useful for convergence analysis, all information on isotropy is lost during the averaging. This is why I believe the polar plots shown above provide a better picture of the structure of the mesh. That is, the polar plot of the path deviation ratio characterizes the structure of the mesh. You can extract quantitative information from it, the most relevant ones being the mean value and the standard deviation.
2-    I am not very familiar with the field of stereology but your points are persuading me to do some literature research in the area to see what can be applied to this field. Thanks!
3-    It is a common practice in computer graphics to look at the Fourier spectrum of a given mesh. This provides useful information about the isotropy of the node distribution. However, this does not provide the full picture for our problem since the way nodes are connected significantly alters the behavior of the path deviation ratio. For example, imagine nodes distributed on a rectangular array. They could be connected in various patterns: square mesh, right triangles with diagonals going from up-left to down-right, right triangles with diagonals going from down-left to up-right, etc. For each of these meshes the polar plot of the path deviation ratio will look different while the power spectrum would be exactly the same.
I hope this clarifies your points. Again, thank you very much for your thoughtful comments!

N. Sukumar's picture


As a follow-up to what is being discussed, couple of questions.  Can convergence in cohesive simulations only be assessed in the statistical sense (eg., Joe's work)?  If quasi-random seed points/nodes are inserted (with some max spacing h between nodes as an input), then in the limit that this spacing h->0, would the crack path converge and in what sense?  Can numerical tests be done to come-up with a converged solution to establish a good benchmark?  It would seem that the number of potential cohesive surfaces (facets of elements) would likely also dictate the crack path, since at every "step" only a finite number of paths are permissible.

Julian J. Rimoli's picture

Dear Suku:

Thank you very much for your follow-up questions. I will do my best to reply to all of them, though altering a little bit the ordering for the sake of clarity.

To the best of my knowledge, the issue of convergence of cohesive models for problems with unknown crack path/paths is still an open problem. I think Joe’s work could be applied to numerically assess the convergence of this particular kind of nonlinear systems. However, even with a converged solution there is no guarantee that it has converged to the right one. As you suggest, performing numerical tests is the only thing we can do for now, and comparing converged solutions to benchmark experimental cases is probably the best way to go.

The fact that cohesive zone modeling provides a limited number of possible crack paths is, from my own perspective, the Achilles’ heel of the cohesive element approach. However, there are not many other options when dealing with problems involving fragmentation or pervasive fracture. The idea of this post was to summarize what people have done (and still are trying to do) to alleviate the undesired consequences of this limitation. For example, your suggestion of using adaptive insertion of nodes with a constraint of minimum distance to existing nodes, is exactly what maximum Poisson sampling does. In this case, saturation of the path deviation ratio to a value larger than one (~1.04 for most meshes) will still occur. That is, there is only so much you can do through this kind of refinement.



Dibakar Datta's picture

 For Beginners in X-FEM ::

I did a small project on X-FEM under the guidance of Prof. Nicolas MOES & Prof. Nicolas Chevaugeon


It may help the beginners in this area. 



Dibakar Datta


sockalsi's picture

Hello all,

I am beginning to explore various methods for modeling crack propagation in a bimaterial interface with no precrack.

In my limited search I have not come across a case where XFEM is used to model a bimaterial interface with no interface. Please let me know if you know any examples/publications on this.



likask's picture


thanks you for the post. This is very good topic, however I am not big fun of this technique when it is applied for homogenous materials. However form my experience it gives great results when its applied for materials with microstructue discretised directly, 


Dog Bone Mesh Dependency 

For more details pleas look here,

Numerical multiscale solution strategy for fracturing heterogeneous materials  


L. Roy Xu's picture

Dear Julian,

When we measured the cohesive (or interfacial) strengths, we used quite large specimens such as width  W=20
mm as suggested by ASTM standards. Indeed, the cohesive zone is quite small in front of a crack.  Do you think the intrinsic  cohesive strength should be larger than the measured one? Thansk for any comments and we plan to conduct some experiments to show the size effect.


rahmani26's picture

dear all;

I need a code matlab for this topic CZM+XFEM, any help 


 I am trying to find out the parameters that can be obtained from XFEM results. By analyzing the XFEM results, what are all possible parameters that can be calculated? I am especially interested in crack propagation of bimaterial interfaces. Any references addressing the same will be helpfull. Thanks.


Subscribe to Comments for "Journal Club Theme of August 2012:  Mesh-Dependence in Cohesive Element Modeling"

Recent comments

More comments


Subscribe to Syndicate