Cohesive Elements for Crack Propagation

I conducted research into finite-element modeling of cracks using cohesive elements as part of my masters thesis. This research area was of interest to Sandia National Labs for the design/maintenance of nuclear reactors. I was advised by Dr. Julian Rimoli of the Daniel Guggenheim School of Aerospace Engineering at the Georgia Institute of Technology. On this page I describe the study, the reasons for undertaking it, the approach taken, and the findings. 


In finite element analysis (FEA), crack propagation is often modeled using cohesive elements, which are well-suited to modeling situations/phenomena such as adhesives, delamination, and fracture (although the modeling strategy depends greatly on the type of effect being modeled). Cohesive elements are placed between the structural elements of the mesh, and propagation of the crack is simulated by sequential damage and subsequent removal of adjoining cohesive elements from the simulation. 

One of the problems with using cohesive elements for crack propagation is that the cohesive elements must lie along the crack path, or conversely, the crack propagates along a path defined by the cohesive elements since the cohesive elements in essence represent the crack. If prior knowledge of the crack path is available (either from experimental testing or other sources), the cohesive elements should be placed along this known path, and the model can then be used to obtain useful information such as whether or not the crack propagates for the given loading conditions, and the distance it travels. However in situations where the crack path is previously unknown, modeling the crack becomes a lot more challenging. One solution is to place cohesive elements between all the structural elements in a region of the mesh through which the crack is likely to propagate, thereby providing the crack with an opportunity to change direction after each element length. However even with the multiple path options this approach provides, the mesh topology still has a significant impact on the path of the simulated crack and can therefore adversely affect the results. 

This type of crack propagation modeling using cohesive elements was of interest to Sandia National Labs for their modeling of fracture and failure in nuclear reactor walls, and this research was conducted at their request. The focus of the study was on how to generate topologies which minimize this dependence of crack path on mesh topology.

Mesh Topologies

At a higher level, there are two types of meshes that can be considered - structured meshes and randomized meshes. 

Structured meshes, which have repeating patterns, provide the same set of choices at each interval. It is possible, if lucky, that one of the choices is aligned with the true direction of the crack. In such a case the model will have a series of consecutive cohesive elements aligned with the crack, and the crack path will be smooth and accurate. However it is more likely that none of the path options will reasonably match the true direction, and the structured geometry will impose a mesh bias on the crack propagation simulation. Refining the mesh is unlikely to change this. 

Randomized meshes on the hand provide a different set of direction choices at each interval. In such a mesh it is unlikely that a series of consecutive cohesive elements will be found all aligned in the same direction. So it will generally not be possible to obtain a smooth crack interface (although mesh refinement will help). Also it is unlikely that the majority of the cohesive elements will be aligned with the true crack direction. However it is equally unlikely that a series of adjacent cohesive elements will all be aligned in the wrong direction, hence, unlike the structured mesh, the biasing is not limited to only one or two directions.

The various mesh toplogies considered as part of this research effort are described below.

1. Default Abaqus mesh

Abaqus was used to generate a triangular mesh using the standard settings of its meshing tool. Mesh refinement was obtained in the region around the crack by seeding the beam boundaries more finely in the vicinity of the crack, while using a larger seed size for the rest of the boundary.

Fig. Default Abaqus mesh: (left) Entire domain, (right) zoomed in on mesh

The meshes generated by Abaqus tend to have uniform patterns over large regions. While it is possible to further modify the mesh within Abaqus/CAE, this default mesh was considered as it is likely to be encountered by engineers using the software, and might give an insight into how accurate such crack propagation simulations are. 

2. 4K mesh

This type of mesh consists of perfectly square shaped cells, divided into 4 triangular elements. Since the cells are all identical, a repeating pattern is obtained all through the mesh domain. Mesh refinement is easily achieved by subdividing the mesh into rows of different sizes, where the cell sizes decrease by ½ compared to the previous row for each mesh refinement. A transition layer separates a coarser row from a finer one. 

Fig. 4k mesh with one level of refinement: (left) entire domain, (right) zoomed in 
Fig. 4k mesh with 4 levels of refinement (zoomed in)

It is evident that the direction choices available to a crack are exactly the same at every interval. Hence the biasing in this mesh to the crack path is expected to be quite heavy.

3. Randomized nodes with K-means and Delaunay triangulation

Nodes are placed at random locations, and then a Delaunay triangulation is applied to create this type of mesh. The k-means algorithm is applied to the nodes so although they are generated at random locations, a certain type of order can be enforced to ensure they are evenly distributed, and the elements are not excessively distorted.

Fig. K-means clustered Delaunay mesh- (left) entire domain, (right) zoomed in

It can be seen that the crack encounters different sets of path choices at each interval. It is generally possible to chart a path in any direction without requiring sharp changes in crack angle at any point.

4. Barycentric mesh

This is an extension of the previous randomized Delaunay mesh. Here each Delaunay triangle is subdivided using a Barycentric technique

Fig. Barycentric subdivided k-means Delaunay mesh (zoomed in)

As with the previous mesh, this barycentric mesh allows a crack to travel in any direction by only making gradual changes in angle as compared with a structured mesh.

Approach to Testing

1. Setup

The setup was modeled on the experiments of Galvez et al. [1] to study mixed mode (Mode I/II) crack propagation in concrete. It consists of a notched beam subject to asymmetric bending which causes the crack to propagate through the beam. The same setup was subsequently used by Areias [2] and Gurses [3]. However instead of cement the material used for this study was alumina, a brittle material with well known material properties. Alumina 99.9% (Al2O3) has a Youngs Modulus of 370GPa, a Poissons Ratio of 0.22, and a fracture toughness of 4.00 MPa-m1/2 [4],[5] & [6]. In addition I set the spring constant k=0, meaning that we do not have a support at that point. Also the applied force was replaced by an applied displacement, because the models with displacement control were found to be much more likely to converge than those with applied force.

Fig. Geometry, loading and boundary conditions of the notched specimen subject to asymmetric bending (Fig.11 from Gurses [6])

It was experimentally determined that the crack trajectory would be an angle of approximately 20 degrees for k=0. This experimental result was used to determine the accuracy of the predicted crack path for each mesh topology is accurate. 

2. Modeling Strategy

In order to generate the test mesh topologies I wrote a MATLAB script which created the coordinate and connectivity matrices for the elements. These meshes could potentially be imported into Abaqus as orphan meshes; however inserting cohesive elements in the correct locations of an orphan mesh is very difficult in Abaqus/CAE (at least when there is no initial crack or it is very thin) hence I wrote an additional MATLAB script to accept a mesh definition and insert cohesive elements in the correct locations.

The way the MATLAB script works is it separates the elements that share an edge, and then cohesive elements (which are rectangular in shape) are placed in this gap. These cohesive elements share two of their edges (the longer sides) with elements that originally shared an edge. The cohesive elements have virtually zero thickness, and therefore the two nodes making up each of the shorter sides have the same coordinates. New coordinate and connectivity matrices are created to account for the presence of both shell elements (triangular) and cohesive elements (rectangular) between them. 

Initial versions of this script placed cohesive elements between all the triangular elements of the mesh. However, due to the large number of degrees of freedom these simulations took a significant amount of time to run, hence the script was later modified to insert cohesive elements only within a specified rectangular region, thus giving the ability to place cohesive elements only in a small region through which the crack is likely to propagate. This reduced simulation time by approximately 70%.

Fig. (left) Sample mesh (right) Same mesh with cohesive elements inserted in rectangular region

Furthermore, once the new coordinate and connectivity matrices are created, the script generates a complete input (.inp) file representing the entire model - mesh, material and section definitions, loads, boundary conditions, analysis steps, field/history output requests, etc. This input file is ready to run using the Abaqus solver.

3. Cohesive Zone Modeling

Abaqus provides different types of cohesive elements, each of which are better suited for a particular type of analysis. The ones chosen for this study were 4-node, two dimensional cohesive element. For the constitutive response of the cohesive elements, traction-separation was chosen as it is appropriate for bonded interfaces where the interface thickness is negligibly small. 

According to the Abaqus/CAE User’s Manual [7], it assumes a linear elastic traction-separation law prior to damage, and the damage process involves progressive degradation of the material stiffness to represent failure of the elements. The failure mechanism consists of a damage initiation criterion, a damage evolution law, and a choice of element removal or deletion upon reaching a completely damaged state [7]. The initial response of the cohesive elements is assumed to be linear.

Fig. Typical traction-separation response (Fig 31.5.6-1 of the Abaqus Analysis User’s Manual v6.12).

This required me to define the initial elasticity, the damage initiation, and the damage evolution in terms of energy dissipated for failure.

The initial stiffness (the linear response slope) does not represent a physically measurable quantity and is a penalty parameter. Ideally it should be infinite, but realistically it should be set to a very high number. This number if set too high causes Abaqus to experience unsuccessful iterations and cut back on step size so a balance must be achieved through trial and error. I used an initial guess for this value of 107 x E x 57/135000 where E is the Young’s modulus of the material of the beam (alumina), and the ratio 57/135000 is the one used by Alfano and Crisfield [8] to calculate a starting value for their model. 

Damage initiation occurs when the stresses and/or strains satisfy a specified damage initiation criteria. Abaqus allows the analyst to specify one of four available damage initiation criteria – maximum nominal stress criterion, maximum nominal strain criterion, quadratic nominal stress criterion, and quadratic nominal strain criterion. I used the quadratic nominal stress criterion (QUADS) in the models, which implies that damage initiates when a quadratic interaction function involving the nominal stress ratios reaches a value of one. In order to figure out the value of the QUADS damage to use, I used a calibration mechanism; I built a calibration model and held the value of the stiffness constant and varied the value of the QUADS damage until the results matched real world expectations (the crack propagated when the stress concentration factor was equal to the known critical stress concentration factor for that material).

With the damage evolution law, the analyst describes the rate at which material stiffness is degraded once the initiation criterion is reached. Abaqus maintains a scalar damage variable D, which represents the overall damage in the material based on the combined effects of all the active mechanisms. Its value is initially 0, and ranges from 0 to 1 upon further loading after the initiation of damage. This value can be viewed during post-processing by requesting SDEG as one of the field output variables before running the simulation. Damage evolution can be defined in two ways – either by specifying the effective displacement at complete failure relative to the effective displacement at the initiation of damage, or by specifying the energy dissipated due to failure. In the models I chose to specify the latter. 

The following animation shows a few cohesive elements being damaged and subsequently removed from the view. The plot is of the field output variable SDEG which ranges from 0 (undamaged) to 1 (damaged).

4. Loading

For the loading I used displacement boundary conditions since force/pressure loads required very small solver increments and even then it was much more difficult to get converged solutions.

5. Generating the various meshes

A different approach was used for each of the meshes generated.

a. Default Abaqus Mesh

For the default Abaqus mesh the notched specimen was modeled in Abaqus/CAE and meshed in the Mesh module. Edge seeds were specified on the boundaries of the specimen in such a way that the mesh is finer in the region of the crack and coarser at the ends of the specimen.

Once the mesh was generated, a job was defined, and Abaqus was instructed to output an input file. This input file (.inp) contained, among other information, the coordinate and connectivity matrices. The coordinate matrix appears on the lines following the keyword “*Node”, and the connectivity matrix appears on the lines following the keyword “*Element”.

b. 4k mesh

A MATLAB script was written to generate the 4K mesh. The script used a number of loops to generate the nodal coordinates and connectivity matrices for each row, accounting for the fact that the number of elements changes with each row refinement.

The script plotted the generated mesh on the screen for display purposes. It also generated a text file with the coordinate and connectivity matrices, which was used as the input to the main script that adds the cohesive elements to the mesh.

c. Randomized nodes with K-means and Delaunay triangulation

A MATLAB script was written to generate this mesh as well. It generated equidistant nodes at points along the boundary of the notched specimen, as well as along the notch. It then then generated nodes at random locations in the interior region of the specimen. The script ensured that nodes were not generated too close to the boundaries, as this results in skewed elements near the edges when the Delaunay triangulation is performed at a later stage. Algorithms determined the distance required based upon the desired mesh size.

Truly random distribution of the nodes results in clustering in some regions. This results in a mesh that has many skewed elements. In order to add some order to the randomness a k-means clustering algorithm was used. This algorithm assigned all  the points/nodes into separate clusters. For each of these clusters the script calculated the centroid, and used this as a new random point. Depending on how many cluster points are desired for finding the centroids, initially a larger number of random nodes needed to be generated. For example, if we would like the notched specimen to have 100 internal nodes, and we wish to create clusters of 40 nodes and find the centroid of each for our k-means clustered nodes, then we need to initially generate 100 x 40 = 4000 random nodes. The higher the number of cluster points, the better the distribution of k-means points.

Initially MATLAB’s built in kmeans() function was used for clustering. However it was found to be too slow especially when a large number of points is used. Hence a parallel k-means data clustering code written in C by Liao [8] was used. Three versions/implementations of this code exist – sequential execution, parallel execution using OpenMP, and parallel execution using MPI. I compiled and implemented the version that uses OpenMP. OpenMP is a specification for parallel programming, and compilers such as GCC on linux, and MinGW (a gcc port) on Windows, possess the libraries required for compiling such a code. The ‘-fopenmp’ option must be used during the compile process to ensure OpenMP support. In addition this software code also calculates the centroids of the clusters so the final clustered points can be obtained directly. My MATLAB script was modified to generate an input file for this software, call it’s executable, and then read the output files produced in order to extract the necessary point coordinates. Execution time depended on the number of random points used. For a mesh size of 1, I needed 100,000 nodes in the specimen, and using 40 points for each cluster necessitated an input of 100,000 x 40 = 4 million nodes. On a PC with an intel i7 920 CPU (2.67 Hz with 4 physical cores, 8 virtual cores), and 8 GB of  RAM, the parallel k-means clustering took approximately 17 hours.

Once the k-means clustered nodes were obtained, a Delaunay triangulation was used to ‘connect the dots’ and form the triangular elements. A Delaunay triangulation of a set of points/vertices is a triangulation with the property that no vertex falls in the interior of the circumcircle of any of the triangles in the triangulation. 

Initially attempts were made to use MATLAB’s built-in delaunay() function. However it was found that this function only works for convex domains. If the domain has a hole or a notch in it, the algorithm does not recognize these boundaries, and places triangular elements inside the hole/notch as well. No method was found to circumvent this problem using MATLAB’s own routines.

Fig. Example - generating a k-means clustered Delaunay mesh with 40 internal nodes, using 50 nodes per cluster. (left) 40 x 50 = 2000 randomly generated internal nodes, (middle) result of k-means clustering, (right) Delaunay triangulation

Instead a software package known as Triangle developed by Shewchuk [9] was used. Triangle is a 2-D mesh generator and Delaunay triangulator. When providing the input to Triangle, it is possible to specify the boundaries of the model, and the algorithm respects these. The MATLAB script generates the input file for Triangle, calls the Triangle executable, and then reads the Triangle output files to extract the necessary information for the new connectivity matrix. The script then displays a plot of the obtained Delaunay triangulation. 

The main MATLAB script which adds cohesive elements to the mesh and generates a ready-to-run input file was then used to complete the mesh generation and model setup process.

d. Barycentric mesh

For this mesh, a k-means clustered Delaunay triangulation was first obtained using the procedure described in the previous section. A barycentric subdivision was then performed on each of the Delaunay triangles. This involved drawing a line from each vertex to the opposite side of the triangle, thus dividing each Delaunay triangle into 6 smaller triangles. 

The mesh size after the barycentric division is smaller than the mesh size after the Delaunay triangulation, hence if one is trying to generate a Delaunay mesh and a Barycentric mesh of the same size for comparison, one must use less nodes for the barycentric case, so coarser Delaunay triangles are first obtained, and these upon subdivision meet the desired size requirement.

6. Testing

Numerous test cases were run for a plate with an edge crack comparable to the setup previously described, and then compared with the simulation results. An example of this is shown in the figure. The plate was fixed at the bottom, and a displacement boundary condition was applied at the top pulling it upward. An initial edge crack is present on the left edge of the plate although this is hard to see in the figure on the left. The figure on the right is the end result of the simulation, where the crack has traveled across the plate cutting it in half.

The animated gif image shows the edge crack propagating across the plate. (This mesh is quite coarse, finer meshes are used for most of our analyses).


Two sets of simulations were run, one for a mesh size of 2 and another for a mesh size of 1. To put this in perspective, the specimen has dimensions of 675 x 150 units. For each mesh size, one simulation was run for the abaqus mesh and k4 mesh. 3 simulations each were set up for the k-means clustered Delaunay and the barycentric meshes, for each mesh size (since these meshes, being random, will look different every time they are generated). For the mesh size of 2, the three Delaunay and three barycentric simulations have been completed. For the mesh size of 1, only one Delaunay and one barycentric simulation has been run, with two more of each ready and scheduled for execution.

Fig. Default Abaqus mesh with size 2.
Fig. Default Abaqus mesh with size 1.
Fig. k4 mesh with size 2.
Fig. k4 mesh with size 1.
Fig. k-means clustered Delaunay mesh with size 2 (first random mesh generated).
Fig. k-means clustered Delaunay mesh with size 2 (second random mesh generated).
Fig. k-means clustered Delaunay mesh with size 2 (third random mesh generated).
Fig. k-means clustered Delaunay mesh with size 1 (first random mesh generated).
Fig. Barycentric mesh with size 2 (first random mesh generated)
Fig. Barycentric mesh with size 2 (second random mesh generated)
Fig. Barycentric mesh with size 2 (third random mesh generated)
Fig. Barycentric mesh with size 1 (first random mesh generated)

The initial propagation angles of the cracks are summarized in the table.

Abaqus mesh, size 265
Abaqus mesh, size 165
K4 mesh, size 245
K4 mesh, size 145
Delaunay, size 2 (#1)35
Delaunay, size 2 (#2)35
Delaunay, size 2 (#3)45
Delaunay, size 1 (#1)35
Barycentric, size 2 (#1)40
Barycentric, size 2 (#2)20
Barycentric, size 2 (#3)25
Barycentric, size 1 (#1)40

As mentioned previously, the experimentally observed crack angle was approximately 20 degrees.


Based on this set of tests, it is observed that results are more accurate for mesh topologies which exhibit a certain level of randomness. Such topologies provide the crack with a different set of directions to choose from at every interval, as compared to more uniform structured meshes with repeating patterns where the same set of path choices are available after every element length. In addition, it is observed that randomized meshes are likely to outperform uniform meshes even if the latter has a higher level of refinement. In fact refinement of a uniform mesh may yield no improvement at all, whereas refinement of a randomized mesh further improves its crack path prediction. 

Keep in mind that this was a short research study, and the number of tests run here is quite small. Further testing (i.e., a larger sample size) would help bolster these results.


  1. Galvez, J. C., Elices, M., Guinea, G. V., & Planas, J. (1999). Mixed mode fracture of concrete under proportional and nonproportional loading.International Journal of Fracture,94(3), 267-284.
  2. Areias, P. M. A. and Belytschko, T. (2005), Analysis of three-dimensional crack initiation and propagation using the extended finite element method. Int. J. Numer. Meth. Engng., 63: 760–788. doi: 10.1002/nme.1305
  3. Ercan Gürses, Christian Miehe, A computational framework of three-dimensional configurational-force-driven brittle crack propagation, Computer Methods in Applied Mechanics and Engineering, Volume 198, Issues 15–16, 15 March 2009, Pages 1413-1428, ISSN 0045-7825, 10.1016/j.cma.2008.12.028.
  4. ASM Engineered Materials Reference Book, Second Edition, Michael Bauccio, Ed. ASM International, Materials Park, OH, 1994.
  5. Stone: Properties, Durability in Man's Environment, Second Edition, E. M. Winkler, Springer-Verlag, New York, 1975.
  6. Oxides and Hydroxides of Aluminum, Alcoa Technical Paper No. 19, Revised, Karl Wefers and Chanakya Misra, Alcoa Laboratories, 1987
  7. Abaqus/CAE User’s Manual – version 6.11 (2011). Dassault Systemes Simulia Corporation, Providence, R.I.
  9. Jonathan Richard Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, in "Applied Computational Geometry: Towards Geometric Engineering'' (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of Lecture Notes in Computer Science, pages 203-222, Springer-Verlag, Berlin, May 1996. (From the First ACM Workshop on Applied Computational Geometry.)


This website uses cookies to deliver services, improve usability, and measure performance. By continuing to use this site you opt-in to receive these cookies. You may disable some of them on the Cookie Settings page. You also acknowledge that you have read and understand our Cookie Policy, Privacy Policy, and Terms of Service.