## 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

 ${\displaystyle {\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

 ${\displaystyle \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

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

and

 ${\displaystyle \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

 ${\displaystyle {\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

 ${\displaystyle 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

 ${\displaystyle {\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

 ${\displaystyle {\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 ${\displaystyle \Gamma }$, under boundary conditions. The right side shows the damage field approximation of the fracture surface.

 ${\displaystyle \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

 ${\displaystyle \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

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

where

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

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

 ${\displaystyle {\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

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle {\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

 ${\displaystyle {\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

 ${\displaystyle {\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

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle {\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

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

to obtain

 ${\displaystyle \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 }$,

 ${\displaystyle \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 }$,

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle 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

 ${\displaystyle 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 =\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

 ${\displaystyle {\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

 ${\displaystyle {\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

 ${\displaystyle \rho {\ddot {\boldsymbol {u}}}-\nabla \cdot {\boldsymbol {\sigma }}=0,}$ (2.33.a) ${\displaystyle 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

 ${\displaystyle {H}-\psi _{e}^{+}\geq 0,}$ (2.34.a) ${\displaystyle {\dot {H}}\geq 0,}$ (2.34.b) ${\displaystyle \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

 ${\displaystyle {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

 ${\displaystyle \nabla \cdot {\boldsymbol {\sigma }}=\rho {\ddot {\boldsymbol {u}}},}$ (2.36.a) ${\displaystyle \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

 ${\displaystyle {\boldsymbol {\sigma }}\mathbf {n} ={\boldsymbol {b}}_{N}(\mathbf {x} ,t),\mathbf {x} \in \Gamma _{N},t\in [0,T],}$ (2.37.a) ${\displaystyle {\boldsymbol {u}}={\boldsymbol {u}}_{D}(\mathbf {x} ,t),\mathbf {x} \in \Gamma _{D},t\in [0,T],}$ (2.37.b) ${\displaystyle \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

 ${\displaystyle {\boldsymbol {u}}(\mathbf {x} ,0)={\boldsymbol {u}}_{0}(\mathbf {x} ),\mathbf {x} \in \Omega ,}$ (2.38.a) ${\displaystyle {\dot {\boldsymbol {u}}}(\mathbf {x} ,0)={\dot {\boldsymbol {u}}}_{0}(\mathbf {x} ),\mathbf {x} \in \Omega ,}$ (2.38.b) ${\displaystyle {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

 ${\displaystyle {\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

 ${\displaystyle 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, ${\displaystyle V_{i}}$, are defined by its radii ${\displaystyle r_{i}}$ and the planes orthogonal to the lines connecting the centroids ${\displaystyle \mathbf {x} _{i}}$ and ${\displaystyle \mathbf {x} _{j}}$ at the point ${\displaystyle {\boldsymbol {q}}_{ij}}$, for all ${\displaystyle 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 ${\displaystyle V_{i}}$ (let ${\displaystyle i=1}$) is bounded by the remaining volumes ${\displaystyle V_{j}}$ (let ${\displaystyle j=2,3,4}$). The right side illustrates the partition, ${\displaystyle P_{h}}$, of the domain, ${\displaystyle \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

 ${\displaystyle 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

 ${\displaystyle \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,

 ${\displaystyle \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

 ${\displaystyle 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}}$,

 ${\displaystyle \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 ${\displaystyle P_{h}}$ is the discretization of the domain ${\displaystyle \Omega }$ into ${\displaystyle N}$ control volumes. The boundary of the control volumes, ${\displaystyle \partial V_{i}}$, is conformed by ${\displaystyle N_{i}}$ flat faces, denoted ${\displaystyle e_{ij}}$. The unit vector ${\displaystyle \mathbf {n} _{ij}}$ is normal to the face ${\displaystyle e_{ij}}$. The faces of the volumes adjacent to the boundary ${\displaystyle \Gamma _{N}}$ are integrated using the condition ${\displaystyle {\boldsymbol {b}}_{N}}$.

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

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

 ${\displaystyle \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,

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle \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}}$,

 ${\displaystyle {\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

 ${\displaystyle {\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

 ${\displaystyle {\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) ${\displaystyle ={\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}}}$.

 ${\displaystyle {\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

 ${\displaystyle {\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

 ${\displaystyle \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

 ${\displaystyle \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}}$,

 ${\displaystyle \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

 ${\displaystyle \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}}$,

 ${\displaystyle \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 ${\displaystyle 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 ${\displaystyle \mathbf {H} _{ij}}$ of the control volume ${\displaystyle V_{i}}$.

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

 ${\displaystyle 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

 ${\displaystyle \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,

 ${\displaystyle \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 ${\displaystyle V_{i}}$ is being integrated. The integral over the subface ${\displaystyle e_{ijk}}$ is calculated using the polynomial approximation of shaded simplex. The integration point must be mapped to (b) Normalized space ${\displaystyle [-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,

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

 ${\displaystyle {P}_{c}\left({\boldsymbol {\xi }}_{[q]}\right),~~~{\hbox{ if}}~~q\in [1,{\hbox{ dim}}],}$ (3.22.a) ${\displaystyle 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

 ${\displaystyle \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

 ${\displaystyle \varphi _{q}({\boldsymbol {\xi }}_{q})=1~~~~{\hbox{ for any }}{\boldsymbol {\xi }}_{q}{\hbox{ forming the simplex,}}}$ (3.24) ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle {\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

 ${\displaystyle {\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) ${\displaystyle =\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

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle {\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,

 ${\displaystyle {\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

 ${\displaystyle \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) ${\displaystyle =\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}}}$ ${\displaystyle =\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

 ${\displaystyle \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),

 ${\displaystyle \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

 ${\displaystyle {\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),}$

 ${\displaystyle 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

 ${\displaystyle \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) ${\displaystyle =\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) ${\displaystyle =\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 ${\displaystyle e_{ij}}$, since its normal is not necessary aligned with ${\displaystyle \mathbf {x} _{\vec {ij}}}$. The integration point is projected into ${\displaystyle \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

 ${\displaystyle {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

 ${\displaystyle 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}}$

 ${\displaystyle 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

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

The first derivative is simply calculated as

 ${\displaystyle {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.

 Smoothness Homeostatic spline ${\textstyle c=0}$ ${\displaystyle {P}_{0}(z)=z}$ ${\displaystyle c=1}$ ${\displaystyle {P}_{1}(z)=3z^{2}-2z^{3}}$ ${\displaystyle c=2}$ ${\displaystyle {P}_{2}(z)=10z^{3}-15z^{4}+6z^{5}}$ ${\displaystyle c=3}$ ${\displaystyle {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 ${\displaystyle c=0}$ to ${\displaystyle 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

 ${\displaystyle {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

 ${\displaystyle {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

 ${\displaystyle 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 ${\displaystyle {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

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle \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),

 ${\displaystyle \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),

 ${\displaystyle {\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 ${\displaystyle \mathbf {x} _{i}}$. (b) The initial mesh is generated from a FEM-like triangular mesh. The calculation points ${\displaystyle \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

 ${\displaystyle \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 ${\displaystyle 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}}$,

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle 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 ${\displaystyle V_{i}}$ is partitioned into pyramids ${\displaystyle V_{ijk}}$, which turns to be triangles in 2D. Pyramids bases correspond to the subfaces ${\displaystyle e_{ijk}}$ resulting from the intersection with the local Delaunay triangulation, and all of them share the calculation point ${\displaystyle \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),

 ${\displaystyle \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)}$ ${\displaystyle =\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,

 ${\displaystyle {\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) ${\displaystyle =\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

 ${\displaystyle \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,

 ${\displaystyle \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,

 ${\displaystyle \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) ${\displaystyle =\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) ${\displaystyle \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

 ${\displaystyle \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) ${\displaystyle =\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) ${\displaystyle =\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) ${\displaystyle =\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

 ${\displaystyle 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}}$ ${\displaystyle =\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

 ${\displaystyle {\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

 ${\displaystyle 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

 ${\displaystyle 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,

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

with the intention of considering internal forces in the approximation,

 ${\displaystyle 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

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

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

 ${\displaystyle 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 ${\displaystyle {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),

 ${\displaystyle {\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

 ${\displaystyle {\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

 ${\displaystyle {\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)~+}$ ${\displaystyle 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),

 ${\displaystyle {\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)~+}$ ${\displaystyle 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]}$,

 ${\displaystyle {\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,

 ${\displaystyle {\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)~+}$ ${\displaystyle 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)~+}$ ${\displaystyle 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,

 ${\displaystyle \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}+}$ ${\displaystyle \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}+}$ ${\displaystyle \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

 ${\displaystyle -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

 ${\displaystyle 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

 ${\displaystyle A={\sqrt {u_{0}+{\frac {m}{k}}{\dot {u}}_{0}}},}$ (5.14) ${\displaystyle \omega ={\sqrt {\frac {k}{m}}}}$ (5.15) ${\displaystyle \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

 ${\displaystyle 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

 ${\displaystyle -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

 ${\displaystyle 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) ${\displaystyle =\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

 ${\displaystyle {\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 ${\displaystyle \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

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