A Facet Strain Tensor Calculation

Discontinuous Cell Method (DCM) for the Simulation of Cohesive Fracture and Fragmentation of Continuous Media


In this paper, the Discontinuous Cell Method (DCM) is formulated with the objective of simulating cohesive fracture propagation and fragmentation in homogeneous solids without issues relevant to excessive mesh deformation typical of available Finite Element formulations. DCM discretizes solids by using the Delaunay triangulation and its associated Voronoi tessellation giving rise to a system of discrete cells interacting through shared facets. For each Voronoi cell, the displacement field is approximated on the basis of rigid body kinematics, which is used to compute a strain vector at the centroid of the Voronoi facets. Such strain vector is demonstrated to be the projection of the strain tensor at that location. At the same point stress tractions are computed through vectorial constitutive equations derived on the basis of classical continuum tensorial theories. Results of analysis of a cantilever beam are used to perform convergence studies and comparison with classical finite element formulations in the elastic regime. Furthermore, cohesive fracture and fragmentation of homogeneous solids are studied under quasi-static and dynamic loading conditions. The mesh dependency problem, typically encountered upon adopting softening constitutive equations, is tackled through the crack band approach. This study demonstrates the capabilities of DCM by solving multiple benchmark problems relevant to cohesive crack propagation. The simulations show that DCM can handle effectively a wide range of problems from the simulation of a single propagating fracture to crack branching and fragmentation.

cohesive fracture, finite elements, discrete models, delaunay triangulation, voronoi tessellation, fragmentation

Center for Sustainable Engineering of Geological and Infrastructure Materials (SEGIM)

[0.1in] Department of Civil and Environmental Engineering

[0.1in] McCormick School of Engineering and Applied Science

[0.1in] Evanston, Illinois 60208, USA


[0.5in] Gianluca Cusatis, Roozbeh Rezakhani, Edward A. Schauffert

[0.75in] SEGIM INTERNAL REPORT No. 16-08/587D


Submitted to Journal of Engineering Fracture Mechanics August 2016

1 Introduction

A quantitative investigation of cohesive fracture propagation necessitates an accurate description of various fracture phenomena including: crack initiation; propagation along complex three-dimensional paths; interaction and coalescence of distributed multi-cracks into localized continuous cracks; and interaction of fractured/unfractured material. The classical Finite Element (FE) method, although it has been used with some success to address some of these aspects, is inherently incapable of modeling the displacement discontinuities associated with fracture. To address this issue, advanced computational technologies have been developed in the recent past. First, the embedded discontinuity methods (EDMs) were proposed to handle displacement discontinuity within finite elements. In these methods the crack is represented by a narrow band of high strain, which is embedded in the element and can be arbitrarily aligned. Many different EDM formulations can be found in the literature and a comprehensive comparative study of these formulations appears in Reference (1). The most common drawbacks of EDM formulations are stress locking (spurious stress transfer between the crack surfaces), inconsistency between the stress acting on the crack surface and the stress in the adjacent material bulk, and mesh sensitivity (crack path depending upon mesh alignment and refinement).

A method that does not experience stress locking and reduces mesh sensitivity is the extended finite element method (X-FEM). X-FEM, first introduced by Belytschko & Black (2), exploits the partition of unity property of FE shape functions. This property enables discontinuous terms to be incorporated locally in the displacement field without the need of topology changes in the initial uncracked mesh. Moës et al. (3) enhanced the primary work of Belytschko et al. (2) through including a discontinuous enrichment function to represent displacement jump across the crack faces away from the crack tip. X-FEM has been successfully applied to a wide variety of problems. Dolbow et. al. (4) applied XFEM to the simulation of growing discontinuity in Mindlin-Reissner plates by employing appropriate asymptotic crack-tip enrichment functions. Belytschko and coworkers (5) modeled evolution of arbitrary discontinuities in classical finite elements, in which discontinuity branching and intersection modeling are handled by the virtue of adding proper terms to the related finite element displacement shape functions. Furthermore, they studied crack initiation and propagation under dynamic loading condition and used a criterion based on the loss of hyperbolicity of the underlying continuum problem (6). Zi et Al. (7) extended X-FEM to the simulation of cohesive crack propagation. The main drawbacks of X-FEM are that the implementation into existing FE codes is not straightforward, the insertion of additional degrees of freedoms is required on-the-fly to describe the discontinuous enrichment, and complex quadrature routines are necessary to integrate discontinuous integrands.

Another approach widely used for the simulation of cohesive fracture is based on the adoption of cohesive zero-thickness finite elements located at the interface between the usual finite elements that discretize the body of interest (10); (11). This method, even if its implementation is relatively simple, tends to be computationally intensive because of the large number of nodes that are needed to allow fracturing at each element interface. Furthermore, in the elastic phase the zero-thickness finite elements require the definition of an artificial penalty stiffness to ensure inter-element compatibility. This stiffness usually deteriorates the accuracy and rate of convergence of the numerical solution and it may cause numerical instability. To avoid this problem, algorithms have been proposed in the literature (12) for the dynamic insertion of cohesive fractures into FE meshes. The dynamic insertion works reasonably well in high speed dynamic applications but is not adequate for quasi-static applications and leads to inaccurate stress calculations along the crack path.

An attractive alternative to the aforementioned approaches is the adoption of discrete models (particle and lattice models), which replace the continuum a priori by a system of rigid particles that interact by means of linear/nonlinear springs or by a grid of beam-type elements. These models were first developed to describe the behavior of particulate materials (13) and to solve elastic problems in the pre-computers era (14). Later, they have been adapted to simulate fracture and failure of quasi-brittle materials in both two (15) and three dimensional problems (16); (17); (18); (19). In this class of models, it is worth mentioning the rigid-body-spring model developed by Bolander and collaborators, which dicretizes the material domain using Voronoi diagrams with random geometry, interconnected by zero-size springs, to simulate cohesive fracture in two and three dimensional problems (20); (21); (22); (23). Various other discrete models, in the form of either lattice or particle models, have been quite successful recently in simulating concrete materials (24); (25); (26); (27); (28); (29).

Discrete models can realistically simulate fracture propagation and fragmentation without suffering from the aforementioned typical drawbacks of other computational technologies. The effectiveness and the robustness of the method are ensured by the fact that: a) their kinematics naturally handle displacement discontinuities; b) the crack opening at a certain point depends upon the displacements of only two nodes of the mesh; c) the constitutive law for the fracturing behavior is vectorial; d) remeshing of the material domain or inclusion of additional degrees of freedom during the fracture propagation process is not necessary. Despite these advantages the general adoption of these methods to simulate fracture propagation in continuous media has been quite limited because of various drawbacks in the uncracked phase, including: 1) the stiffness of the springs is defined through a heuristic (trial-and-error) characterization; 2) various elastic phenomena, e.g. Poisson’s effect, cannot be reproduced exactly; 3) the convergence of the numerical scheme to the continuum solution cannot be proved; 4) amalgamation with classical tensorial constitutive laws is not possible; and 5) spurious numerical heterogeneity of the response (not related to the internal structure of the material) is inherently associated with these methods if simply used as discretization techniques for continuum problems.

The Discontinuous Cell Method (DCM) presented in this paper provides a framework unifying discrete models and continuum based methods. The Delaunay triangulation is employed to discretize the solid domain into triangular elements, the Voronoi tessellation is then used to build a set of discrete polyhedral cells whose kinematics is described through rigid body motion typical of discrete models. Tonti (30) presented a somewhat similar approach to discretize the material domain and to compute the finite element nodal forces using dual cell geometries. Furthermore, the DCM formulation is similar to that of the discontinuous Galerkin method which has primarily been applied in the past to the solution of fluid dynamics problems, but has also been extended to the study of elasticity (31). Recently, discontinuous Galerkin approaches have also been used for the study of fracture mechanics (32) and cohesive fracture propagation (33). The DCM formulation can be considered as a discontinuous Galerkin approach which utilizes piecewise constant shape functions. Another interesting feature of DCM is that the formulation includes rotational degrees of freedom. Researchers have attempted to introduce rotational degrees of freedom to classical finite elements by considering special form of displacement functions along each element edge to improve their performance in bending problems (34); (35). This strategy leads often to zero energy deformation modes and to singular element stiffness matrix even if the rigid body motions are constrained. DCM formulation simply incorporates nodal rotational degrees of freedom, without suffering from the aforementioned problem.

2 Governing Equations

Equilibrium, compatibility, and constitutive laws for Cauchy continua can be formulated as limit case of the governing equations for Cosserat continua in which displacements and rotations are assumed to be independent fields but the couple stress tensor is identically zero (36). For small deformations and for any position vector in the material domain , one has




In the above equations, the summation rule of repeated indices applies; and are displacement and rotation fields, respectively. is the strain tensor; is the stress tensor; are body forces per unit volume; is the mass density. Subscripts , , and represent components of Cartesian coordinate system which can be in three dimensional problems; is the Levi-Civita permutation symbol. Considering any position dependent field variable such as , represent partial derivative of with respect to the th component of the coordinate system, while is the time derivative of the variable. The partial differential equations above need to be complemented by appropriate boundary conditions that can either involve displacements, on (essential boundary conditions); or tractions, on (natural boundary conditions); where is the boundary of the solid volume .

In the elastic regime, the constitutive equations can be written as


where is the volumetric strain; and are the volumetric and deviatoric moduli that can be expressed through Young’s modulus and Poisson’s ratio : ; . It is worth observing that since the solution of the problem formulated above requires the stress tensor to be symmetric (see second equation in Equation 2), the constitutive equations imply the symmetry of the strain tensor as well, which, in turn leads to displacements and rotations to be related through the following expression:

The weak form of the equilibrium equations can be obtained through the Principle of Virtual Work (PVW) as


where , and are arbitrary strains, displacements, and rotations, respectively, satisfying compatibility equations with homogeneous essential boundary conditions. It must be observed here that the PVW in Equation 4 is the weak formulation of both linear and angular momentum balances. Hence, the symmetry of the stress tensor and, consequently, the symmetry of the strain tensor are imposed in average sense. This is a significant difference compared to classical formulations for Cauchy continua in which the symmetry of the stress tensor is assumed “a priori”.

3 Discontinuous Cell Method Approximation

3.1 Domain Discretization

Let us consider a three-dimensional primal cell complex, which, according to the customary terminology in algebraic topology (37), is a subdivision of the three-dimensional space through sets of vertices (0-cells), edges (1-cells), faces (2-cells), and volumes (3-cells). Next let us construct a dual cell complex anchored to the primal. This can be achieved, for example, by associating a primal 3-cell with a dual 0-cell , a primal 2-cell with a dual 1-cell, etc. The primal/dual complex obtained through the Delaunay triangulation of a set of points and its associated Voronoi tessellation is a very popular choice in many fields of study for its ability to discretize complex geometry and it is adopted in this study.

Let us consider a material domain and discretize it into tetrahedral elements by using the centroidal Delaunay tetrahedralization, and the associated Voronoi tessellation which leads to a system of polyhedral cells (38). Figure 1a, shows a typical tetrahedral element with the volume , external boundary , and oriented surfaces located within the volume. The interior oriented surfaces are derived from the Voronoi tessellation and are hereinafter called “facets”. In 3D, the facets are triangular areas of contact between adjacent polyhedral cells. In the Voronoi tessellation procedure, the triangular facets are perpendicular to the element edges , which is a crucial feature of DCM formulation as it will be shown later, and their geometry is such that one node of each facet is placed in the middle of the tetrahedral element edge, one is located on one of the triangular faces of the tetrahedral element, and one is located inside the tetrahedral element. As a result, each tetrahedral element contains twelve facets in a 3D setting, Figure 1a. Figure 1b illustrates a portion of the tetrahedral element associated with one of its four nodes and the corresponding facets. Combining such portions from all the tetrahedral elements connected to the same node, one obtains the corresponding Voronoi cell. Each node in the 3D DCM formulation has six degrees of freedom, three translational and three rotational, which are shown in Figure 1b. The same figure depicts, for a generic facet, three unit vectors, one normal and two tangential ones and , defining a local system of reference. In the rest of the paper, the facet index is dropped when possible to simplify notation.

Figure 1: (a) Three dimensional Delaunay tetrahedralization and Voronoi tesselation. (b) Tetrahedron portion associated with node . c) Voronoi cell.

3.2 Discretized Kinematics

The DCM approximation is based on the assumptions that displacement and rotation fields can be approximated by the rigid body kinematics of each Voronoi cell, that is


where , are displacements and rotations of node ; is the volume of the cell associated with node . Obviously with this approximating displacement and rotation functions, strain versus displacement/rotation relationships in Equations 1 cannot be enforced locally – as typically done in displacement based finite element formulation.

Let us consider a generic node of spatial coordinates and the adjacent nodes of spatial coordinates , where is the vector connecting the two nodes. One can write in which is a unit vector in the direction of . Note that, to simplify notation the two indices and are dropped when supposed to appear together. Without loss of generality, let us also assume that node is located on the negative side of the facets whereas node is located on the positive side of the associated facet oriented through its normal unit vector . Moreover, it is useful to introduce the vector connecting node to a generic point P on the facet. The displacement jump at point P reads


where and are the values of displacements on the positive and negative side of the facet, respectively; ; ; ; and , are magnitude and direction of vector . Equation 6 can be rewritten in tensorial notation as


By expanding the displacement and the rotation fields in Taylor series around and by truncating the displacement to the second order and the rotation to the first, one obtains


By projecting the displacement jump in the direction orthogonal to the facet and dividing it by the element edge length , one can write


where is the strain tensor. At convergence the discretization size tends to zero () and one can write . Furthermore, if the dual complex adopted for the volume discretization is such that the facets are orthogonal to the element edges – this condition is verified, for example, by the Delaunay-Voronoi complex – then, , and one has : the normal component of the displacement jump normalized with the element edge length represents the projection of the strain tensor onto the facet. Similarly, it can be shown that the components of the displacement jump tangential to the facets can be expressed as and .

Before moving forward, a few observations are in order. Since the facet are flat the unit vector is the same for any point belonging to a given facet and the projection of the strain tensor is uniform over each facet for a uniform strain field. The variation of the displacement jump over the facet is due to the curvature and it is an high order effect that can be neglected (see the last term in Equation 9). Based on the previous observation one can conclude that the analysis of the interaction of two adjacent nodes can be based on the average displacement jump which, given the linear distribution of the jump, can be calculated as the displacement jump, , at the centroid of the facet C. This leads naturally to the following definition of “facet strains”:




Equations 10 and 11 show that the “facet strains” correspond to the projection of the strain tensor onto the facet local system of reference.

Let us now consider a uniform hydrostatic stress/strain state in one element and . In this case the tractions on each facet must correspond to the volumetric stress and energetic consistency requires that which gives


where is the set of facets belonging to one element. By using Equation 12 and Equations 10, 11, facet deviatoric strains can also be calculated as .

By introducing Equation 5 into Equations 10 and 11 one can also write


where , , , , , , and , are vectors connecting the facet centroid C with nodes , , respectively.

3.3 Discretized Equilibrium Equations

For the adopted discretized kinematics in which all deformability is concentrated at the facets, the PVW in Equation 4 can be rewritten as


where is the set of all facets in the domain.

By introducing Equations 5, 13, 14, and 15 into Equation 16 and considering displacement and traction continuity at the inter-element interfaces, one can write


where , is the generic Voronoi cell, and are the volume and the set of facets of the cell , respectively. Note that the first term on the LHS of Equation 17 also includes the contribution of external tractions for cells located on the domain boundary.

Since Equation 17 must be satisfied for any virtual variation and , it is equivalent to the following system of algebraic equations (, = total number of Voronoi cells):


where = external force resultant, = mass, = first-order mass moments, = external moment resultant, = second-order mass moments, of cell . Equations 18 and 19 coincides with the force and moment equilibrium equations for each Voronoi cell.

Note that and (for uniform body force), if the vertex of the nodes of the Delaunay discretization coincide with the mass centroid of the Voronoi cells. This is the case for all the cells in the interior of the mesh if a centroidal Voronoi tessellation is adopted. Also, in general, for , and this leads to a non-diagonal mass matrix. A diagonalized mass matrix can be obtained simply by discarding the non-diagonal terms.

3.4 Discretized Constitutive Equations

In the DCM framework, the constitutive equations are imposed at the facet level where the facet tractions need to be expressed as function of the facet strains. For elasticity, by projecting the tensorial constitutive equations reported in Equation 3 in the local system of reference of each facet, one has


4 Two-Dimensional Implementation

4.1 Three-Node Triangular Element

In order to pursue a two-dimensional implementation of DCM, let us consider a 2D Delaunay-Voronoi discretization as shown in Figures 2a and b. A generic triangle can be considered as the triangular base of a prismatic volume as shown in Figure 2c and characterized by 6 vertexes, 9 edges, 2 triangular faces, and 3 rectangular faces. By considering the Voronoi vertexes on the two parallel triangular faces and face/edge points located at mid-thickness (see, e.g., points and ) a complete tessellation of the volume in six sub-volumes, one per vertex, can be obtained by triangular facets. Of these facets, are orthogonal to the triangular faces, set (see, e.g., the one connecting points , , in Figure 2c) and are parallel to the triangular faces, set , (see, e.g., the one connecting points , , in Figure 2c).

One can write


where is the out-of-plane thickness.

For plane strain conditions for the facet set and simply the second term in Equation 23 is zero. For plane stress, instead, for the facet set . Therefore, for the facet set , and , also, . Using these relations, one has . Substituting this relation in Equation 23 leads to . In addition the facet set is composed by 3 sets of 4 planar triangular facets. For each set, strains and tractions are the same on the 4 facets because the response is uniform through the thickness. Consequently the 4 facets can be grouped into one rectangular facet of area where is the in length of the facet (see Figure 2b)

By taking everything into account the volumetric strain for 2D problem can be written as


where is the area of the triangular element and for plane stress and for plane strain.

The triangular DCM element has 9 degrees of freedom, two displacements and one rotation for each node, which can be collected in one vector in which , , and are the element node indexes, and 1, 2, and 3 represent the three Cartesian coordinate axes. By using Equations 13 to 15, one can write , , , where


and ,

Figure 2: (a) Delaunay triangulation and Voronoi tesselation. (b) Facet geometry and associated edge length. (c) Triangular element geometry through the thickness.

The elastic stiffness matrix for a single DCM triangle can be computed similarly to classical FEM elements by writing the internal energy as function of the nodal degrees of freedom. By using Equations 20 and 21, one has