The Discrete Element Method has been used to simulate fracture dynamics beacuse its inherent capacity to reproduce multibody interaction, but in the case of elasticity mechanics the microparameters of the numerical model, required to replicate the properties of the material, are difficult to calibrate. On the other hand, damage models based on finite element strategies can easily reproduce the properties of the media but they can not simulate the dynamics of multiple fractures.
We propose a numerical approach, the Discrete Volume Method, to simulate fracture of brittle materials without the disadvantages mentioned, by combining the benefits of variational formulations and the numerical convenience of discrete element method to capture the dynamics of cracks. The Discrete Volume Method does not have microparameters, since the displacements are computed using the material properties and the fracture mechanism is controlled by an auxiliary damage field.
Within this work we discuss a numerical strategy to solve the elasticity problem upon unstructured and non conforming meshes, allowing all kinds of flatfaced elements (polygons in 2D and polyhedra in 3D). The core of the formulation relies on two numerical procedures the Control Volume Function Approximation (CVFA), and the polynomial interpolation in the neighborhood of the control volumes, which is used to solve the surface integrals resulting from applying the divergence theorem. By comparing the estimated stress against the analytical stress field of the well known test of an infinite plate with a hole, we show that this conservative approach is robust and accurate. A similar strategy is used to get the damage field solution.
In order to coupling both fields, displacement and damage, we use a finite increment arrangement for reducing the resdidual of elastic equation within each time step.
We develop a numerical formulation for time discretization based on the analytical solution of the differential equation resulting from assuming a continuous variation of internal forces of the system between time steps.
Finally, we show the effectiveness of the methodology by performing numerical experiments and comparing the solutions with published results.
The present investigation was sponsored by a CONACYT scolarship from the Mexican government and the TCAiNMaND project, an IRSES Marie Curie initiative under the European Union 7th Framework Programme.
In addition, the author want to express his gratitude to friends and mentors at CIMNE for all his support and shared wisdom, to friends at CIMAT for being always available for discussing the topics of this work and for his insightful commments and illuminating explanations about mathematical concepts, to Dr. Arturo Hernández for his support in promoting and divulging our discoveries in several conferences, and to Dr. Rafael Herrera for his priceless comments about numerical procedures and observations about the splines used here.
And last but not least, I want to thank to my beloved wife, Jimena, for all his support, tremendous patience and unconditional love since the beginning of this project, and to my Champion for teaching me every day how valuable is life. If angels exist, I already have a pair in my life.
One of the main aims in engineering is creating tools, structures and systems to enhance the quality of life in our society. In the course of the creation process, the design stage is critical for the final outcome. During this stage the engineer have to predict the prototype response when interacting with the physical world. Many of the observed phenomena in the physical world, such as solid mechanics, fluid dynamics, heat diffusion, and others, can be described with Partial Differential Equations (PDEs) by assuming time and space as a continuum.
Computational Continuum Mechanics (CCM) is the area dedicated to develop numerical methods and heuristics to solve these PDEs. Most of the methods can be classified into these two families: weighted residual and conservative methods. The Galerkin formulations are popular and widely used weighted residual methods, such as the Finite Element Method (FEM), which is a well established technique in Computational Solid Mechanics (CSM). Alternatively, the Finite Volume Method (FV) and the Control Volume Function Approximation (CVFA) are common approaches of conservative methods. The main difference between both families is that weighted residuals methods do not conserve quantities locally, but globally instead, resulting in linear systems with commendable numerical properties (symmetrical and wellconditioned matrices, for example). Nevertheless, due to its conservative nature, the second group is more attractive for fluid structure interaction ([1,2]) and multiphysics simulations ([3,4]), where several PDEsolvers must be coupled. For that reason, in recent years FV has been subject of interest for solving CSM problems.
Most of the CSM nonlinear strategies depend on the accuracy of the estimated stress field for the elasticity problem, such as those for plasticity and damage (see [5,6]). Hereafter we refer as elasticitysolver to the numerical computation that calculates the displacement and stress fields for a given domain and boundary conditions.
In industrial design it is critical to predict the cracks on materials in order to prevent a major failure on the whole system, especially in automotive, aeronautic and civil structures, where human lives can be lost. The three most important features that should be predicted with accuracy are the crack's morphology, tip's nucleation and evolution of the existing tips.
There are two main approaches to predict these cracks' features, the variational formulation which assumes a continuum where the crack is approximated by means of a function, and the multibody system where the cracks emerges naturally by the separation of the rigid bodies. The first approach estimates the internal mechanics of materials with high accuracy, and the second approach is more suitable to capture the dynamics of systems where the initial continuum is broken apart into several subdomains.
The main objective of this work is to describe a numerical method to predict cracks by combining the accuracy and efficency of variational formulations and the ability to capture the dynamics of multibody systems.
The prediction and analysis of brittle fracture is an intense research area with applicability to a wide range of industrial problems, such as the failure mechanism of structures, the fracking process, the detonations impact upon structures and the rock cutting. Moreover, the prevention of cracks is a main requirement in structural designs.
In his influential papers, Oñate et al [7,8], propose a FV format for structural mechanics based on triangular meshes, discussing the cell vertex scheme, the cell centred finite volume scheme and its corresponding mixed formulations, showing that the cell centred strategy produces the same symmetrical global stiffness matrix that FEM using linear triangular elements. Analogously, Bailey et al [9,10], develop a similar approach, but using quadrilateral elements to produce cellcentred volumes. Even though, the shapes of the volumes in both formulations are completely defined by the FEMlike mesh (triangular or quadrilateral) and it is not possible to handle arbitrary polygonal shapes, as we might expect when the mesh elements are produced by cracks.
Slone et al [11] extends the investigation of [7] by developing a dynamic solver based on an implicit Newmark scheme for the temporal discretization, with the motivation of coupling an elasticitysolver with his multiphysics modelling software framework, for later application to fluid structure interaction.
Another remarkable algorithm is the proposed by Demirdzic et al [12,13,14,15,16] The numerical procedure consists in decoupling the strain term into the displacement Jacobian and its transpose in a cellcentred scheme. The Jacobian is implicitly estimated by approximating the normal component of each face as the finite difference with respect to the adjacent nodes, while the Jacobian transpose is an explicit average of Taylor approximations around the same adjacent nodes. This decoupling produces a smaller memory footprint than FEM because the stiffness matrix is the same for all the components. The solution is found by solving one component each iteration in a coordinate descent minimization. This line of work has shaped most of the state of the art techniques in FV for coupling elasticitysolvers to Computational Fluid Dynamics (CFD) via finite volume practices (usually associated to CFD), such as the schemes proposed by [17,18,19]. However, this segregated algorithm may lead to slow convergence rates when processing nonlinear formulations, for example, when it is required to remove the positive principal components of the stress tensor in phasefield damage formulations [20]. In addition, if some nonlinear strategy requires multiple iterations of the linear elasticitysolver, such as finite increments in damage models, the nested iterations will increase the processing requirements for simple problems. To circumvent this drawback Cardiff et al [21] presents a fully block coupled direct solution procedure, which does not require multiple iterations, at expense of decomposing the displacement Jacobian of any arbitrary face into a) the Jacobian of the displacement normal component, b) the Jacobian of the displacement tangential component, c) the tangential derivative of the displacement normal component and d) the tangential derivative of the displacement tangential component. This decomposition complicates the treatment of the stress tensor in the iterative nonlinear solvers mentioned before for plasticity and damage.
A generalized finite volume framework for elasticity problems on rectangular domains is proposed by Cavalcante et al [22]. They use higher order displacement approximations at the expense of fixed axisaligned grids for discretization.
Nordbotten [23] proposes a generalization of the multipoint flux approximation (MPFA), which he names multipoint stress approximation (MPSA). The MPSA assembles unique linear expressions for the face average stress with more than two points in order to capture the tangential derivatives. The stress field calculated with this procedure is piecewise constant.
In this work, we propose an elasticitysolver based on CVFA techniques (see [24,25]), using piecewise polynomial interpolators for solving the surface integrals on the volumes boundaries. The polynomials degree can be increased without incrementing the system degrees of freedom, which make this method more suitable for nonlinear models and dynamic computations. Furthermore, this algorithm can handle polygonal/polyhedral, unstructured and non conforming meshes, and does not require the decomposition of the stress tensor.
There are remarkable methodologies to solve the nonlinear behaviour of brittle fracture using FEM, such as the damage models proposed in [26,27,28,29,30,31], the phasefield approaches to estimate the fracture surface described in [20,32,33] and the models of Extended FEM (XFEM) explained in [34,35,36]. However, these methods can not easily handle large displacements of the resulting subbodies after the fracture, such as the fragments blown up by a detonation. The Element Deletion Method could deal with these large displacements (see [37,38]), but none of these techniques can manage the collision between multiples bodies and the selfcollision of boundaries.
The Discrete Element Method (DEM) has been used to solve problems involving granular material with success (as presented in [39,40,41]), since it can handle discontinuities in the domain without special considerations. DEM defines the interaction mechanism of multiple rigidspheres (disks in 2D), such interaction is characterized by a set of microparameters which pretend to emulate the material properties. In order to approximate a continuum behaviour, the discrete elements are linked with cohesive bonds to its adjacent neighbours in the initial discretization. The fracture emerges when the cohesive bonds are broken systematically, this occurs if the force applied to them is superior to some threshold (which is a microparameter), a complete review of DEM is in [42,43,44,45,46].
There are two major challenges when we are working with the continuum using DEM. The first challenge is the approximation of the material properties with the microparameters, there are techniques to calculate these from a given material, as the proposed in [42], but none of these proposals proofs that the resulting behaviour of the body corresponds to the material properties. The second challenge is the computation of the system, due to the huge quantity of discrete elements (billions for some engineering problems) and the tiny time steps to maintain the numerical stability (a large time step could produce overlapping discrete elements and the wrong evolution of the displacements).
To handle these challenges, Oñate [45] proposes a DEM/FEM formulation with an underlying DEM discretization which is enabled when the finite elements are completely damaged, but this approach is expensive almost as much as the simple DEM. Zárate [47] proposes a FEM/DEM coupling scheme for fast computing simulations, but it requires the same microparameters than DEM. In the literature exists similar schemes to couple atomistic and continuum models [48,49,50,51,52], but all of them need microparameters to fix the interface between the discrete and the continuum model, and require small enough time steps to make the computation slow.
The Discrete Volume Method (DVM) aims to reduce the computational effort to perform a simulation of brittle fracture without the need of microparameters. The strategy is to solve the elasticity problem using the Control Volume Function Approximation method (CVFA), introduced in [53,24,25], on a coarse mesh and utilize an auxiliary damage field to refine the mesh in the damaged zones, separating the control volumes adjacent to completely damaged faces during the fracture process. The control volumes are named discrete volumes because them can be isolated from the domain.
DVM exploits the accuracy and robustness of CVFA and the ability to create cracks and handle multiple collisions of DEM.
We consider an arbitrary body, , with boundary . The displacement of a point at time is denoted by . We assume small deformations and deformation gradients, and the infinitesimal strain tensor, denoted , is given by

(2.1) 
Since we assume isotropic linear elasticity, the elastic energy density is defined

(2.2) 
where and are the Lamé parameters characterizing the material. These parameters are related with Young's modulus, , and Poisson's ratio, , by the following equivalences

(2.3) 
and

(2.4) 
where for plane stress analysis, and for plane strain and 3D cases.
The stress components are given by the partial derivative of the elastic energy density with respect to the corresponding strain component

(2.5) 
to simplify the notation, we use the fourth order tensor, , to map the strain field to the stress field

(2.6) 
where is the Kronecker delta. This tensor is symmetric, (major symmetry), (minor symmetry), and positive definite. The equation (2.5) is equivalent to

(2.7) 
where using the summation convention over repeated indices. Furthermore, since the strain tensor is symmetric, we can simplify the tensorial product to

(2.8) 
where is the identity matrix, defined in tensorial notation.
To model the loss of stiffness and the rupture of the material we use the damage field, denoted , which goes to one in the failure zones and it is equal to zero in the rest of the domain, as illustrated in Figure 1. We redefine the elastic energy density, , to consider the damage field effectsFigure 1: The left side shows the body with an internal fracture, denoted , under boundary conditions. The right side shows the damage field approximation of the fracture surface. 

(2.9) 
where is the energy contribution due to tension, calculated with the positive part of the principal strains, denoted , and is the energy contribution due to compression, calculated with the negative part of the principal strains, denoted . To simplify the notation we use and . The principal strains are calculated through a spectral decomposition of the strain tensor

(2.10) 
where is the diagonal matrix containing the principal strains, denoted , and is conformed by their orthonormal eigenvectors. The positive and negative contributions are defined by

where

with . The equation (2.14) implies

(2.15) 
Observe that if there is not damage, , the energy density of the equation (2.9) is equivalent to the elastic energy density of the equation (2.2). The energy contribution due to tension is obtained from

(2.16) 
using the equation (2.15), the contribution due to compression is given by

(2.17) 
The stress of equation (2.5) is now calculated as

(2.18) 
developing the derivatives, the stress is expressed as

(2.19) 
and rearranging the terms we obtain

(2.20) 
From here, we are going to use the symbol to refer the linear elastic stress.
Observe that for the equation (2.20) is equal to (2.8), however for we have only the compression contribution.
According to Griffith's theory of brittle fracture (see [20]), the energy required to create a unit area of fracture surface, , is equal to the critical fracture energy density, denoted , also known as critical energy release rate. The potential energy of the body, , is given by the sum of the elastic energy and the fracture energy

(2.21) 
Since we do not know the fracture surface, we use a crack surface density function, , to estimate the contribution of such surface in terms of the damage

(2.22) 
The damage field decays exponentially when goes away from the crack surface (see the work of Miehe [32,33]), this behaviour is given by the following differential equation

(2.23) 
where is a length scale parameter to control the smooth approximation of the crack. We take (2.23) as the Euler equation of the general form of the variational calculus problem

(2.24) 
to obtain

(2.25) 
By substituting (2.25) into (2.22) we approximate the fracture energy without a priori knowledge of the fracture surface, , with an integral over the entire domain, ,

(2.26) 
Replacing (2.26) into (2.21) we get the potential energy using only integrals over the domain ,

(2.27) 
The kinetic energy of the body is given by

(2.28) 
where is the density and is the velocity. Observe that the kinetic energy is unaffected by the damage field, resulting in a mass conservative scheme. The potential and kinetic energies defines the Lagrangian of the discrete fracture problem as

(2.29) 
Expanding the terms we have

(2.30) 
According to the principle of least action (see [33]), the displacement field is obtained from the following minimization

(2.31) 
and the damage field is given in a similar calculation

(2.32) 
Using the EulerLagrange equations to solve the minimization problems we get the strong form equations of motion

These equations of motion should be solved to found the displacement and damage fields.
The cracking process is irreversible, , this condition is enforced introducing a strain history field, , in the strong form equations of motion, which satisfies the KuhnTucker conditions for loading and unloading

In this work the strain history field is defined as the maximum elastic energy density due to tension from to current time

(2.35) 
where is the dummy time variable.
Replacing the elastic energy density due to tension, , by the strain history field, , in (2.33.b) we get the system to be solved

The displacement field satisfies the timedependent Neumann conditions given by upon the boundary and Dirichlet conditions given by upon the boundary , where . The damage gradient must be zero along the external boundary, . These conditions could be imposed by means of

The initial state of the system is characterized by

The strain history field, , could be used to model initial fracture surfaces (see appendix A of [20]).
For a given set of centroids, denoted , the discrete volumes are spheres (disks in 2D) with radii truncated by planes orthogonal to the line connecting the centroids and at the following point

(2.39) 
the point is in the middle of and if is equal to . Formally, the discrete volumes of the partition are defined by

(2.40) 
Figure 2 helps to visualize the discrete volume defined by the equation (2.40). The left side of Figure 3 illustrates the domain of the discrete volume with respect to the remaining volumes , and the right side shows the discrete volumes forming a continuum in the domain, .
Figure 2: The discrete volumes, , are defined by its radii and the planes orthogonal to the lines connecting the centroids and at the point , for all . 
The mass of the volumes is timeinvariant and its center of mass is assumed to be the centroid. To enforce these assumptions, we associate a mass, denoted , an initial density, , and an initial volume, , to the discrete volumes, such quantities are calculated as

(2.41) 
Then, the density associated to the discrete volumes at any time, denoted , is given by

(2.42) 
Figure 4 shows the density of calculated from (2.42) for three cases.
Figure 4: The density is updated depending on the current volume of the sphere (disk in 2D) in order to conserve the mass. 
The integrals over the faces of the discrete volumes requires the normal of their surface, , but only the shared faces have a constant normal, the integrals on the curved faces are considered with a Neumann condition equal to zero, since such faces are not interacting with other discrete volumes,

(2.43) 
We want to remark that the elastic energy is transferred from one volume to its neighbours through the shared faces and the size of such faces has a nonlinear behaviour with respect to the distance between its adjacent centroids. Most of the methodologies dealing with discrete bodies, such as the Discrete Element Method, assumes that this behaviour is linear. Figure 5 shows the surface area of the face shared by two discrete volumes with the same radius as a function of the distance between their centroids.
Figure 5: The curves shows the surface area of the face shared by two discrete volumes with the same radius as a function of their distance, also referred as penetration. 
On this chapter we go into the details of the numerical procedure by discussing the discretization with CVFA, the control volumes integration, the subfaces integrals, the simplexwise polynomial approximation, the pairwise polynomial approximation, the homeostatic splines used within the shape functions, the linear system assembling, how to impose boundary conditions, and two special cases of the formulation.
For the sake of legibility, in some parts of the text we unfold the matrices only for the bidimensional case, but the very same procedures must be followed for the 3D case.
The domain is discretized into control volumes, denoted , using the Control Volume Function Approximation (CVFA) proposed by Li et al [24,25]. The partition of is defined by

(3.1) 
where the boundary of each control volume, , is composed by flat faces, denoted ,

(3.2) 
Figure 7: The boundary of the three dimensional control volume is subdivided into flat faces, denoted . The unit vector is normal to the face . 
Every control volume must have a calculation point

(3.3) 
which is used to estimate the displacement field. Such a point is the base location to calculate the stiffness of the volume. In the volumes adjacent to the boundary , it is convenient to establish the calculation point over the corresponding boundary face,

(3.4) 
in order to set the Dirichlet condition directly on the point.
In this chapter we will focus our attention on the spatial discretization and numerical treatment of the stress term in first equation of motion (2.36.a), for simplicity assume , later we will remove this assumption.
We begin by integrating the stress divergence over the control volume

(3.5) 
using the divergence theorem we transform the volume integral into a surface integral over the volume boundary

(3.6) 
The evaluation of the surface integrals is based on the approximation of the displacement field inside the neighborhood of the volume, denoted ,

(3.7) 
making use of a group of piecewise polynomial interpolators, denoted . We are going to discuss these interpolators later in this section.
For that reason, the displacement field is decoupled from the stress tensor by using the strain (2.1) and stress (2.8) definitions. Taking advantage of the stress tensor symmetry , we rewrite the stress normal to the boundary as

(3.8) 
where is the face orientation matrix and is the engineering stress vector. Developing the stress definition (2.8) componentwise we can decompose it into the constitutive matrix, denoted , and the engineering strain vector, denoted , as follows

then the components of the strain vector are retrieved from the equation (2.1), and it is decomposed into the matrix differential operator and the displacement function .

(3.11) 
Summarizing the equations (3.8), (3.10) and (3.11) we have

(3.12) 
where is the stiffness of the volume boundary.
Once the displacement field is decoupled, we rewrite the equation (3.6) as

(3.13) 
Using the fact that the control volume boundary is divided into flat faces, as in equation (3.2), we split the integral (3.13) into the sum of the flat faces integrals

(3.14) 
Notice that the face orientation along the flat face, denoted , is constant. Furthermore, if the control volumes are considered to be made of a unique material and the flat faces to be formed by pairs of adjacent volumes, then the constitutive matrix along the flat face, denoted , is also considered constant. The matrix is estimated from the harmonic average of the Lamé parameters assigned to the adjacent volumes, where and are the properties of the volume ,

(3.15) 
With and we simplify the equation (3.14) as

(3.16) 
where is the strain integral along the flat face ,

(3.17) 
The accuracy of the method depends on the correct evaluation of this integral.
The surface integrals along the flat faces are calculated using an auxiliary piecewise polynomial approximation of the displacement field. This approximation is based on the simplices (triangles in 2D or tetrahedra in 3D) resulting from the Delaunay triangulation of the calculation points from the neighborhood of . The Delaunay triangulation is the best triangulation for numerical interpolation, since it maximizes the minimum angle of the simplices, which means that its quality is maximized as well. We define the neighborhood of volume as the minimum set of calculation points such that the simplices intersecting does not change if we add another calculation point to the set. Observe that the neighborhood does not always coincide with the set of calculation points in volumes adjacent to , as in most of the FV formulations. Once the neighborhood is triangulated, we ignore the simplices with angles smaller than degrees, and the simplices formed outside the domain, which commonly appear in concavities of the boundary . The local set of simplices resulting from the neighborhood of is denoted . Figure 8 illustrates the difference between (a) the simplices resulting from the triangulation of the calculation points in adjacent volumes and (b) those resulting from the triangulation of the proposed neighborhood .
The face is subdivided into subfaces, denoted ,

(3.18) 
these subfaces result from the intersection between and the control volume . Figure 8.b illustrates six key points of this approach, 1) the simplices are used to create a polynomial interpolation of over the boundary of the control volume, 2) most of the faces are intersected by several simplices, such faces must be divided into subfaces to be integrated, 3) some few faces are inside a single simplex, as illustrated in the face formed by and , 4) there are volumes that require information of nonadjacent volumes to calculate its face integrals, such as requires , 5) the dependency between volumes is not always symmetric, which means that if requires does not implies that requires , and 6) non conforming meshes are supported, as shown in the faces formed by , , , and .
The integral (3.17) is now rewritten in terms of the subfaces

(3.19) 
Each subface is bounded by a simplex, where the displacement , and it derivatives, , are a polynomial interpolation. Hence the integrals in equation (3.19) are solved exactly by using the GaussLegendre quadrature with the required number of integration points, denoted , depending on the polynomial degree,

(3.20) 
Most of the cases, the displacement is interpolated inside the simplices, but in some geometrical locations these can not be created, in consequence, the displacement is interpolated pairwise using the volumes adjacent to the subface . We discuss both strategies in the following subsections.
In the general case, the simplices are formed by points. The points forming the simplex that is bounding the subface are denoted , and its displacements .
The shape functions used for the polynomial interpolation are defined into the normalized space. A point in such space is denoted , its component is denoted , and the point forming the simplex is expressed . The nodes of the normalized simplex are given by the origin, , and the standard basis vectors,

(3.21) 
Figure 10: (a) The simplex formed by the points , and in the original space contains an interior point that is mapped to (b) into the normalized 2Dsimplex formed by the points , and . 
Figure 11: (a) The 3Dsimplex formed by the points , , and in the original space contains an interior point that is mapped to (b) into the normalized 3Dsimplex formed by the points , , and . 
The shape functions, denoted , are used to interpolate the displacement field inside the normalized simplex. Such functions are nonnegative and are given by

where is the homeostatic spline, which is the simplest polynomial defined in the interval that have derivatives equal to zero in the endpoints of the interval. We will discuss this spline later.
The set of shape functions is a partition of unity, which means that the sum of the functions in the set is equal to one into the interpolated domain

(3.23) 
furthermore, these functions are equal to one in its corresponding node, which implies that

The gradients of the shape functions with respect to the normalized space are denoted . The norm of the sum of such gradients is zero

(3.26) 
which means that there are not numerical artifacts into the strain field.
Any point inside the simplex can be formulated as a function of a point in the normalized space, , by using the shape functions and the points forming the simplex

(3.27) 
In order to calculate the normalized point, denoted , associated to the integration point , we use the shape functions definitions to rewrite the equation (3.27) in matrix form

where is the vector resulting from evaluating the spline for componentwise, and is the distortion matrix given by the concatenation of the following column vector differences

(3.30) 
Now, from equation (3.29) we retrieve the point as

(3.31) 
and solving for we obtain

(3.32) 
where is the inverse function of applied componentwise to the product of the matrixvector operation.
Similar to the approximation in equation (3.27), within the simplex enclosing the subface , the displacement field evaluated at is defined as,

(3.33) 
Hence, when calculating the quadrature of equation (3.20), the strain evaluated at the integration point is given by

where captures the deformation at , and is the vector with the concatenated displacement components of the points forming the simplex.
In order to calculate the deformation matrix , we require the derivatives of the shape functions with respect to , denoted . These derivatives are calculated by solving the linear systems resulting from the chain rule

(3.36) 
where is the geometric jacobian evaluated at . This jacobian relates both spaces, captures the distortion of the simplex, and is derivated from equation (3.27),

(3.37) 
The gradients of the shape functions with respect to inside the sum are obtained straightforward once we have the spline first derivative . Notice that the geometric jacobian is equivalent to the distortion matrix if and only if the homeostatic spline is .
Since we are not making any assumption about the volumes distribution through the mesh, neither about the internal location of its calculation points, then we have to deal with portions of the mesh that are no covered by any simplex. Figure 12 illustrates the two most common cases. The first case takes place in meshes where the calculation points of volumes contiguous to the boundary are in the interior of such volumes, producing subfaces not intersected by any simplex, and the second case occurs when elongated sections of the domain are discretized with a queue of aligned volumes, where each volume has only two neighbors on opposite faces and no simplex can be formed.
In such cases, the displacement field within the subface is a pairwise polynomial approximation between the adjacent volumes, and , regardless the dimension

(3.38) 
where and are the shape functions, and is the normalized projection of the integration point into the vector which goes from to , denoted

(3.39) 
When calculating the quadrature of equation (3.20), the pairwise strain is given by

This pairwise approximation must be used only when necessary because it can not capture the deformation orthogonal to .
The homeostatic spline is a function of a single variable defined from to , denoted , and curved by the parameter , which indicates the level of smoothness. This spline is the simplest polynomial with derivatives equal to zero at the endpoints and . The polynomial degree is given by , and such a polynomial requires integration points to calculate the exact integral in equation (3.20) using the GaussLegendre quadrature.
When designing this spline, we wanted to gain accuracy by building a piecewise bellshaped interpolation function around the calculation points, inspired on the infinitely smooth kernels used in other numerical techniques. Therefore, we force the derivatives of the polynomial to be zero over such points in order to homogenize the function. For that reason, we use the term homeostatic spline when referring to this spline.
To fulfill the smoothness requisites commented before, we solved a linear system for calculating the coefficients of the polynomial. The equations of this system were obtained by forcing the derivatives to be zero at the endpoints. Once we solved the coefficients for the first twenty polynomials, from to , we found out that the first half of such coefficients are null, and the entire polynomial can be calculated directly as

(3.43) 
where is the not null coefficient

(3.44) 
is the number of factors needed to calculate

(3.45) 
and is accumulation of the coefficients for normalizing the spline

(3.46) 
The first derivative is simply calculated as

(3.47) 
Table 1 shows the polynomials resulting from low values of and Figure 14 depicts (a) the evolution of the spline as we increase the smoothness parameter from to , and (b) the evolution of it first derivative.
Smoothness  Homeostatic spline 
Figure 14: (a) The evolution of the homeostatic spline from to illustrates the smoothness requirements at the endpoints of each spline and its (b) first derivatives. 
Since the derivatives of the homeostatic spline (3.43) are zero at the endpoints of the interval , the inverse function is not defined in that points. However, we estimate a pseudoinverse within this interval, , by finding the coefficients of a polynomial of the same degree, , such that the endpoints coincide with the spline and the first derivative at the midpoint is equivalent to the inverse of the spline first derivative, that is

(3.48) 
The higher derivatives in the midpoint are forced to be zero. Once we calculated the coefficients for the first twenty polynomials, from to , we found out that the pseudoinverse can be approximated directly from the following formulae

(3.49) 
where is the coefficient for , and is the factor that distinguish higher order coefficients. Such terms are calculated as

(3.50) 
respectively. Figure 15 exhibits the curves for the first seven levels of smoothness. The null higher derivatives requirement is noticeable at the midpoint.
Figure 15: Curves of the pseudoinverse for the first seven levels of smoothness. The slope at the midpoint exposes the null higher derivatives requirement when increasing the polynomial order. 
By using the simplexwise (3.35) or the pairwise (3.42) approximation, the strain face integral (3.19) is reformulated as

(3.51) 
then, the volume equilibrium equation (3.16) is

(3.52) 
reordering the terms we obtain

(3.53) 
where the matrix

(3.54) 
is the stiffness contribution at the integration point , within the subface when integrating the volume. Observe that the stiffness matrix is rectangular and the resulting global stiffness matrix is asymmetric.
The Neumann boundary conditions are imposed over the volume faces intersecting the boundary, by replacing the corresponding term in the sum of equation (3.14) with the integral of the function provided in (2.37.a),

(3.55) 
The Dirichlet conditions are imposed over the volumes calculation points by fixing the displacement as it is evaluated in the function given in (2.37.b),

(3.56) 
Since the Dirichlet conditions are imposed directly on the calculation points, these points must be located along the face which intersects the boundary with the condition .
By making some considerations, we identify two special cases where the calculations can be simplified, in order to increase the performance of the total computation, at the expense of losing control over the volumes shape. These cases are 1) the Voronoi mesh assumption and 2) the FVFEM correlation.
In the first case, we assume that the initial mesh is equivalent to the Voronoi diagram and that the Voronoi centres correspond to the calculation points . Hence, the subdivision of the neighborhood is already given by the Delaunay triangulation which is dual to the Voronoi mesh, as illustrated in Figure 18.a. Moreover, the integrals of subfaces using pairwise approximations can be exactly integrated with the midpoint rule, since the faces are orthogonal to the vector joining the calculation points , and the derivatives along the subface are constants.
In the second case, the initial mesh is generated from a FEMlike triangular mesh and the approximations are assumed to be linear. In such a case, the calculation points are defined to be the nodes of the triangular mesh, and the volume faces are created by joining the centroids of the triangles with the midpoint of the segments, as presented in Figure 18.b. This particular version is equivalent to the cellcentred finite volume scheme introduced by Oñate et al [7], who proved that the global linear system produced by this FV scheme is identical to that produced by FEM if the same mesh is used.
In this chapter we will focus on the numerical treatment of the second equation of motion (2.36.b), this equation describes the damage mechanics within the physical system by considering the potential energy produced by tensile stress.
As discussed in the mathematical formulation, the damage field is a smooth approximation of the fracture surface, a benefit of this approach is that fracture morphology is completely defined by the solution of this equation and we do not have to track the crack propagation with auxiliary checking procedures neither to check for new crack nucleations. However, it is important to be aware about the effects over the stress field produced by the scale length parameter which controls the smoothness of damage field solution. We observe that a length parameter proportional to the average size of control volumes, denoted , produces accurate results, these mesh size is taken as

(4.1) 
Figure 19: a) Damage field above control volumes shows how the crack arises along volumes boundaries. b) Scale length parameter controls the smoothness of the damage field. 
For assembling the system of equations We will follow a similar path to that used in the first equation of motion by using the same partition and interpolators, simplewise and pairwise approximations also apply for the damage field.
We start by integrating the strong form equation of motion (2.36.b) over the control volumes of the partition ,

(4.2) 
using the divergence theorem on the second integral we get

(4.3) 
Since is a material property, we assume that it is constant along the control volume, and dividing the first integral in two terms we obtain

(4.4) 
In order to solve volume integrals involving the strain history field, we use the following partition of the control volume

(4.5) 
The surface integral is solved along subfaces defined in (3.18), and the remaining volume integrals are solved using partition (4.5),

(4.6) 
The damage field is estimated using the same shape functions, (3.22.a) and (3.22.b), that we use for the displacement field,

where is the point corresponding to in the normalized space, is the vector containing the shape functions evaluated at , and is the vector containing the estimation of the damage field at the nodes forming the simplex. The gradient of the damage field is given by

(4.9) 
where is calculated from the chain rule in (3.36). Now the equation is fully discretized, the next step is to solve the integrals.
The first integral in equation (4.6) is approximated using the midpoint rule,

(4.10) 
where is the damage estimated at calculation point . Due to the simple nature of polygonal subvolumes , we can always reduce them to simplices in order to use the GaussLegendre quadrature to solve the volume integrals, in 2D is straightforward,

where is the number of points in the quadrature, and is the strain history field evaluated with the strain information of the simplex corresponding to subface . Last but not least, we solve the surface integral that unfolds the damage gradient defined in (4.9), using again GaussLegendre quadrature

where is the matrix containing the derivatives of the shape functions evaluated at . For simplicity, the matrix notation in previous equation shows only values for 2D case.
Substituting, equations (4.10), (4.12), (4.13) and (4.17) into (4.6) we get

(4.19) 
The damage of the face produced by excesive loads is captured by vector

(4.20) 
on the other hand, the potential energy to create new crack surfaces is captured by

(4.21) 
Now we can rewrite the damage equation (4.19) for the control volume as follows

(4.22) 
Since damage is not a physical quantity, there is no damage flow between the system and the exterior, for that reason all the Neumann conditions are null, and Dirichlet conditions can be numerically set, but in our formulation these are defined in the initial strain history field .
In this chapter we will remove the assumption done in equation (3.5) about null acceleration, , and we will discuss in detail the discretization of time derivatives.
A common approach to approximate these derivatives in dynamic stress analysis is a staggered scheme by means of Finite Differences (FD), such as in [11], [34] and [52]. The simplicity of FD makes easy the incorporation of spatial nonlinear phenomena, for instance fracture and damage, nevertheless FD does not consider the stress state within its approximation and we are forced to use tiny time steps to diminish spurious stress waves that produce undesired artifical internal forces.
In this work we build a customized numerical scheme considering the time variation of internal forces in order to get an approximation capable of performing bigger and more accurate time steps.
In order to analyze the dynamic component of elasticity equation (2.36.a) we define the stress state of control volume as a function of time,

(5.1) 
with the intention of considering internal forces in the approximation,

(5.2) 
Equation (5.2) is an ordinary differential equation that can be solved analytically for a time step with the following Dirichlet conditions

We assume that temporal variation of the internal forces is given by

(5.4) 
where is the stress state calculated at , is the stress state which will be estimated at time , and is the shape function modelling time variation between concecutive stress states. This shape function has only two constraints and , for that reason we use as a normalizer in equation (5.4). In the discussion of this chapter we utilize “stress state” and “internal forces” as synonyms to refer the same term in equation (5.1).
Figure 21 illustrates the variation of the stress state in terms of the shape function that is used as interpolator between the value at two contiguous time steps.
Using the asumption in (5.4), we get the analytical solution of the equation (5.2) for the interval by means of the Laplace transform (see appendix A for details),

(5.5) 
where is the velocity at time , and is the convolution between the spline and the function , as defined in appendix A.
By setting the second boundary condition, , into the analytical solution (5.5), we can find the velocity required to fulfill both Dirichlet conditions

(5.6) 
where is the convolution evaluated at . Thus, we replace equation (5.6) into (5.5) to get the analytical solution in terms of the known Dirichlet conditions

(5.7) 
now we can obtain the analytical time derivative (velocity),

(5.8) 
where is the time derivative of . Since the analytical solution (5.5) requires the initial conditions (displacement and velocity), we calculate the initial velocity by using equation (5.8) for a previous time interval ,

(5.9) 
where and . Finally, we replace equation (5.9) into (5.5) in order to get an analytical solution for as a function of two history system states,

(5.10) 
evaluating such an equation at , denoted , and rearranging the terms we get a numerical approximation dependent of the convolution of choosen spline,

(5.11) 
observe that even in the simplest case this approximation is more accurate than simple central finite differences applied directly on equation (5.2), because it takes into account variable time steps and the time variation of the system internal forces.
The analytical solution (5.11) of the ordinary differential equation (5.2) can be used to generate a family of numerical approximations, these approximations has a similar structure but distinct coefficients that depend on the shape function used for time variation of stress state. In this work we explore distinct families of functions in order to get a continuous stress state in contiguous time steps.
In order to select a good shape function for stress time variation we used the harmonic oscillator to measure the sensibility of the numerical scheme to distinct shape functions. The harmonic oscilator differential equation is

(5.12) 
where is the stiffness of the system, is the mass of the body and is the onedimensional displacement. The analytical solution of equation (5.12) is

(5.13) 
where is the oscillation amplitude, the oscillation frequency and the phase, such constants are calculated in terms of material properties

with and as initial displacement and initial velocity respectively. In our numerical tests, the one dimensional stress state, denoted , is assumed to be

(5.17) 
For simplicity, in this sensibility analysis we used a constant time step .
By using a central finite differences scheme, equation (5.12) can be rewritten as

(5.18) 
and the solution for next time step is calculated from

To measure the relative error with respect to analytical solution, we used (5.20) to compute the solution in the interval . To make evident the numerical drawbacks of FD we utilized a big enough . In Figure 22 we show the experiment results in four plots, the first one shows the displacement against time with a solid line for analytical solution and a dashed line for the numerical one, in this plot is clear that the system is loosing energy through time, no matter how small is the system will always loose energy proportionally to the time step. The second plot shows the phase space (solid line is analytical solution), which is velocity against displacement, in this plot we see the closing spiral when displacement and velocity decreases. The third plot shows the total energy in the system to emphasize that it is loosing energy, while total energy of analytical solution (solid line) remains constant. The fourth plot shows the cumulative relative error for distinct , such an error remains almost consant for since the numerical system looses all its energy in the first few time steps. In this plot we compute the comulative error as

(5.21) 
If we choose a linear shape function, , in order to set a constant variation of the internal forces in the interval , the convolution and its time derivative are given by
