Virtual element formulation for isotropic damage

Abstract In the paper we present a low-order virtual element formulation for modelling the strain-softening response of quasi-brittle materials. For this purpose, a formulation in two-dimensions is considered, with virtual elements having arbitrary shape. The method is based on minimization of an incremental energy expression, with a novel construction of the stabilization energy for isotropic elasto-damage. A set of numerical examples, illustrating the efficiency of the proposed method, complements the paper.


Introduction
The Virtual Element Method (VEM) has been recently developed in keeping with Mimetic Finite Difference [12], characterizing it as a Galerkin finite element-type reformulation. The basic principles have been presented in the seminal work by [5]. Subsequent significant contributions in explaining the theoretical basis and providing examples of implementations can be found in [7,13,19]. The VEM permits the numerical solution of boundary value problems on arbitrary polyhedral meshes, including convex and non-convex elements, very stretched elements, hanging nodes and collapsing nodes. The great flexibility in dealing with very general geometries and the robust mathematical basis of the method pave the way for possible use in very general cases. They are ranging from crack propagation in fractured solids, to modelling the texture evolution in polycrystalline materials, up to reproducing the complex behaviour of structured materials. Applications have been devoted so far to linear elastic two-and threedimensional problems [1,6,19], discrete fracture network simulations [10], eigenvalue problems [29], contact problems [45], topology optimization [20] and nonlinear problems [2,8,14]. Stabilization procedures for the virtual element method, which are well known from the work of [9] for finite elements, are described in [13] for linear Poisson problems. In the VEM formulation a stabilization term is mandatory. The structure of the VEM, indeed, typically comprises a term in the weak formulation or energy functional in which the quantity φ v , here deformation, is replaced by its projection Πφ v onto a polynomial space.
This results in a rank-deficient structure, so that it is necessary to add a stabilization term to the formulation, see [5,6,14], where in the latter the scalar stabilization parameter of the linear case being replaced by one that depends on the fourth-order elasticity tensor. In the framework of compressible and incompressible nonlinear elasticity, a novel stabilization technique has been proposed in [44], inspired by an idea first proposed in [30], generalized in [11] and simplified in [24]. The key innovative aspect is to add to the positive semidefinite mean strain energy, Ψ, a positive-definite energy, Ψ, evaluated using full quadrature. Moreover, for consistency a term involving Ψ as a function of the mean strain is subtracted. The resulting strain energy is the sum of the original energy as a function of the projected displacement, Ψ(Πu u u v ), and the term Ψ(Πu u u v ) to which a positive definite stabilization energy as a function of the displacement and its projection are respectively added and subtracted. In the terms involving Ψ(u u u v ), the quadrature is carried out by constructing a triangular mesh in the element, without introducing additional degrees of freedom, since the nodal points are those of the original element. The same stabilization has been successfully applied also to problems of finite strain plasticity [43]. In the present paper, this stabilization technique is modified and exploited in the framework of a virtual element formulation for 2D scalar damage problems. The main idea is to test the effectiveness of the VEM in dealing with highly localized strains due to material instabilities, typically exhibited by quasi-brittle materials undergoing severe loading conditions. In the last decades, many efforts have been focused on the formulation of continuum damage models able to predict the irreversible phenomena related to the onset and evolution of void formation, micro-cracking and strain-softening in quasi-brittle materials. Based on Kachanov's pioneered work, [22], different phenomenological models have been proposed both accounting for isotropic [16,25,27,33,40,41], or anisotropic [21,26,34,41] damage response. Special attention has been devoted to develop strategies able to overcome the well-known spurious mesh sensitivity problems, occurring in finite elements computations when, in the presence of softening behaviour, the governing differential equations may loose ellipticity and resort to ill-posed boundary value problems. Successful remedies have been proposed by recourse to non-local continuum theories, [17,18], which possess intrinsic regularization properties due to the natural introduction of characteristic lengths, or either recourse to viscous regularization [31,42]. Alternative solutions to the pathological mesh dependency in damage models are non-local constitutive models, see e.g. [3] for a comprehensive survey, or local models properly enriched by a dependence of the material properties on the element size, as in [32,33,40]. We focus here on an isotropic damage law, in the framework of linearised kinematics, with different thresholds in tensions and compression, as proposed in [33], and we adopt as regularization technique, alternatively, the last two aforementioned approaches to which we will refer as "nonlocal" and "local". The paper is organized as follows: in Section 2 the governing equations of the local and non-local scalar damage model are briefly recalled; section 3 is devoted to the construction of the linear ansatz functions of the proposed virtual element approach; the virtual element formulation is, then, fully developed in Section 4; a number of illustrative applications is proposed in Section 5, and, finally, in Section 6 some concluding remarks are reported.

Governing equations for isotropic damage model
Consider an initially elastic body that occupies the bounded domain Ω ⊂ R 2 . Let Γ = Γ D ∪ Γ N be the boundary of Ω, with Γ D the Dirichlet and Γ N the Neumann boundaries, such that Γ D ∩ Γ N = ∅. Each material point x is characterized by the displacement field u that is the primary unknown variable. The symmetric Cauchy stress σ satisfies the linear momentum balance with f being the body force. The Dirichlet and Neumann boundary conditions hold, respectively, as: with u the prescribed displacement, n the outward unit normal vector and t the surface traction.
The strain-displacement relation is given by For isotropic damage a scalar variable d is introduced, satisfying 0 ≤ d ≤ 1. At a given point of the damaged material, the free energy is expressed as a function approaching zero as the damage d increases: Here Ψ 0 is the initial undamaged elastic energy that reads as where σ EL = Dε = λtr(ε)I + 2µε is the stress for an initially homogeneous isotropic elastic material, with D being the elasticity tensor, and λ and µ the Lamé constants. The energetic consistency of the constitutive model is ensured by the fulfilment of the Clausius-Planck inequality, [28], whereḊ is the rate of the mechanical energy dissipation defined for arbitrary infinitesimal variationsε. Using the Coleman's method [15] it follows that the constitutive relation is Following [33,41], an equivalent effective stress τ is defined as a suitable energy norm of the undamaged stress tensor σ 0 and used to compare different material states. Here we adopt a damage model with different thresholds in tension and compression. Hence we choose τ as where ζ is a weight factor depending on the elastic principal stresses σ 0 i , • is the Macaulay bracket, and n = f c / f t is the ratio between compressive f c and tensile f t strength of the material. The damage criterion is defined via the limit damage surface, i.e. a function F(τ t , r t ) that splits the admissible stress space into the elastic domain (when F < 0) and the damage domain (when F = 0). It depends both on the equivalent effective stress τ t and on a material parameter representing the damage threshold r t at the current time t. The most general form of F(τ, r) among different possibilities is where G(•) is a suitable monotonic scalar function. The damage is governed by the following evolution equationṡ withγ being the damage multiplier adopted to define the loading-unloading conditions (equivalent to the plastic multiplier in the rate independent plasticity). The following Kuhn-Tucker relations have also to be satisfied It is possible to directly integrate the evolution of the internal variables, as in [41,33], and obtain where r 0 is a characteristic of the material, i.e. the initial damage threshold for the virgin material and max (τ s ) is the maximum equivalent effective stress attained until current time t.
Among different possible choices of explicit functions for the scalar damage d, two options are frequently used in literature, e.g. [33,40]: 1. Linear damage law where , E 0 is the elastic modulus, and G f the fracture energy per unit area.

Exponential damage law
where

Regularization technique
In order to overcome the well known problems of strain localization and mesh dependency, characterizing the numerical response in the case of strain softening materials, we adopt two alternative regularization techniques. The first approach consists in introducing an explicit dependence of the constitutive model on the size of the numerical discretization element. Following [32,33,40] among others, it is possible to modify the expressions of H 1 in (14) and H 2 in (15) as ≥ 0, with l i c being the characteristic length that in 2D has been defined either as l 1 c = √ Ω e , [33,40], or as l 2 c =Ω e / A (∂φ/∂x)dA, [32], where Ω e is the area of the element, φ is a non-dimensional, continuous and differentiable function describing the displacement in the localization band and x is a linear coordinate, perpendicular to the localization band. The second approach is the well-established non-local damage theory proposed by [36] and widely adopted in later works, see [3]. The key points of the non-local damage theory are then briefly recalled. The non-local equivalent stress τ can be defined as the weighted average of the local equivalent effective stress τ over a representative spherical volume surrounding each material point x, that is where Ω is the area of the considered body, Ω r (x) = Ω ψ(x − s) ds is the representative volume and ψ(x − s) is a properly chosen weighted function that usually, see [4,40], takes the form The internal length l 3 c is introduced, acting as a parameter to control the localization zone. As pointed out in [4], its values depend on the material considered and, i.e. in the case of standard concrete this parameter can be taken as 3d a , that is directly related to the maximum size d a of aggregates. The non-local equivalent stress τ replaces its local counterpart τ in the damage laws (14), for the linear case, and (15), for the exponential case.
As an alternative to the weak form, it is possible to construct a scalar damage virtual element starting from the pseudo potential energy This form is particularly suitable for codes automatically generated using the software ACEGEN, see [23].

The VEM formulation
In a VEM discretization, the domain Ω is partitioned into non-overlapping polygonal elements of general shape and not necessarily convex. Linear ansatz functions are here considered. It follows that element nodes are placed only at the vertices of the polygonal elements. We denote the discrete space of test functions on Ω by V v , and for a conforming approach we require that V v ⊂ V (V being the space of test functions that satisfy Dirichlet boundary conditions of Γ D ). This means that the adopted shape (or basis) functions in V v comprise continuous functions whose restriction to an element Ω e includes (but is larger than) linear functions. Furthermore, as a consequence of the choice of linear ansatz functions, the restriction of the element shape functions to the element boundaries are linear functions. Referring to a generic element, the nodal displacement vector lists the displacement components u i = {u xi , u yi }, i=1,..n V of the n V vertices. As usual in the virtual element method, we split the ansatz space into a linear part, the projection Πu v , and a remainder where we choose with the unknowns a being introduced.
In consistency with with [8,44], the projection in the element Ω e satisfies where N is the outward normal to the boundary Γ e of the domain Ω e . Once Πu v | e is known, the strain tensor can be computed as (21) it follows that the projected gradient is constant and has the form The evaluation of the boundary integral in (22) can be performed exactly using the trapezoidal Gauss-Lobatto rule. To this purpose, a linear ansatz for the deformation along the element edges is introduced, that is with ξ k = x k /L k being a local dimensionless coordinate on the boundary segment k of length L k of the virtual element. Moreover, 1-2 are the local nodes of the boundary edges. It is, thus, possible to directly relate the unknowns a 3 to a 6 to normal vectors of the boundary segments of the elements and to nodal displacement d.o.f., see [44], as where summations are taken over the n V number of segments forming the element edges. Also, (N x ) k and (N y ) k are the components of the vector N k , normal to the k-th segment defined as with {X k , Y k }, k=1,2 being the local coordinates of the two vertices of the k-th segment. The remaining unknowns a 1 and a 2 can be finally evaluated as in [8], by ensuring that the sums of the displacements u v and the projections Πu v evaluated at each node, are the same. It follows that

Construction of the Virtual Element
The pseudo potential energy in (18) can be rewritten by exploiting the split in (20). By summing up all the element contributions for the N e virtual elements of area Ω e , it gives where U c is a consistency term related only to the projection of the displacements, and thus involving a constant strain field over the elements, and U stab is a stabilization term, required to avoid element rank deficiencies in the case of elements with n V >3, as discussed in [44]. Let's start from the first term in (28), that takes the form where the integrals are evaluated, respectively, over the element area Ω e and the associated boundary Γ e = Γ N e ∪ Γ D e . The pseudo strain energy is then and both the damage variable d and the components of the strain tensor ε are constant in each element Ω e . As a consequence, the following relation holds with Ψ(ε(Πu v | e ), d) a non-linear function of the displacement nodal d.o.f. u v and of the scalar damage variable d.
The stabilization term has now to be derived for the case of scalar damage. This derivation is based on the non-linear stabilization procedure proposed in [44] for compressible and incompressible finite deformations. Consistently with this approach, we define U stab on the virtual element Ω e as the sum of two contributions, depending either on the displacement field or its projection and defined on the basis of a modified positive definite strain pseudo energy Ψ that will be discussed later on. The second term in (32), that is required to ensure the consistency of the total pseudo potential energy in (28), can be evaluated similarly to the consistency term (31), while for the first term, depending on the displacement field u v | e at the damage state d, the novel procedure proposed in [44] is here exploited.
The key idea for computing the term U(u v ) is to consider a triangular finite element mesh, embedded in the virtual element one, made up of n int triangles that share all the nodes with the virtual elements. Each i-th virtual element with n nodes is divided in n − 2 triangles as schematically shown in Figure 1. is to avoid computing the solution at interior points of the element also for the stabilization term involving full quadrature. The linear ansatz adopted for the displacement field implies a constant deformation measure within each triangle Ω i T of the inscribed mesh. The pseudo potential U(u v | e , d), related to the virtual element e, can be thus obtained by assembling the contributions of all the internal triangular elements where u v | T is the part of the nodal displacements in a virtual element pertaining to the triangle Ω i T . For damaged materials the modified pseudo potential is formulated as where the strain tensor (•) is in turn u v | T or Πu v | e , and D(d) is a modified damaged isotropic constitutive tensor characterized by the Lamé parameters λ and µ, that are evaluated on the basis of the following modified values of Young's modulus E and Poisson ratio ν with E being the Young's modulus of the actual material and d the damage variable. The value of the Poisson ratio is kept constant since it does not influence the convergence behaviour of the element.
Concerning the algorithmic treatment of the damage model, the simple solution algorithm reported in [33] is here extended to the VEM. Assuming that the integration of the evolution equation is performed within the time step ∆t = t − t n , the constitutive equations, briefly summarized in Box 1, are solved at the element level. In what follows quantities without an index represent the current time and quantities with index n are related to the time t n .
In the particular case of nonlocal regularization, the value of τ e is replaced by its nonlocal counterpartτ e evaluated, according to equations (16) and (17), by averaging the stress values pertaining to virtual elements that fall within the representative circular volume surrounding the given element. The actual value of the element damage d e is, thus, evaluated accordingly. Note that the current element damage is used both for the calculation of the consistency terms and for the stabilization terms of the pseudo potential energy.

Box 1: Damage evaluation at the element level
Starting from (28), the global equilibrium is achieved from its minimization, i.e. ∂U(u v , d)/∂u = 0. Adopting a Newton-Raphson algorithm, the following standard expression for the linearised system is obtained The current element residual and tangent matrix are R(u v | e , d e ) = R c (u v | e , d e )+R stab (u v | e , d e ) and K(u v | e , d e ) = K c (u v | e , d e ) + K stab (u v | e , d e ). The residual vector and the tangent matrix, pertaining both to the consistency and to the stabilization term, are derived by exploiting the symbolic tool ACEGEN. The element residual vector for the consistency term results as where the formalism d e =C means that, within the process of automatic differentiation performed by ACEGEN, the value of damage d e is kept constant, and Πu v | e depends on u v | e . Referring to the stabilization term, in turn, it reads where the pseudo potential U(u v | e , d e ) pertains to the triangular finite element mesh (equation (33)), while Ψ(ε(Πu v | e ), d e ) is the modified positive definite pseudo energy directly evaluated on the virtual element mesh. Analogously, the contributions to the tangent matrix result as

Numerical Examples
This Section is devoted to some illustrative examples aimed at investigating the performances of the virtual elements with the isotropic damage law, adopting either local or non-local regularization techniques. Comparison with existing finite element formulations will help to highlight the specific features of the VEM.
Two simple examples are first proposed. The first reproduces a constant stress field, and no regularization is required. In the second example a localization is induced by a nonuniform uni-axial tensile stress field. The last two examples reproduce, instead, experimental results available in literature, related to a three-point bending test on a concrete beam and to a splitting test of a granite cylindrical specimen, respectively. All computations are performed by using a Newton-Raphson algorithm with load stepping. The meshes used were automatically generated by the meshing tools in Mathematica. Finally, note that the mesh refinement, when adopted and unless otherwise specified, is uniform in the sense that the finer meshes are included into the coarser meshes both for regular and distorted meshes.

Uniform damage
In the first example we consider a square 1x1mm 2 specimen undergoing an ideal constant stress tension-compression path. Plane strain conditions are assumed. Table 1 collects the material parameter adopted. In this case no localization appears and the solution is exact and independent on the mesh size, see [32]. The capability of the Virtual Elements for correctly reproducing the constitutive response

Tension specimen test
As second example, a classical test proposed by [32] is considered. The specimen, 20 cm thick, is shown in Figure 4 and the material parameters are reported in Table 2. Due to the specimen shape, the central neck induces localization. Also in this case different meshes are adopted, including distorted eight node quadrilaterals and Voronoi elements.
In this example we adopt both the regularization approaches recalled in Section 2.1. First, the local model assuming a linear damage law with l 2 c is exploited. In this case, the response shown in [32] is recovered first adopting an 8 nodes VEM element, with the regular meshes B1, B2, B3 and B4 , shown in Figure 5. As reported in the reference paper [32], the damage is found to localize within the first row of elements, close to the neck. The convergence trend of the global curves is perfectly matched and the finest mesh correctly reproduce the expected behaviour, as reported in Figure 6. As a second step, distorted meshesB5 and B6, including non convex shapes, depicted in Figure 7 and corresponding to meshes B3 and B4, are tested. The global curves are very close to the corresponding ones in Figure 6 and are almost undistinguishable in the plot. The second regularization approach is the non local damage for which we adopt an internal length l 3 c = 1 cm. In this case, in order to exploit the nonlocal regularization technique, finer meshes are required. We here consider centroidal Voronoi tessellations, in which the seeds coincide with centroids of the resulting Voronoi cells. In Figure 8 the global curves obtained adopting, respectively, meshes of 2000 (mesh V1) and 8000 (mesh V2) elements are shown. A nearly perfect match is observed between the two curved paths, that deviate from the linear curves obtained with local regularization. A reasonable explanation of this is suggested by the damage distributions within the specimen. It emerges that, in this case, the damage onsets in the central part of the specimen, as expected, and develops within a curved region, as shown in Figure 9 where damage plots, both for V1 and V2, correspond to the final state of the loading process. The two damage distributions are spread in the same region of the specimen, irrespective of the mesh density, thus remarkably confirming the effectiveness of the regularization technique.

Three-point bending test
As a first structural example, we simulate a three-point bending test of a concrete beam that is not reinforced. The experimental test proposed in [35] generates mode I fracture. A 2000x200x50 mm 3 single-edged notched-beam having a central notch 10 mm wide and 100 mm deep is considered. The material properties, available from the experimental work [35] are adopted in the numerical simulations, see Table 3. Plane strain conditions are assumed. In what follows, the global curves (total vertical reaction versus vertical displacement) ob-  larization technique, recalled in Section 2.1 with l 1 c , with experimental results. The same discretization with 5831 elements, shown in Figure 10, is adopted both for finite elements and virtual elements. In Figure 11 the global response curves are reported. The comparison with the experimental result (gray filled area) shows a good agreement. For both FEM and VEM approaches, the initial stiffness is well reproduced, together with the peak value and the post-peak softening trend. The two numeric curves show slight discrepancies in the pre-peak branch, while more pronounced differences take place in the post-peak branch. These circumstances are likely attributable to the different interpolation orders of the elements.  As a second investigation, we adopt both distorted and Voronoi meshes, shown in Figure  12, and the non local regularization technique, recalled in Section 2.1 assuming l 3 c = 6 mm. Figure 13 collects global response curves referred to two centroidal Voronoi tessellations with 3127 virtual elements (mesh V1) and 7216 virtual elements(mesh V2), respectively. A good agreement is found between the numerical and experimental results (gray filled area). In this case, both numerical solutions slightly underestimate the initial stiffness, but satisfactorily reproduce the peak value and the post peak behaviour. A slightly oscillating behaviour is observed in the post peak zone. Such response is reasonably related to the Voronoi discretization resulting in a more unstructured mesh. Finally, in Figure 14 zoomed views around the notch of the scalar damage distributions are shown corresponding to assigned displacement A, B and C of the global curves ( Figure 13) for the meshes V1 and V2. As expected, the damage onsets on the notch and gradually spreads symmetrically around it.  Figure 13 for Voronoi meshes V1 and V2.

Brazilian Test
As last experimental example, a splitting test of a cylindrical granite specimen described in [37,39,38] is considered. In Figure 15a a scheme of the geometry and the load configuration is presented, and in Table 4 the material parameters used in the numerical simulations are listed. The vertical load is applied through two bearing strips and two strain gauges (L 0 = 0.8 D, with D being the diameter of the specimen and L 0 the gauge length of the extensometer) are fixed on the opposite faces of the cylinder to measure the diametric deformations along the axis orthogonal to the loading direction. The diameter of the specimen is D = 120mm, the length of the bearing strips is b = 19.2 mm and the thickness is s = 30 mm. The outcomes of the experimental work show that the cracking process onsets in the center of the specimen and a first single crack vertically grows along the loading direction. At the end of this first stage the specimen is split in two halves that start to resist the load separately. A subsequent mechanism is then observed, involving the formation of secondary cracks. In general this type of failure mechanism can induce nearly vertical crack development at both sides of the bearing strips.
In the actual experimental case, a single crack takes place at the right bottom side of the bearing strips, propagating in subvertical direction. See Figure 15b where the observed cracks are highlighted in red. The direct consequence of this two-steps rupture mechanism results in the presence of two peak loads in the global response curve. The numerical tests are performed in plane strain state, using centroidal Voronoi virtual elements. The geometry is discretized with about 2000 Voronoi elements. The non-local damage model is adopted with a characteristic length l 3 c = 3 mm, that is 3d g with d g the estimated maximum size of grains in granite. In Figure 16 the comparison between the experimental and the numerical global curves load vs. diametral deformation shows that the highly fragile behaviour is well predicted by the numerical simulations. The vertical load P is normalized to the first peak load P max , as well as the relative horizontal displacement ∆w is normalized to the respective value ∆w( P max ), as proposed in [37]. The peak load estimated in the numerical analysis is 64.5 kN that is in good agreement with the value specified by the authors of the experimental work, i.e. P u = f t πsD 2[1−(b/D) 2 ] 3/2 65 kN. Also the post peak behaviour is predicted quite satisfactorily: the magnitude of the load drop is comparable with the experimental data. In Figure 17   Since the Voronoi elements are randomly generated, the analysis has been repeated for different meshes, characterized by about the same number of elements, to make sure that the damage path and the associated global response are not affected by mesh dependency. The same mechanism of rupture in the splitting test is recovered in all the simulations and differences of about 5% are found for the peak value.

Final Remarks
A virtual element method for the modelling of isotropic damage, in the framework of linearised kinematics, has been presented. The focus is on 2D problems and a loworder formulation is adopted, based on the novel definition of a stabilization term obtained through an "ad-hoc" modification of an idea previously proposed in [44] for nonlinear elasticity and extended in [43] to finite strain plasticity. Both local and nonlocal regularization techniques of the constitutive model are taken into account to avoid spurious mesh dependencies in the presence of softening behaviour. Numerical simulations reproducing either theoretical or experimental examples provide encouraging results. The VEM method confirms its effectiveness in dealing with arbitrary geometries of the elements also in the considered nonlinear regime. Extensions to higher order formulations will enable a richer description of the damaging within each virtual element.