Discrete Volume Method. A variational approach for brittle fracture

Abstract

The Discrete Element Method has been used to simulate fracture dynamics beacuse its inherent capacity to reproduce multi-body 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 flat-faced 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.

1 Introduction

1.1 Problem definition

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 well-conditioned 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 PDE-solvers must be coupled. For that reason, in recent years FV has been subject of interest for solving CSM problems.

Most of the CSM non-linear 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 elasticity-solver 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 multi-body 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.

1.2 State of the art

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 cell-centred volumes. Even though, the shapes of the volumes in both formulations are completely defined by the FEM-like 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 elasticity-solver with his multi-physics 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 cell-centred 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 elasticity-solvers 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 non-linear formulations, for example, when it is required to remove the positive principal components of the stress tensor in phase-field damage formulations [20]. In addition, if some non-linear strategy requires multiple iterations of the linear elasticity-solver, 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 non-linear 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 axis-aligned grids for discretization.

Nordbotten [23] proposes a generalization of the multi-point flux approximation (MPFA), which he names multi-point 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 piece-wise constant.

In this work, we propose an elasticity-solver based on CVFA techniques (see [24,25]), using piece-wise 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 non-linear 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 non-linear behaviour of brittle fracture using FEM, such as the damage models proposed in [26,27,28,29,30,31], the phase-field 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 sub-bodies 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 self-collision 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 rigid-spheres (disks in 2D), such interaction is characterized by a set of micro-parameters 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 micro-parameter), 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.

2 Mathematical formulation

2.1 Continuum mechanics

We consider an arbitrary body, ${\textstyle \Omega \in \mathbb {R} ^{\hbox{dim}}}$, with boundary ${\textstyle \partial \Omega }$. The displacement of a point ${\textstyle \mathbf {x} \in \Omega }$ at time ${\textstyle t\in [0,T]}$ is denoted by ${\textstyle {\boldsymbol {u}}(\mathbf {x} ,t)\in \mathbb {R} ^{\hbox{dim}}}$. We assume small deformations and deformation gradients, and the infinitesimal strain tensor, denoted ${\textstyle {\boldsymbol {\varepsilon }}(\mathbf {x} ,t)\in \mathbb {R} ^{{\hbox{dim}}\times {\hbox{dim}}}}$, is given by

${\boldsymbol {\varepsilon }}(\mathbf {x} ,t)=\nabla _{s}{\boldsymbol {u}}={\frac {1}{2}}\left(\nabla {\boldsymbol {u}}+\left[\nabla {\boldsymbol {u}}\right]^{T}\right).$

(2.1)

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

$\psi _{e}({\boldsymbol {\varepsilon }})={\frac {1}{2}}\lambda ~tr({\boldsymbol {\varepsilon }})^{2}+\mu ~tr({\boldsymbol {\varepsilon }}^{T}{\boldsymbol {\varepsilon }}),$

(2.2)

where ${\textstyle \lambda }$ and ${\textstyle \mu }$ are the Lamé parameters characterizing the material. These parameters are related with Young's modulus, ${\textstyle E}$, and Poisson's ratio, ${\textstyle \nu }$, by the following equivalences

$\mu ={\frac {E}{2(1+\nu )}},$

(2.3)

and

$\lambda ={\frac {\nu E}{(1+\nu )(1-\lambda _{\nu })}}$

(2.4)

where ${\textstyle \lambda _{\nu }=\nu }$ for plane stress analysis, and ${\textstyle \lambda _{\nu }=2\nu }$ 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

${\boldsymbol {\sigma }}_{ij}={\frac {\partial \psi _{e}}{\partial {\boldsymbol {\varepsilon }}_{ij}}},$

(2.5)

to simplify the notation, we use the fourth order tensor, ${\textstyle C}$, to map the strain field to the stress field

$C_{ijkl}=\lambda ~\delta _{ij}\delta _{kl}+\mu ~(\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}),$

(2.6)

where ${\textstyle \delta _{ij}}$ is the Kronecker delta. This tensor is symmetric, ${\textstyle C_{ijkl}=C_{klij}}$ (major symmetry), ${\textstyle C_{ijkl}=C_{ijlk}}$ (minor symmetry), and positive definite. The equation (2.5) is equivalent to

${\boldsymbol {\sigma }}=C:{\boldsymbol {\varepsilon }},$

(2.7)

where ${\textstyle C:{\boldsymbol {\varepsilon }}=C_{ijkl}{\boldsymbol {\varepsilon }}_{kl}}$ using the summation convention over repeated indices. Furthermore, since the strain tensor is symmetric, we can simplify the tensorial product to

${\boldsymbol {\sigma }}=2\mu ~{\boldsymbol {\varepsilon }}+\lambda ~tr({\boldsymbol {\varepsilon }})\mathbf {I} ,$

(2.8)

where ${\textstyle \mathbf {I} }$ is the identity matrix, defined ${\textstyle \mathbf {I} _{ij}=\delta _{ij}}$ in tensorial notation.

To model the loss of stiffness and the rupture of the material we use the damage field, denoted ${\textstyle {\boldsymbol {d}}(\mathbf {x} ,t)\in [0,1]}$, 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, ${\textstyle \psi _{e}}$, to consider the damage field effects

Figure 1: The left side shows the body with an internal fracture, denoted $\Gamma$, under boundary conditions. The right side shows the damage field approximation of the fracture surface.

$\psi _{d}({\boldsymbol {\varepsilon }},{\boldsymbol {d}})=(1-{\boldsymbol {d}})^{2}\psi _{e}({\boldsymbol {\varepsilon }}^{+})+\psi _{e}({\boldsymbol {\varepsilon }}^{-}),$

(2.9)

where ${\textstyle \psi _{e}({\boldsymbol {\varepsilon }}^{+})}$ is the energy contribution due to tension, calculated with the positive part of the principal strains, denoted ${\textstyle {\boldsymbol {\varepsilon }}^{+}}$, and ${\textstyle \psi _{e}({\boldsymbol {\varepsilon }}^{-})}$ is the energy contribution due to compression, calculated with the negative part of the principal strains, denoted ${\textstyle {\boldsymbol {\varepsilon }}^{-}}$. To simplify the notation we use ${\textstyle \psi _{e}^{+}=\psi _{e}({\boldsymbol {\varepsilon }}^{+})}$ and ${\textstyle \psi _{e}^{-}=\psi _{e}({\boldsymbol {\varepsilon }}^{-})}$. The principal strains are calculated through a spectral decomposition of the strain tensor

$\mathbf {P} {\boldsymbol {\Lambda }}\mathbf {P} ^{T}\leftarrow {\boldsymbol {\varepsilon }},$

(2.10)

where ${\textstyle {\boldsymbol {\Lambda }}}$ is the diagonal matrix containing the principal strains, denoted ${\textstyle \lambda _{i}}$, and ${\textstyle \mathbf {P} }$ is conformed by their orthonormal eigenvectors. The positive and negative contributions are defined by

${\boldsymbol {\varepsilon }}^{+}=\mathbf {P} {\boldsymbol {\Lambda }}^{+}\mathbf {P} ^{T},$ (2.11) ${\boldsymbol {\varepsilon }}^{-}=\mathbf {P} {\boldsymbol {\Lambda }}^{-}\mathbf {P} ^{T},$ (2.12)

where

${\boldsymbol {\Lambda }}^{+}=diag(\langle \lambda _{1}\rangle ,\langle \lambda _{2}\rangle ,\langle \lambda _{3}\rangle ),$ (2.13) ${\boldsymbol {\Lambda }}^{-}={\boldsymbol {\Lambda }}-{\boldsymbol {\Lambda }}^{+},$ (2.14)

with ${\textstyle \langle \lambda \rangle =\max(\lambda _{i},0)}$. The equation (2.14) implies

${\boldsymbol {\varepsilon }}^{-}={\boldsymbol {\varepsilon }}-{\boldsymbol {\varepsilon }}^{+}.$

(2.15)

Observe that if there is not damage, ${\textstyle {\boldsymbol {d}}=0}$, 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

$\psi _{e}^{+}={\frac {1}{2}}\lambda ~\langle tr({\boldsymbol {\varepsilon }})\rangle ^{2}+\mu ~tr\left(\left[{\boldsymbol {\varepsilon }}^{+}\right]^{2}\right),$

(2.16)

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

$\psi _{e}^{-}={\frac {1}{2}}\lambda ~\left(tr({\boldsymbol {\varepsilon }})-\langle tr({\boldsymbol {\varepsilon }})\rangle \right)^{2}+\mu ~tr\left(\left[{\boldsymbol {\varepsilon }}-{\boldsymbol {\varepsilon }}^{+}\right]^{2}\right)$

(2.17)

The stress of equation (2.5) is now calculated as

${\boldsymbol {\sigma }}_{ij}=(1-{\boldsymbol {d}})^{2}{\frac {\partial \psi _{e}^{+}}{\partial {\boldsymbol {\varepsilon }}_{ij}}}+{\frac {\partial \psi _{e}^{-}}{\partial {\boldsymbol {\varepsilon }}_{ij}}},$

(2.18)

developing the derivatives, the stress is expressed as

${\boldsymbol {\sigma }}=(1-{\boldsymbol {d}})^{2}\underbrace {\left(\lambda ~\langle tr({\boldsymbol {\varepsilon }})\rangle \mathbf {I} +2\mu ~{\boldsymbol {\varepsilon }}^{+}\right)} _{{\boldsymbol {\sigma }}^{+}{\hbox{(Stress due to tension)}}}+\underbrace {\left(\lambda \left(tr({\boldsymbol {\varepsilon }})-\langle tr({\boldsymbol {\varepsilon }})\rangle \right)\mathbf {I} +2\mu \left({\boldsymbol {\varepsilon }}-{\boldsymbol {\varepsilon }}^{+}\right)\right)} _{{\boldsymbol {\sigma }}^{-}{\hbox{(Stress due to compression)}}},$

(2.19)

and rearranging the terms we obtain

${\boldsymbol {\sigma }}=\underbrace {2\mu ~{\boldsymbol {\varepsilon }}+\lambda ~tr({\boldsymbol {\varepsilon }})\mathbf {I} } _{{\boldsymbol {\sigma }}^{e}{\hbox{(Linear elastic stress)}}}+\left({\boldsymbol {d}}^{2}-2{\boldsymbol {d}}\right){\boldsymbol {\sigma }}^{+}.$

(2.20)

From here, we are going to use the symbol ${\textstyle {\boldsymbol {\sigma }}^{e}}$ to refer the linear elastic stress.

Observe that for ${\textstyle {\boldsymbol {d}}=0}$ the equation (2.20) is equal to (2.8), however for ${\textstyle {\boldsymbol {d}}=1}$ we have only the compression contribution.

2.2 Fracture mechanics

According to Griffith's theory of brittle fracture (see [20]), the energy required to create a unit area of fracture surface, ${\textstyle \Gamma }$, is equal to the critical fracture energy density, denoted ${\textstyle {G}}$, also known as critical energy release rate. The potential energy of the body, ${\textstyle \Psi _{P}}$, is given by the sum of the elastic energy and the fracture energy

$\Psi _{P}({\boldsymbol {u}},\Gamma )=\int _{\Omega }\psi _{d}(\nabla {\boldsymbol {u}})dV+\int _{\Gamma }{G}dS.$

(2.21)

Since we do not know the fracture surface, we use a crack surface density function, ${\textstyle \gamma ({\boldsymbol {d}})}$, to estimate the contribution of such surface in terms of the damage

$\Gamma ({\boldsymbol {d}})=\int _{\Omega }\gamma ({\boldsymbol {d}})dV.$

(2.22)

The damage field decays exponentially when ${\textstyle \mathbf {x} }$ goes away from the crack surface (see the work of Miehe [32,33]), this behaviour is given by the following differential equation

${\boldsymbol {d}}-h^{2}\nabla ^{2}{\boldsymbol {d}}=0,$

(2.23)

where ${\textstyle h}$ 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

${\boldsymbol {d}}(\mathbf {x} )=arg\,min_{\boldsymbol {d}}\left\{\Gamma ({\boldsymbol {d}})\right\},$

(2.24)

to obtain

$\gamma ({\boldsymbol {d}})={\frac {{\boldsymbol {d}}^{2}}{2h}}+{\frac {h}{2}}(\nabla {\boldsymbol {d}}\cdot \nabla {\boldsymbol {d}}).$

(2.25)

By substituting (2.25) into (2.22) we approximate the fracture energy without a priori knowledge of the fracture surface, ${\textstyle \Gamma }$, with an integral over the entire domain, ${\textstyle \Omega }$,

$\int _{\Gamma }{G}dS\approx \int _{\Omega }{G}\left({\frac {{\boldsymbol {d}}^{2}}{2h}}+{\frac {h}{2}}\nabla {\boldsymbol {d}}\cdot \nabla {\boldsymbol {d}}\right)dV.$

(2.26)

2.3 Formulation of the equations of motion

Replacing (2.26) into (2.21) we get the potential energy using only integrals over the domain ${\textstyle \Omega }$,

$\Psi _{P}({\boldsymbol {u}},{\boldsymbol {d}},\nabla {\boldsymbol {d}})=\int _{\Omega }\psi _{d}(\nabla {\boldsymbol {u}})dV+\int _{\Omega }{G}\left({\frac {{\boldsymbol {d}}^{2}}{2h}}+{\frac {h}{2}}\nabla {\boldsymbol {d}}\cdot \nabla {\boldsymbol {d}}\right)dV.$

(2.27)

The kinetic energy of the body is given by

$\Psi _{K}({\dot {\boldsymbol {u}}})={\frac {1}{2}}\int _{\Omega }\rho \left({\dot {\boldsymbol {u}}}\cdot {\dot {\boldsymbol {u}}}\right)dV,$

(2.28)

where ${\textstyle \rho (\mathbf {x} ,t)\in \mathbb {R} }$ is the density and ${\textstyle {\dot {\boldsymbol {u}}}(\mathbf {x} ,t)\in \mathbb {R} ^{\hbox{dim}}}$ 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

$L({\boldsymbol {u}},{\dot {\boldsymbol {u}}},{\boldsymbol {d}},\nabla {\boldsymbol {d}})=\Psi _{K}({\dot {\boldsymbol {u}}})-\Psi _{P}({\boldsymbol {u}},{\boldsymbol {d}},\nabla {\boldsymbol {d}}).$

(2.29)

Expanding the terms we have

$L=\displaystyle \int _{\Omega }\left[{\frac {\rho }{2}}\left({\dot {\boldsymbol {u}}}\cdot {\dot {\boldsymbol {u}}}\right)-\psi _{d}(\nabla {\boldsymbol {u}})-{G}\left({\frac {{\boldsymbol {d}}^{2}}{2h}}+{\frac {h}{2}}\nabla {\boldsymbol {d}}\cdot \nabla {\boldsymbol {d}}\right)\right]dV,$ $=\displaystyle \int _{\Omega }\left[{\frac {\rho }{2}}\left({\dot {\boldsymbol {u}}}\cdot {\dot {\boldsymbol {u}}}\right)-(1-{\boldsymbol {d}})^{2}\psi _{e}^{+}-\psi _{e}^{-}-{G}\left({\frac {{\boldsymbol {d}}^{2}}{2h}}+{\frac {h}{2}}\nabla {\boldsymbol {d}}\cdot \nabla {\boldsymbol {d}}\right)\right]dV.$

(2.30)

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

${\boldsymbol {u}}=arg\,min_{\boldsymbol {u}}\left\{\int _{0}^{T}L({\boldsymbol {u}},{\dot {\boldsymbol {u}}},{\boldsymbol {d}},\nabla {\boldsymbol {d}})dt\right\},$

(2.31)

and the damage field is given in a similar calculation

${\boldsymbol {d}}=arg\,min_{\boldsymbol {d}}\left\{\int _{\Omega }L({\boldsymbol {u}},{\dot {\boldsymbol {u}}},{\boldsymbol {d}},\nabla {\boldsymbol {d}})dV\right\}.$

(2.32)

Using the Euler-Lagrange equations to solve the minimization problems we get the strong form equations of motion

$\rho {\ddot {\boldsymbol {u}}}-\nabla \cdot {\boldsymbol {\sigma }}=0,$ (2.33.a) $2(1-{\boldsymbol {d}})\psi _{e}^{+}-{\frac {G}{h}}{\boldsymbol {d}}+{G}h\nabla ^{2}{\boldsymbol {d}}=0.$ (2.33.b)

These equations of motion should be solved to found the displacement and damage fields.

The cracking process is irreversible, ${\textstyle \Gamma (\mathbf {x} ,t)\subseteq \Gamma (\mathbf {x} ,t+\Delta t)}$, this condition is enforced introducing a strain history field, ${\textstyle {H}}$, in the strong form equations of motion, which satisfies the Kuhn-Tucker conditions for loading and unloading

${H}-\psi _{e}^{+}\geq 0,$ (2.34.a) ${\dot {H}}\geq 0,$ (2.34.b) $\left({H}-\psi _{e}^{+}\right){\dot {H}}=0.$ (2.34.c)

In this work the strain history field is defined as the maximum elastic energy density due to tension from ${\textstyle t=0}$ to current time

${H}(\mathbf {x} ,t)=\max _{\tau }\left\{\psi _{e}^{+}(\mathbf {x} ,\tau )\right\},\tau \in [0,t],$

(2.35)

where ${\textstyle \tau }$ is the dummy time variable.

Replacing the elastic energy density due to tension, ${\textstyle \psi _{e}^{+}}$, by the strain history field, ${\textstyle {H}}$, in (2.33.b) we get the system to be solved

$\nabla \cdot {\boldsymbol {\sigma }}=\rho {\ddot {\boldsymbol {u}}},$ (2.36.a) $\left(1+{\frac {2h{H}}{G}}\right){\boldsymbol {d}}-h^{2}\nabla ^{2}{\boldsymbol {d}}={\frac {2h{H}}{G}}.$ (2.36.b)

The displacement field satisfies the time-dependent Neumann conditions given by ${\textstyle {\boldsymbol {b}}_{N}}$ upon the boundary ${\textstyle \Gamma _{N}}$ and Dirichlet conditions given by ${\textstyle {\boldsymbol {u}}_{D}}$ upon the boundary ${\textstyle \Gamma _{D}}$, where ${\textstyle \partial \Omega =\Gamma _{N}\cup \Gamma _{D}}$. The damage gradient must be zero along the external boundary, ${\textstyle \partial \Omega }$. These conditions could be imposed by means of

${\boldsymbol {\sigma }}\mathbf {n} ={\boldsymbol {b}}_{N}(\mathbf {x} ,t),\mathbf {x} \in \Gamma _{N},t\in [0,T],$ (2.37.a) ${\boldsymbol {u}}={\boldsymbol {u}}_{D}(\mathbf {x} ,t),\mathbf {x} \in \Gamma _{D},t\in [0,T],$ (2.37.b) $\nabla {\boldsymbol {d}}\cdot \mathbf {n} =0,\mathbf {x} \in \partial \Omega ,t\in [0,T].$ (2.37.c)

The initial state of the system is characterized by

${\boldsymbol {u}}(\mathbf {x} ,0)={\boldsymbol {u}}_{0}(\mathbf {x} ),\mathbf {x} \in \Omega ,$ (2.38.a) ${\dot {\boldsymbol {u}}}(\mathbf {x} ,0)={\dot {\boldsymbol {u}}}_{0}(\mathbf {x} ),\mathbf {x} \in \Omega ,$ (2.38.b) ${H}(\mathbf {x} ,0)={H}_{0}(\mathbf {x} ),\mathbf {x} \in \Omega {.}$ (2.38.c)

The strain history field, ${\textstyle {H}}$, could be used to model initial fracture surfaces (see appendix A of [20]).

2.4 Volumes definition

For a given set of centroids, denoted ${\textstyle \mathbf {x} _{i}}$, the discrete volumes are spheres (disks in 2D) with radii ${\textstyle r_{i}}$ truncated by planes orthogonal to the line connecting the centroids ${\textstyle \mathbf {x} _{i}}$ and ${\textstyle \mathbf {x} _{i}}$ at the following point

${\boldsymbol {q}}_{ij}=\mathbf {x} _{i}+{\frac {1}{2}}\left(1+{\frac {r_{i}^{2}-r_{j}^{2}}{||\mathbf {x} _{i}-\mathbf {x} _{j}||^{2}}}\right)(\mathbf {x} _{i}-\mathbf {x} _{j}),$

(2.39)

the point ${\textstyle {\boldsymbol {q}}_{ij}}$ is in the middle of ${\textstyle \mathbf {x} _{i}}$ and ${\textstyle \mathbf {x} _{j}}$ if ${\textstyle r_{i}}$ is equal to ${\textstyle r_{j}}$. Formally, the discrete volumes of the partition ${\textstyle P_{h}}$ are defined by

$V_{i}=\left\{\mathbf {x} \in \Omega {\Big |}||\mathbf {x} -\mathbf {x} _{i}||\leq r_{i},(\mathbf {x} -{\boldsymbol {q}}_{ij})\cdot (\mathbf {x} _{i}-\mathbf {x} _{j})\leq 0,\forall i\neq j\right\}.$

(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 ${\textstyle V_{i}}$ with respect to the remaining volumes ${\textstyle V_{j}}$, and the right side shows the discrete volumes forming a continuum in the domain, ${\textstyle \Omega }$.

Figure 2: The discrete volumes, $V_{i}$, are defined by its radii $r_{i}$ and the planes orthogonal to the lines connecting the centroids $\mathbf {x} _{i}$ and $\mathbf {x} _{j}$ at the point ${\boldsymbol {q}}_{ij}$, for all $i\neq j$.

Figure 3: The left side shows four discrete volumes colliding, the elastic response between two volumes is proportional to the size of the face formed, the domain of the volume $V_{i}$ (let $i=1$) is bounded by the remaining volumes $V_{j}$ (let $j=2,3,4$). The right side illustrates the partition, $P_{h}$, of the domain, $\Omega$, with discrete volumes (forming a continuum).

The mass of the volumes is time-invariant and its center of mass is assumed to be the centroid. To enforce these assumptions, we associate a mass, denoted ${\textstyle m_{i}}$, an initial density, ${\textstyle \rho _{i}^{o}}$, and an initial volume, ${\textstyle V_{i}^{o}}$, to the discrete volumes, such quantities are calculated as

$m_{i}=\int _{V_{i}^{o}}\rho dV,\rho _{p}^{o}={\frac {m_{i}}{V_{i}^{o}}}.$

(2.41)

Then, the density associated to the discrete volumes at any time, denoted ${\textstyle \rho _{i}}$, is given by

$\rho _{p}=\left(\rho _{i}^{o}\right){\frac {V_{p}^{o}}{V_{p}}}.$

(2.42)

Figure 4 shows the density of ${\textstyle V_{i}}$ 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, ${\textstyle \mathbf {n} _{ij}}$, 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,

$\oint _{e_{ij}}{\boldsymbol {\sigma }}\cdot \mathbf {n} _{ij}dS=0{\hbox{if }}e_{ij}{\hbox{ is curved.}}$

(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 non-linear 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.

3 First equation of motion

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 simplex-wise polynomial approximation, the pair-wise 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.

3.1 Discretization of domain into control volumes

The domain ${\textstyle \Omega }$ is discretized into ${\textstyle N}$ control volumes, denoted ${\textstyle V_{i}}$, using the Control Volume Function Approximation (CVFA) proposed by Li et al [24,25]. The partition ${\textstyle P_{h}}$ of ${\textstyle \Omega }$ is defined by

$P_{h}=\bigcup _{i=1}^{N}V_{i},~~~{\hbox{with}}~~~V_{i}\cap V_{j}=\emptyset ,~~~i\neq j,$

(3.1)

where the boundary of each control volume, ${\textstyle \partial V_{i}}$, is composed by ${\textstyle N_{i}}$ flat faces, denoted ${\textstyle e_{ij}}$,

$\partial V_{i}=\bigcup _{j=1}^{N_{i}}e_{ij},~~~{\hbox{with}}~~~e_{ij}\cap e_{ik}=\emptyset ,~~~j\neq k.$

(3.2)

Figure 6 illustrates the partition ${\textstyle P_{h}}$ of ${\textstyle \Omega }$ into ${\textstyle N}$ control volumes defined in the equations (3.1) and (3.2).

Figure 6: The partition $P_{h}$ is the discretization of the domain $\Omega$ into $N$ control volumes. The boundary of the control volumes, $\partial V_{i}$, is conformed by $N_{i}$ flat faces, denoted $e_{ij}$. The unit vector $\mathbf {n} _{ij}$ is normal to the face $e_{ij}$. The faces of the volumes adjacent to the boundary $\Gamma _{N}$ are integrated using the condition ${\boldsymbol {b}}_{N}$.

Figure 7 shows a three dimensional control volume.

Figure 7: The boundary $\partial V_{i}$ of the three dimensional control volume $V_{i}$ is subdivided into $N_{i}$ flat faces, denoted $e_{ij}$. The unit vector $\mathbf {n} _{ij}$ is normal to the face $e_{ij}$.

Every control volume ${\textstyle V_{i}}$ must have a calculation point

$\mathbf {x} _{i}\in V_{i}\cup \partial V_{i},$

(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 ${\textstyle \Gamma _{D}}$, it is convenient to establish the calculation point over the corresponding boundary face,

$\mathbf {x} _{i}\in \partial V_{i}\cap \Gamma _{D},$

(3.4)

in order to set the Dirichlet condition directly on the point.

3.2 Control volumes integration

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 ${\textstyle {\ddot {\boldsymbol {u}}}=0}$, later we will remove this assumption.

We begin by integrating the stress divergence over the control volume

$\int _{V_{i}}\nabla \cdot {\boldsymbol {\sigma }}~dV=0,$

(3.5)

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

$\int _{\partial V_{i}}{\boldsymbol {\sigma }}\mathbf {n} ~dS=0.$

(3.6)

The evaluation of the surface integrals is based on the approximation of the displacement field inside the neighborhood of the volume, denoted ${\textstyle {B}_{i}}$,

${\boldsymbol {u}}_{i}(\mathbf {x} )=\sum \limits _{q\in {B}_{i}}\varphi _{q}\mathbf {x} _{q},$

(3.7)

making use of a group of piece-wise polynomial interpolators, denoted ${\textstyle \varphi _{q}}$. 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 ${\textstyle ~{\boldsymbol {\sigma }}}$, we rewrite the stress normal to the boundary as

${\boldsymbol {\sigma }}\mathbf {n} ={\begin{bmatrix}{\boldsymbol {\sigma }}_{[11]}&{\boldsymbol {\sigma }}_{[12]}\\{\boldsymbol {\sigma }}_{[12]}&{\boldsymbol {\sigma }}_{[22]}\end{bmatrix}}{\begin{bmatrix}\mathbf {n} _{[1]}\\\mathbf {n} _{[2]}\end{bmatrix}}={\begin{bmatrix}\mathbf {n} _{[1]}&&\mathbf {n} _{[2]}\\&\mathbf {n} _{[2]}&\mathbf {n} _{[1]}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {\sigma }}_{[11]}\\{\boldsymbol {\sigma }}_{[22]}\\{\boldsymbol {\sigma }}_{[12]}\end{bmatrix}}=\mathbf {T} {\vec {\boldsymbol {\sigma }}},$

(3.8)

where ${\textstyle \mathbf {T} }$ is the face orientation matrix and ${\textstyle ~{\vec {\boldsymbol {\sigma }}}}$ is the engineering stress vector. Developing the stress definition (2.8) component-wise we can decompose it into the constitutive matrix, denoted ${\textstyle \mathbf {D} }$, and the engineering strain vector, denoted ${\textstyle {\vec {\boldsymbol {\varepsilon }}}}$, as follows

${\vec {\boldsymbol {\sigma }}}={\begin{bmatrix}{\boldsymbol {\sigma }}_{[11]}\\{\boldsymbol {\sigma }}_{[22]}\\{\boldsymbol {\sigma }}_{[12]}\end{bmatrix}}={\begin{bmatrix}2\mu ~{\boldsymbol {\varepsilon }}_{[11]}+\lambda \left({\boldsymbol {\varepsilon }}_{[11]}+{\boldsymbol {\varepsilon }}_{[22]}\right)\\2\mu ~{\boldsymbol {\varepsilon }}_{[22]}+\lambda \left({\boldsymbol {\varepsilon }}_{[11]}+{\boldsymbol {\varepsilon }}_{[22]}\right)\\2\mu ~{\boldsymbol {\varepsilon }}_{[12]}\end{bmatrix}}$ (3.9) $={\begin{bmatrix}\left(2\mu +\lambda \right)&\lambda &\\\lambda &\left(2\mu +\lambda \right)&\\&&\mu \end{bmatrix}}{\begin{bmatrix}{\boldsymbol {\varepsilon }}_{[11]}\\{\boldsymbol {\varepsilon }}_{[22]}\\2{\boldsymbol {\varepsilon }}_{[12]}\end{bmatrix}}=\mathbf {D} {\vec {\boldsymbol {\varepsilon }}},$ (3.10)

then the components of the strain vector are retrieved from the equation (2.1), and it is decomposed into the matrix differential operator ${\textstyle \mathbf {S} }$ and the displacement function ${\textstyle {\boldsymbol {u}}}$.

${\vec {\boldsymbol {\varepsilon }}}={\begin{bmatrix}{\boldsymbol {\varepsilon }}_{[11]}\\{\boldsymbol {\varepsilon }}_{[22]}\\2{\boldsymbol {\varepsilon }}_{[12]}\end{bmatrix}}={\begin{bmatrix}\displaystyle {\frac {\partial {\boldsymbol {u}}_{[1]}}{\partial \mathbf {x} _{[1]}}}\\\displaystyle {\frac {\partial {\boldsymbol {u}}_{[2]}}{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial {\boldsymbol {u}}_{[1]}}{\partial \mathbf {x} _{[2]}}}+\displaystyle {\frac {\partial {\boldsymbol {u}}_{[2]}}{\partial \mathbf {x} _{[1]}}}\end{bmatrix}}={\begin{bmatrix}\displaystyle {\frac {\partial }{\partial \mathbf {x} _{[1]}}}&\\&\displaystyle {\frac {\partial }{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial }{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial }{\partial \mathbf {x} _{[1]}}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {u}}_{[1]}\\{\boldsymbol {u}}_{[2]}\end{bmatrix}}=\mathbf {S} {\boldsymbol {u}},$

(3.11)

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

${\boldsymbol {\sigma }}\mathbf {n} =\mathbf {T} {\vec {\boldsymbol {\sigma }}}=\mathbf {T} \mathbf {D} {\vec {\boldsymbol {\varepsilon }}}=\mathbf {T} \mathbf {D} \mathbf {S} {\boldsymbol {u}},$

(3.12)

where ${\textstyle \mathbf {T} \mathbf {D} \mathbf {S} }$ is the stiffness of the volume boundary.

Once the displacement field is decoupled, we rewrite the equation (3.6) as

$\int _{\partial V_{i}}\mathbf {T} \mathbf {D} \mathbf {S} {\boldsymbol {u}}~dS=0.$

(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

$\sum \limits _{j=1}^{N_{i}}\int _{e_{ij}}\mathbf {T} \mathbf {D} \mathbf {S} {\boldsymbol {u}}~dS=0.$

(3.14)

Notice that the face orientation ${\textstyle \mathbf {T} }$ along the flat face, denoted ${\textstyle \mathbf {T} _{ij}}$, 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 ${\textstyle \mathbf {D} }$ along the flat face, denoted ${\textstyle \mathbf {D} _{ij}}$, is also considered constant. The matrix ${\textstyle \mathbf {D} _{ij}}$ is estimated from the harmonic average of the Lamé parameters assigned to the adjacent volumes, where ${\textstyle \lambda _{i}}$ and ${\textstyle \mu _{i}}$ are the properties of the volume ${\textstyle V_{i}}$,

$\mu _{ij}={\frac {2\mu _{i}\mu _{j}}{\mu _{i}+\mu _{j}}}~~~~{\hbox{and}}~~~~\lambda _{ij}={\frac {2\lambda _{i}\lambda _{j}}{\lambda _{i}+\lambda _{j}}},$

(3.15)

With ${\textstyle \mathbf {T} _{ij}}$ and ${\textstyle \mathbf {D} _{ij}}$ we simplify the equation (3.14) as

$\sum \limits _{j=1}^{N_{i}}\mathbf {T} _{ij}\mathbf {D} _{ij}\mathbf {H} _{ij}=0,$

(3.16)

where ${\textstyle \mathbf {H} _{ij}}$ is the strain integral along the flat face ${\textstyle e_{ij}}$,

$\mathbf {H} _{ij}=\int _{e_{ij}}\mathbf {S} {\boldsymbol {u}}~dS.$

(3.17)

The accuracy of the method depends on the correct evaluation of this integral.

3.3 Calculating face integrals

The surface integrals ${\textstyle \mathbf {H} _{ij}}$ along the flat faces ${\textstyle e_{ij}}$ are calculated using an auxiliary piece-wise 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 ${\textstyle \mathbf {x} _{i}}$ from the neighborhood of ${\textstyle V_{i}}$. 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 ${\textstyle {B}_{i}}$ of volume ${\textstyle V_{i}}$ as the minimum set of calculation points ${\textstyle \mathbf {x} _{j}}$ such that the simplices intersecting ${\textstyle V_{i}}$ does not change if we add another calculation point to the set. Observe that the neighborhood ${\textstyle {B}_{i}}$ does not always coincide with the set of calculation points in volumes adjacent to ${\textstyle V_{i}}$, as in most of the FV formulations. Once the neighborhood ${\textstyle {B}_{i}}$ is triangulated, we ignore the simplices with angles smaller than ${\textstyle 10}$ degrees, and the simplices formed outside the domain, which commonly appear in concavities of the boundary ${\textstyle \partial \Omega }$. The local set of simplices resulting from the neighborhood of ${\textstyle V_{i}}$ is denoted ${\textstyle P_{\alpha }}$. 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 ${\textstyle {B}_{i}}$.

Figure 8: (a) The dotted line illustrates the triangulation of the calculation points of adjacent volumes to $V_{i}$, used by most of the FV methods. (b) The dotted line shows the simplices forming the piece-wise approximation used to solve the integrals $\mathbf {H} _{ij}$ of the control volume $V_{i}$.

The face ${\textstyle e_{ij}}$ is subdivided into ${\textstyle N_{ij}}$ subfaces, denoted ${\textstyle e_{ijk}}$,

$e_{ij}=\bigcup \limits _{k=1}^{N_{ij}}e_{ijk},~~~{\hbox{with}}~~~e_{ijk}\cap e_{ijl}=\emptyset ,~~~l\neq k,$

(3.18)

these subfaces result from the intersection between ${\textstyle P_{\alpha }}$ and the control volume ${\textstyle V_{i}}$. Figure 8.b illustrates six key points of this approach, 1) the simplices are used to create a polynomial interpolation of ${\textstyle {\boldsymbol {u}}(\mathbf {x} )}$ 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 ${\textstyle V_{i}}$ and ${\textstyle V_{k}}$, 4) there are volumes that require information of non-adjacent volumes to calculate its face integrals, such as ${\textstyle V_{i}}$ requires ${\textstyle V_{k}}$, 5) the dependency between volumes is not always symmetric, which means that if ${\textstyle V_{i}}$ requires ${\textstyle V_{k}}$ does not implies that ${\textstyle V_{k}}$ requires ${\textstyle V_{i}}$, and 6) non conforming meshes are supported, as shown in the faces formed by ${\textstyle V_{a}}$, ${\textstyle V_{b}}$, ${\textstyle V_{c}}$, ${\textstyle V_{d}}$ and ${\textstyle V_{j}}$.

The integral (3.17) is now rewritten in terms of the subfaces

$\mathbf {H} _{ij}=\sum \limits _{k=1}^{N_{ij}}~\int _{e_{ijk}}\mathbf {S} {\boldsymbol {u}}~dS,$

(3.19)

Each subface ${\textstyle e_{ijk}}$ is bounded by a simplex, where the displacement ${\textstyle {\boldsymbol {u}}_{ijk}}$, and it derivatives, ${\textstyle ~(\mathbf {S} {\boldsymbol {u}})_{ijk}}$, are a polynomial interpolation. Hence the integrals in equation (3.19) are solved exactly by using the Gauss-Legendre quadrature with the required number of integration points, denoted ${\textstyle N_{g}}$, depending on the polynomial degree,

$\int _{e_{ijk}}\mathbf {S} {\boldsymbol {u}}~dS=\sum \limits _{g=1}^{N_{g}}~~w_{g}(\mathbf {S} {\boldsymbol {u}})_{ijk}|_{\mathbf {x} _{g}}.$

(3.20)

where ${\textstyle w_{g}}$ is the corresponding quadrature weight and ${\textstyle (\mathbf {S} {\boldsymbol {u}})_{ijk}|_{\mathbf {x} _{g}}}$ is the strain evaluation of the Gauss point with the proper change of interval, denoted ${\textstyle \mathbf {x} _{g}}$. Figure 9 shows the change of interval required for a 2D face. A 3D face (a polygon) must be subdivided to be integrated with a triangular quadrature.

Figure 9: (a) Blue shaded volume $V_{i}$ is being integrated. The integral over the subface $e_{ijk}$ is calculated using the polynomial approximation of shaded simplex. The integration point must be mapped to (b) Normalized space $[-1,1]$ in order to use the Gauss- Legendre quadrature.

Most of the cases, the displacement ${\textstyle ~{\boldsymbol {u}}_{ijk}}$ is interpolated inside the simplices, but in some geometrical locations these can not be created, in consequence, the displacement ${\textstyle ~{\boldsymbol {u}}_{ijk}}$ is interpolated pair-wise using the volumes adjacent to the subface ${\textstyle ~e_{ijk}}$. We discuss both strategies in the following subsections.

3.4 Simplex-wise polynomial approximation

In the general case, the simplices are formed by ${\textstyle ~({\hbox{dim}}+1)~}$ points. The points forming the simplex that is bounding the subface ${\textstyle e_{ijk}}$ are denoted ${\textstyle \mathbf {x} _{q}}$, and its displacements ${\textstyle {\boldsymbol {u}}_{q}}$.

The shape functions used for the polynomial interpolation are defined into the normalized space. A point in such space is denoted ${\textstyle {\boldsymbol {\xi }}}$, its ${\textstyle d^{th}}$ component is denoted ${\textstyle {\boldsymbol {\xi }}_{[d]}}$, and the ${\textstyle q^{th}}$ point forming the simplex is expressed ${\textstyle {\boldsymbol {\xi }}_{q}}$. The nodes of the normalized simplex are given by the origin, ${\textstyle {\boldsymbol {0}}}$, and the standard basis vectors,

${\boldsymbol {\xi }}_{q}=\left\{{\begin{array}{rl}\mathbf {e} _{q},&{\hbox{ for}}~~q\in [1,{\hbox{ dim}}],\\{\boldsymbol {0}},&{\hbox{ if}}~~q={\hbox{ dim}}+1\end{array}}\right.$

(3.21)

where ${\textstyle \mathbf {e} _{q}}$ is the ${\textstyle q^{th}}$ standard basis vector. Figures 10 and 11 illustrate the original and the normalized simplices with the corresponding node numeration for 2D and 3D respectively.

Figure 10: (a) The simplex formed by the points $\mathbf {x} _{1}$, $\mathbf {x} _{2}$ and $\mathbf {x} _{3}$ in the original space contains an interior point $\mathbf {x} _{g}$ that is mapped to (b) ${\boldsymbol {\xi }}_{g}$ into the normalized 2D-simplex formed by the points ${\boldsymbol {\xi }}_{1}$, ${\boldsymbol {\xi }}_{2}$ and ${\boldsymbol {\xi }}_{3}$.

Figure 11: (a) The 3D-simplex formed by the points $\mathbf {x} _{1}$, $\mathbf {x} _{2}$, $\mathbf {x} _{3}$ and $\mathbf {x} _{4}$ in the original space contains an interior point $\mathbf {x} _{g}$ that is mapped to (b) ${\boldsymbol {\xi }}_{g}$ into the normalized 3D-simplex formed by the points ${\boldsymbol {\xi }}_{1}$, ${\boldsymbol {\xi }}_{2}$, ${\boldsymbol {\xi }}_{3}$ and ${\boldsymbol {\xi }}_{4}$.

The shape functions, denoted ${\textstyle \varphi _{q}}$, are used to interpolate the displacement field inside the normalized simplex. Such functions are non-negative and are given by

${P}_{c}\left({\boldsymbol {\xi }}_{[q]}\right),~~~{\hbox{ if}}~~q\in [1,{\hbox{ dim}}],$ (3.22.a) $1-\sum \limits _{d=1}^{\hbox{dim}}{P}_{c}\left({\boldsymbol {\xi }}_{[d]}\right),~~~{\hbox{ for}}~~q={\hbox{ dim}}+1,$ (3.22.b)

where ${\textstyle {P}_{c}(\cdot )}$ is the homeostatic spline, which is the simplest polynomial defined in the interval ${\textstyle [0,1]}$ that have ${\textstyle ~c~}$ 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

$\sum \limits _{q=1}^{{\hbox{dim}}+1}\varphi _{q}({\boldsymbol {\xi }})=1~~~~{\hbox{for any }}{\boldsymbol {\xi }}{\hbox{ inside the simplex,}}$

(3.23)

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

$\varphi _{q}({\boldsymbol {\xi }}_{q})=1~~~~{\hbox{ for any }}{\boldsymbol {\xi }}_{q}{\hbox{ forming the simplex,}}$ (3.24) $\varphi _{q}({\boldsymbol {\xi }}_{p})=0~~~~{\hbox{ for any }}{\boldsymbol {\xi }}_{p}\neq {\boldsymbol {\xi }}_{q}{\hbox{ forming the simplex,}}$ (3.25)

The gradients of the shape functions with respect to the normalized space are denoted ${\textstyle \nabla _{\boldsymbol {\xi }}\varphi _{q}}$. The norm of the sum of such gradients is zero

$\left|\left|\sum \limits _{q=1}^{{\hbox{dim}}+1}\nabla _{\boldsymbol {\xi }}\varphi _{q}\left({\boldsymbol {\xi }}\right)\right|\right|=0~~~~{\hbox{for any }}{\boldsymbol {\xi }}{\hbox{ inside the simplex,}}$

(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, ${\textstyle {\boldsymbol {p}}\left({\boldsymbol {\xi }}\right)}$, by using the shape functions and the points forming the simplex

${\boldsymbol {p}}\left({\boldsymbol {\xi }}\right)=\sum \limits _{q=1}^{{\hbox{dim}}+1}\varphi _{q}\left({\boldsymbol {\xi }}\right)\mathbf {x} _{q},$

(3.27)

In order to calculate the normalized point, denoted ${\textstyle {\boldsymbol {\xi }}_{g}}$, associated to the integration point ${\textstyle \mathbf {x} _{g}={\boldsymbol {p}}\left({\boldsymbol {\xi }}_{g}\right)}$, we use the shape functions definitions to rewrite the equation (3.27) in matrix form

${\boldsymbol {p}}\left({\boldsymbol {\xi }}\right)=\underbrace {{\begin{bmatrix}\mathbf {x} _{3[1]}\\\mathbf {x} _{3[2]}\end{bmatrix}}+\left[{\begin{matrix}\left(\mathbf {x} _{1[1]}-\mathbf {x} _{3[1]}\right)&\left(\mathbf {x} _{2[1]}-\mathbf {x} _{3[1]}\right)\\\left(\mathbf {x} _{1[2]}-\mathbf {x} _{3[2]}\right)&\left(\mathbf {x} _{2[2]}-\mathbf {x} _{3[2]}\right)\end{matrix}}\right]{\begin{bmatrix}{P}_{c}\left({\boldsymbol {\xi }}_{[1]}\right)\\{P}_{c}\left({\boldsymbol {\xi }}_{[2]}\right)\end{bmatrix}}} _{\hbox{2D case (triangle)}}$ (3.28) $=\mathbf {x} _{({\hbox{dim}}+1)}+\mathbf {J} _{\Delta }~{P}_{c}\left({\boldsymbol {\xi }}\right),$ (3.29)

where ${\textstyle {P}_{c}\left({\boldsymbol {\xi }}\right)}$ is the vector resulting from evaluating the spline for ${\textstyle {\boldsymbol {\xi }}}$ component-wise, and ${\textstyle \mathbf {J} _{\Delta }~}$ is the distortion matrix given by the concatenation of the following column vector differences

$\mathbf {J} _{\Delta }~=\left[(\mathbf {x} _{1}-\mathbf {x} _{({\hbox{dim}}+1)}),...,(\mathbf {x} _{\hbox{dim}}-\mathbf {x} _{({\hbox{dim}}+1)})\right]$

(3.30)

Now, from equation (3.29) we retrieve the point ${\textstyle \mathbf {x} _{g}}$ as

$\mathbf {x} _{g}={\boldsymbol {p}}\left({\boldsymbol {\xi }}_{g}\right)=\mathbf {x} _{({\hbox{dim}}+1)}+\mathbf {J} _{\Delta }~{P}_{c}\left({\boldsymbol {\xi }}_{g}\right),$

(3.31)

and solving for ${\textstyle {\boldsymbol {\xi }}_{g}}$ we obtain

${\boldsymbol {\xi }}_{g}={Q}_{c}\left(\left(\mathbf {J} _{\Delta }~\right)^{-1}\left(\mathbf {x} _{g}-\mathbf {x} _{({\hbox{dim}}+1)}\right)\right),$

(3.32)

where ${\textstyle {Q}_{c}}$ is the inverse function of ${\textstyle {P}_{c}}$ applied component-wise to the product of the matrix-vector operation.

Similar to the approximation in equation (3.27), within the simplex enclosing the subface ${\textstyle e_{ijk}}$, the displacement field evaluated at ${\textstyle \mathbf {x} _{g}}$ is defined as,

${\boldsymbol {u}}_{ijk}|_{\mathbf {x} _{g}}=\sum \limits _{q=1}^{{\hbox{dim}}+1}\varphi _{q}\left({\boldsymbol {\xi }}_{g}\right)~{\boldsymbol {u}}_{q}$

(3.33)

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

$\left(\mathbf {S} {\boldsymbol {u}}\right)_{ijk}|_{\mathbf {x} _{g}}~~~=\sum \limits _{q=1}^{{\hbox{dim}}+1}\mathbf {S} \varphi _{q}({\boldsymbol {\xi }}_{g})~~{\boldsymbol {u}}_{q},$ (3.34) $=\left[{\begin{matrix}\displaystyle {\frac {\partial \varphi _{1}}{\partial \mathbf {x} _{[1]}}}&\\&\displaystyle {\frac {\partial \varphi _{1}}{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial \varphi _{1}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{1}}{\partial \mathbf {x} _{[1]}}}\end{matrix}}~~~~~{\begin{matrix}\displaystyle {\frac {\partial \varphi _{2}}{\partial \mathbf {x} _{[1]}}}&\\&\displaystyle {\frac {\partial \varphi _{2}}{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial \varphi _{2}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{2}}{\partial \mathbf {x} _{[1]}}}\end{matrix}}~~~~~{\begin{matrix}\displaystyle {\frac {\partial \varphi _{3}}{\partial \mathbf {x} _{[1]}}}&\\&\displaystyle {\frac {\partial \varphi _{3}}{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial \varphi _{3}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{3}}{\partial \mathbf {x} _{[1]}}}\end{matrix}}\right]_{|\mathbf {x} _{g}}{\begin{bmatrix}{\boldsymbol {u}}_{1[1]}\\{\boldsymbol {u}}_{1[2]}\\{\boldsymbol {u}}_{2[1]}\\{\boldsymbol {u}}_{2[2]}\\{\boldsymbol {u}}_{3[1]}\\{\boldsymbol {u}}_{3[2]}\end{bmatrix}}$ $=\mathbf {B} _{ijk}|_{\mathbf {x} _{g}}~~~{\vec {\boldsymbol {u}}}_{ijk},$ (3.35)

where ${\textstyle \mathbf {B} _{ijk}|_{\mathbf {x} _{g}}}$ captures the deformation at ${\textstyle \mathbf {x} _{g}}$, and ${\textstyle {\vec {\boldsymbol {u}}}_{ijk}}$ is the vector with the concatenated displacement components of the points forming the simplex.

In order to calculate the deformation matrix ${\textstyle \mathbf {B} _{ijk}}$, we require the derivatives of the shape functions with respect to ${\textstyle \mathbf {x} }$, denoted ${\textstyle \nabla \varphi _{q}}$. These derivatives are calculated by solving the linear systems resulting from the chain rule

$\nabla _{\boldsymbol {\xi }}\varphi _{q}={\begin{bmatrix}\displaystyle {\frac {\partial \varphi _{q}}{\partial {\boldsymbol {\xi }}_{[1]}}}\\\displaystyle {\frac {\partial \varphi _{q}}{\partial {\boldsymbol {\xi }}_{[2]}}}\end{bmatrix}}={\begin{bmatrix}\displaystyle {\frac {\partial \mathbf {x} _{[1]}}{\partial {\boldsymbol {\xi }}_{[1]}}}&\displaystyle {\frac {\partial \mathbf {x} _{[2]}}{\partial {\boldsymbol {\xi }}_{[1]}}}\\\displaystyle {\frac {\partial \mathbf {x} _{[1]}}{\partial {\boldsymbol {\xi }}_{[2]}}}&\displaystyle {\frac {\partial \mathbf {x} _{[2]}}{\partial {\boldsymbol {\xi }}_{[2]}}}\end{bmatrix}}{\begin{bmatrix}\displaystyle {\frac {\partial \varphi _{q}}{\partial \mathbf {x} _{[1]}}}\\\displaystyle {\frac {\partial \varphi _{q}}{\partial \mathbf {x} _{[2]}}}\end{bmatrix}}=\left(\nabla _{\boldsymbol {\xi }}{\boldsymbol {p}}\right)^{T}~\nabla \varphi _{q},$

(3.36)

where ${\textstyle \nabla _{\boldsymbol {\xi }}{\boldsymbol {p}}}$ is the geometric jacobian evaluated at ${\textstyle {\boldsymbol {\xi }}}$. This jacobian relates both spaces, captures the distortion of the simplex, and is derivated from equation (3.27),

$\nabla _{\boldsymbol {\xi }}{\boldsymbol {p}}=\sum \limits _{q=1}^{{\hbox{dim}}+1}\mathbf {x} _{q}\left(\nabla _{\boldsymbol {\xi }}\varphi _{q}\right)^{T},$

(3.37)

The gradients of the shape functions with respect to ${\textstyle {\boldsymbol {\xi }}}$ inside the sum are obtained straightforward once we have the spline first derivative ${\textstyle {P}_{c}'}$. Notice that the geometric jacobian is equivalent to the distortion matrix ${\textstyle \mathbf {J} _{\Delta }~}$ if and only if the homeostatic spline is ${\textstyle {P}_{c}(z)=z}$.

3.5 Pair-wise polynomial approximation

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.

Figure 12: (a) When the calculation points of volumes contiguous to the boundary are in the interior of such volumes, there will arise subfaces next to the boundary that can not be covered by any simplex. (b) Portions of the mesh formed by a queue of aligned volumes do not allow the formation of simplices through that queue and there will be whole faces not covered by any simplex.

In such cases, the displacement field within the subface ${\textstyle e_{ijk}}$ is a pair-wise polynomial approximation between the adjacent volumes, ${\textstyle \mathbf {x} _{i}}$ and ${\textstyle \mathbf {x} _{j}}$, regardless the dimension

${\boldsymbol {u}}_{ijk}(\mathbf {x} _{g})=\underbrace {\left(1-{P}_{c}\left(z_{g}\right)\right)} _{\varphi _{i}}~{\boldsymbol {u}}_{i}+\underbrace {{P}_{c}\left(z_{g}\right)} _{\varphi _{j}}~{\boldsymbol {u}}_{j},$

(3.38)

where ${\textstyle \varphi _{i}}$ and ${\textstyle \varphi _{j}}$ are the shape functions, and ${\textstyle z_{g}}$ is the normalized projection of the integration point ${\textstyle \mathbf {x} _{g}}$ into the vector which goes from ${\textstyle \mathbf {x} _{i}}$ to ${\textstyle \mathbf {x} _{j}}$, denoted ${\textstyle \mathbf {x} _{\vec {ij}}=\left(\mathbf {x} _{j}-\mathbf {x} _{i}\right),}$

$z_{g}={\frac {\left(\mathbf {x} _{g}-\mathbf {x} _{i}\right)^{T}\mathbf {x} _{\vec {ij}}}{||\mathbf {x} _{\vec {ij}}||^{2}}}.$

(3.39)

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

$\left(\mathbf {S} {\boldsymbol {u}}\right)_{ijk}|_{\mathbf {x} _{g}}~~~=\mathbf {S} \varphi _{i}({\boldsymbol {\xi }}_{g})~{\boldsymbol {u}}_{i}~~+~~\mathbf {S} \varphi _{j}({\boldsymbol {\xi }}_{g})~{\boldsymbol {u}}_{j},$ (3.40) $=\left[{\begin{matrix}\displaystyle {\frac {\partial \varphi _{i}}{\partial \mathbf {x} _{[1]}}}&\\&\displaystyle {\frac {\partial \varphi _{i}}{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial \varphi _{i}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{i}}{\partial \mathbf {x} _{[1]}}}\end{matrix}}~~~~~{\begin{matrix}\displaystyle {\frac {\partial \varphi _{j}}{\partial \mathbf {x} _{[1]}}}&\\&\displaystyle {\frac {\partial \varphi _{j}}{\partial \mathbf {x} _{[2]}}}\\\displaystyle {\frac {\partial \varphi _{j}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{j}}{\partial \mathbf {x} _{[1]}}}\end{matrix}}\right]_{|\mathbf {x} _{g}}{\begin{bmatrix}{\boldsymbol {u}}_{i[1]}\\{\boldsymbol {u}}_{i[2]}\\{\boldsymbol {u}}_{j[1]}\\{\boldsymbol {u}}_{j[2]}\end{bmatrix}}$ (3.41) $=\mathbf {B} _{ijk}|_{\mathbf {x} _{g}}~~~{\vec {\boldsymbol {u}}}_{ijk},$ (3.42)

In the general case, the gradient is not constant along the face ${\textstyle e_{ij}}$, since its normal is not necessary aligned with ${\textstyle \mathbf {x} _{\vec {ij}}}$, as illustrated in Figure 13.

Figure 13: The gradient of the pairwise approximation is not constant along the face $e_{ij}$, since its normal is not necessary aligned with $\mathbf {x} _{\vec {ij}}$. The integration point is projected into $\mathbf {x} _{\vec {ij}}$ to evaluate the deformation matrix.

This pairwise approximation must be used only when necessary because it can not capture the deformation orthogonal to ${\textstyle \mathbf {x} _{\vec {ij}}}$.

3.6 Homeostatic spline

The homeostatic spline is a function of a single variable defined from ${\textstyle z=0}$ to ${\textstyle z=1}$, denoted ${\textstyle {P}_{c}(z)}$, and curved by the parameter ${\textstyle ~c}$, which indicates the level of smoothness. This spline is the simplest polynomial with ${\textstyle ~c~}$ derivatives equal to zero at the endpoints ${\textstyle z=0}$ and ${\textstyle z=1}$. The polynomial degree is given by ${\textstyle 2c+1}$, and such a polynomial requires ${\textstyle N_{g}=c+1}$ integration points to calculate the exact integral in equation (3.20) using the Gauss-Legendre quadrature.

When designing this spline, we wanted to gain accuracy by building a piece-wise bell-shaped 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 ${\textstyle 2c+2}$ coefficients of the polynomial. The equations of this system were obtained by forcing the ${\textstyle ~c}$ derivatives to be zero at the endpoints. Once we solved the coefficients for the first twenty polynomials, from ${\textstyle c=0}$ to ${\textstyle c=19}$, we found out that the first half of such coefficients are null, and the entire polynomial can be calculated directly as

${P}_{c}(z)={\frac {1}{b_{c}}}\sum \limits _{k=1}^{1+c}(-1)^{k}~b_{k}~~z^{(k+c)},$

(3.43)

where ${\textstyle b_{k}}$ is the ${\textstyle k^{th}}$ not null coefficient

$b_{k}={\frac {1}{k+c}}\left(\prod \limits _{l=1}^{C_{k}}{\frac {\left(1+c\right)}{l}}-1\right),$

(3.44)

${\textstyle C_{k}}$ is the number of factors needed to calculate ${\textstyle b_{k}}$

$C_{k}=(c/2)-{\big |}1+(c/2)-k{\big |},$

(3.45)

and ${\textstyle b_{c}}$ is accumulation of the coefficients for normalizing the spline

$b_{c}=\sum \limits _{k=1}^{1+c}(-1)^{k}b_{k},$

(3.46)

The first derivative is simply calculated as

${P}_{c}'(z)={\frac {1}{b_{c}}}\sum \limits _{k=1}^{1+c}(-1)^{k}~b_{k}~~(k+c)~~z^{(k+c-1)}$

(3.47)

Table 1 shows the polynomials resulting from low values of ${\textstyle c}$ and Figure 14 depicts (a) the evolution of the spline as we increase the smoothness parameter from ${\textstyle c=0}$ to ${\textstyle c=6}$, and (b) the evolution of it first derivative.

Table. 1 Coefficients for the first few levels of smoothness of homeostatic spline

Smoothness	Homeostatic spline
${\textstyle c=0}$	${P}_{0}(z)=z$
$c=1$	${P}_{1}(z)=3z^{2}-2z^{3}$
$c=2$	${P}_{2}(z)=10z^{3}-15z^{4}+6z^{5}$
$c=3$	${P}_{3}(z)=35z^{4}-84z^{5}+70z^{6}-20z^{7}$

Smoother splines produces higher order polynomials which increases the accuracy of the stress field approximation. This feature is specially important when solving non-linear problems sensibles to the stress field.

Figure 14: (a) The evolution of the homeostatic spline from $c=0$ to $c=6$ 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 ${\textstyle [0,1]}$, the inverse function is not defined in that points. However, we estimate a pseudo-inverse within this interval, ${\textstyle {Q}_{c}\approx {P}_{c}^{-1}}$, by finding the coefficients of a polynomial of the same degree, ${\textstyle 2c+1}$, 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

${Q}_{c}(0)={P}_{c}(0)=0,~~~{Q}_{c}(1)={P}_{c}(1)=1,~~~{\hbox{and}}~~~{Q}_{c}'(0.5)={\frac {1}{{P}_{c}'(0.5)}}$

(3.48)

The higher derivatives in the midpoint are forced to be zero. Once we calculated the coefficients for the first twenty polynomials, from ${\textstyle c=0}$ to ${\textstyle c=19}$, we found out that the pseudo-inverse can be approximated directly from the following formulae

${Q}(z)=~a_{1}~z~~+~~(a_{1}-1)(2c+1)~\sum \limits _{k=1}^{2c}(-1)^{k}~a_{k}~z^{(k+1)}$

(3.49)

where ${\textstyle a_{1}}$ is the coefficient for ${\textstyle z}$, and ${\textstyle a_{k}}$ is the factor that distinguish higher order coefficients. Such terms are calculated as

$a_{1}=\left({\frac {c}{2{\sqrt {2}}}}+1\right)^{2},~~~{\hbox{and}}~~~a_{k}=2^{(k-1)}~\prod _{l=1}^{k-1}\left({\frac {2c-l}{2+l}}\right),$

(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 pseudo-inverse ${Q}_{c}$ for the first seven levels of smoothness. The slope at the midpoint exposes the null higher derivatives requirement when increasing the polynomial order.

Figure 16 shows the shape functions for the 2D case. The top displays the last node function and the bottom the first node function, the function of the second node is equivalent to that of the first one. The columns separate the first three levels of smoothness. Top and bottom functions coincides at the edges in order to create a continuous field, but only the bottom functions decay uniformly from the node to the opposite edge. The shape functions with ${\textstyle c=0}$ are the only case where all the shape functions are indistinguishable, these are planes.

Figure 16: For the bidimensional case, the top displays the last node function and the bottom the first node function, the function of the second node is equivalent to that of the first one. The columns separate the first three levels of smoothness.

Figure 17 shows the magnitude of the gradient with respect to the normalized space. With the same tabular configuration of Figure 16, the columns separate the first three levels of smoothness, the top displays the last node gradient and the bottom the first node gradient, the gradient of the second node is equivalent to that of the first one. Only the gradient magnitude at the bottom has a uniform variation from the node to the opposite face, and the value of the node does not contribute to the value of such a face. On the contrary, in the top can be observed that the value of the node contributes to the gradient at the opposite face, which means that using ${\textstyle ~c>0~}$ the continuity on the stress field is only guaranteed at the calculation points, but not in the simplices edges.

Figure 17: For the bidimensional case, the top displays the last node gradient magnitudes and the bottom the first node gradient magnitudes, the gradient magnitudes of the second node is equivalent to that of the first one. The columns separate the first three levels of smoothness.

3.7 Assembling volume's equation

By using the simplex-wise (3.35) or the pair-wise (3.42) approximation, the strain face integral (3.19) is reformulated as

$\mathbf {H} _{ij}=\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {B} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk},$

(3.51)

then, the volume equilibrium equation (3.16) is

$\sum \limits _{j=1}^{N_{i}}\mathbf {T} _{ij}\mathbf {D} _{ij}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {B} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk}=0,$

(3.52)

reordering the terms we obtain

$\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {K} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk}=0,$

(3.53)

where the matrix

$\mathbf {K} _{ijk}|_{\mathbf {x} _{g}}=\mathbf {T} _{ij}\mathbf {D} _{ij}\mathbf {B} _{ijk}|_{\mathbf {x} _{g}},$

(3.54)

is the stiffness contribution at the integration point ${\textstyle \mathbf {x} _{g}}$, within the subface ${\textstyle e_{ijk}}$ when integrating the ${\textstyle i^{th}}$ volume. Observe that the stiffness matrix ${\textstyle \mathbf {K} _{ijk}}$ is rectangular and the resulting global stiffness matrix is asymmetric.

3.8 Boundary conditions

The Neumann boundary conditions are imposed over the volume faces ${\textstyle e_{ij}}$ 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),

$\int _{e_{ij}}\mathbf {T} \mathbf {D} \mathbf {S} {\boldsymbol {u}}~dS=\int _{e_{ij}}{\boldsymbol {b}}_{N}(\mathbf {x} )~dS$

(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),

${\boldsymbol {u}}_{i}={\boldsymbol {u}}_{D}(\mathbf {x} _{i}),$

(3.56)

Since the Dirichlet conditions are imposed directly on the calculation points, these points must be located along the face ${\textstyle e_{ij}}$ which intersects the boundary with the condition ${\textstyle \Gamma _{D}}$.

3.9 Special cases

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 FV-FEM correlation.

Figure 18: (a) The initial mesh is equivalent to the Voronoi diagram and the Voronoi centres correspond to the calculation points $\mathbf {x} _{i}$. (b) The initial mesh is generated from a FEM-like triangular mesh. The calculation points $\mathbf {x} _{i}$ 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.

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 ${\textstyle \mathbf {x} _{i}}$. Hence, the subdivision of the neighborhood ${\textstyle {B}_{i}}$ 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 ${\textstyle e_{ijk}}$ using pair-wise approximations can be exactly integrated with the midpoint rule, since the faces are orthogonal to the vector joining the calculation points ${\textstyle \mathbf {x} _{\vec {ij}}}$, and the derivatives along the subface are constants.

In the second case, the initial mesh is generated from a FEM-like triangular mesh and the approximations are assumed to be linear. In such a case, the calculation points ${\textstyle \mathbf {x} _{i}}$ 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 cell-centred 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.

4 Second equation of motion

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 ${\textstyle h}$ which controls the smoothness of damage field solution. We observe that a length parameter proportional to the average size of control volumes, denoted ${\textstyle \Delta \mathbf {x} }$, produces accurate results, these mesh size is taken as

$\Delta \mathbf {x} =\left({\frac {1}{N}}\sum \limits _{i=1}^{N}V_{i}\right)^{\frac {1}{\hbox{dim}}}$

(4.1)

Figure 19 illustrates the graphical meaning of scale length parameter.

Figure 19: a) Damage field above control volumes shows how the crack arises along volumes boundaries. b) Scale length parameter $h$ 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 ${\textstyle P_{h}}$ and interpolators, simple-wise and pair-wise approximations also apply for the damage field.

4.1 Discretization

We start by integrating the strong form equation of motion (2.36.b) over the control volumes of the partition ${\textstyle P_{h}}$,

$\int _{V_{i}}\left(1+{\frac {2h}{G}}{H}\right){\boldsymbol {d}}~~dV-\int _{V_{i}}h^{2}\nabla ^{2}{\boldsymbol {d}}~~dV=\int _{V_{i}}{\frac {2h}{G}}{H}~~dV,$

(4.2)

using the divergence theorem on the second integral we get

$\int _{V_{i}}\left(1+{\frac {2h}{G}}{H}\right){\boldsymbol {d}}~~dV-h^{2}\int _{\partial V_{i}}\nabla {\boldsymbol {d}}\cdot \mathbf {n} ~dS=\int _{V_{i}}{\frac {2h}{G}}{H}~~dV$

(4.3)

Since ${\textstyle {G}}$ is a material property, we assume that it is constant along the control volume, and dividing the first integral in two terms we obtain

$\int _{V_{i}}{\boldsymbol {d}}~~dV+{\frac {2h}{G}}\int _{V_{i}}{H}{\boldsymbol {d}}~~dV-h^{2}\int _{\partial V_{i}}\nabla {\boldsymbol {d}}\cdot \mathbf {n} ~dS={\frac {2h}{G}}\int _{V_{i}}{H}~~dV$

(4.4)

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

$V_{i}=\bigcup \limits _{j=1}^{N_{i}}\bigcup \limits _{k=1}^{N_{ij}}V_{ijk},~~~{\hbox{with no subvolume intersections}},$

(4.5)

where ${\textstyle V_{ijk}}$ are the pyramids (triangles in 2D) which base corresponds to the subfaces ${\textstyle e_{ijk}}$ and its apex is the calculation point ${\textstyle \mathbf {x} _{i}}$ as illustrated in figure 20.

Figure 20: The control volume $V_{i}$ is partitioned into pyramids $V_{ijk}$, which turns to be triangles in 2D. Pyramids bases correspond to the subfaces $e_{ijk}$ resulting from the intersection with the local Delaunay triangulation, and all of them share the calculation point $\mathbf {x} _{i}$ as its apex.

The surface integral is solved along subfaces ${\textstyle e_{ijk}}$ defined in (3.18), and the remaining volume integrals are solved using partition (4.5),

$\int _{V_{i}}{\boldsymbol {d}}~~dV~~+\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\left({\frac {2h}{G}}\int _{V_{ijk}}{H}{\boldsymbol {d}}~~dV-h^{2}\int _{e_{ijk}}\nabla {\boldsymbol {d}}\cdot \mathbf {n} ~dS\right)$ $=\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}~~{\frac {2h}{G}}\int _{V_{ijk}}{H}~~dV$

(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,

${\boldsymbol {d}}_{ijk}|_{\mathbf {x} _{g}}=\sum \limits _{q=1}^{{\hbox{dim}}+1}\varphi _{q}\left({\boldsymbol {\xi }}_{g}\right)~{\boldsymbol {d}}_{q},$ (4.7) $=\left({\vec {\varphi }}|_{{\boldsymbol {\xi }}_{g}}\right)^{T}~{\vec {\boldsymbol {d}}}_{ijk},$ (4.8)

where ${\textstyle {\boldsymbol {\xi }}_{g}}$ is the point corresponding to ${\textstyle \mathbf {x} _{g}}$ in the normalized space, ${\textstyle {\vec {\varphi }}|_{{\boldsymbol {\xi }}_{g}}}$ is the vector containing the shape functions evaluated at ${\textstyle {\boldsymbol {\xi }}_{g}}$, and ${\textstyle {\vec {\boldsymbol {d}}}_{ijk}}$ 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

$\nabla {\boldsymbol {d}}_{ijk}|_{\mathbf {x} _{g}}=\sum \limits _{q=1}^{{\hbox{dim}}+1}\nabla \varphi _{q}\left({\boldsymbol {\xi }}_{g}\right)~{\boldsymbol {d}}_{q},$

(4.9)

where ${\textstyle \nabla \varphi _{q}}$ is calculated from the chain rule in (3.36). Now the equation is fully discretized, the next step is to solve the integrals.

4.2 Assembling system of equations

The first integral in equation (4.6) is approximated using the midpoint rule,

$\int _{V_{i}}{\boldsymbol {d}}~~dV=V_{i}~{\boldsymbol {d}}_{i}$

(4.10)

where ${\textstyle {\boldsymbol {d}}_{i}}$ is the damage estimated at calculation point ${\textstyle \mathbf {x} _{i}}$. Due to the simple nature of polygonal subvolumes ${\textstyle V_{ijk}}$, we can always reduce them to simplices in order to use the Gauss-Legendre quadrature to solve the volume integrals, in 2D is straightforward,

$\int _{V_{ijk}}{H}{\boldsymbol {d}}~~dV=\sum \limits _{p=1}^{N_{p}}w_{p}\left({H}_{ijk}|_{\mathbf {x} _{p}}\right)\left({\boldsymbol {d}}_{ijk}|_{\mathbf {x} _{p}}\right),$ (4.11) $=\sum \limits _{p=1}^{N_{p}}w_{p}\left({H}_{ijk}|_{\mathbf {x} _{p}}\right)\left({\vec {\varphi }}|_{{\boldsymbol {\xi }}_{p}}\right)^{T}~{\vec {\boldsymbol {d}}}_{ijk},$ (4.12) $\int _{V_{ijk}}{H}~~dV=\sum \limits _{p=1}^{N_{p}}w_{p}{H}_{ijk}|_{\mathbf {x} _{p}},$ (4.13)

where ${\textstyle N_{p}}$ is the number of points in the quadrature, and ${\textstyle {H}_{ijk}|_{\mathbf {x} _{p}}}$ is the strain history field evaluated with the strain information of the simplex corresponding to subface ${\textstyle e_{ijk}}$. Last but not least, we solve the surface integral that unfolds the damage gradient defined in (4.9), using again Gauss-Legendre quadrature

$\int _{e_{ijk}}\nabla {\boldsymbol {d}}\cdot \mathbf {n} ~dS=\sum \limits _{g=1}^{N_{g}}w_{g}~\mathbf {n} _{ijk}~\cdot ~\nabla {\boldsymbol {d}}|_{\mathbf {x} _{g}},$ (4.14) $=\sum \limits _{g=1}^{N_{g}}w_{g}~\mathbf {n} _{ijk}\cdot \left(\sum \limits _{q=1}^{{\hbox{dim}}+1}\nabla \varphi _{q}\left({\boldsymbol {\xi }}_{g}\right)~{\boldsymbol {d}}_{q}\right),$ (4.15) $=\sum \limits _{g=1}^{N_{g}}w_{g}~\left[\mathbf {n} _{[1]}~~\mathbf {n} _{[2]}\right]_{ijk}\underbrace {{\begin{bmatrix}\displaystyle {\frac {\partial \varphi _{1}}{\partial \mathbf {x} _{[1]}}}&\displaystyle {\frac {\partial \varphi _{2}}{\partial \mathbf {x} _{[1]}}}&\displaystyle {\frac {\partial \varphi _{3}}{\partial \mathbf {x} _{[1]}}}\\\displaystyle {\frac {\partial \varphi _{1}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{2}}{\partial \mathbf {x} _{[2]}}}&\displaystyle {\frac {\partial \varphi _{3}}{\partial \mathbf {x} _{[2]}}}\end{bmatrix}}_{|_{\mathbf {x} _{g}}}} _{\mathbf {Z} _{ijk}|_{\mathbf {x} _{g}}}{\begin{bmatrix}{\boldsymbol {d}}_{1}\\{\boldsymbol {d}}_{2}\\{\boldsymbol {d}}_{3}\end{bmatrix}},$ (4.16) $=\sum \limits _{g=1}^{N_{g}}w_{g}~\left(\mathbf {n} _{ijk}\right)^{T}~\mathbf {Z} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {d}}}_{ijk}$ (4.17) (4.18)

where ${\textstyle \mathbf {Z} _{ijk}|_{\mathbf {x} _{g}}}$ is the matrix containing the derivatives of the shape functions evaluated at ${\textstyle \mathbf {x} _{g}}$. 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

$V_{i}~{\boldsymbol {d}}_{i}~~+\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\left({\frac {2h}{G}}\sum \limits _{p=1}^{N_{p}}w_{p}\left({H}_{ijk}|_{\mathbf {x} _{p}}\right)\left({\vec {\varphi }}|_{{\boldsymbol {\xi }}_{p}}\right)^{T}-h^{2}\sum \limits _{g=1}^{N_{g}}w_{g}~(\mathbf {n} _{ijk})^{T}~\mathbf {Z} _{ijk}|_{\mathbf {x} _{g}}\right)~{\vec {\boldsymbol {d}}}_{ijk}$ $=\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}~~{\frac {2h}{G}}\sum \limits _{p=1}^{N_{p}}w_{p}~{H}_{ijk}|_{\mathbf {x} _{p}}$

(4.19)

The damage of the ${\textstyle j^{th}}$ face produced by excesive loads is captured by vector

${\vec {\mathbf {D} }}_{ijk}=\left({\frac {2h}{G}}\sum \limits _{p=1}^{N_{p}}w_{p}\left({H}_{ijk}|_{\mathbf {x} _{p}}\right){\vec {\varphi }}|_{{\boldsymbol {\xi }}_{p}}-h^{2}\sum \limits _{g=1}^{N_{g}}w_{g}~\left(\mathbf {Z} _{ijk}|_{\mathbf {x} _{g}}\right)^{T}~\mathbf {n} _{ijk}\right),$

(4.20)

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

$W_{i}={\frac {2h}{G}}\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{p=1}^{N_{p}}w_{p}~{H}_{ijk}|_{\mathbf {x} _{p}}$

(4.21)

Now we can rewrite the damage equation (4.19) for the ${\textstyle i^{th}}$ control volume as follows

$V_{i}~{\boldsymbol {d}}_{i}~~+\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\left({\vec {\mathbf {D} }}_{ijk}\right)^{T}~{\vec {\boldsymbol {d}}}_{ijk}=~~W_{i}$

(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 ${\textstyle {H}}$.

5 Time discretization

In this chapter we will remove the assumption done in equation (3.5) about null acceleration, ${\textstyle {\ddot {\boldsymbol {u}}}=0}$, 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 non-linear 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.

5.1 Time variation

In order to analyze the dynamic component of elasticity equation (2.36.a) we define the stress state of control volume ${\textstyle i}$ as a function of time,

$S_{i}(t)=\nabla \cdot {\boldsymbol {\sigma }}_{i},$

(5.1)

with the intention of considering internal forces in the approximation,

$S_{i}(t)=\rho {\ddot {\boldsymbol {u}}}_{i}$

(5.2)

Equation (5.2) is an ordinary differential equation that can be solved analytically for a time step ${\textstyle t\in [0,\Delta t]}$ with the following Dirichlet conditions

${\boldsymbol {u}}_{i}(0)={\boldsymbol {u}}_{i}^{0},$ (5.3.a) ${\boldsymbol {u}}_{i}(\Delta t)={\boldsymbol {u}}_{i},$ (5.3.b)

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

$S_{i}(t)=\left(1-{P}\left({\frac {t}{\Delta t}}\right)\right)S_{i}^{0}+{P}\left({\frac {t}{\Delta t}}\right)~S_{i},~~~~~~t\in [0,\Delta t],$

(5.4)

Figure 21: The time variation of the stress state is defined by the shape function ${P}$ that interpolates the stress states of two contiguous time steps. During this work we found that continuous functions like the shown here produces more accurate approximations in the stress field than the numerical schemes that does not consider this variation.

where ${\textstyle S_{i}^{0}}$ is the stress state calculated at ${\textstyle t=0}$, ${\textstyle S_{i}}$ is the stress state which will be estimated at time ${\textstyle t=\Delta t}$, and ${\textstyle {P}(\cdot )}$ is the shape function modelling time variation between concecutive stress states. This shape function has only two constraints ${\textstyle {P}(0)=0}$ and ${\textstyle {P}(1)=1}$, for that reason we use ${\textstyle \Delta t}$ 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 ${\textstyle {P}}$ that is used as interpolator between the value at two contiguous time steps.

5.2 Analytical solution

Using the asumption in (5.4), we get the analytical solution of the equation (5.2) for the interval ${\textstyle t\in [0,\Delta t]}$ by means of the Laplace transform (see appendix A for details),

${\boldsymbol {u}}_{i}(t)={\boldsymbol {u}}_{i}^{0}+t~{\dot {\boldsymbol {u}}}_{i}^{0}+{\frac {1}{2\rho }}~t^{2}~S_{i}^{0}+{\frac {1}{\rho }}{C}_{P}\left(t\right)\left(S_{i}-S_{i}^{0}\right),$

(5.5)

where ${\textstyle {\dot {\boldsymbol {u}}}_{i}^{0}}$ is the velocity at time ${\textstyle t=0}$, and ${\textstyle {C}_{P}(t)}$ is the convolution between the spline ${\textstyle {P}(t/\Delta t)}$ and the function ${\textstyle t}$, as defined in appendix A.

By setting the second boundary condition, ${\textstyle {\boldsymbol {u}}_{i}(\Delta t)={\boldsymbol {u}}_{i}}$, into the analytical solution (5.5), we can find the velocity required to fulfill both Dirichlet conditions

${\dot {\boldsymbol {u}}}_{i}^{0}=\left({\frac {{\boldsymbol {u}}_{i}-{\boldsymbol {u}}_{i}^{0}}{\Delta t}}\right)-{\frac {1}{2\rho }}~\Delta t~S_{i}^{0}-{\frac {1}{\rho ~\Delta t}}~{C}_{P}^{\Delta t}~\left(S_{i}-S_{i}^{0}\right),$

(5.6)

where ${\textstyle {C}_{P}^{\Delta t}}$ is the convolution evaluated at ${\textstyle \Delta t}$. Thus, we replace equation (5.6) into (5.5) to get the analytical solution in terms of the known Dirichlet conditions

${\boldsymbol {u}}_{i}(t)=~~{\boldsymbol {u}}_{i}^{0}+t~\left(\left({\frac {{\boldsymbol {u}}_{i}-{\boldsymbol {u}}_{i}^{0}}{\Delta t}}\right)-{\frac {1}{2\rho }}~\Delta t~S_{i}^{0}-{\frac {1}{\rho ~\Delta t}}~{C}_{P}^{\Delta t}~\left(S_{i}-S_{i}^{0}\right)\right)~+$ $t^{2}~{\frac {1}{2\rho }}S_{i}^{0}~+{C}_{P}\left(t\right){\frac {1}{\rho }}\left(S_{i}-S_{i}^{0}\right),$

(5.7)

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

${\dot {\boldsymbol {u}}}_{i}(t)=~\left({\frac {{\boldsymbol {u}}_{i}-{\boldsymbol {u}}_{i}^{0}}{\Delta t}}\right)-{\frac {1}{2\rho }}~\Delta t~S_{i}^{0}-{\frac {1}{\rho ~\Delta t}}~{C}_{P}^{\Delta t}~\left(S_{i}-S_{i}^{0}\right)~+$ $t~{\frac {1}{\rho }}S_{i}^{0}+~{\frac {1}{\rho }}~~{\dot {C}}_{P}\left(t\right)\left(S_{i}-S_{i}^{0}\right)$

(5.8)

where ${\textstyle {\dot {C}}_{P}}$ is the time derivative of ${\textstyle {C}_{P}}$. 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 ${\textstyle t\in [-\Delta h,0]}$,

${\dot {\boldsymbol {u}}}_{i}^{0}=~\left({\frac {{\boldsymbol {u}}_{i}^{0}-{\boldsymbol {u}}_{i}^{00}}{\Delta h}}\right)-{\frac {\Delta h}{\rho }}\left({\frac {1}{2}}~S_{i}^{0}-S_{i}^{00}\right)+{\frac {1}{\rho }}\left({\dot {C}}_{P}^{\Delta h}-{\frac {1}{\Delta h}}{C}_{P}^{\Delta h}\right)\left(S_{i}^{0}-S_{i}^{00}\right),$

(5.9)

where ${\textstyle {\boldsymbol {u}}_{i}^{00}={\boldsymbol {u}}_{i}(-\Delta h)}$ and ${\textstyle S_{i}^{00}=S_{i}(-\Delta h)}$. Finally, we replace equation (5.9) into (5.5) in order to get an analytical solution for ${\textstyle t\in [0,\Delta t]}$ as a function of two history system states,

${\boldsymbol {u}}_{i}(t)=~~{\boldsymbol {u}}_{i}^{0}+t~\left({\frac {{\boldsymbol {u}}_{i}^{0}-{\boldsymbol {u}}_{i}^{00}}{\Delta h}}\right)-t~{\frac {\Delta h}{\rho }}\left({\frac {1}{2}}~S_{i}^{0}-S_{i}^{00}\right)~+$ $t~{\frac {1}{\rho }}\left({\dot {C}}_{P}^{\Delta h}-{\frac {1}{\Delta h}}{C}_{P}^{\Delta h}\right)\left(S_{i}^{0}-S_{i}^{00}\right)~+$ $t^{2}~{\frac {1}{2\rho }}S_{i}^{0}~+{C}_{P}\left(t\right){\frac {1}{\rho }}\left(S_{i}-S_{i}^{0}\right),$

(5.10)

evaluating such an equation at ${\textstyle t=\Delta t}$, denoted ${\textstyle {\boldsymbol {u}}_{i}={\boldsymbol {u}}_{i}(\Delta t)}$, and rearranging the terms we get a numerical approximation dependent of the convolution of choosen spline,

$\rho \left({\boldsymbol {u}}_{i}-\left(1+{\frac {\Delta t}{\Delta h}}\right){\boldsymbol {u}}_{i}^{0}+{\frac {\Delta t}{\Delta h}}~{\boldsymbol {u}}_{i}^{00}\right)={C}_{P}^{\Delta t}~S_{i}+$ $\left({\frac {1}{2}}(\Delta t^{2}-\Delta t\Delta h)+\Delta t\left({\dot {C}}_{P}^{\Delta h}-{\frac {1}{\Delta h}}{C}_{P}^{\Delta h}\right)-{C}_{P}^{\Delta t}\right)S_{i}^{0}+$ $\left(\Delta t\Delta h-\Delta t\left({\dot {C}}_{P}^{\Delta h}-{\frac {1}{\Delta h}}{C}_{P}^{\Delta h}\right)\right)S_{i}^{00},$

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

5.3 Numerical scheme

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 ${\textstyle {P}}$ 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.

5.3.1 Harmonic oscillator sensibility

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

$-ku=m{\ddot {u}},$

(5.12)

where ${\textstyle k}$ is the stiffness of the system, ${\textstyle m}$ is the mass of the body and ${\textstyle u}$ is the one-dimensional displacement. The analytical solution of equation (5.12) is

$u(t)=A~cos(\omega t+\gamma ),$

(5.13)

where ${\textstyle A}$ is the oscillation amplitude, ${\textstyle \omega }$ the oscillation frequency and ${\textstyle \gamma }$ the phase, such constants are calculated in terms of material properties

$A={\sqrt {u_{0}+{\frac {m}{k}}{\dot {u}}_{0}}},$ (5.14) $\omega ={\sqrt {\frac {k}{m}}}$ (5.15) $\gamma =-arctan\left({\sqrt {\frac {m}{k}}}{\frac {{\dot {u}}_{0}}{u_{o}}}\right),$ (5.16)

with ${\textstyle u_{0}}$ and ${\textstyle {\dot {u}}_{0}}$ as initial displacement and initial velocity respectively. In our numerical tests, the one dimensional stress state, denoted ${\textstyle s}$, is assumed to be

$s(t)=-k~u(t)$

(5.17)

For simplicity, in this sensibility analysis we used a constant time step ${\textstyle \Delta h=\Delta t}$.

5.3.1.1 Central Finite Differences

By using a central finite differences scheme, equation (5.12) can be rewritten as

$-ku={\frac {m}{\Delta t^{2}}}\left(u-2u^{0}+u^{00}\right),$

(5.18)

and the solution for next time step is calculated from

$u={\frac {m}{k\Delta t^{2}}}\left(2u^{0}-u^{00}\right)\left(1+{\frac {m}{k\Delta t^{2}}}\right)^{-1}$ (5.19) $=\left({\frac {k\Delta t^{2}}{m}}+1\right)^{-1}\left(2u^{0}-u^{00}\right)$ (5.20)

To measure the relative error with respect to analytical solution, we used (5.20) to compute the solution in the interval ${\textstyle t\in [0,7]}$. To make evident the numerical drawbacks of FD we utilized a big enough ${\textstyle \Delta t=0.1}$. 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 ${\textstyle \Delta t}$ 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 ${\textstyle \Delta t}$, such an error remains almost consant for ${\textstyle \Delta t>0.06}$ since the numerical system looses all its energy in the first few time steps. In this plot we compute the comulative error as

${\hbox{Error}}(T)=\left(\int \limits _{0}^{T}{\frac {(U_{n}-U_{a})^{2}}{U_{a}}}dt\right)^{1/2},$

(5.21)

where ${\textstyle T}$ indicates the simulation duration, ${\textstyle U_{a}(t)}$ is the analytical total energy and ${\textstyle U_{n}(t)}$ is the numerical total energy.

Figure 22: Dashed lines show numerical results, while solid lines are for analytical solution. Top-left plot is the direct solution, the curve of displacement vs time. Top-right plot is the phase space, velocity vs displacement, the closing spiral tell us that numerical system is loosing energy. Bottom-left plot is the total energy in the system, the analytical solution is constant and the numerical energy decreases to zero. Bottom-right plot is the cumulative relative error for distinct values of $\Delta t$.

5.3.1.2 Linear spline

If we choose a linear shape function, ${\textstyle {P}(t)=t}$, in order to set a constant variation of the internal forces in the interval ${\textstyle [0,\Delta t]}$, the convolution and its time derivative are given by

${C}_{P}(t)={\frac {t^{3}}{6\Delta t}}~~~{\hbox{and}}~~~{\dot {C}}_{P}(t)={\frac {t^{2}}{2\Delta t}},$

(5.22)

respectively, and the resulting numerical scheme 5.11 is

${\frac {1}{6}}\left({\frac {\Delta t}{\Delta h}}\right)~S_{i}+\left({\frac {1}{3}}\left({\frac {\Delta t}{\Delta h}}+{\frac {\Delta h}{\Delta t}}\right)-{\frac {1}{2}}\right)~S_{i}^{0}+\left(1-{\frac {1}{3}}~{\frac {\Delta h}{\Delta t}}\right)~S_{i}^{00}=$ ${\frac {\rho }{\Delta t\Delta h}}\left({\boldsymbol {u}}_{i}-\left(1+{\frac {\Delta t}{\Delta h}}\right){\boldsymbol {u}}_{i}^{0}+{\frac {\Delta t}{\Delta h}}~{\boldsymbol {u}}_{i}^{00}\right),$

(5.23)

by applying the assumption of constant time steps, we reduce previous equation to

${\frac {1}{6}}~S_{i}+{\frac {1}{6}}~S_{i}^{0}+{\frac {4}{6}}~S_{i}^{00}={\frac {\rho }{\Delta t^{2}}}\left({\boldsymbol {u}}_{i}-2{\boldsymbol {u}}_{i}^{0}+~{\boldsymbol {u}}_{i}^{00}\right),$

(5.24)

then we use this numerical approximation to solve the harmonic oscillator and we get

${\frac {1}{6}}~(-k~u)+{\frac {1}{6}}~(-k~u^{0})+{\frac {4}{6}}~(-k~u^{00})={\frac {m}{\Delta t^{2}}}\left(u-2u^{0}+~u^{00}\right),$

(5.25)

and the numerical solution is given by

$u=\left(2u^{0}-~u^{00}\right)\left(1+{\frac {k\Delta t^{2}}{6m}}\right)^{-1}-\left(u^{0}+4~u^{00}\right)\left(1+{\frac {6m}{k\Delta t^{2}}}\right)^{-1},$

(5.26)

for displacement and

${\dot {u}}=\left({\frac {u-u^{0}}{\Delta t}}\right)-{\frac {k\Delta t}{6m}}\left(2u+u^{0}\right),$

(5.27)

for velocity.

In our experiments we used the same ${\textstyle \Delta t=0.1}$ than with Finite Differences. Figure 23 shows the experimental results in four plots, analytical solution is the solid line and numerical results are depicted with a dashed line. In the first plot we show the direct numerical solution, displacement vs time, and we see how the system gains energy through time, reducing time step alleviates the problem but it does not solve it, since the artificial energy increasing is proportional to the time step. The second plot shows the phase space, which is velocity against displacement, here we observe how the artificial generated energy creates an opening spiral producing greater waves as the simulation moves in time. The third plot reflects how the total energy in the system grows with respect to time. The fourth plot shows the cumulative relative error (5.21) in the interval ${\textstyle t\in [0,7]}$ with respect to ${\textstyle \Delta t}$. From here we noticed that for ${\textstyle \Delta t<0.05}$ this scheme is slightly better than Finite Differences, and for ${\textstyle \Delta t>0.05}$ both schemes are useless in long term simulations, at least that we use a numerical mechanism to rebalance the energy (dampers for instance).

Figure 23: Dashed lines show numerical results, while solid lines are for analytical solution. Top-left plot is the direct solution, the curve of displacement vs time. Top-right plot is the phase space, velocity vs displacement, the opening spiral indicates that artifical energy is being generated. Bottom-left plot is the total energy in the system, the analytical solution is constant and the numerical energy increases as simulation moves on time. Bottom-right plot is the cumulative relative error for distinct values of $\Delta t$.

5.3.1.3 Numerical equilibrium

Using the general scheme in (5.11) and considering constant time steps, ${\textstyle \Delta h=\Delta t}$, we define a numerical approximation for harmonic oscillator in terms of the convolution and its derivative

${\frac {-m}{k}}\left(u-2u^{0}+u^{00}\right)={C}_{P}^{\Delta t}~u+\left(\Delta t{\dot {C}}_{P}^{\Delta t}-2~{C}_{P}^{\Delta t}\right)u^{0}+\left(\Delta t^{2}-\Delta t{\dot {C}}_{P}^{\Delta t}+{C}_{P}^{\Delta t}\right)u^{00},$

(5.28)

then the solution for displacement is

$u=\left({\frac {2u^{0}-u^{00}}{{\frac {k}{m}}{C}_{P}^{\Delta t}+1}}\right)+\left({\frac {2~{C}_{P}^{\Delta t}-\Delta t{\dot {C}}_{P}^{\Delta t}}{{C}_{P}^{\Delta t}+{\frac {m}{k}}}}\right)u^{0}+\left({\frac {\Delta t{\dot {C}}_{P}^{\Delta t}-{C}_{P}^{\Delta t}-\Delta t^{2}}{{C}_{P}^{\Delta t}+{\frac {m}{k}}}}\right)u^{00}$

(5.29)

and the velocity is given by

${\dot {u}}=~\left({\frac {u-u^{0}}{\Delta t}}\right)-\left({\frac {k\Delta t}{2m}}\right)~u^{0}+{\frac {k}{m}}\left({\frac {{C}_{P}^{\Delta t}}{~\Delta t}}-{\dot {C}}_{P}^{\Delta t}\right)~\left(u-u^{0}\right)$

(5.30)

In equation (5.29), observe that by choosing ${\textstyle {P}=2}$ , we get a convolution of ${\textstyle {C}(t)=t^{2}}$ and a convolution derivative of ${\textstyle {\dot {C}}(t)=2t}$, which produces the very same numerical scheme that finite differences in (5.20).

Notice that no matter which shape function we choose, the convolution evaluated at ${\textstyle \Delta t}$ always have the form of ${\textstyle {C}_{P}^{\Delta t}=\beta \Delta t^{2}}$ and its derivative the form of ${\textstyle {\dot {{C}_{P}^{\Delta t}}}=\alpha \Delta t}$, where ${\textstyle \alpha }$ and ${\textstyle \beta }$ are variations of ${\textstyle \Delta t}$ and ${\textstyle \Delta t^{2}}$ respectively. From this fact we will use ${\textstyle \alpha }$ as an optimization variable for minimizing the error, and we set ${\textstyle \beta =1}$ for simplificate the formula (since ${\textstyle \beta }$ is not involved in the minimization). Now we simplify equation (5.29) as

$u=\left({\frac {2u^{0}-u^{00}}{{\frac {k\Delta t^{2}}{m}}+1}}\right)+\left({\frac {2-\alpha }{1+{\frac {m}{k\Delta t^{2}}}}}\right)\left(u^{0}-u^{00}\right)$

(5.31)

With previous equation and using the same ${\textstyle \Delta t=0.1}$ that we use in previous numerical tests, we calculate ${\textstyle \alpha }$ as

$\alpha =arg\,min_{\alpha }\left\{{\hbox{Error}}(7,\alpha )\right\}$

(5.32)

where Error() is the function defined in (5.21). Figure 24 shows the Error as a function of ${\textstyle \alpha }$, in this plot is evident that the optimal value is ${\textstyle \alpha =1}$.

Figure 24: The plot shows the cumulative relative error of equation (5.21) as a function of the optimization variable $\alpha$. It is clear that the minimum is in $\alpha =1$. The error curve is asymptotic to zero in the left and converges to some constant to the right.

Since ${\textstyle \alpha }$ is the proportion of ${\textstyle \Delta t}$ in convolution derivative, ${\textstyle {\dot {{C}_{P}^{\Delta t}}}=\alpha \Delta t}$, from this experiment we found out that

${\hbox{if}}~~~{\dot {C}}_{P}^{\Delta t}<\Delta t~~~{\hbox{ the scheme increases energy}}$ (5.33) ${\hbox{ if}}~~~{\dot {C}}_{P}^{\Delta t}=\Delta t~~~{\hbox{ the energy is stable}}$ (5.34) ${\hbox{ if}}~~~{\dot {C}}_{P}^{\Delta t}>\Delta t~~~{\hbox{ the scheme decreases energy}}$ (5.35)

Examples of energy state of conditions (5.33) and (5.35) can be observed in Figures 22 and 23 respectively.

In our experiments we notice that the variation of ${\textstyle \beta }$ has a little impact on the results, but values of ${\textstyle \beta \leq 1}$ produce smaller oscillations of total energy than ${\textstyle \beta >1}$. For that reason, we constrain our search of shape functions ${\textstyle {P}}$ to those functions that produce ${\textstyle {\dot {C}}_{P}^{\Delta t}=\Delta t}$ and ${\textstyle {C}_{P}^{\Delta t}<\Delta t^{2}}$.

Figure 25 shows experimental results with same ${\textstyle \Delta t=0.1}$ of the numerical scheme resulting from taking ${\textstyle \alpha =1}$ and ${\textstyle \beta =2/5}$ (value of ${\textstyle \beta =2/5}$ is arbitrary chosen only for plotting purposes), which implies that ${\textstyle {C}_{P}^{\Delta t}=\Delta t^{2}(2/5)}$ and ${\textstyle {\dot {C}}_{P}^{\Delta t}=\Delta t}$. The results are displayed in the same format that previous experiments of harmonic oscillator, the analytical solution is the solid line, the numerical solution is the dashed line, and we have 4 plots to show the curves of displacement vs time, the phase space (velocity vs displacement), the total energy and the error when moving ${\textstyle \Delta t}$. In this plots we can appreciate the stability of the system, which have little oscillations of the total energy.

Figure 25: Experimental results of the numerical scheme resulting of setting ${\dot {C}}_{P}^{\Delta t}=\Delta t$. Dashed lines show numerical results and solid lines are for analytical solution. Top-left plot is the direct solution, the curve of displacement vs time. Top-right plot is the phase space, velocity vs displacement. Bottom-left plot is the total energy in the system, the analytical solution is constant and the numerical energy has little oscillations around such constant. Bottom-right plot is the cumulative relative error for distinct values of $\Delta t$.

5.3.2 Trigonometric shape function

In order to build our numerical time discretization, we propose this trigonometric shape function shown in Figure 21

${P}(t)={\frac {1}{2}}\left(1-\cos(\pi ~t)\right)$

(5.36)

The appendix B discuss another proposals based on polynomials that produces accurate results, nevertheless this trigonometric function introduces less artificial energy that impact results in long term simulations.

The convolution corresponding to (5.36) is

${C}_{P}(t)={\frac {\Delta t^{2}}{2\pi ^{2}}}\left(\cos \left({\frac {\pi ~t}{\Delta t}}\right)-1\right)+{\frac {1}{2}}~t^{2},$

(5.37)

and its derivative is

${\dot {C}}_{P}(t)=t-{\frac {\Delta t}{2\pi }}\sin \left({\frac {\pi ~t}{\Delta t}}\right),$

(5.38)

although when we evaluate at ${\textstyle {C}_{P}(\Delta t)}$ and ${\textstyle {\dot {C}}_{P}(\Delta t)}$ we obtain

${C}_{P}^{\Delta t}=\left({\frac {\pi ^{2}-2}{2~\pi ^{2}}}\right)\Delta t^{2}~~~~~{\hbox{and}}~~~~~~{\dot {C}}_{P}^{\Delta t}=\Delta t$

(5.39)

as expected by previous discussion in this section. Notice that plots in Figure 25 were produced with ${\textstyle \beta =2/5}$, and for this shape function we have

$\beta ={\frac {\pi ^{2}-2}{2\pi ^{2}}}\approx {\frac {2}{5}}$

(5.40)

that produces almost identical plots.

5.3.3 Dynamic stress analysis

Considering discretized stress equation (3.53) into this numerical scheme, we have that

$S_{i}=\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {K} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk},$

(5.41)

replacing (5.36) into such stress equation we get our final numerical system for the first equation of motion (2.36.a)

${\boldsymbol {u}}_{i}-\left({\frac {\pi ^{2}}{2}}-1\right){\frac {\Delta t^{2}}{\rho ~\pi ^{2}}}\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {K} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk}=$ $\left(1+{\frac {\Delta t}{\Delta h}}\right){\boldsymbol {u}}_{i}^{0}~-~{\frac {\Delta t}{\Delta h}}~{\boldsymbol {u}}_{i}^{00}~+$ ${\frac {\Delta t\Delta h}{\rho ~\pi ^{2}}}\left(\left({\frac {\pi ^{2}-2}{2}}\right)~S_{i}^{00}+\left({\frac {\Delta t}{\Delta h}}+1\right)S_{i}^{0}\right),$

(5.42)

If we choose a constant time step, ${\textstyle \Delta h=\Delta t}$, we can simplify the equation to

${\boldsymbol {u}}_{i}-\left(\pi ^{2}-2\right){\frac {\Delta t^{2}}{2\rho ~\pi ^{2}}}\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {K} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk}=$ $2~{\boldsymbol {u}}_{i}^{0}-{\boldsymbol {u}}_{i}^{00}+{\frac {\Delta t^{2}}{2\rho ~\pi ^{2}}}\left(\left(\pi ^{2}-2\right)~S_{i}^{00}+4~S_{i}^{0}\right),$

(5.43)

The results shown in this work use the trigonometric time variation with fixed length time steps, although our software supports automatic (variable) time steps.

5.3.4 Time step calculation

Numerically speaking, the time step ${\textstyle \Delta t}$ is bounded by the maximum propagation speed of stress waves, denoted ${\textstyle \mathbf {v} _{\sigma }\in \mathbb {R} }$, which means that ${\textstyle \Delta t}$ should be small enough to capture the stress state. Using the solution of the one dimensional wave equation, ${\textstyle \mathbf {v} _{\sigma }}$ can be estimated as

$\mathbf {v} _{\sigma }^{2}={\frac {E}{\rho }},$

(5.44)

where ${\textstyle E}$ and ${\textstyle \rho }$ are material properties. This equation implies that in order to reproduce stress waves numerically, the following relation should be satisfied

$\Delta t~=~C_{o}~\Delta \mathbf {x} ~{\sqrt {\frac {\rho }{E}}}$

(5.45)

where ${\textstyle C_{o}}$ is the Courant number and ${\textstyle \Delta \mathbf {x} }$ is the average area of control volumes. In this work we use Courant numbers ${\textstyle C_{o}\in (0.05,0.5)}$ proposed by [17] for stress analysis. At the beginning of simulation we use ${\textstyle C_{o}\approx 0.3}$, but when damage field starts producing failures (interfering with stress analysis), we use a smaller Courant number ${\textstyle C_{o}\approx 0.05}$ to produce more accurate results and better conditioned sytems of equations.

6 Coupled system

In this chapter we will discuss in detail how the discretized versions of both equations of motion, (5.42) and (4.22), are coupled in a segregated approach. We will take first equation of motion as a cornerstone of the numerical scheme, since displacements and internal forces are first class fields in the physical system, while the damage field is an auxiliary abstraction to approximate fracture surface.

6.1 Residual minimization

Due to the fact that a single time step can release enough energy to create wide crack surfaces, the numerical structure can become unstable. In order to make small changes in stress field to produce numerically manageable systems we dose internal forces in equation (5.42) into ${\textstyle N_{\boldsymbol {k}}}$ finite increments,

${\boldsymbol {f_{k}}}={\frac {1}{N_{\boldsymbol {k}}}}\left(\left(1+{\frac {\Delta t}{\Delta h}}\right){\boldsymbol {u}}_{i}^{0}-{\frac {\Delta t}{\Delta h}}~{\boldsymbol {u}}_{i}^{00}~+k_{\pi }{\frac {\Delta h}{\Delta t}}~S_{i}^{00}+{\frac {\Delta t^{2}}{\rho ~\pi ^{2}}}\left({\frac {\Delta h}{\Delta t}}+1\right)S_{i}^{0}\right)$

(6.1)

where coefficient ${\textstyle k_{\pi }}$ is given by

$k_{\pi }=\left({\frac {\pi ^{2}}{2}}-1\right){\frac {\Delta t^{2}}{\rho ~\pi ^{2}}},$

(6.2)

thus displacement and damage fields at time ${\textstyle t+\Delta t}$ are calculated by computing ${\textstyle N_{\boldsymbol {k}}}$ finite increments. Each finite increment ${\textstyle ~{\boldsymbol {k}}\in [1,N_{\boldsymbol {k}}]~}$ is solved by minimizing residual

$R_{i}^{\boldsymbol {k}}={\boldsymbol {u}}_{i}^{\boldsymbol {k}}-k_{\pi }\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {K} _{ijk}|_{\mathbf {x} _{g}}~{\vec {\boldsymbol {u}}}_{ijk}^{\boldsymbol {k}}-{\boldsymbol {k}}~{\boldsymbol {f_{k}}},$

(6.3)

The displacement corresponding to finite increment ${\textstyle {\boldsymbol {k}}}$, referred as ${\textstyle {\boldsymbol {u}}_{i}^{\boldsymbol {k}}}$, is the accumulation of ${\textstyle n}$ displacement increments in residual minimization, denoted ${\textstyle \Delta {\boldsymbol {u}}^{[n]}}$. Residual is minimized using Newton-Raphson iteration [41], starting with

$R_{i}^{{\boldsymbol {=}}0}=-{\boldsymbol {k~f_{k}}}$

(6.4)

Such numerical scheme is composed by the linear systems in equations (5.42) and (4.22),

$[left=n^{th}{\hbox{ iteration}}\{\,]\Delta {\boldsymbol {u}}_{i}^{[n]}-k_{\pi }\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\sum \limits _{g=1}^{N_{g}}~~w_{g}~\mathbf {K} _{ijk}^{[n]}|_{\mathbf {x} _{g}}~\Delta {\vec {\boldsymbol {u}}}_{ijk}^{[n]}=R_{i}^{{\boldsymbol {=}}n-1},$ (6.5) ${\boldsymbol {u}}_{i}^{\boldsymbol {k}}=\sum _{r=1}^{n}\Delta {\boldsymbol {u}}_{i}^{[r]}$ (6.6) $V_{i}~{\boldsymbol {d}}_{i}^{[n]}~~+\sum \limits _{j=1}^{N_{i}}\sum \limits _{k=1}^{N_{ij}}\left({\vec {\mathbf {D} }}_{ijk}^{[n]}\right)^{T}~{\vec {\boldsymbol {d}}}_{ijk}^{[n]}=~~W_{i}^{[n]}$ (6.7)

where superindex ${\textstyle ()^{[n]}}$ indicates the iteration within the residual minimization for solving finite increment ${\textstyle {\boldsymbol {k}}}$. Finally, minimization finish when residual converges. We have two criteria for detecting convergence, the first one consist in checking if residual norm is small enough

$||R^{\boldsymbol {k}}||<\epsilon ,~~~~{\hbox{with}}~~~~\epsilon <<1,$

(6.8)

and the second criterion consist in checking that the derivative of the norm of the residual with respect to the iterations is approximately zero

${\frac {d||R^{\boldsymbol {k}}||}{dn}}\approx 0~~~~~~{\hbox{if convergence is reached}},$

(6.9)

which means that solution is not changing any more in concecutive iterations. In order to estimate this derivative we calculate a linear regression of ${\textstyle ||R^{\boldsymbol {k}}||}$ (with respect to ${\textstyle n}$) of the fifty previous iterations, and we take the slope coefficient as the derivative.

6.2 Discrete Fracture

The fracture between the ${\textstyle i^{th}}$ control volume and the volume adjacent to its ${\textstyle j^{th}}$ face is produced when such face satisfies

$\int _{e_{ij}}{\boldsymbol {d}}(\mathbf {x} )dS=\int _{e_{ij}}dS,\mathbf {x} \in e_{ij},$

(6.10)

The condition is reached when the damage field is enable along all the face, this means that we consider a fracture when the face is completely damaged. Thus the control volumes could be separated from their adjacent control volumes due to this fracture mechanism, and for that reason, the title of the monograph is the Discrete Volume Method.

The discrete volumes are integrated with the CVFA method, considering discrete faces as continuum faces completely damaged, ${\textstyle {\boldsymbol {d}}=1}$, when calculating the damaged strain equation (3.17). All the faces which not appear in the initial discretization are treated as discrete faces, this allows the collision of separated bodies and the self-collision of the body boundaries. New discrete faces arise from the fracture process.

7 Results

In this chapter we will discuss the following numerical experiments: plate with a hole to compare elasticity solver against analytical solution, stress wave in a bar to compare dynamic solver against published results, perfored strip under tension for fracture due to direct tension, three point bending bar test to compare our solver with lab experiments, brazilian test to analyze fracture due to indirect tensile strain, compressive test to verify cracking patterns against those produced with similar numerical methods, four point notched bar test for analyzing fracture mode II, Three point bending bar with asymmetric perforations test to evaluate sensibility of crack morphology, notched plate under shear to contrast our results with those of other authors, dynamic shear loading, and dynamic crack branching to demonstrate multicrack generation.

For each test we review the material properties, the geometrical description of the body, the configuration of boundary conditions, the numerical parameters and other considerations. We also provide figures to portray the specification of each experiment.

We use LU decomposition to solve both systems of equations (elasticity and damage), the displacement field requires memory to store a vector per discrete volume, whereas that damage field only requres a floating point number (double precision) per node. The discrete volumes are numbered using a nested dissection technique [54] before assembling the systems in order to reduce the fill-in in LU decomposition. In our experiments, initial meshes were generated with a central Voronoi-cells procedure.

7.1 Plate with a hole

In order to test the numerical performance of the proposed method for solving the first equation of motion, we use the well known analytical experiment of an infinite plate with a circular hole in the origin (see [55]). In such a experiment, the plate is stretched along the horizontal axis with a uniform tension of ${\textstyle \mathbf {f} _{[1]}=10~~kPa}$ from each side, as is shown in Figure 26. The material is characterized by the Poisson's ratio, ${\textstyle \nu =0.3}$, and Young modulus, ${\textstyle E=10~~MPa}$. Plane stress is assumed with thickness equivalent to the unity. The dimensions of the computational domain are ${\textstyle a=0.5m}$ and ${\textstyle b=2m}$.

Figure 26: (a) Infinite plate with a hole being stretched along the horizontal axis with a force of $\mathbf {f} _{[1]}=10~~kPa$ from each side. (b) Computational domain, $a=0.5m$ and $b=2m$, with axysymmetrical assumptions used to test the numerical method. The polar coordinates, $r$ and $\theta$, for calculating the analytical stress field.

The analytical solution is given by the following formulae

${\boldsymbol {\sigma }}_{[11]}=\mathbf {f} _{[1]}\left[1-{\frac {a^{2}}{r^{2}}}\left({\frac {3}{2}}\cos(2\theta )+\cos(4\theta )\right)+{\frac {3a^{4}}{2r^{4}}}\cos(4\theta )\right],$ (7.1) ${\boldsymbol {\sigma }}_{[22]}=\mathbf {f} _{[1]}\left[-{\frac {a^{2}}{r^{2}}}\left({\frac {1}{2}}\cos(2\theta )-\cos(4\theta )\right)-{\frac {3a^{4}}{2r^{4}}}\cos(4\theta )\right],$ (7.2) ${\boldsymbol {\sigma }}_{[12]}=\mathbf {f} _{[1]}\left[-{\frac {a^{2}}{r^{2}}}\left({\frac {1}{2}}\sin(2\theta )+\sin(4\theta )\right)+{\frac {3a^{4}}{2r^{4}}}\sin(4\theta )\right],$ (7.3)

where the polar coordinates,

$r={\sqrt {\mathbf {x} _{[1]}^{2}+\mathbf {x} _{[2]}^{2}}},~~~{\hbox{and}}~~~\theta =\tan ^{-1}\left({\frac {\mathbf {x} _{[2]}}{\mathbf {x} _{[1]}}}\right),$

(7.4)

are used within the calculus. Figure 27 exhibits (a) the discretization of the computational domain into 2411 polygonal volumes used to compare the numerical results (using smoothness ${\textstyle c=0}$) against the analytical stress field. This mesh is not equivalent to the Voronoi diagram. (b) Level sets of ${\textstyle {\boldsymbol {\sigma }}_{[11]}}$ between ${\textstyle 0}$ to ${\textstyle 30~~kPa}$, with steps of ${\textstyle 1~~kPa}$. (c) Level sets of ${\textstyle {\boldsymbol {\sigma }}_{[22]}}$ between ${\textstyle -10}$ and ${\textstyle 6~~kPa}$, with steps of ${\textstyle 0.8~~kPa}$. (d) Level sets of ${\textstyle {\boldsymbol {\sigma }}_{[12]}}$ between ${\textstyle -10}$ and ${\textstyle 2~~kPa}$, with steps of ${\textstyle 0.6~~kPa}$.

Figure 27: (a) Polygonal mesh used for comparison of numerical results. (b) Level sets of ${\boldsymbol {\sigma }}_{[11]}$ between $0$ to $30~~kPa$. (c) Level sets of ${\boldsymbol {\sigma }}_{[22]}$ between $-10$ and $6~~kPa$. (d) Level sets of ${\boldsymbol {\sigma }}_{[12]}$ between $-10$ and $2~~kPa$.

The Dirichlet conditions are imposed on the bottom and left side of the computational domain as is shown in Figure 27.b. Next in order, the analytic stress of equations (7.1), (7.2) and (7.3) is imposed as Neumann condition over the top and right side of the computational domain.

Figure 28.a presents the averaged error as a function of mesh size, denoted ${\textstyle \Delta \mathbf {x} }$, as we might expect, the error is proportional to the mesh refinement. For a mesh of 628 volumes, Figure 28.b shows the percentage of the error with respect to the error of ${\textstyle c=0}$, for different smoothing levels, ${\textstyle c=0}$ correspond the linear interpolator. Observe that the error of the stress field does not decreases significantly in the first three levels of smoothness, this is because we do not increase the degrees of freedom of the linear system (is the same mesh), although we built a better field description, which can be useful when solving non-linear formulations. The increasing error after ${\textstyle c=2}$ is produced by floating point truncation, since ${\textstyle c>2}$ implies computing integrals for polynomials of ${\textstyle 7^{th}}$ order or higher.

Figure 28: (a) The averaged error decreasing as a function of mesh size, denoted $\Delta \mathbf {x}$. (b) Using a mesh of 628 volumes, percentage of error for different smoothing levels with respect to the error of $c=0$, which is the linear interpolator, error increases after $c=2$.

7.2 Stress wave in a bar

This experiment consists in analyzing the pattern of a stress wave in a long bar to assess performance of time discretization developed in this work.

We assume plane stress with a thickness of ${\textstyle 1~m}$, and we choose a smoothness ${\textstyle c=1}$. The material properties are those of steel, elasticity modulus ${\textstyle E=200~GPa}$, Poisson's ratio ${\textstyle \nu =0.3}$ and density ${\textstyle \rho =7854~kg/m^{3}}$. The size of the bar is ${\textstyle 10~m}$ long and ${\textstyle 1~m}$ wide. The initial displacement and velocity at the interior of the bar are null. Figure 29 illustrates the geometrical distribution and numerical parameters.

Figure 29: Large bar used for propagating a stress wave, with a size of $10~m$ long and $1~m$ wide. We assume Plane stress with a thickness of 1, and properties of steel are used. The bar is fixed to right side and is pushed $\delta =1~mm$ from left side at the beginning.

The bar is fixed from right side and is pushed ${\textstyle \delta =1~mm}$ from left side at time ${\textstyle t=0}$, this produces a propagation of the stress wave at speed of sound through the bar, that for steel is

${\hbox{Sound speed in steel }}={\sqrt {\frac {E}{\rho }}}=5.0462~km/s$

(7.5)

when the stress wave arrives to the fixed side, it is reflected in reverse direction, at this moment perpendicular waves resulting from Poisson's effect start interfering with our frontal wave. Figure 30 shows the bubble formed by stress wave after time ${\textstyle t=1.92~ms}$, when frontal wave is reflected, as predicted by [17]. This bubble is the contour of ${\textstyle C_{\boldsymbol {\sigma }}=0.95}$, where

$C_{\boldsymbol {\sigma }}={\frac {1-\exp \left(-{\boldsymbol {\sigma }}_{[11]}\right)}{1+\exp \left(-{\boldsymbol {\sigma }}_{[11]}\right)}}$

(7.6)

is an auxiliary field used for rescaling the horizontal stress component to filter smaller and negative stress waves.

Figure 30: Stress wave produced by initial imposed displacement. The top shows partition $P_{h}$ of bar into discrete volumes. Then a sequence of images from top to bottom illustrates the moment when the wave is being reflected, the contour of the bubble is produced by $C_{\boldsymbol {\sigma }}=0.95$, which is an auxiliary field to rescale horizontal component of stress tensor, this rescaling is performed to filter negative and small waves.

7.3 Perfored strip under tension

This test involves a perfored strip under tension, the stress field induces a fracture around the perforation because its tensions are greater there than in the rest of the domain. The goal of this experiment is producing a pure mode I failure. The strip is a perfect square of ${\textstyle 40~cm}$ each side and perforation has a radius of ${\textstyle 1~cm}$. We assume plane stress with thikness = ${\textstyle 1~m}$, and perform a quasi-static analysis with 100 finite increments, using a smoothness ${\textstyle c=1}$. Taking advantage of the symmetry we analyze only the right half of the body by imposing symmetry conditions. In the experiment, the strip is pulled apart vertically with an equivalent displacement from top and bottom ${\textstyle \delta =0.1~mm}$. The material properties are Young modulus ${\textstyle E=30~GPa}$, Poisson's ratio ${\textstyle \nu =0.2}$ and energy release rate ${\textstyle G=100~J/m^{2}}$. In order to compare our results to those published by [26], we discretize the strip in two distinct meshes with average size of ${\textstyle \Delta \mathbf {x} =2.5~mm}$ (12277 volumes) and ${\textstyle \Delta \mathbf {x} =5~mm}$ (3281 volumes). Figure 31 depicts a) the geometrical specifications and material properties, and b) the portion used for numerical analysis, where symmetry conditions were imposed.

Figure 31: a) Geometry of the strip, a perfect square of $40~cm\times 40~cm$ with a hole in the middle with radius $1~cm$. Assumption of plane strain is considered. b) Due to symmetry only the right half analyzed numerically, the figure shows the corresponding boundary conditions.

The numerical experiment is performed in meshes exposed in Figure 32. Fracture location coincides in both meshes, such fracture corresponds to an horizontal line in the middle of the strip. Curves shown in Figure 32 are similar to those published by [26]. The area under curve obtained with mesh size ${\textstyle \Delta \mathbf {x} =5~mm}$ is equal to ${\textstyle 8.93~J}$, whereas that theoretical energy released by fracture is

${\hbox{Energy due to crack}}={\hbox{ Crack length }}\times G\times {\hbox{ thickness}}$ $=9~cm\times 100~J/m^{2}\times 1~m=9~J$

(7.7)

Area under curve obtained with mesh size ${\textstyle \Delta \mathbf {x} =2.5~mm}$ is almost theoretical value ${\textstyle 9~J}$.

Figure 32: a) Vertical reaction vs Total vertical displacement, area under the curve is close to theoretical $9~J$ required to generate the crack. Solid line indicates results of mesh size $\Delta \mathbf {x} =5~mm$ and dashed line indicates results of mesh $\Delta \mathbf {x} =2.5~mm$. b) Right side depicts mesh size $\Delta \mathbf {x} =5~mm$, and lef side shows a reflected version of mesh size $\Delta \mathbf {x} =2.5~mm$.

Figure 33 depicts damage field for three distincts displacements in both meshes, with a fracture arising next to the perforation.

Figure 33: First column illustrates damage field and discretization for mesh $\Delta \mathbf {x} =5$ mm with x50 deformation factor, and second column shows same damage field but in discretization of mesh $\Delta \mathbf {x} =2.5$ mm

7.4 Three point bending bar

In this example we have a bar with size of ${\textstyle 45~cm}$ long and ${\textstyle 10~cm}$ wide, a vertical notch from center to bottom. The material has elasticity modulus ${\textstyle E=20~GPa}$, Poisson's modulus ${\textstyle \nu =0.2}$ and energy release rate ${\textstyle G=113~J/m^{2}}$. The bar is vertically displaced, ${\textstyle \delta =1mm}$, from top to bottom. Figure 34 depicts geometrical specification, material properties, numerical parameters and boundary conditions.

Figure 34: Three point bending bar. Geometrical specification, material properties, numerical parameters and boundary conditions.

We assume plane strain in a quasi-static analysis using 100 finite increments and smoothness ${\textstyle c=1}$. The goal of this test is comparing experimental results published by [56] with our numerical approximation. We use a domain partition of 5297 discrete volumes with an average size ${\textstyle \Delta \mathbf {x} =2.9~mm}$. Figure 35 shows such partition in the top, and below it depicts b) the damage field calculated and c) the displacement with a deformation factor of x35.

Figure 35: Three point bending bar. a) Partition used in numerical analysis with average size $\Delta \mathbf {x} =2.9~mm$ (5297 discrete volumes). b) Damage field and c) Deformation scaled with a factor of x35.

Figure 36 shows curve of reaction (load) agains displacement, gray area corresponds to experimental results obtained by [56] and dashed line with black dots is related to our numerical results.

Figure 36: Three point bending bar. Reaction (load) vs displacement, gray area corresponds to experimental results, whereas that dashed line and black dots are related to numerical analysis.

7.5 Brazilian test

The brazilian tensile strength (BTS) test was designed to assess strength of brittle materials [57]. This experiment consists in compressing a disk to generate, by Poisson's effect, indirect tensions to produce a vertical fracture. The disk has a radius of ${\textstyle 10~cm}$ and is fixed to the bottom from a plain side (circular chord) of ${\textstyle 2~cm}$ and is pushed from top to bottom using another plain side in the top of ${\textstyle 2~cm}$. The material has elasticity modulus ${\textstyle E=21~GPa}$, Poisson's ratio ${\textstyle \nu =0.2}$, and energy release rate ${\textstyle G=1~mJ/m^{2}}$. We assume plane stress with thickness of ${\textstyle 10~cm}$ in a quasi-static analysis using 100 finite increments and smoothness ${\textstyle c=1}$. Within this experiment (see [57]), vertical stress is given by

${\boldsymbol {\sigma }}_{[22]}={\frac {2{\boldsymbol {p}}}{\pi \times {\hbox{diameter}}\times {\hbox{thickness}}}}={\frac {\boldsymbol {p}}{100\pi }}$

(7.8)

where ${\textstyle {\boldsymbol {p}}}$ is the applied load. Figure 37 shows geometrical specification, material properties and boundary conditions.

Figure 37: Brazilian test. Geometry is described by a disk with radius and thickness of $10~cm$, assuming plane stress. Material properties and boundary conditions for quasi-static analysis are displayed.

We analyze three distinct discretization, first one has an average size ${\textstyle \Delta \mathbf {x} =4~mm}$ (1926 discrete volumes), the second one has ${\textstyle \Delta \mathbf {x} =2.8~mm}$ (3896 discrete volumes), and the third has ${\textstyle \Delta \mathbf {x} =2~mm}$ (7553 discrete volumes). Figure 38 illustrates damage field obtained for three discretizations and meshes are deformed with a factor of x5000 (last one x1000). The result of finest mesh produces the theoretical vertical fracture. All three numerical calculations fail close to predicted by formula (7.8), which is at ${\textstyle 314.16N}$.

Figure 38

7.6 Four point notched bar

This test is intended to produce a mode II failure and it is selected because its geometry has two initial crack tips at the end of the notches, where stress field has its highest values. The experiment consists in bar with size of ${\textstyle 134~cm}$ long and ${\textstyle 30.6~cm}$ wide with two vertical notches in the middle, onte from top and one from bottom. The bar is fixed to two boxes of ${\textstyle 4~cm\times 4~cm}$ that work as main support, the first one is close to the bottom-right corner, and the second is close to the top-left corner. The bar is pushed from top to bottom by displacing ${\textstyle \delta =0.1~mm}$ a small box next to the right of the centered upper notch, whereas that a second displaced box (same ${\textstyle \delta }$) is pushing from bottom to top next to the left of the centered lower-notch. The material properties are elasticity modulus ${\textstyle E=30GPa}$, Poisson's ratio ${\textstyle \nu =0.2}$, and energy release rate ${\textstyle G=100J/m^{2}}$. Figure 39 depicts geometrical specification, position of notches and supports, material properties, and boundary conditions.

Figure 39: Four point bending bar. Geometry specification, material properties, numerical parameters and boundary conditions.

Plane strain is assumed in a quasi-static analysis using 100 finite increments and smoothness ${\textstyle c=1}$. We perform two separate analysis in distinct meshes, the first one has an average size ${\textstyle \Delta \mathbf {x} =12~mm}$ (2849 discrete volumes) and the second one has ${\textstyle \Delta \mathbf {x} =7.5~mm}$ (7255 discrete volumes). Figure 40 illustrates damage field obtained with our numerical approach, and its corresponding displacement deformated with a factor x100. The differences between both discretizations can be appreciated in this figure. In contrast with most solutions obtained with damage models based on Finite Element Method, we get asymmetrical crack morphology induced by discrete volumes shape.

Figure 40: At top is shown discretization, damage field and displacement computed with mesh size $\Delta \mathbf {x} =12~mm$ (2849 discrete volumes) with a deformation factor x100. At bottom we can appreciate discretization, damage field and displacement using a mesh size $\Delta \mathbf {x} =7.5~mm$ (7255 discrete volumes).

In Figure 41 we can observe vertical reaction against vertical displacement for both numerical experiments, and the curve obtained by Cervera et al [28] using a FEM damage model with triangular elements and average size ${\textstyle \Delta \mathbf {x} =5~mm}$ (5909 FEM nodes). This graph exposes how coarse discretizations increase brittleness in material, which is an expected behaviour in formulations where continuum is dislocated toproduce new crack surfaces. In this experiment we observe that discrete volumes shapes interfere with damage field computation, which is a numerical artifact since our mathematical formulation assumes a continuum, however this effect occurs in non-continuum fractures and in continuum but not homogenous materials, which are more likely to fail in regions with lower density for example.

Figure 41: Four point bending bar test. Reaction vs displacement curves, solid line corresponds to published results, dashed line shows results for mesh size $\Delta \mathbf {x} =8~mm$, and dotted line depicts results for mesh size $\Delta \mathbf {x} =12~mm$.

7.7 Three point bending bar with asymmetric perforations

The three point bending bar with asymmetric perforations is the classical example of fracture mechanics in brittle materials, this experiment was proposed by [58]. This test consists in a bar with size of ${\textstyle 20~in}$ long and ${\textstyle 8~in}$ wide, it has three perforations in left half with ${\textstyle 0.5~in}$ radius, these perforations are horizontally aligned and have a vertical separation of ${\textstyle 2~in}$. The bar has a vertical notch ${\textstyle 1~in}$ long in the bottom-left quadrant. Such a bar is fixed from one point close to the bottom-left corner, vertical displacement is forbidden in a point close to the bottom-right corner, and a vertical displacement ${\textstyle \delta }$. Plane stress is assumed, thickness = 1, in a quasi-static analysis using 100 finite increments and smoothness ${\textstyle c=1}$. The material properties are those of PMMA, Polymethyl-methacrylate (also known as acrylic glass or plexiglas). We choose the average Poisson's ratio ${\textstyle \nu =0.375}$ and according to [58] Young modulus is ${\textstyle E=474\times 10^{3}}$ psi, which corresponds to ${\textstyle E=3.27~GPa}$. For numerical experiments, material properties must be transformed from units based on meters to units based on inches, or when defining the geometry and boundary conditions inches must be transformed to meters. Figure 42 depicts geometrical specifications, material properties and boundary conditions.

Figure 42: Three point bending bar with asymmetric perforations. Geometrical specification, material properties and boundary conditions.

Figure 43 illustrates damage field obtained.

Figure 43: Three point bending bar with asymmetric perforations. Damage field obtained.

Figure 44 shows experimental results obtained by Bittencourt et al [58]. We can observe that crack trajectory is similar to that obtained numerically in this calculation.

Figure 44: Experimental results obtained by Bittencourt et al [58].

7.8 Notched plate under shear

This test is intended to demonstrate rupture by shear displacement, it consists in plate with an horizontal notch from left side to the center. The plate is fixed from bottom and is displaced ${\textstyle \delta =1.5\times 10^{-2}~mm}$ from top to the right. The material properties are Young modulus ${\textstyle E=210~GPa}$, Poisson's ratio ${\textstyle \nu =0.3}$ and energy release rate ${\textstyle G=2.7kJ/m^{2}}$. Numerical analysis is performed assuming plane strain and using 100 finite increments and smoothness ${\textstyle c=1}$. Figure 45 illustrates a) the geometrical specification, material properties and boundary conditions, whereas that b) shows the damage field obtained with displacements scaled 1000.

Figure 45: Notched plate under shear, a) the geometrical specification, material properties and boundary conditions, and b) damage field obtained with displacements scaled 1000.

7.9 Dynamic shear loading

This test consists in a bar with size of ${\textstyle 20~cm}$ height and ${\textstyle 10~cm}$ width, it has to horizontal notches symmetrical to the horizontal line that splits geometry in two halves, these notches goes from left side to central vertical line. Material properties are Young's modulus ${\textstyle E=190~GPa}$, Poisson's ratio ${\textstyle 0.3}$, density ${\textstyle \rho =8000~kg/m^{3}}$ and energy release rate ${\textstyle G=22.13~kJ/m^{3}}$. We simulate a projectile impacting the left side in between the notches at a velocity of

$\mathbf {v} _{c}=\left\{{\begin{matrix}{\frac {t}{t_{0}}}\mathbf {v} _{0}&t<t_{0}\\\mathbf {v} _{0}&t\geq t_{0}\end{matrix}}\right.$

(7.9)

where ${\textstyle \mathbf {v} _{0}=16.5~m/s}$ and ${\textstyle t_{0}=~\mu s}$. Since the geometry and boundary conditions are vertically symmetric, we take the superior half to perform our numerical analysis, by assuming plane strain and using 100 finite increments and smoothness ${\textstyle c=1}$. Figure 46 illustrates the geometrical specifications, material properties, initial conditions and other considerations for numerical analysis.

Figure 46: a) Dynamic shear loading test simulates a bar being impacted by a projectile from left side in between two notches. b) Assuming symmetrical conditions we analyze the superior half of the geometry.

Figure 47 depicts damage field obtained.

Figure 47: Damage field obtained in dynamic shear loading test.

7.10 Dynamic crack branching

This experiment was proposed by [33] to generate a dynamic crack branching. It consists of a bar with size of ${\textstyle 10~cm}$ long and ${\textstyle 4~cm}$ wide. A horizontal notch is inserted from left side to center to produce an initial crack. Material properties are Young's modulus ${\textstyle E=32~GPa}$, Poisson's ratio ${\textstyle 0.2}$, density ${\textstyle \rho =2450~kg/m^{3}}$ and energy release rate ${\textstyle G=3~J/m^{3}}$. The bar is being pulled apart from top and bottom with a pressure of ${\textstyle {\boldsymbol {\sigma }}=1~MPa}$. The right side can not be displaced horizontally but it can over the vertical, the middle point of this side is completely fixed. For numerical analysis we assume plane strain using 100 finite increments and smoothness ${\textstyle c=1}$. Although the geometry and boundary conditions are vertically symmetric, we perform the analysis over the whole domain. Figure 48 shows numerical and geometrical specification, material properties and boundary conditions.

Figure 48: Geometrical specification, numerical parameters, material properties and boundary conditions for dynamic crack branching experiment.

Figure 49 depicts damage field obtained.

Figure 49: Damage field obtained int dynamic crack branching experiment. Displacements are scaled up a factor of 1000

8 Conclusions

In this work we proposed a numerical technique for simulating the mechanics of brittle fracture by using an alternative definition of the elastic potential energy that involves a damage field to decrease energy due to tensile strain over the fractured surface. Such damage field is a smooth approximation of the crack morphology, with a value of one to describe fractured surfaces and zero for the rest of the elastic body. In the mathematical formulation the total potential energy is determined by the contribution of the elastic potential energy and the potential energy of the body to nucleate new cracks. The elastic potential energy is charaterized by the sum of the elastic energy due to compression plus the elastic energy due to tension, the second term is scaled by a quadratic expression of the damage field, nullifying it when damage is equal to one. The potential energy to generate new cracks is related with the length of the existing cracks and the energy release rate, a material property. A bigger crack increases the potential energy to propagate it. The equations of motion are obtained from the solution of the variational problem for minimizing the Lagrangian of our system, that is, finding the optimal displacement and damage fields for reducing the difference between the potential and the kinetic energy of the body.

The solution of the system is calculated by applying finite increments, with an inner loop within each time step until reaching equilibrium of elasticity equation. That is solving elasticity equation and using the computed displacement field to solve the damage equation, in the next iteration we use damage field estimation to solve again elasticity equation, repeating the process until the residual norm is zero in first equation of motion.

In order to solve the partial differential equations we employ a numerical technique to discretize the body using unstructured and non conforming meshes formed by elements of any arbitrary polygonal/polyhedral shape. The elastic solver is based on a finite volume formulation that, using the divergence theorem, represent the volume integral of the stress divergence in terms of the surface integral of the stress over the volume boundary. Since the stress term is calculated directly on the boundary of the control volumes, this strategy can be used in our fracture formulation where volumes are treated as indivisible components and the rupture occurs across the volumes boundaries. The damage solver follows a similar approach, but considering volume integrals of damage field apart of the surface integral resulting of applying divergence theorem to damage diffusive term. Control volume boundary is divided into flat faces for considering the normal unit vector as a constant. Conforming and non-conforming meshes are processed without distinction. Both fields, displacement and damage, are a piece-wise polynomial approximation surrounding the volumes, built on the top of the simplices resulting from the Delaunay triangulation of the volume neighborhood. A pair-wise polynomial interpolation is used for neighborhoods where the simplices are exceedinlgy distorted or it can not be formed.

On the other hand, time discretization is based on the analytical solution, obtained by means of Laplace transform, of the ordinary differential equation resulting from assuming a continuous variation in time of the stress state.

In spatial discretization, we introduced the homeostatic splines and its pseudo-inverses for higher order polynomial interpolations without the necessity of increasing the discretization points, but adding a computational cost for numerical integration. The rate of increasing computational cost is greater than the rate of decreasing numerical error when choosing high degree homeostatic splines. This situation makes computationally expensive smoothness of solution at nodes.

In time discretization, we propose a trigonometric shape function to describe time variation of stress state, which produces an energy-stable numerical scheme and tolerates bigger time steps than methods based on simply finite differences. Leaving aside the stability of the method, choosing big time steps will increase the number of iterations of the finite increments strategy performed in every time step. In the results presented here we use a Courant number of 0.05 as a reasonable trade-off.

Finally we present numerical experiments for the well known plate with a hole to compare our elasticity solver against analytical solution, stress wave in a bar to compare our dynamic solver against published results, perfored strip under tension for checking fracture due to direct tension, three point bending bar test to compare our solver with lab experiments of fracture due to tensile strain, brazilian test to analyze fracture due to indirect tensile strain comparing our results against published by other authors, a compressive test to verify cracking patterns against those produced with similar numerical methods, four point notched bar test for analyzing fracture mode II, three point bending bar with asymmetric perforations to evaluate sensibility of crack morphology, notched plate under shear to contrast our results with those of other authors, dynamic shear loading, and dynamic crack branching to demonstrate multicrack generation.

In future work, we would like to analize numerical results of 3D tests, to explore adaptable meshes in order refine elements if damage is likely to occur, and to investigate the response of the system to shock wave impacts, such as those produced by detonations. Furthermore, it is possible to track position and stress state of each discrete volume and we would like to develop a contact interface for interacting with classical Discrete Element formulations. We also would like to develop a similar mathematical formulation for topology optimization problems, by redefining the elastic potential energy in terms of a solidification field that predicts the optimal shape of an elastic body to the given boundary conditions.

Appendix A. Analytical solution for time

In order to get an accurate approximation of the temporal derivative in the interval ${\textstyle t\in [0,\Delta t]}$, we replace the stress state function of time ${\textstyle S_{i}(t):\mathbb {R} \rightarrow \mathbb {R} ^{\hbox{dim}}}$ into the differential equation (5.2),

$\left(1-{P}\left({\frac {t}{\Delta t}}\right)\right)~S_{i}^{0}+{P}\left({\frac {t}{\Delta t}}\right)~S_{i}=\rho {\ddot {\boldsymbol {u}}}_{i},$

(A.1)

stress function is defined in equation (5.4). Reordering terms we have

$S_{i}^{0}+{P}\left({\frac {t}{\Delta t}}\right)~\left(S_{i}-S_{i}^{0}\right)=\rho {\ddot {\boldsymbol {u}}}_{i},$

(A.2)

with initial conditions

${\boldsymbol {u}}_{i}(0)={\boldsymbol {u}}_{i}^{0}~~~~{\hbox{and}}~~~~{\dot {\boldsymbol {u}}}_{i}(0)={\dot {\boldsymbol {u}}}_{i}^{0}$

(A.3)

In order to solve (A.2) by means of Laplace Transform

$\mathbf {F} _{i}(s)={\mathcal {L}}\{{\boldsymbol {u}}_{i}(t)\}=\int \limits _{0}^{\infty }{\boldsymbol {u}}_{i}(t)~e^{-st}dt,~~~~~~~~\mathbf {F} _{i}(s):\mathbb {R} \rightarrow \mathbb {R} ^{\hbox{dim}},$

(A.4)

we change from time domain in equation (A.2) to frequency domain

$\left({\frac {1}{s}}\right)~S_{i}^{0}+{\mathcal {L}}_{P}(s)~\left(S_{i}-S_{i}^{0}\right)=\rho \left(s^{2}~\mathbf {F} _{i}(s)-s~{\boldsymbol {u}}_{i}^{0}-{\dot {\boldsymbol {u}}}_{i}^{0}\right),$

(A.5)

where ${\textstyle s}$ is the frequency variable, ${\textstyle {\mathcal {L}}_{P}}$ is the Laplace transform of ${\textstyle {P}(t/\Delta t)}$,

${\mathcal {L}}_{P}(s)={\mathcal {L}}\left\{{P}\left({\frac {t}{\Delta t}}\right)\right\},~~~~{\mathcal {L}}_{P}(s):\mathbb {R} \rightarrow \mathbb {R} ,$

(A.6)

and the Laplace transform of acceleration term includes initial conditions (A.3)

${\mathcal {L}}\{{\ddot {\boldsymbol {u}}}_{i}\}=\left(s^{2}\mathbf {F} _{i}(s)-s~{\boldsymbol {u}}_{i}^{0}-{\dot {\boldsymbol {u}}}_{i}^{0}\right)$

(A.7)

We can rewrite equation (A.5) as

$\mathbf {F} _{i}(s)=\left({\frac {1}{s}}\right){\boldsymbol {u}}_{i}^{0}+\left({\frac {1}{s^{2}}}\right){\dot {\boldsymbol {u}}}_{i}^{0}+\left({\frac {2}{s^{3}}}\right)\left({\frac {1}{2\rho }}\right)S_{i}^{0}+\left({\frac {1}{s^{2}}}\right){\mathcal {L}}_{P}(s)~\left({\frac {1}{\rho }}\right)\left(S_{i}-S_{i}^{0}\right)$

(A.8)

and applying the inverse Laplace transform, ${\textstyle {\boldsymbol {u}}_{i}(t)={\mathcal {L}}^{-1}\{\mathbf {F} _{i}(s)\}}$, we obtain

${\boldsymbol {u}}_{i}(t)={\boldsymbol {u}}_{i}^{0}+t~{\dot {\boldsymbol {u}}}_{i}^{0}+{\frac {1}{2\rho }}~t^{2}~S_{i}^{0}+\left({\frac {1}{\rho }}\right){C}_{P}(t)~\left(S_{i}-S_{i}^{0}\right),$

(A.9)

where ${\textstyle {C}_{P}(t)}$ is a convolution. Such convolution is defined as

${C}_{P}(t)=\left({P}\left({\frac {t}{\Delta t}}\right)*t\right)(t)=\int \limits _{0}^{t}{P}\left({\frac {\tau }{\Delta t}}\right)(t-\tau )d\tau ,$

(A.10)

and it derivative is

${\dot {C}}_{P}(t)=\left({P}\left({\frac {t}{\Delta t}}\right)*1\right)(t)=\int \limits _{0}^{t}{P}\left({\frac {\tau }{\Delta t}}\right)d\tau ,$

(A.11)

where ${\textstyle \tau }$ is the integration variable. Developing previous definitions we get

${C}_{P}(t)=t~{\dot {C}}_{P}(t)-{\frac {t^{2}}{2}}{P}\left({\frac {t}{\Delta t}}\right)+{\frac {t^{3}}{6}}{\dot {P}}_{c}\left({\frac {t}{\Delta t}}\right)+O\left(t^{4}~{\ddot {P}}_{c}\left({\frac {t}{\Delta t}}\right)~\right)$ (A.12) $=t~{\dot {C}}_{P}(t)+\sum \limits _{n=0}^{\infty }{\frac {(-1)^{(n+1)}}{(n+2)!}}~~t^{(n+2)}~{\frac {d^{n}}{dt^{n}}}{P}\left({\frac {t}{\Delta t}}\right),$ (A.13)

Finally, analytical solution of (A.2) is equation (A.9) and it is completely dependent of the shape function ${\textstyle {P}}$. This solution is used for building an accurate numerical squeme for discretizing time.

Appendix B. Polynomial shape functions for time

The analytical solution (5.11) of equation (5.2) is used to generate numerical schemes for time discretization, these approximations has a similar structure but distinct coefficients that depend on the shape function ${\textstyle {P}}$ used for time variation of stress state. In this appendix we propose two families of polynomial functions in order to get a continuous stress state in contiguous time steps, such polynomial functions meet the condition of ${\textstyle {\dot {C}}_{P}^{\Delta t}=\Delta t}$ for producing stable schemes in terms of total energy.

The first polynomial family is defined by

${P}(t)=\left(\sum \limits _{i=0}^{p}(1+2~i)\right)\left(t^{p}-t^{(p+1)}\right)+t^{(p+1)},$

(B.1)

where ${\textstyle p+1}$ is the polynomial order. The first three functions generated by this equation are shown in Table 2 and Figure 50 depicts the curves.

Table. 2 First few polynomials generated with equation (B.1) and its respective convolutions ${C}_{P}(t)$.

${\textstyle p}$	Shape function	Convolution
1	${P}(t)=4t-3t^{2}$	${C}_{P}(t)=(8\Delta t~t^{3}-3t^{4})/(12\Delta t^{2})$
2	${P}(t)=9t^{2}-8t^{3}$	${C}_{P}(t)=(3\Delta t~t^{4}-8t^{5})/(20\Delta t^{3})$
3	${P}(t)=16t^{3}-15t^{4}$	${C}_{P}(t)=(8\Delta t~t^{5}-5t^{6})/(10\Delta t^{4})$

Figure 50: Curves for first few polynomials generated with equation (B.1)

The second polynomial family is given by

${P}(t)=\left(1+\sum \limits _{i=1}^{p}(3~i+1)\right)t^{p}-\left(\sum \limits _{i=1}^{p}2(2~i+1)\right)t^{(p+1)}+\left(\sum \limits _{i=1}^{p}(i+1)\right)t^{(p+2)},$

(B.2)

where ${\textstyle p+2}$ is the polynomial order. The first few functions generated by this equation are shown in Table 3 and Figure 51 depicts the curves.

Table. 3 First few polynomials generated with equation (B.2) and its respective convolutions ${C}_{P}(t)$.

${\textstyle p}$	Shape function	Convolution
1	${P}(t)=5t-6t^{2}+2t^{3}$	${C}_{P}(t)=(3t^{5}-15\Delta t~t^{4}+25\Delta t^{2}t^{3})/(30\Delta t^{3})$
2	${P}(t)=12t^{2}-16t^{3}+5t^{4}$	${C}_{P}(t)=(5t^{6}-24\Delta t~t^{5}+30\Delta t^{2}t^{4})/(30\Delta t^{4})$

Figure 51: Curves for first few polynomials generated with equation (B.2)

Bibliography

[1] A.K. Slone and K. Pericleous and C. Bailey and M. Cross. (2002) Dynamic fluid-structure interaction using finite volume unstructured mesh procedures, Volume 80. Elsevier. Computers and Structures 371-390

[2] A.K. Slone and K. Pericleous and C. Bailey and M. Cross and C. Bennett. (2004) A finite volume unstructured mesh approach to dynamic fluid-structure interaction: an assessment of the challenge of predicting the onset flutter, Volume 28. Elsevier. Applied Mathematical Modeling 211-239

[3] C. Bailey and G. A. Taylor and M. Cross and P. Chow. (1999) Discretisation procedures for multi-physics phenomena, Volume 103. Elsevier. Journal of Computational and Applied Mathematics 3-17

[4] M. Souli and A. Ouahsine and L. Lewin. (2000) ALE formulation for fluid-structure interaction problems, Volume 190. Elsevier. Computer Methods in Applied Mechanics and Engineering 659-675

[5] J. Lubliner and J. Oliver and S. Oller and E. Oñate. (1989) A Plastic-Damage model for concrete, Volume 25. Elsevier. International Journal of Solids and Structures 299-326

[6] Francisco Armero and S. Oller. (2000) A general framework for continuum damage models. I. Infinitesimal plastic damage models in stress space, Volume 37. Elsevier. International Journal of Solids and Structures 7409-7436

[7] E. Oñate and M. Cervera and C. Zienkiewicz. (1994) A finite volume format for structural mechanics, Volume 37. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 181-201

[8] S. R. Idelsohn and E. Oñate. (1994) Finite volumes and finite elements: Two 'Good Friends', Volume 37. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 3323-3341

[9] Y. D. Fryer and C. Bailey and M. Cross and C. H. Lai. (1991) A control volume procedure for solving the elastic stress-strain equations on an unstructured mesh, Volume 15. Elsevier. Applied Mathematical Modeling 639-645

[10] C. Bailey and M. Cross. (1995) A Finite Volume procedure to solve elastic solid mechanics problems in three dimensions on an unstructured mesh, Volume 38. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 1757-1776

[11] A. K. Slone and C. Bailey and M. Cross. (2003) Dynamic solid mechanics using finite volume methods. Elsevier. Applied Mathematical Modeling 69-87

[12] I. Demirdzic and D. Martinovic. (1993) Finite volume method for thermo-elasto-plastic stress analysis, Volume 109. Elsevier. Computer Methods in Applied Mechanics and Engineering 331-349

[13] I. Demirdzic and S. Muzaferija. (1994) Finite volume methods for stress analysis in complex domains, Volume 137. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 3751-3766

[14] I. Demirdzic and S. Muzaferija and M. Peric. (1997) Benchmark solutions of some structural analysis problems using finite-volume method and multigrid acceleration, Volume 40. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 1893-1908

[15] I. Bijelonja and I Demirdzic and S. Muzaferija. (2005) A finite volume method for large strain analysis of incompressible hyperelastic materials, Volume 64. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 1594-1609

[16] I. Bijelonja and I. Demirdzic and S. Muzaferija. (2006) A finite volume method for incompressible linear elasticity, Volume 195. Elsevier. Computer Methods in Applied Mechanics and Engineering 6378-6390

[17] H. Jasak and H. G. Weller. (2000) Application of the finite volume method and unstructured meshes to linear elasticity, Volume 48. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 267-287

[18] Z. Tukovic and A. Ivankovic and A. Karac. (2012) Finite-volume stress analysis in multi-material linear elastic body. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering

[19] Tian Tang and Ole Hededal and Philip Cardiff. (2015) On finite volume method implementation of poro-elasto-plasticity soil model. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering

[20] M. J. Borden and C. V. Verhoosel and M. Scott and T. J. R. Hughes and C. M. Landis. (2012) A phase-field description of dynamic brittle fracture. Elsevier. Computer Methods in Applied Mechanics and Engineering 217-220

[21] P. Cardiff and Z. Tukovic and H. Jasak and A. Ivankovic. (2016) A block-coupled Finite Volume methodology for linear elasticity and unstructured meshes. Elsevier. Computers and Structures 100-122

[22] Marcio A. A. Cavalcante and Marek-Jerzy Pindera. (2012) Generalized Finite-Volume Theory for Elastic Stress Analysis in Solid Mechanics - Part I: Framework. The American Society of Mechanical Engineers. Journal of Applied Mechanics

[23] Jan Martin Nordbotten. (2014) Cell-centered finite volume discretizations for deformable porous media. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering

[24] B. Li and Z. Chen and G. Huan. (2003) Control volume function approximation methods and their applications to modeling porous media flow. Elsevier. Advances in water resources 435-444

[25] B. Li and Z. Chen and G. Huan. (2004) Control volume function approximation methods and their applications to modeling porous media flow II: the black oil model. Elsevier. Advances in water resources 99-120

[26] M. Cervera and M. Chiumenti. (2006) Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique. Elsevier. Computer Methods in Applied Mechanics and Engineering

[27] M. Cervera and L. Pela and R. Clemente and P. Roca. (2010) A crack-tracking technique for localized damage in quasi-brittle materials. Elsevier. Engineering Fracture Mechanics

[28] M. Cervera and M. Chiumenti and R. Codina. (2011) Mesh objective modeling of cracks using continuous linear strain and displacement interpolations. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering

[29] A. D. Hanganu and E. Oñate and A. H. Barbat. (2002) A finite element methodology for local/global damage evaluation in civil engineering structures. Elsevier. Computers and Structures 1667-1687

[30] E. Oñate. (2000) Desarrollos y aplicaciones de modelos de fractura en la escuela de ingenieros de caminos de barcelona. Escuela Técnica Superior de Ingenieros de Caminos, Canales y Puertos. Universidad Politécnica de Catalunya

[31] E. Oñate and J. Rojek. (2004) Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems. Elsevier. Computer Methods in Applied Mechanics and Engineering 3087-3128

[32] C. Miehe and M. Hofacker and F. Welschinger. (2010) A phase field model for rate independent crack propagation: robust algorithmic implementation based on operator splits. Elsevier. Computer Methods in Applied Mechanics and Engineering

[33] C. Miehe and M. Hofacker and F. Welschinger. (2010) Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. Elsevier. Computer Methods in Applied Mechanics and Engineering

[34] J. Hoon Song and H. Wang and T. Belytschko. (2007) A comparative study on finite element methods for dynamic fracture. Computational Mechanics

[35] N. Moes and T. Belytschko. (2002) Extended finite element method for cohesive crack growth. Elsevier. Engineering Fracture Mechanics 813-833

[36] G. Zi and T. Belytschko. (2003) New crack-tip elements for XFEM and applications to cohesive cracks. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering 2221-2240

[37] L. Mishnaevsky and S. Shmauder. (1998) FE-simulation of crack growth using a damage parameter and the cohesive zone concept. In ECF 12- Fracture from defects, Proc. 12th European Conference on Fracture. London, EMAS 2, 1053-1059.

[38] L. Mishnaevsky and N. Lippmann and S. Shmauder. (2003) Computational modeling of crack propagation in real microstructures of steels and virtual testing of artificially designed materials. International Journal of Fracture

[39] P. Cundall and O. Strack. (1979) A discrete numerical method for granular assemblies, Volume 1. Geotechniques 47-65

[40] J. Rojek and F. Zarate and C. Agelet and C. Gilbourne and P. Verdot. (2004) Discrete element modelling and simulation of sand mould manufacture for the lost foam process. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering

[41] O. C. Zienkiewicz and R. L. Taylor. (2005) The Finite Element Method. For Solid and Structural Mechanics. Elsevier

[42] J.P. Morris and M.B. Rubin and G.I. Blocka and M.P. Bonnera. (2006) Simulations of fracture and fragmentation of geologic materials using combined FEM/DEM analysis. International Journal of Impact Engineering 463-473

[43] A. Munjiza. (2004) The combined finite-discrete element-method. John Wiley & sons

[44] D. Potyondy and P. Cundall. (2004) A bonded-particle model for rock. Elsevier. International Journal of rock mechanics and mining sciences 1329-1364

[45] J. Rojek and E. Oñate. (2007) Multiscale analysis using a coupled discrete/finite element model. Interaction and Multiscale Mechanics

[46] E. Oñate and C. Labra and F. Zárate and J. Rojek. (2012) Modelling and simulation of the effect of blast loading on structures using an adaptive blending of discrete and finite element methods. Taylor & Francis group. Risk Analysis, Dam Safety, Dam Security and Critical Infrastructure Management

[47] Francisco Zárate and Eugenio Oñate. (2015) A simple FEM-DEM technique for fracture prediction in materials and structures, Volume 2. Springer. Computational Particle Mechanics 301-314

[48] E. B. Tadmor and M. Ortiz and R. Phillips. (1996) Quasicontinuum analysis of defects in solids, Volume 73. Philosophical Magazine A. 1529-1563

[49] E. B. Tadmor and Rob Phillips. (1996) Mixed Atomistic and Continuum Models of Deformation in Solids, Volume 12. Americal Chemical Society 4529-4534

[50] Gregory J. Wagner and W. K. Liu. (2003) Coupling of Atomistic and Continuum Simulations using a Bridging Scale Decomposition. Journal of Computational Physics

[51] S.P. Xiao and T. Belytschko. (2004) A bridging domain method for coupling continua with molecular dynamics. Computer Methods in Applied Mechanics and Engineering

[52] Mei Xu and Ted Belytschko. (2004) Conservation properties of the bridging domain method for coupled molecular/continuum dynamics. John Wiley & Sons, Ltd. International Journal for Numerical Methods in Engineering

[53] Z. Chen and G. Huan and B. Li. (2003) Modeling 2D and 3D horizontal wells using CVFA. Comm. Math. Sci.

[54] Richard Lipton and Donald Rose and Robert Endre Tarjan. (1979) Generalized nested dissection, Volume 16. SIAM. Journal on Numerical Analysis 346-358

[55] S. P. Timoshenko and J. N. Goodier. (1970) Theory of Elasticity. McGraw-Hill, 3 Edition

[56] Jan Gerrit Rots. (1988) Computational Modeling of Concrete Fracture. Technische Universiteit Delft

[57] Jose Rafael Capua Proveti and Gerard Michot. (2006) The Brazilian test: a tool for measuring the toughness of a material and its brittle to ductile transition. International Journal of Fracture 455-460

[58] T. N. Bittencourt. (1996) Quasi-automatic simulation of crack propagation for 2D LEFM problems, Volume 55. Elsevier. Engineering Fracture Mechanics 321-334