3D mixed virtual element formulation for dynamic elasto-plastic analysis

The virtual element method (VEM) for dynamic analyses of nonlinear elasto-plastic problems undergoing large deformations is outlined within this work. VEM has been applied to various problems in engineering, considering elasto-plasticity, multiphysics, damage, elastodynamics, contact- and fracture mechanics. This work focuses on the extension of VEM formulations towards dynamic elasto-plastic applications. Hereby low-order ansatz functions are employed in three dimensions with elements having arbitrary convex or concave polygonal shapes. The formulations presented in this study are based on minimization of potential function for both the static as well as the dynamic behavior. Additionally, to overcome the volumetric locking phenomena due to elastic and plastic incompressibility conditions, a mixed formulation based on a Hu-Washizu functional is adopted. For the implicit time integration scheme, Newmark method is used. To show the model performance, various numerical examples in 3D are presented.


Introduction
Numerical methods for solving differential equations have already been extensively investigated in recent years. The further development of these methods represents an important point to increase their efficiency. Several method can be used for the spatial discretization of the domain. Among them, the finite element method (FEM) is one of the most used one. While FEM is restricted to the usage of regular shaped elements, the recently developed virtual element method (VEM) represents one further step towards a generalization of the finite element method, see [1,2]. The virtual ele-  1 Institute for Continuum Mechanics, Leibniz Universität Hannover, Hannover, Germany ment method (VEM) allows the usage of meshes with highly irregular shaped elements, including non-convex shapes, as outlined in [2]. 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 employed to linear elastic deformations in [2,3], contact problems in [4,5], elasto-plastic deformations in [6][7][8], anisotropic materials in [9][10][11], curvilinear virtual elements for 2D solid mechanics applications in [12], hyperelastic materials at finite deformations in [13,14], crack-propagation for 2D elastic solids at small strains in [15], phase-field modeling of brittle and ductile fracture in [16,17]. Dynamic behavior has a strong influence on the mechanical behaviour of solids, failure of components and the prediction of their real response. Thus a large amount of work has been devoted to these class of problems, as well from the theoretical side as from the numerical point of view, see e.g. [18][19][20][21][22][23] and [24]. In the range of virtual element methods most of the investigations are related to static problems so far. First investigations in dynamics can be found in [25] who proposed a virtual element method for linear elastodynamics problems. However their formulations are restricted to small strain settings. This has motivated the authors into enlarge the application of VEM from the static to the dynamic case for finite deformations, see [26]. The presented contribution extends now the application of the virtual element method to 3D finite strain dynamic elasto-plastic problems, using a mixed formulation to correctly account for the case of plastic incompressibility. 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 a rankdeficient 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. Various possibilities exist to evaluate this stabilization term. To this end, Da Veiga et al. [27] proposed a stabilization term, where integrations take place at the element boundaries. Wriggers et al. presented in [14] a novel stabilization technique, based on a triangulated sub-mesh, which uses the same nodes as the original mesh. This formulation however needed an integration within the volume of the virtual element. The stabilization parameters for the latter stabilization were based on an approach first described for finite elements in Nadler and Rubin [28], generalized in Boerner et al. [29] and simplified in Krysl [30] for the stabilization of a reduced order mean-strain hexahedron. The stabilzation method described in [14] is used in this paper as well. The elasto-plastic dynamic behavior of the solid is modelled on the basis of a numerical integration scheme. Here we utilize the implicit Newmark method as documented in [31,32]. It has the advantage that mass and tangent stiffness matrix of the virtual element are combined for the solution. Thus a rank deficient mass matrix does not provide any problem once the tangent stiffness matrix is stabilized as described above. Based on this observation it is sufficient to compute the mass matrix only for the consistency part, using the projectionφ , see [1] and [26]. This approach provides some advantage in comparison with [25] who need stabilization of the mass matrix within an explicit time integration scheme. The structure of the presented work is as follows. In Sect. 2.1 the governing equations for nonlinear dynamic elasto-plasticity are outlined. In addition to the pure displacement formulation, a mixed formulation based on Hu-Washizu functional is introduced. Section 3 summarizes the virtual element method. It includes details on the computation of the element mass-matrix and the algorithmic treatment of finite strain plasticity. To verify the proposed virtual element formulations, a various number of examples are demonstrated and discussed in Sect. 4. Section 5 briefly summarizes the work and gives some concluding remarks. In this section the kinematical relations and balance laws are described together with the boundary and initial conditions for three-dimensionale solids. Variational forms are provided for the pure displacement and mixed forms as well as constitutive equations for finite plasticity.

Pure displacement formulation
In this section we summarize the finite strain elasto-plastic formulation (see e.g. [6,7]) and supplement it by the dynamic behavior. For that we first consider an elastic Body ⊂ R 3 with boundary . This boundary is decomposed into a nonoverlapping Dirichlet D and Neumann N boundary, such that D ∪ N = , see Fig. 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 and t the time.
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). A key quantity for the description of finite deformations is the deformation gradient F, defined by the Fréchet derivative as where the gradient is evaluated with respect to the initial configuration X. For elasto-plastic material behavior, the deformation gradient is decomposed into an elastic F e and a plastic part F p as which is also known as the multiplicative split of the deformation gradient for finite strain plasticity. We further consider the isochoric J 2 -plasticity theory which leads to the fact that the volume change due to plasticity can be neglected where J e and J p are the elastic and plastic part of the Jacobian J . Next, the elastic part of the left Cauchy-Green tensor b e is defined as whereb e is the deviatoric part of the elastic left Cauchy-Green tensor, with detb e = 1 and C p is the plastic part of the right Cauchy-Green tensor. The solid has to satisfy the balance of linear momentum, with the body forces f , where P, S are the 1st and 2nd Piola-Kirchhoff stresses, respectively. The right side of the equation (6) 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 Fig. 1. Furthermore we have the initial conditions at time t = 0 With the introduction of a strain energy function e for the elastic part of the deformation, the quantities such as the 1st Piola-Kirchhoff stress P and 2nd Piola-Kirchhoff stress S can be derived as follows: ∂ e ∂C e and P = FS (11) The push forward of the 2nd Piola-Kirchhoff stress to the current configuration yields the Cauchy σ and Kirchhoff stress τ : F T e and τ = J e σ (12) 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. The elastic part of the Jacobian J e is computed as Next, we use the potential energy function as a starting point for the development of a discretization method for the elastodynamics problem in (6). The static part of the potential is defined for elastic materials 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.
With the above set of equations, the finite strain elastodynamic problem is well formulated. Next the model for an elasto-plastic material behavior requires additionally the formulation of a yield function, a hardening law and an evolution equation for the plastic variables. We adopt the following yield function, which follows from the assumption of J 2plasticity with nonlinear isotropic hardening, with Here, σ v M is the von Mises stress, Y 0 the initial yield limit, Y ∞ the infinite yield stress, δ the saturation parameter, H the hardening modulus and α the hardening variable. The stress s represents the deviatoric part of the Kichhoff stress. To account for phenomenological hardening/softening response, the equivalent plastic strain α is defined as a local internal variable. For the formulation of an elasto-plastic problem an equation for the evolution of the plastic variable is needed, which can be written as, see e.g. [33], where L v b e denotes the Lie derivative in time of the left Cauchy-Green Tensor,γ ≥ 0 is the plastic Lagrange multiplier, and n is the flow direction. For the algorithmic treatment of plasticity, we adopt the Kuhn-Tucker conditions for our model: Equation (19) can be written in an equivalent form as, see [33]: The evolution equation (21) will be used for the algorithmic treatment of plasticity within the numerical solution algorithm.

Mixed Hu-Washizu formulation
To obtain a locking free behavior in the framework of J 2plasticity, we adopt the Hu-Washizu functional, see [34], for the virtual element formulation. Several mixed formulations were already discussed in the framework of virtual elements in [14] applied to finite strain hyperelastic solids. Here, based on the classical formulations in the context of finite element methods for finite strain plasticity, see [35], the following Hu-Washizu functional will be employed with In this formulation, the energy is split into an isochoric and volumetric part. In addition to that, a constraint associated with the volumetric deformation is added to the potential. Within this framework, the pressure p and the volume dilatation occur as additional independent variables.

Formulation of the virtual element method
The virtual element method is using a Galerkin projection, which maps the primary variables to a specific polynomial ansatz space. Unlike the finite element method, the isoparametric mapping for VEM is not simply obtainable. Thus simple polynomial ansatz functions can be given in terms of the coordinates X in the initial configuration. Since this contributions deals with low order VEM, the virtual element is based on linear functions and therefore contains nodes, which are placed at the element vertices.

VEM ansatz
In general, for finite strains the deformation map ϕ = X + u has to be discretized (1). But since the coordinates X in the initial configuration are exactly known, the discretization of the displacement field u = u i E i is sufficient. Here, E i are the basis vectors with respect to the initial configuration in the three-dimensional space i ∈ {1, 2, 3}. The main concept of the virtual element method is 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 j which have to be determined. The goal is to express the projection u within a virtual element v in terms of the element degrees of freedom u v : where˜ ∇ is a projection operator that will be computed next.
While the VEM ansatz for the displacements is linear, the pressure p and the dilatation are considered to be constant over the entire element, see Fig. 2. To determine the virtual parameters, u has to fulfil the orthogonality condition (29), see [1]. Thus ∇u is computed through the Galerkin projection as where p is a polynomial function which has been chosen similarly to the projection u , see (27). Since linear ansatz functions are used, ∇p and ∇u are constant at element level and can be shifted out of the integral as Applying integration by parts to (30), integral over the element volume can be transformed to an equivalent integral over its boundary, obtaining Here N denotes the normal vector on the reference boundary e of the domain e , which belongs to a virtual element e. By employing a low order ansatz, the ansatz space is linear and thus the left hand side of (31) takes the simple form To be able to determine the virtual parameters, (31) needs to be computed on the element boundary. For the 2D case, the right hand side of (31) is evaluated along the straight edges. As the displacements are known at the boundary, which are straight line segments, a linear ansatz for the displacements is used, see [14]. However, in the 3D case, the element boundary consists of polygonal faces. Therefore the evaluation of the integral in (31) is not straight forward, unless an appropriate ansatz is found. For the evaluation, there are many different methods available [3]. One option is to subdivide the element faces into 3 noded triangles, see Fig. 3. 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 [6].: Here u T h denotes the linear ansatz for the displacements at each triangle of the polygonal faces. u T is a list which contains the three nodal displacement vectors of the triangle T . ξ and η are the local reference ξ, η ∈ [0, 1] coordinates at the triangle level. The local nodes of T and the outward normal vector N i are visualized in Fig. 4.
The right hand side of (31) can now be computed. Using (34), the integral in (31) takes the form: Here n f is the number of element faces. For an integration over triangles with linear shape functions (33) one point quadrature with n g = 1 Gauss point and w g = 1/2 Gauss weight is sufficient. N ζ is the Jacobian of the 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 (32) and (35) the unknown virtual parameters a = a i j i∈ (1,...,3)∧ j∈ (2,...,4) can be linked to the nodal displacements Since only the gradient of the projection ∇u is needed to define the strain energy function of the static part, the virtual parameters which are related to the constant parts due not have to be computed in general. However, for the construction of the mass-matrix, the complete ansatz for the projected displacements u is necessary which includes the virtual parameters that are related to the constant parts. Thus our formulation has to be supplemented by a further condition to obtain the constants a = a i j i∈ (1,...,3)∧ j=1 . This condition (see for example [2]) relates the sum of the nodal values of u h to the values of the projection u . This leads for each virtual element e to the following condition where n V is the number of boundary nodes and X I are the initial coordinates of the nodal point I . The sum includes all boundary nodes n V . By substituting (27) in (40) the missing three parameters can expressed in terms of the nodal displacements and the already known virtual parameters a = a i j i∈ (1,...,3)∧ j∈ (2,...,4) as: Finally with equation (35) and (41) the ansatz function u of the virtual element is completely defined in terms of the element nodal displacements u e which can be written with the projector introduced in (28) as We note that a matrix formulation of the projector˜ ∇ can be derived, but it is not necessary since we use the symbolic tool AceGen which does this automatically.

Construction of the element mass-matrix for VEM
The inertia term in the balance equation (6) leads in a weak sense for a virtual element v to v ρü · δu d (43) where δu is the test function. Discretization of this term yields a mass-matrix which has to be constructed for the virtual element method. Our starting point is the same split for the accelerationsü h as we already used for the displacement in (26): For a linear ansatz, the projected accelerationsü at element level take for three-dimensional elements by using (42) the form: where H is the same ansatz as for the displacements (27) and u e are the accelerations of the nodal degrees of freedom. For the construction of the elastodynamic virtual element, we employ the software tool AceGen, see [36]. It generates code automatically and provides the most efficient element routines when a potential formulation is used. Thus the part of the weak form related to the mass matrix (43) will be expressed by the specific pseudo-potential whereü needs to be hold constant during the first variation to obtain exactly (43). Inserting both equations (26) and (44) for the displacements and the accelerations in (46) yields for the virtual The coupled terms are zero due to the orthogonality condition, see e.g. [1]. The first term in (47) is the consistency part, whereas the second term is the stabilization part. For the construction of the mass-matrix, it is sufficient to use the consistency term without any stabilization. This is valid, when the problem is not reaction dominated, as stated in [1,37] and shown in [26].
Using the relationship between the projected values and the unknown values for the displacement (42) and the accelerations (45), the following expression for the consistency part of the pseudo dynamic potential results for one element e Note, that the projector˜ ∇ is constant over the element.
Therefore they can be shifted out of the integral. The argument of the integral is a polynomial function up to second order for constant density For the evaluation of the integral in (48), several ways can be utilized, see [26]. As in finite element procedures, the mass-matrix for a virtual element e can now be defined

Construction of the virtual element
As already introduced in Sect. 3.1, 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. Since the nodal degrees of freedom are in each element approximated with one interpolation function per coordinate direction, the consistency part does not lead to a stable formulation. Thus an appropriate 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 [30]. The starting point for the construction of the virtual element method is the potential function The variation of (51) yields exactly the weak form of (6) when considering the nonlinear dependency of the 2nd Piola-Kirchhoff stress S on the displacement u.
Assembling of all element contributions for the n v virtual elements yields the following expression: where h v denote the plastic history variables at element level.

Consistency part
For the consistency part the projection u , introduced in Sect. 3.1, is used and therefore the first part of equation (52) for each element is given by The gradient of the projection ∇u is constant on the entire domain e . Therefore, all kinematic quantities, that steam from ∇u , such as F, b e and C −1 p are constant as well: Thus the integration of the strain energy function can be simplified as: which is still nonlinear with respect to the unknown nodal degrees of freedom and plastic history variables. As already mentioned in Sect. 3.2, the third integral in (53) related to the dynamic part can be computed in different ways.
In [26] it has been shown, that the evaluation of the third integral in (53) at the centroid X c of the virtual element yields sufficient results and needs less computational effort when compared to the other evaluation schemes. Thus the massmatrix (50) takes the form The element residual related to the inertia term is obtained by the first derivative of the dynamic potential U dyn = u T v M vüv , see (48) and (50), holding the accelerationü v constant Before computing the second derivative, the Newmark method, see e.g. [31], is used for the implicit time integration which leads to a approximation of the accelerationü e in time within a time increment t = t n+1 − t n where the Newmark parameters are chosen as ζ = 1/4 and γ = 1/2. Thus, the residual for the dynamic part follows as The derivative of R dyn v with respect to the current displacement u v,n+1 leads then to the dynamic part of the tangent for an element e

Algorithmic treatment of finite strain plasticity
The discretized form of (61) follows from [33,38] and yields together with (20) the local residual: Here, C −1 p and α n are the converged history variables from the previous step and therefore given. Equation (61) contains one equation for each of the six unique components of C −1 p,n and one additional equation for the hardening variable α. For < 0, a pure elastic step follows and therefore the history variables, h e , will remain the same as from the previous time step, i.e. C −1 p = C −1 p,n and α = α n . If > 0, the set of equations (61) needs to be solved locally at the centroid of the element v which yields to an updated history field array

Box 1 Summary of the finite strain elasto-plastic material model
The resulting equations, which need to be solved at the centroid of each virtual element v , are the residual Q v in (61) which stem from the plastic routine and the residual R resulting from the first variation of the pseudo-potential (52): The above equations are solved in a nested algorithm, where first (62) needs to be solved locally at the element level in a inner Newton-Raphson loop for a fixed u v to update the The summary of the finite strain plasticity model, that leads to Q v is given in Box 1. Thus the local tangent matrix for the inner loop yields: Next, the outer Newton-Raphson loop is solved globally by using standard Newton-Rhaphson iteration procedure: K u = R. The residual R c v and tangent matrix K v at each virtual element v (63) are obtained by utilizing AceFEMs automatic differentiation techniques, which will yield to the residual and tangent of the virtual element: Note that residual R c v is obtained by holding history variables h v constant during differentiation procedure. Additionally, when deriving the tangent K c e with respect to the primary variables u e , providing the dependency Dh v DF is necessary to ensure a consistent linearization. For further details see [33].

Stabilization part
Using only the consistency term yields a rank deficient stiffness matrix and thus needs to be stabilized. In [14] Wriggers et al. a new positive definite energyÛ was introduced, with the help of which the stabilization term is redefined as: Furthermore, the positive definite energyÛ can be defined in terms of a stabilization parameter β ∈ [0, 1] and the U c : Thus the stabilization term takes the form:   Fig. 7 Necking problem-force-displacement response for two different meshes The computation of the first term of equation (69) can be done as explained in Sect. 3.3.1. The second term U c (u h , h v ) needs an approximation. An approach how to compute this part is introduced in [14]. An additional internal Mesh of triangles in 2D and tetrahedrons in 3D is introduced. The approximation of the displacement field is done by standard linear finite element ansatz functions. The nodes of the generated submesh belong to the set of nodes in the virtual element, such that no additional nodes are introduced. The potential U c (u h , h v ) can then be calculated on internal/embedded tetrahedron mesh: Based on the triangulated submesh, the displacement gradient is computed as ∇u h = Du T h DX T by employing standard FEM shape functions N T (analogous to (34) and (36)) for linear tetrahedron. The stabilization term U c (u h , h v ) contains both the static U stat c (u h , h v ) and dynamic part U dyn c (u h ). 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 tetrahedron n T , as sketched in Fig. 4. The plastic history variables need to be computed once h e and than be used in both parts of the strain energies in (69) and (72). Therefore, the computation of the left Cauchy Green tensor is performed in a approximative way by using the contact plastic strains C −1 p (h e ) from the consistency part. By doing so, this approximation yields in a non-symmetric tangent.
Since β ∈ [0, 1], it can be seen as a ratio parameter. The stabilization parameter β can be chosen freely. For β = 1 the total energy is computed using only the stabilization part. Thus the tangent results from the internal FEM-submesh with three noded triangles in 2D or four noded tetrahedron in 3D. Using β = 0 yields in a tangent, which is solely calculated from the projection part. Thus for β = 0 the computation results in a rank deficient tangent. The choice for the stabilization parameter β for hyperelasticity was analyzed in [6,7] and it has been shown that the optimal value is in the range β ∈ [0.2, 0.6]. In [6,7] the stabilization parameter β was chosen as a function of the accumulated plastic strains. Thus, with increasing amount of plastic deformation, the stabilization parameter decreases. For our investigations, we choose the approach from [6,7].
In both, the pure displacement as well as the Hu-Washizu based element framework, the following equation for β is used. Here η = 10 −3 denotes the minimum amount of stabilization, see Fig. 5. Without η, the stabilization parameter would decrease during the simulation and tend to be zero, which would result in a rank deficient tangent.
Due to the combination of both VEM and FEM, the outlined stabilization procedure is called mixed VEM-FEM-Stabilization. The total element residual vector R e and tangent matrix K e are obtained as the first and second derivative of the element total energy U (u h , h) with respect to the global unknowns u e , analogous to (65), keeping the same internal variables and its dependencies as for projected part.

Hu-Washizu VEM
For the Hu-Washizu formulation, the same potential as in (69) is used as well, But it is only necessary to use the mixed formulation only for the consistency part of the virtual element. Thus for the consistency part is exchanged by the Hu-Washizu potential (22). By inserting the projected quantities in to the consistency part, the total energy yields: Fig. 11 Taylor Anvil test-a geometry and boundary value problem, b deformed state (schematic illustration) Note that the derivative of (72) needs to be taken with respect to all variables. Looking at (72), it is clear that the Hu-Washizu virtual element is based only on the projection part of the mixed terms. The stabilization is constructed purely on the displacement potential. Test computations have demonstrated that a stabilization of the mixed part is not necessary and additionally such a reduced formulation leads to a more efficient element since the mixed variables need not to be treated within the internal triangular mesh. This leads to an element tangent K v with the structure Since p and are constant over the entire virtual element, static condensation can be applied to eliminate these constant variables at element level.

Numerical examples
In this section the performance of the proposed mixed virtual element formulation will be investigated. 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 provided in Table 1, unless it is otherwise specified.
In this contribution, the following mesh types for first order virtual element discretizations are introduced: The stabilization parameter of the static part β stat is chosen in all the simulations using (71), unless it is otherwise specified. For the dynamic part the mass-matrix is computed according to (56) without any stabilization.

Necking problem
In the first numerical example the proposed element formulations will be tested and compared for the quasi static case. Necking of cylindrical bar due to prescribed displacements along axial direction is considered, see [6]. This example serves to illustrate the robustness of the mixed virtual element method for localization of plastic strains in the necking area.   The geometrical setup and the boundary conditions of the cylindrical bar with diameter d = 1 mm and length L = 10 mm is depicted in Fig. 6. The material parameters can be taken from Table 1. Figure 7 depicts the load-displacement curves for two different mesh discretization. The prescribed displacement is applied at the center of the cross section. It can be observed that all elements give nearly the same force response until the necking appears. Thereafter at aboutū = 0.7 mm, the FEM H1 element demonstrates stiffer results compared to all other elements due to an expected locking behaviour. However, the similar, displacement based virtual elements (i.e. VEM H1 and VEM VO) perform way better but still not as good as the reference FEM H2 element with a quadratic ansatz. In this regard, the newly developed mixed VEM formulation produces very good results that compare with the higher order FEM H2 element. The results are even better than the ones using the mixed FEM H1JP element as shown in Fig. 7 for both, coarser and finer meshes.

3D beam
In the second example a three-dimensional beam is dynamically loaded by a surface load p(t) at the end of the beam, as illustrated in Fig. 8. The load is applied as a half sine function with the time period T 0 = 0.0008 and an amplitude of 45 N / mm 2 . Thereafter, the force is released and the beam is oscillating around its new position of rest. The time increment is set to t = 1μs.
The key goal of this test is to demonstrate the performance of  Table 1. Figure 9 shows the time history of the displacement at the tip of the beam. For a compressible material (i.e. ν = 0.3 outlined in Fig. 9a-c), the bending of the beam converges with increasing nodes number to nearly u = 70 mm, see Fig.  9c).
Nevertheless, the finite and virtual elements, which are based on a pure displacement formulation H1/VO, tend to provide stiffer responses after getting into the plastic regime. Such an observation is in line with the artificial stiffening due to volumetric locking. Since this example is bending dominated bending locking can also appear. By increasing the Poisson's ratio up to a nearly incompressible material (i.e. ν = 0.499999 outlined in Fig. 9d-f), a strongly stiffer response is observed for the pure displacement elements in comparison with the stable and robust mixed finite and virtual element formulations. Thus the Hu-Washizu based finite and virtual elements produce a much softer response and hence can handle incompressible material behaviour well.
For a representative comparison between all elements, the relative error of the maximum displacement (related to Fig. 9) is plotted in Fig.10 for different elements and Poisson ratios. Hereby, the error is computed with respect to an overkill solution, that is obtained from the mixed finite element FEM H1JP using 100000 elements. In Fig. 10 a (for ν = 0.3), the error is remarkably reduced by increasing the number of element for all types. In this regard, the pure displacement elements H1/VO demonstrate a high error in comparison with mixed FEM and VEM formulations in the case of coarse meshes. When increasing the Poisson's ratio, the error of the pure displacement elements is further increased, reaching its maximum for ν = 0.499999. The mixed finite and virtual elements stay nearly constant and are not effected by any kind of locking phenomena. This illustrates the importance of using a mixed formulation for virtual element, when it comes to elastic and plastic incompressibility.

Taylor anvil test
The next example presents the Taylor-Anvil problem, which is widely used to test the dynamical behaviour of metals but it is also a validation test for discretization schemes that simulate finite strains elasto-plasticity undergoing dynamic loadings, see [20,39]. Within this framework, a rod impacts at high velocity a rigid plate. This is modelled by fixing in longitudinal direction one side of the rod and by prescribing an initial velocity to all other parts of the body, as depicted in Fig. 11a.
The material parameters for the simulations are taken from the literature, see [40][41][42][43][44], and summarized in Table  2. Hereby, the saturation parameter is set to zero (δ = 0), hence the exponential term in (17) disappears and the model is reduced to linear hardening. The initial velocity is set to v 0 = 227 m/s. The time increment for the dynamic simulation is t = 0.01μs. During the impact a plastic front develops and moves upwards leading to a deformed state as shown in Fig. 11 b).
The equivalent plastic strain at the final deformation state for all element formulations is depicted in Fig. 12 and is obtained with 10000 elements. As expected large plastic deformations are observed at the end of the rod, as well documented in the literature [39,44,45]. This is due to the influence of the kinetic energy resulting in high stresses at the front of the rod where the essential boundary condition is applied. When a certain energy is dissipated, the stresses are not reaching the yield stress anymore. Therefore some elastic energy is still stored in the upper part of the rod as shown in Fig. 12. From the contour plots, all element yield similar results except the stiffer FEM H1. Figure 13 demonstrates the length change over time for different element formulations and two mesh discretization. All element types show nearly the same displacement curves over time. Locking effects for this impact test occur only for the FEM H1 discretization. For all formulations a small oscillation with low frequency can be observed which is due to the elastic response at the upper part of the rod. Figure 14 depicts the development of the mushroom radius at the lower part (z = 0). Good agreement between all element formulations -besides the FEM H1 -is also observed. Again the mixed formulation converges for finite and virtual elements (FEM/VEM H1JP) the best fast and results in the best coarse mesh accuracy, see Fig. 14a. Next, we illustrate in Table 3 the results obtained by different authors with different methods and compare them with the values from the current work depicted in Table 4. Good agreement is achieved for the proposed mixed virtual element formulation.

Punch problem
The last test example is concerned with the capability of the proposed mixed VEM formulations to solve dynamic elastic-plastic problems. For this purpose, a punch problem is selected which is subjected to high compression loading. A surface load is applied on one quarter of the block with the geometrical properties H = B = L = 50 mm. The boundary and loading conditions can be taken from Fig. 15. For a comprehensive comparison, a convergence study is performed where the number of elements is set to N = {64, 512, 1728, 4096, 10648, 21952, 39304}. The load is applied in time as a half-sine function with a time period of T 0 = 0.4 ms and an amplitude of 2.5 kN/mm 2 . Similar to the previous examples, the force is released after a half sine. The material parameters used for the numerical simulations are same as in the previous examples, see Table  1. The time increment used in this example is t = 1μs. Figure 16 illustrates the accumulated plastic strain α at the end of the simulation. As expected the pure displacement formulations H1 of FEM and VEM underestimate the large deformation behavior due to locking phenomena, resulting to small maximum values of α, see Fig. 16a, b. For VEM VO with voronoi cells, the locking phenomena is even more significant as depicted in Fig. 16c. This nonphysical behavior is overcome for both, FEM and VEM elements, by the mixed form based on the Hu-Washizu formulation as shown in Fig.  16d-f. Figure 17a depicts the time history of the displacement at the corner of the block, where the maximum displacement appears. The presented curves are obtained with 40000 elements. It can be seen, that the mixed finite element FEM H1JP leads to the largest deformation followed by the mixed virtual element VEM H1JP. Those elements illustrate a much softer response compared with the pure displacement finite and virtual elements which is related to their locking free behaviour. The same can be seen in Fig. 17b, c, which is showing the maximum displacement at the corner and its relative error for different numbers of elements. The reference solution for the error analyses is computed with the mixed finite element FEM H1JP, using around 100000 elements. A closer look reveals, that the mixed finite and virtual elements are providing a much softer response, compared to the pure displacement elements. Especially for course mesh, the mixed elements behave softer and thus are not affected by volumetric locking phenomena.

Summary and conclusions
A mixed low order virtual element formulation for threedimensional dynamic elasto-plasticity was developed in this work. The mixed approach is based on a three field Hu-Washizu potential function, which leads to a softer response of the body undergoing large deformations. This yields also for virtual elements a superior caorse grid accuracy in comparision with pure displacement elements. The presented formulation is based on a minimization of a specific pseudopotential, considering the dynamic behavior of the solid. The treatment of VEM for elasto-plasticity in this contribution is in line with the authors previous works, see [6,7]. The extension towards dynamic problems was performed using a fast and simple computation of the mass-matrix, see [26]. It has been shown that the mixed formulation for virtual elements, can prevent volumetric locking under elastic and plastic incompressibility conditions, especially for Voronoi shaped virtual elements.
Funding Information Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.