Virtual Element Formulation For Finite Strain Elastodynamics

This work provides an efficient virtual element scheme for the modeling of nonlinear elastodynamics undergoing large deformations. The virtual element method (VEM) has been applied to various engineering problems such as elasto-plasticity, multiphysics, damage and fracture mechanics. This work focuses on the extension of VEM towards dynamic applications. Within this framework, we employ low-order ansatz functions in one, two and three dimensions that having arbitrary convex or concave polygonal elements. The formulations considered in this contribution are based on minimization of potential function for both the static and the dynamic behavior. While the stiffness-matrix needs a suitable stabilization, the mass-matrix can be calculated using only the projection part. For the implicit time integration scheme, Newmark-Method is used. To show the performance of the method, various numerical examples in 1D, 2D and 3D are presented.


Introduction
The virtual element method (VEM) can be seen as an extension of the classical finite element method (FEM) based on Galerkin projection. It allows meshes with highly irregular shaped elements, including non-convex shapes, as outlined in [1,2]. This gives more flexibility and new possibilities to geometry discretization in solid-and fluid-mechanics. The large number of positive properties of VEM increases the variety of possible applications in engineering and science. Recent works on virtual elements have been devoted to linear elastic deformations in [3,4,5], contact problems in [6], elasto-plastic deformations in [7,8,9], anisotropic materials in [10,11,12], curvilinear virtual elements for 2D solid mechanics applications in [13], hyperelastic materials at finite deformations in [14,15], crack-propagation for 2D elastic solids at small strains in [16] and phase-field modeling of brittle and ductile fracture in [17,18]. Despite the fact that dynamic behavior has a strong influence on the mechanical properties and the prediction of their real response, most of the investigations introduced above are only done for static problems so far. Thus the element mass-matrix is needed to be calculated. In this regard, [19] proposed a virtual element method for linear elastodynamics problems. However their formulations are restricted to small strain setting, hence it is not appropriate for large deformations. This has motivated the presented contribution to extend the application of VEM from the static to the dynamic case in the finite deformation range. Typically the construction of a virtual element is divided into a projection step and a stabilization step. Within the projection step, a quantity ϕ h is replaced by its projection ϕ Π onto a polynomial space. Using this projected quantity in the weak formulation or energy functional yields to a rank-deficient structure which needs to be stabilized. In the second step, the stabilization term, which is a function of the difference ϕ h − ϕ Π between the original variable and the projected quantity needs to be evaluated. There are various possibilities to evaluate this stabilization term. To this end, Da Veiga et al. [20] proposed a stabilization term, where all integrations take place on the element boundaries. Wriggers et al. presented in [6] a novel stabilization technique, which was first described for finite elements in Nadler and Rubin [21], generalized in Boerner et al. [22] and simplified in Krysl [23] for the stabilization procedure of a mean-strain hexahedron. In this framework, the integration is carried out over a triangulated sub-mesh, which uses the same nodes as the original mesh. The method presented in this contribution based on the stabilization technique of [6]. In order to model the dynamic behavior of the body, we define a specific potential function, where the second derivative of it with respect to the global unknowns yields the mass-matrix. As a key advantage of this approach in comparison to [19], only the stiffness-matrix needs to be stabilized, whereas the mass-matrix is only computed using the projection part. Hence no stabilization is required for the mass-matrix as in the case of [19]. For the time integration scheme, we utilized the implicit Newmark method as documented in [24,25]. The structure of the presented work is as follows. In Section 2 the governing equations for nonlinear elastodynamics are outlined. Section 3 summarizes the virtual element formulation. It includes details on the computation of the element mass-matrix. To verify the proposed virtual element formulations, a various number of examples are demonstrated and discussed in Section 4. Section 5 briefly summarizes the work and gives some concluding remarks.

Governing equations for Nonlinear Elastodynamics
In this section we summarize the finite strain elasto-static formulation (see e.g. [26,27,28]) and supplement it by the dynamic effect. For that consider an elastic Body Ω ⊂ R 3 with boundary Γ. This boundary is decomposed into a non-overlapping Dirichlet Γ D and Neumann Γ N boundary conditions such that Γ D ∪ Γ N = Γ, see Figure 1. The position x of a material point in the current configuration is given by the deformation map where X is the position of a material point in the initial configuration and u(X, t) is the displacement. In the further course of this work we will skip the explicit specification of the dependence of variables on the initial configuration and time thus we will write: u = u(X, t). In order to transform quantities which are defined with respect to the deformed configuration to the reference configuration and vice versa, we define the deformation gradient F as where the gradient is evaluated with respect to the initial configuration X. We further define the right Cauchy-Green tensor C(u) with C = F T F as a strain measure and the Jacobian J(u) with J = det F as a volume map.
The solid Ω has to satisfy the balance of linear momentum where f are the body forces and P, S are the 1 st and 2 nd Piola-Kirchhoff stresses, respectively. The right side of the equation (3) 1 is taking the dynamic effects ρü into consideration. The Dirichlet and Neumann boundary conditions are defined by here N is the outward unit normal vector related to the initial configuration,ū represents the prescribed displacement on the Dirichlet boundary Γ D , andt depicts the surface traction at the Neumann boundary Γ N , as illustrated in Figure 1. The weak formulation of the elastodynamics problem in (3) then takes the form where δu is the test function of the displacement u. A homogeneous compressible isotropic elastic material is considered, here we use the Neo-Hookean strain energy function in terms of the bulk κ and shear µ modulus and the invariants of the right Cauchy-Green tensor I 1 = tr C, I 3 = det C.
With the above set of equations, the finite strain elastodynamic problem is well formulated. Next, we use the potential function as a starting point for the development of a discretization method. The static part of the potential is defined as whereas the dynamic part of the potential is the kinetic energy that describes inertial effects takes the form where ρ is the density of the solid.

Formulation of the Virtual Element Method
The main idea of the virtual element method is to use a Galerkin projection, which maps the primary fields (displacement in this work) to a specific polynomial ansatz space. Thus, the domain Ω can be discretized with non-overlapping polygon in 2D or polyhedral elements in 3D which do not need to have convex shapes [3]. Since VEM has no isoparametric mapping, the ansatz functions are given in terms of the coordinates X in the initial configuration. Here the ansatz functions for the virtual element are based on linear functions, therefore the nodes are placed at the element vertices.

VEM Ansatz
In general, for finite strains the deformation map ϕ = X + u has to be discretized. But as the coordinates X in the initial configuration are exactly known, we can reduce the discretization to the displacement field u = u i E i where E i are the basis vectors with respect to the initial configuration in the three-dimensional space i ∈ {1, 2, 3}. The central concept of the virtual element method relies on the split of the ansatz space u h into a projected part u Π and a remainder u h − u Π as For a linear ansatz, the projection u Π at element level takes for three-dimensional elements the form where a represents the twelve unknown virtual parameter a = a i which have to be determined. Instead of using the polynomial N Π in equation (11) as the interpolation function, a scaled ansatz can be used, for details see e.g. [5]. Furthermore, the projection u Π has to fulfill the orthogonality condition, as discussed in the work of [29]. Hence ∇u Π is computed through the Galerkin projection as where p is a polynomial function which has been chosen similarly to the projection u Π , see (11). Since we use linear ansatz functions, ∇p and ∇u Π are constant and can be shifted out of the integral as Applying integration by parts to (13), we obtain at the element level. Here N denotes the normal vector on the reference boundary Γ e of the domain Ω e , which belongs to a virtual element e. Element quantities, which have constant values within the entire element e, are denoted by | e . With this simplification the projection u Π is defined as outlined in [15]. By employing the linear ansatz space, the left hand side of (14) takes the simple form In the 2D case, the right hand side of (14) is evaluated along the edges. As the displacements are known at the boundary, which are straight line segments, a linear ansatz for the displacements is used, see [15]. However, in the 3D case, the element boundary consists of polygonal faces. Therefore the evaluation of the integral in (14) is not straight forward, unless an appropriate ansatz is found. For the evaluation, there are two possible methods available. The first one is presented in [4], where the faces are subdivided in quadrilateral elements where the corners of the quadrilateral elements have certain positions. Finally the evaluation of the integral is carried out on those quadrilateral elements. An alternative option is to subdivide the element faces into 3 noded triangles. The integration is then carried out over the triangles of the polygonal faces by using the standard ansatz function for a linear triangle and Gauss integration: as outlined in [2]. Here u T h denotes the linear ansatz for the displacements at each triangle of the polygonal faces. u I is a list which contains the three nodal displacement vectors of the triangle T . ξ and η are the local dimensionless coordinates at the element level. The local nodes of T and the outward normal vector N i are visualized in Figure 2. Finally we are able to compute the right hand side of (14). Using (17), the integral in (14) takes the form: Here n f is the number of element faces. For an integration over triangles with linear shape functions (16) one point quadrature with n g = 1 Gauss point and w g = 1/2 Gauss weight is sufficient. N ζ is the Jacobian of transformation from the reference to the initial configuration. g denotes quantities which are evaluated at the Gauss point with the local coordinates ξ = 1/3 and η = 1/3. The normal vector N and the Jacobian of the isoparametric mapping N ζ are evaluated as follows: All quantities are related to the initial configuration.
Comparing (15) and (18) we obtain the unknown virtual parameters a| i∈(4,...,12) by inspection, for further details see e.g. [15]. Since only the gradient of the projection ∇u Π | e is needed to define the strain energy function of the static part, we actually do not have to compute the virtual parameters which are related to the constant parts. However, we will later see that the projected displacements u Π , including the virtual parameters which are related to the constant parts, are needed to construct the mass-matrix. Thus we need to supplement our formulation by a further condition to obtain the constants a| i∈ (1,2,3) . For this purpose we adopt the condition (see for example [3]) that the sum of the nodal values of u h and of its projection u Π are equal over the entire domain. This yields for each virtual element Ω e to the following condition where n V is the number of boundary nodes and X I the initial coordinates of the nodal point I. The sum includes all boundary nodes n V . By substituting (11) in (22) we can express the missing three parameters in terms of the nodal displacements and the already known virtual parameters a| i∈(4,...,12) : Finally with equation (14) and (23) the ansatz function u Π of the virtual element is completely defined in terms of the element nodal displacements u e :

Implicit time integration
For the time integration scheme, we use the implicit Newmark method outlined in [24,25]. The equations for the velocityu = v and accelerationü at time step t n+1 are given as where the Newmark parameters are chosen as ζ = 1/4 and γ = 1/2.

Construction of the element mass-matrix for VEM
For the calculation of the mass-matrix, we start from the dynamic potential in equation (9). From (14) and (24) we obtain the virtual parameters a in terms of the unknown displacements u e as a =Π ∇ u e , whereΠ ∇ is a constant matrix. Hence the projection u Π can be rewritten as where H is the matrix representation of the ansatz functions N Π . It is defined in the three-dimensional case as The matrix representation of the projection operatorΠ ∇ in (28) can be derived with the help of (18) and (23), as discussed in [15].
To compute a good approximation for the mass-matrix, we use the projection u Π which is now defined by the orthogonality condition: Usually one has to compute a new projection which is not the same as for the construction of the stiffness matrix. Nevertheless, it has been shown in [29,30] that it is possible to use the same projection operatorΠ ∇ for k ≤ 2 where k is the order of the ansatz functions and it does not lead to an error. Thus we use the same projection operatorΠ ∇ , which has been used for the static part asü For the unknown accelerationsü h we are using the same split as for the unknown displacements in equation (10), yields For the construction of the elastodynamic virtual element, it is computationally advantageous to employ the software tool AceGen, see [27]. It provides the most efficient element routines when a potential formulation is used. Thus we construct a specific pseudo-potential for inertia term where the first variation has to be performed for fixedü. The total potential function is now defined as which depends on both the static and inertia parts. The variation of (34) yields exactly the weak form of (3) when considering the nonlinear dependency of the 2 nd Piola-Kirchhoff stress S on the displacement u. Therefore the usage of the pseudo potential is equivalent to using the weak form (6) directly.
If we insert both equations (10) and (32) for the displacements and the accelerations in (33), we obtain: where coupled terms vanish due to the orthogonality condition (30). The first term in (35) is the consistency part, whereas the second term is the stabilization part. It is sufficient to use the consistency term for the construction of the mass-matrix, without any stabilization, when the problem is without any reaction term, as shown in [29,30]. To compute the mass-matrix from U dyn , we need to compute the first and second derivative with respect to the global unknowns. By utilizing the relationship between the projected values and the unknown values for the displacement and the accelerations, the following expression for the modified dynamic potential function results Hereby the virtual accelerations and displacements are constant over the entire domain of the element therefore they can be shifted out of the integral. Thus the integral yields This integral can be evaluated in different ways which will be explained in Section 3.4.1.
The first derivative of the dynamic potential U dyn is computed holding the acceleration u e constant: Before computing the second derivative, the Newmark method is used for the implicit time integration. With (26), the residual for the dynamic part follows as: The second derivative of U dyn leads then to the dynamic part of the tangent

Construction of the virtual element
As introduced in Section (3), the formulation of a virtual element undergoing large deformations is based on a split of the energy into a constant part and an associated stabilization term. The nodal degrees of freedom of an element are approximated with one interpolation function per coordinate direction in each element. Thus the consistency part does not lead to a stable formulation and a stabilization term is required. The idea of stabilizing the formulation is analogous to the stabilization of the classical finite elements with reduced integration, developed by [23]. For the construction of the virtual element method we start from the potential function (34). After summing up all element contributions for the n e virtual elements we obtain the following expression:

Consistency part
For the consistency part we use the projection u Π as introduced in Section 3.1 and therefore the first part of equation (41) for each element is given by The gradient of the projection ∇u Π | e is constant on the entire domain Ω e thus all kinematic quantities such as F| e = 1 + ∇u Π | e are constant as well. Hence the integration of the strain energy function can be simplified as: which is still nonlinear with respect to the unknown nodal degrees of freedom.
As already mentioned in Section 3.3, we can compute the third integral in (42) related to the dynamic part in different ways: 1. First possibility is to evaluate the integral at the centroid X c of the polygon in 2D and of the polyhedra in 3D. The displacement and the acceleration are then evaluated at the centroid and multiplied by the area of the element: 2. Another possibility is to introduce a sub-triangulation of the polygon or polyhedra and again use one point Gauss integration which yields an evaluation at the centroid X c | T of each triangle: Since the integral contains quadratic terms of X and Y , the two ways introduced above approximate the integral.
3. As a third option, we can compute the integral (37) exactly, using the nodal coordinates at the boundary, see [31,32]: where x 0 = x n V and y 0 = y n V .

Stabilization part
The consistency term is computable but yields to a rank deficient stiffness matrix and thus needs to be stabilized. The idea is to introduce a new positive definite energyÛ , with the help of which the stabilization term is redefined, as introduced in [15]: We further define a stabilization parameter β ∈ [0, 1] for the definition of the positive definite energyÛ as: Thus the stabilization term takes the form: Applying equation (49) in equation (41), the final form of the total potential energy function takes the form: The computation of the first term of equation (50) can be done as explained in Section 3.4.1. The second term U c (u h | e ) needs an approximation. An approach how to compute this part is introduced in [15]. The displacement field is approximated by introducing an internal mesh of 3 noded triangles in 2D or 4 noded tetrahedra in 3D with linear ansatz functions. The nodes of the generated submesh belong to the set of nodes in the virtual element, such that no additional nodes are needed. Based on the triangulated submesh, the displacement gradient is computed. The stabilization term U c (u h | e ) contains both the static U stat c (u h | e ) and dynamic part U dyn c (u h | e ). As the ansatz is linear, the gradient is constant and thus the integral for the static part can be simply evaluated at the centroid X c | T of each triangle n T , as sketched in Figure 2. For the dynamic part, the integral can also be evaluated at the centroid of each triangle as shown with elements VEM H2S-I and H2S-II in Section 4. For this case, we would use (45) but replace the projected quantities with the nodal values for the displacement u h and the acceleration u h . The nodal accelerations can be simply computed with the use of (25) and (26). The stabilization parameter β can be chosen freely. For β = 1 the total energy is calculated using only the stabilization part and thus the solution is purely related to the FEM results with three noded triangles in 2D or four noded tetrahedron in 3D. β = 0 yields to a rank deficient stiffness matrix. The choice for the stabilization parameter β was analyzed in [7,17] and it has been shown that the optimal value is in the range β ∈ [0.2, 0.6]. For our investigations we choose β = 0.4. Because of the coupling of VEM and FEM, we call this stabilization procedure mixed VEM-FEM-Stabilization. To obtain the total residual vector R e and total element tangent matrix K e , we compute the first and second derivative of the total energy U (u h ) with respect to the global unknowns u e as: For the separate computation of the mass-and the stiffness matrix, we split β into β stat and β dyn . Equation (51) and (52) become: (53) Here, β stat and β dyn are the stabilization parameters for the static and dynamic part, where β stat , β dyn ∈ [0, 1]. For the computation of the derivatives in equation (53) and (54), the Mathematica Package AceGen is used, see [27].

Numerical Examples
In this section we demonstrate the performance of the presented 3D virtual element formulation for dynamic problems at finite deformations. For comparison purposes results of the standard finite element method (FEM) are also included. The material parameters used in this work are the same for all examples and are given in Table 1.
In this contribution, the following mesh types with first order virtual element discretizations are introduced. Herein, the different ways to evaluate the integral in (37) are employed in the following element types: • VEM Q2S: A regular 2D mesh with 8 noded quadrilateral elements. The argument of the integral is evaluated at the centroid of the polygon and is multiplied by the area of the element. β dyn = 0, which means that the mass-matrix is computed using only the projection part. This represents the classical way as introduced in (44).
• VEM Q2S BI: A regular 2D mesh with 8 noded quadrilateral elements. The argument of the integral is computed exactly on the boundary with the moments of area (46). β dyn = 0, which means that the mass-matrix is computed using only the projection part.
• VEM VO: 2D/3D Voronoi cell mesh with arbitrary number of element nodes. The argument of the integral is evaluated at the centroid of the polygon and is multiplied by the area of the element. β dyn = 0, which means that the mass-matrix is computed using only the projection part.
• VEM VO BI: 2D/3D Voronoi cell mesh with arbitrary number of element nodes. The argument of the integral is computed exactly on the boundary with the moments of area (46). β dyn = 0, which means that the mass-matrix is computed using only the projection part.
• VEM H2S: A regular 3D mesh with 20 noded hexahedral elements. The argument of the integral is evaluated at the centroid of the polygon and is multiplied by the area of the element. β dyn = 0, which means that the mass-matrix is computed using only the projection part.
• VEM H2S-I: A regular 3D mesh with 20 noded hexahedral elements. The polyhedra is subdivided into tetrahedrons and the argument of the integral is evaluated at the centroid of each tetrahedron and is multiplied by the volume of the element. β dyn = 0, which means that the mass-matrix is computed using only the projection part. This procedure is same as in (45).
• VEM H2S-II: A regular 3D mesh with 20 noded hexahedral elements. The polyhedra is subdivided into tetrahedrons and the argument of the integral is evaluated at the centroid of each tetrahedron and is multiplied by the volume of the element. β dyn = 1, which means that the mass-matrix is computed using only the stabilization part.
For a representative comparison, the following finite element formulations were selected: • FEM T1: A regular 2D mesh with 3 noded triangular first order finite elements.
• FEM H1: A regular 3D mesh with 8 noded first order finite elements.
• FEM H2: A regular 3D mesh with 27 noded second order finite elements.
For the stabilization parameter of the static part we choose β stat = 0.4 for all the simulations. Unless otherwise specified, all the results for the dynamic part are obtained with β dyn = 0. In this case the mass-matrix is computed using only the projection part having no stabilization. Therefore, equations (53) and (54) simplify to:

Wave propagation in longitudinal beams
The first model problem is concerned with analyzing the wave propagation in longitudinal beams. The geometric setup and the loading conditions of the specimen are depicted in Figure 3. The height of the beam is chosen to be h = 0.3 mm and the length = 30 mm, where the degrees of freedom are fixed in longitudinal direction on the right side. The height of the specimen is much smaller than the length thus this problem is considered to be a one dimensional problem. The initial velocity is set to v 0 = 20 m/s. After a certain time, we observe the wave propagation through the elastic body. Next, the virtual element method will be compared with the finite element method and the analytical solution. The analytical solution in (58) is obtained by solving the wave equation (57). Figure 4a and 2v 0 c ω n 2 sin w n x c sin(w n t) with w n = 1 0, 95 In VEM Q2S the integral for the dynamic part in equation (37) is evaluated at the centroid of the element, hence this method seems to be sufficient.

Transversal Beam Vibration
The second benchmark test is concerned with analyzing the transversal vibration in beams. The geometric setup and the loading conditions of the beam are depicted in Figure 5a. The length of the bar is set to = 30 mm and the hight is h = 5 mm.
The force is applied transversal at the end of the specimen as shown in Figure 5b. The temporal course of the force is given by a half sine, where the maximum of the force is set to P max = 100 kN . The period T of the applied force is adjusted to the bending stiffness of the beam and defined as   of the dynamic part, we used different type of meshes which can be seen in Figure 6. The "animal" mesh ( Figure 6b) includes non convex elements. To see the effect of using nonconvex elements where the centroid of the element is outside of the element domain, we use a special mesh with elements like C's, where the centroid of the element is outside of the element domain (Figure 6c). Figure 7a and 7b show the displacement over time response in the center at x = /2 and at the end of the beam at x = . The finite element solution is computed for 1000 elements, whereas VEM results are obtained with 100 virtual elements. The comparison of the virtual elements Q2S and Q2S BI shows that it makes no difference for regular shaped meshes if the integral in equation (37) is evaluated approximately on the centroid of the element or exactly on the boundary using the moments of area. Furthermore, we can see that the displacements in the center of the beam are slightly higher than the finite element results. However the period fits very well compared with FEM results. In general the virtual element results are in a good agreement with the compared finite element results.
The comparison of the different meshes shows, that the C-mesh yields a higher deflection, compared to the other results. Nevertheless qualitatively the shape of the displacement over time response fits very well the finite element Q2 results and the virtual element Q2S results. Again, the evaluation of the integral at the centroid of the element compared to computing the integral at the boundary exactly using the moments of area does not affect the results.

Cook's membrane problem
The next example is the Cook's membrane problem in 2D. Here as well the virtual element performance will be compared with the finite element results. The geometrical setup and boundary conditions are demonstrated in Figure 9b. In this test a force driven scenario is applied at the right edge as a line load as depicted in Figure 9b. The force is applied as shown in 5b with P max = 10000 kN/mm. VEM VO mesh and regular VEM Q2S mesh are also plotted in 9a and 9c, respectively. The contour plots of the vertical displacement evolution for different deformation states {t = 0 s, 0.0001 s, 0.0002 s, 0.00035 s, 0.00055 s, 0.00065 s} are sketched in Figure 10. The   nonlinear behavior is clearly observed in the deformation process due to the dynamic effects at finite strains. Figure 11 shows a mesh refinement study with the element division of 2 N for N=1,2,3,4. For N=3 and higher the solution converges. A comparison with FEM depicts that the results are in a very good agreement.
This study shows that again, that the evaluation of the integral in (37) at the element centroid is absolutely sufficient to compute the mass-matrix.

3D Boundary value problems 4.2.1 Wave propagation in a bar
The previously introduced 2D model of a bar is here extended to the third dimension. The length of the bar is set to = 30mm and the height is equal to the width h = b = 5mm. We apply an initial velocity of v 0 = 20m/s in longitudinal direction. The virtual element results are obtained using 400 elements, were the finite element results were obtained with 4320 elements. In this example we compare the virtual elements H2S, H2S-I and H2S-II with the finite element H1 and the analytical solution which was obtained for the 1D case in equation (58). As already introduced before, the variable β dyn indicates how the mass-matrix is going to be evaluated. For β dyn = 0, the mass-matrix is calculated using only the projection part. Whereas for β dyn = 1 the computation of the mass-matrix is carried out using the stabilization part. Figure

Transversal vibration of a thick beam
In this benchmark test a 3D cantilever beam is investigated. The geometric setup and the loading conditions of the specimen are depicted in Figure 13. Here a line load is applied at the end of the beam with P max = 6 kN/mm. The temporal course of the force is again given by a half sine, as shown in 5b. Furthermore, similar to the 2D case, we set the beam length as = 30 mm with equal height and width as h = b = 5 mm. We compare the virtual element H1 with the finite elements H1 and H2. For this purpose a mesh refinement The maximum deformation state is sketched in Figure 14b, representing the deflection w. Here the nonlinear behavior is clearly observed due to the dynamic effects at finite strains. Figure 15 illustrates the displacement over time response at x = for the mesh refinement study. This response is plotted for the center of the cross section. We observe that both, VEM and FEM results are converging to the reference solution for increasing number of elements. Still there is a shift with increasing time, this is due to the less accuracy of VEM/FEM H1 element compared with the FEM H2 quadratic ansatz function.
Additionaly, we employ the virtual element VEM H2S computed with 256 elements in Figure 14a and compare it with the reference solution (FEM H2 with 3200 elements). It is interesting to note that despite the use of linear ansatz functions VEM H2S produces nearly the same results as the reference solution. This is due to the fact that the stabilization uses the bending modes. In conclusion, the presented formulation depicts very good results also in the 3D case.
This test also confirms that the evaluation of the integral in (37) for the computation of the mass-matrix only at the centroid of the polygon/polyhedra is absolutely enough to get satisfying results. (d) N=4 Figure 15 -Convergence Study -Displacement over time response for 3D Example -Transversal vibration of a thick beam. Element division 2 N , where N increases from (a) to (d).

Vibration of a thick plate
The last example is related to the vibration of a thick plate which is discretized using three-dimensional elements. The plate has a length = 30 mm, a thickness h = 5 mm and a width b = 30 mm as shown in Figure 16. We further set an initial condition such that the initial velocity of the plate is set to v 0 = 200m/s in the z-direction, see Figure  16. Figure 17 demonstrates the evolution of the displacement in the z-direction for different deformation states using the VEM-VO element. Herein a nonlinear response undergoing large deformation is observed due to the elastodynamic behavior. In 18 the vertical displacement over time response is plotted at the center of the plate at the thickness z = h/2. We observe a good match between the virtual element results and the finite element re- sults. Here we used in addition to regular shaped elements, Voronoi elements which have an arbitrary number of nodes and element shapes. The computation is performed with 1024 virtual elements of type H1, H2S and VO and the finite element H1. We used 6400 elements for FEM-H2 as a reference solution. Again one can observe that the computation of the mass-matrix using only the projection part and evaluating the integral in (37) at the centroid of the element yields accurate results.

Summary and Conclusions
An efficient low order virtual element formulation was developed for nonlinear elastodynamics within this work. The virtual element results show a very good match with finite element results and analytical results. Arbitrary shaped elements with a various number of nodes could be used successfully for the simulations. The computation of the mass-matrix was performed using the projection part and does not need any stabilization. This is valid only for problems, where the governing equations have no reaction terms [29,30], which was the case in this work. To compute the integral in equation (37), the argument can be evaluated at the element centroid. This is sufficient accurate as shown in the examples. Hence, there is no need to perform any sub-triangulation of the element or use the moment of areas in (46) for the computation of the mass-matrix. The extension of VEM to other applications open a wide range of new research directions such as dynamic elasto-plasticity, contact or impact.