A Virtual Element Method for elastic and inelastic problems on polytope meshes

A Virtual Element Method for elastic and inelastic problems on polytope meshes

L. Beirão da Veiga Dipartimento di Matematica, Università degli Studi di Milano, Italy, E-mail: lourenco.beirao@unimi.it C. Lovadina Dipartimento di Matematica, Università di Pavia, Italy, E-mail: carlo.lovadina@unipv.it D. Mora GIMNAP, Departamento de Matemática, Universidad del Bío-Bío, Casilla 5-C, Concepción, Chile, and Centro de Investigación en Ingeniería Matemática (CI2MA), Universidad de Concepción, Concepción, Chile, E-mail: dmora@ubiobio.cl

We present a Virtual Element Method (VEM) for possibly nonlinear elastic and inelastic problems, mainly focusing on a small deformation regime. The numerical scheme is based on a low-order approximation of the displacement field, as well as a suitable treatment of the displacement gradient. The proposed method allows for general polygonal and polyhedral meshes, it is efficient in terms of number of applications of the constitutive law, and it can make use of any standard black-box constitutive law algorithm. Some theoretical results have been developed for the elastic case. Several numerical results within the 2D setting are presented, and a brief discussion on the extension to large deformation problems is included.

1 Introduction

The Virtual Element Method (VEM), introduced in [2], is a recent generalization of the Finite Element Method which is characterized by the capability of dealing with very general polygonal/polyhedral meshes and the possibility to easily implement highly regular discrete spaces [11, 6]. Indeed, by avoiding the explicit construction of the local basis functions, the VEM can easily handle general polygons/polyhedrons without complex integrations on the element (see [4] for details on the coding aspects of the method). The interest in numerical methods that can make use of general polytopal meshes has recently undergone a significant growth in the mathematical and engineering literature. Among the large number of papers, we cite as a minimal sample [10, 5, 28, 19, 12, 30, 31, 17, 18, 16]. Indeed, polytopal meshes can be very useful for a wide range of reasons, including meshing of the domain (such as cracks) and data (such as inclusions) features, automatic use of hanging nodes, use of moving meshes, adaptivity.

In the framework of Structural Mechanics, recent applications of Polygonal Finite Element Methods, which is a different technology employing direct integration of complex non-polynomial functions, have shed light on some very interesting advantages of using general polygons to mesh the computational domain. This include, for instance, the greater robustness to mesh distortion [13], a reduced mesh sensitivity of solutions in topology optimization [12, 20], better handling of contact problems [7] and crack propagation [24]. Unfortunately, Polygonal Finite Elements suffer from some serious drawbacks, such as the strong difficulties in the three dimensional case (polyhedrons) and in the use of non convex elements. On the contrary, the VEM is free from the above-mentioned troubles, and thus it represents a very promising approach for Computational Structural Mechanics problems.

Aim of the present paper is to initiate the investigation on the VEM when applied to non-linear elastic and inelastic problems in small deformations. More precisely, we mainly focus on the following cases: 1) non-linear elastic constitutive laws in a small deformation regime which, however, pertain to stable materials; 2) inelastic constitutive laws in a small deformation regime as they arise, for instance, in classical plasticity problems. We remark that we are not going to consider here situations with internal constraints, such as incompressibility, which require additional peculiar numerical treatment. Virtual elements for the linear elasticity problem where introduced in [3, 21]. The scheme in the present paper is one of the very first developements of the VEM technology for nonlinear problems, and it is structured in such a way that a general non linear constitutive law can be automatically included. Indeed, on every element of the mesh the constitutive law needs only to be applied once (similarly to what happens in one-point Gauss quadrature scheme) and the constitutive law algorithm can be independently embedded as a self-standing black-box, as in common engineering FEM schemes. Therefore, in addition to the advantage of handling general polygons/polyhedra, the present method is computationally efficient, in the sense that the constitutive law needs to be applied only once per element at every iteration step. The risk of ensuing hourglass modes is avoided by using an evolution of the standard VEM stabilization procedure used in linear problems. However, we highlight that the proposed method is described for general -dimensional problems (), but the performed numerical experiments are confined to the two dimensional setting.

A brief outline of the paper is as follows. In Section 2 we describe the continuous problems we are interested in. In particular, we distinguish between the elastic, possibly non-linear, case (Section 2.1), and the general inelastic case (Section 2.2). Section 3 deals with the VEM discretization. After having introduced the approximation spaces and the necessary projection operators (Section 3.1), we detail the discrete problems for the elastic case in Section 3.2, and for the inelastic case in Section 3.3. In Section 4, combining ideas and techniques from [14] and [2], we provide some theoretical results concerning the convergence of the proposed scheme in the elastic situation. We remark that our analysis is confined to cases where the non-linear costitutive law fulfills suitable continuity and stability properties, as stated at the beginning of the Section. Section 5 presents several numerical examples which asses the actual behaviour of the proposed scheme. In Sections 5.1 and 5.2 we consider non-linear elastic cases, while in Section 5.3 a von Mises plasticity problem with hardening is detailed. Furthermore, an initial brief discussion about a possible extension to large deformation problem is included (Section 5.4). Finally, we draw some conclusion in Section 6.

Throughout the paper, we will make use of standard notations regarding Sobolev spaces, norms and seminorms, see [8], for instance. In addition, will denote a constant independent of the meshsize, not necessarily the same at each occurrence. Finally, given two real quantities and , we will write to mean that there exists such that .

2 The continuous problems

In the present section we describe the problem considered in this paper. Although the elastic case could be considered as a particular instance of the inelastic case, we prefer to keep the presentation of the two problems separate. This will allow us a clearer presentation of the ideas of the virtual element scheme in the following section.

2.1 The elastic case

We consider an elastic body () clamped on part of the boundary and subjected to a body load . We are interested, assuming a regime of small deformations, in finding the displacement of the deformed body.

We are given a constitutive law for the material at every point , relating strains to stresses , through the function


where represents the gradient of the displacement .

Given the law (1), the deformation problem reads


where denotes the unit outward normal to .

Let now denote the space of admissible displacements and the space of its variations; both spaces will, in particular, satisfy the homogeneous Dirichlet boundary condition on . The variational formulation of the elastic deformation problem reads

Remark 2.1.

The generalization of the results of the present paper to other type of loadings (for instance in the presence of boundary forces) and boundary conditions (for instance in the presence of enforced displacements) is trivial. Our choice in (2) allows to keep the exposition shorter.

2.2 The inelastic case

We assume a small deformation regime and restrict ourselves to rate independent inelasticity. We consider a material body () clamped on part of the boundary and subjected to a body load depending also on a pseudo-time variable . The interested reader can find more details in [23, 27], for instance. We are interested in finding the displacement of the deformed body at a given final time .

We are given an inelastic constitutive law for the material, relating strains to stresses , through the function


where the vector contains all history variables at the point .

The above rule is to be coupled with an evolution law for the history variables in time


where, as usual, a dot above a function stands for a pseudo-time derivative. Since we consider a quasi-static problem, at each time instant the stresses and displacements must satisfy the equilibrium and boundary conditions in (2).

We here avoid to write a rigorous variational formulation for the problem above, and limit ourselves to the minimal setting that will be needed to introduce the associated discrete problem. As in the elastic case, let denote the space of admissible displacements and the space of its variations. Then, assuming an initial value for the history variables, the quasi-static inelastic deformation problem can be written as


where the displacements and history variables are sufficiently regular in time and must satisfy the evolution law (5) almost everywhere.

3 The virtual element approximation

In the present section we introduce the virtual element discretization of problems (3) and (6). In what follows, given any subset of () and , we denote by (respectively ) the scalar (respectively vector with components) polynomials of degree up to on .

3.1 The virtual spaces and operators

We consider a mesh for the domain , made of general polygonal/polyhedral conforming elements. For the time being, we only assume that such mesh is compatible with the boundary conditions, i.e. that is union of faces (edges) of the mesh. We denote by the generic element of the mesh and by the generic face (or edge if ). The symbols and will represent, respectively, diameter and volume (or area) of the element . As usual, will indicate the maximum element size.

We start by introducing the discrete virtual space for displacements, that is essentially the same as in [3]. We first consider the two dimensional case. Given any , let the local virtual space


where denotes the component-wise Laplace operator. The space is a space of harmonic functions that on the boundary of the element are piecewise linear (edge by edge) and continuous. Such space is virtual in the sense that is well defined but not known explicitly inside the element.

Note that ; in the case of a triangular element, we recover exactly the standard space. It is easy to check [3] that a set of degrees of freedom for the space is simply given by the collection of the vertex values:

Once the above degrees of freedom values are given, since is linear on each edge, the value of on the boundary is completely determined. Therefore, an integration by parts allows to compute the integral average of the gradient


with indicating the outward unit normal at each edge .

We now define the virtual local spaces for the three dimensional case. Given a polyhedron , any face is now a polygon. We denote by the virtual bi-dimensional space (7) on the polygon adjusted with three components:


where the symbol represents the generic edge of the polyhedron and denotes the planar laplacian on . We then define


The space is a space of harmonic functions that on the boundary of the element are continuous and, on each face, functions of . Note that, as a consequence, the functions of are linear on each edge of the polyhedron.

Again we note that ; in the case of a tetrahedral element, we recover exactly the standard space. It is easy to check that a set of degrees of freedom for the space is again given by

An integration by parts exactly as in (8) allows to compute, for all the integral average of the gradient, provided one is able to compute the face integrals for all and . Such face integrals can be easily computed by introducing the virtual space modification proposed in [1], that we do not detail here. The result is

where the scalars are the weights of any integration rule on the face that is exact for linear functions.

Once the local virtual spaces are defined, all that follows holds identically in two and three dimensions. We can now present the global virtual space

A set of degrees of freedom for is given by all pointwise values of on all vertices of , excluding the vertices on (where the value vanishes).

In the following, we will denote by the tensor valued projection operator on the space of piecewise constant functions and by its restriction to the generic element . More precisely, for any , we have with the local operators defined as


We have the following important remark, which is a direct consequence of (8).

Remark 3.1.

For all functions and all elements , the operators are explicitly computable.

We moreover introduce a second projection operator , defined on as follows. For any , we have with the local operators defined as


for all in . Note that, by definition, is a (discontinuous) piecewise linear function on . On each element , is the unique linear function such that:

  1. its (constant) gradient equals the mean value over of the function ;

  2. its vertex value average equals the vertex value average of .

We notice that the second condition in (12) is only to fix the constant part of on each element. Recalling Remark 3.1, it is immediate to check that the operator is explicitly computable.

3.2 The elastic case

The main missing step is to introduce the local forms that will be used in the discrete variational formulation. We assume that the constitutive law (1) is piecewise constant with respect to the mesh . Therefore, instead of , we will write to represent the constitutive law on , and . In addition, foe every pair and , we introduce the bilinear forms and as:


Therefore, it holds


and, recalling (3), the elastic problem can be written as


We now consider, for all and all , the following preliminary form


where the identity above follows since all the involved functions are constant on the element. The above form is -consistent, in the sense that it recovers exactly the original form whenever the first entry is a linear polynomial. Indeed, it follows from (11) and (12) that


However, unless the elements are triangular/tetrahedral, the form has a non-physical kernel that may lead to spurious modes in the solution. We therefore follow the idea proposed initially in [2] and introduce the discrete bilinear form


As discussed in [3, 2], under suitable mesh regularity assumptions detailed in Section 4, there exist positive constants independent of the element such that


for all with . In other words, on the orthogonal complement of with respect to , the bilinear form behaves as the local energy of a linearly elastic body with unitary material constants and is thus suitable to stabilize form in such case. In order to take into account different material constants and also nonlinear materials, the form needs to be multiplied by a positive constant that may depend on the discrete solution.

We therefore introduce the following local virtual form on . For all and all


where the stabilizing parameter depends on the additional entry . We remark that the bilinear form is still -consistent. This follows from (17) and the observation that for every . The choice that we here propose for the parameter is


with representing any norm on the fourth order tensor space, for instance the maximum of the absolute values of all the entries, see Remark 3.2.

We present also the global form


Given , a possible virtual discretization of Problem (3) is


Above, the load approximation term

is a vertex-based quadrature rule with weights chosen to provide the exact integral on when applied to linear functions. Furthermore, a reasonable choice for could be .

We instead propose a modification of (23), that is more practical from the implementation viewpoint. We assume the usual incremental loading procedure for the solution of the nonlinear discrete problem: given a positive integer , let the partial loadings for all . Then, given the initial displacement (for instance the zero function), one applies for the iterative procedure


The final solution is then . The nonlinear problems above can be solved with the Newton scheme. Note that, since the stability constants (see (20)) are computed by using , the tangent matrix in the Newton iterations turns out to be simpler. Since is typically taken large (at least 10, but often much more) the effect of such modification is not detrimental for the discrete approximation; the constants are only used as scaling parameters and do not enter the accuracy of the algorithm.

We close the section with some observations regarding the local forms used in the scheme. First, we recall that the proposed forms are -consistent, in the sense that for all , we have:


Identity (25) is a fundamental condition for approximation and, in particular, guarantees the satisfaction of the patch test. Moreover, such forms are explicitly computable for any polygonal/polyhedral element (even non-convex). Finally, the constitutive law needs to be computed only once per element and thus the method, from this point of view, is as cheap as finite elements with one point gauss integration rule. This observation has an even bigger impact in the inelastic case, where the constitutive laws are typically more expensive to compute.

Remark 3.2.

The motivation for choice (21) and (24) is to better mimic the stability properties of the material for the current displacement. For materials in which the stress-strain incremental relation does not depend too strongly on the value of the current displacement, the constants can be taken as independent of . For instance, a scaling directly proportional to the local material constants could be used. On the other hand, the choice proposed in (21) and (24) give good results for a wider range of materials. Examples and investigations in this direction can be found in Section 5.

3.3 The inelastic case

We start by introducing a sub-division of the “time” interval into smaller intervals for , where for simplicity we assume that . We will denote the partial loadings by for all .

We assume, as in standard engineering procedures, a constitutive algorithm that is an approximation of the constitutive and evolution laws (4), (5). In Finite Element analysis, this pointwise algorithm can be coded independently from the global FE construction and can be regarded as a “black-box” procedure that is applied at every Gauss point and at every iteration step. In the present Virtual Element method, we want to keep the same approach; in other words, our scheme will be compatible with any black-box constitutive algorithm that follows in the general setting below and that can be imported from other independent sources.

We assume that the constitutive law is piecewise constant with respect to the mesh . Let represent the constitutive algorithm for the element . For any , given a value for the displacement gradient at time , a value for the history variables at time and a tentative value for the displacement gradient at time , the algorithm computes the stresses (and updates the history variables) at time . We thus write the computed stress as

As part of the approximation procedure of our method, we assume that the history variables are piecewise constant with respect to the mesh, and therefore write to represent the value assumed on the element at time . Consistently, will represent the collection of all .

In our scheme, instead of applying the constitutive algorithm at Gauss points, we make use of the projections introduced in the previous sections and of the same stabilization as in the elastic case. The Virtual Element scheme reads, for :


where the form

with, for all ,

Here above, the bilinear form and the scalar are calculated as already shown in (18) and (21), respectively. Note that, as already mentioned in Section 3.2, the constitutive algorithm needs to be applied only once per element.

4 Theoretical results

We here develop an error analysis for the method described in Section 3.2, under some additional hypotheses on the function . More precisely, we assume that the following properties are satisfied.

Hypotheses (RPC)

  • The function belongs to for every ;

  • for every , the differential satisfies

    1. there exists such that

    2. there exists such that


We moreover explicit here the shape regularity conditions that are needed for the theoretical results of the present paper. We assume that there exists a positive constant such that all the elements of the mesh sequence are star shaped with respect to a ball of radius and that all the edges of have length .

Lemma 4.1.

Let the bilinear forms , , and be defined by (13), (20) and (22). Suppose that the Hypotheses (RPC) introduced above are satisified. Then, it holds


We first note that (27) and (28), together with (21), imply the existence of positive constants and such that


Step (i): proof of (29). From (27), we deduce that


Therefore, for every we have


by which


For every , we have (see (20))


We now notice that (see (16))


First using (33) with and , then recalling (12) we get


In addition, we have, using (32) and (19):


Combining (36) with (38) and (39), we infer


Summing up over all the elements, we get (29):


Step (ii): proof of (30) and (31). From (28), we deduce that


from which we easily get (30):


We now notice that (see (16))


Using (42), identity (44) yields


To continue, since is a bilinear form and using continuity arguments, we have for every (see (19))


From (20), using (45) and (46), we deduce (31). ∎

Theorem 4.1.

Let be the solution of Problem (3). Given any , let be the solution of Problem (23):


For any and such that , it holds:


where denotes the -scalar product.


Given , we set . For every such that , using (29) we have


Since (25) implies , from (49) we get