Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.016851

ARTICLE

Virtual Element Formulation for Finite Strain Elastodynamics

Dedicated to Professor Karl Stark Pister for his 95th birthday

Mertcan Cihan, Blaž Hudobivnik, Fadi Aldakheel and Peter Wriggers*

Institute for Continuum Mechanics, Leibniz Universität Hannover, Hannover, 30419, Germany
*Corresponding Author: Peter Wriggers. Email: Wriggers@ikm.uni-hannover.de
Received: 1 April 2021; Accepted: 12 May 2021

1  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 concave shapes, as outlined in [1]. 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 [2–4], contact problems in [5], elasto-plastic deformations in [6–8], anisotropic materials in [9–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] and phase-field modeling of fracture in [16,17].

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. In this regard, Park et al. [18] proposed a virtual element method for linear elastodynamics problems which is restricted to small strains. This has motivated the presented contribution in which we extend the application of VEM from the static to the dynamic case in the finite deformation range. Damping effects will not be considered in this investigation.

Typically the construction of a virtual element is divided into a projection part and a stabilization part. Within the projection part, a quantity φh, exactly know only at the edges, is first projected onto a fully known 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 part, 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, Beião da Veiga et al. [19] proposed a stabilization term, where all integrations take place on the element boundaries. Wriggers et al. [5] presented in a novel stabilization technique, which was first described for finite elements in Nadler et al. [20], generalized in Boerner et al. [21] and simplified in Krysl [22] 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 is based on the stabilization technique of [5]. 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 [18], only the stiffness matrix needs to be stabilized, whereas the mass matrix is only computed using the it projection part. For the time integration scheme, we utilized the implicit Newmark method as documented in [23,24].

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. Section 4 includes details on the computation of the element mass matrix. To verify the proposed virtual element formulations, various examples are demonstrated and discussed in Section 5. Section 6 briefly summarizes the work and gives some concluding remarks.

2  Governing Equations for Nonlinear Elastodynamics

In this section we summarize the finite strain elasto-static formulation (see e.g., [25 –27]) and supplement it by the dynamic effect. For that consider an elastic Body Ω⊂R3 with boundary Γ. This boundary is decomposed into a non-overlapping Dirichlet ΓD and Neumann ΓN boundary conditions such that ΓD∪ΓN=Γ, see Fig. 1.

Figure 1: Solid with boundary conditions

The position x of a material point in the current configuration is given by the deformation map

x=φ(X,t)=X+u(X,t), (1)

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

F=Gradφ=∇Xφ, (2)

where the gradient is evaluated with respect to the initial configuration X. We further define the right Cauchy-Green tensor C(u) with C=FTF as a strain measure and the Jacobian J(u) with J=detF as a volume map.

The solid Ω has to satisfy the balance of linear momentum

DivP+f¯=ρu¨with P=FS, (3)

where f¯ are the body forces and P, S are the 1st and 2nd Piola-Kirchhoff stresses, respectively. The right side of the Eq. (3) is taking the dynamic effects ρu¨ into consideration. The Dirichlet and Neumann boundary conditions are defined by

u=u¯on ΓD, (4)

PN=t¯on ΓN, (5)

here N is the outward unit normal vector related to the initial configuration. Eq. (4) represents Dirichlet boundary condition, where u¯ is the prescribed displacement at the constrained boundary ΓD and Eq. (5) represents Neumann boundary condition, with t¯ depicting the surface traction applied over the loaded boundary ΓN, as illustrated in Fig. 1. The weak formulation of the elastodynamics problem in (3) then takes the form

G(u,δu)=∫Ω[12S(u): C(δu)−f¯⋅δu+ρu¨⋅δu]dΩ−∫ΓNt¯⋅δudΓ, (6)

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

Ψ=κ4(I3−1−ln⁡I3)+μ2(I3−13I1−3), (7)

in terms of the bulk κ and shear μ modulus and the invariants of the right Cauchy-Green tensor I1=trC, I3=detC.

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 method1. The static part of the potential is defined as

Ustat(u)=∫Ω[Ψ(u)−f¯⋅u]dΩ−∫ΓNt¯⋅udΓ. (8)

Generally it would be necessary to start from a Hamiltonian to obtain a potential form in elastodynamics. However for the purpose of generating residuals and tangents it is sufficient to formulate a dynamic pseudo potential that describes inertial effects

Udyn(u)=∫Ωρu¨⋅udΩ, (9)

and to keep the inertia term ρu¨ fixed in the first variation (which is equivalent to using the classical weak form). The density of the solid is denoted by ρ.

The total potential is then defined as the sum of static and inertia parts

U(u)=Ustat(u)+Udyn(u). (10)

3  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 (displacements in this work) to a specific polynomial ansatz space. Thus, the domain Ω can be discretized with non-overlapping polygonal in 2D or polyhedral elements in 3D which do not need to have convex shapes [2]. 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.

3.1 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=uiEi where Ei are the orthonormal 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 uh into a projected part uΠ and a remainder uh−uΠ as

uh=uΠ+(uh−uΠ) (11)

For a linear ansatz, the projection uΠ at element level takes for three-dimensional elements the form

uΠ=(NΠ⋅ai)Ei,i∈{1,2,3},NΠ=(1,X,Y,Z),ai=(ai1,ai2,ai3,ai4), (12)

where a represents the twelve unknown virtual parameter a=⋃ai which have to be determined. Instead of using the polynomial NΠ in Eq. (12) as the interpolation function, a scaled ansatz can be used, for details see e.g., [4].

Furthermore, the projection uΠ has to fulfill an orthogonality condition, as discussed in the work of [28]. Hence the projection parameters can be computed by using a L2- and a Galerkin projection

∫Ωep⋅(uΠ−uh)dΩ=0,∫Ωe∇p⋅(∇uΠ−∇uh)dΩ=0, (13)

where p is a polynomial function which has been chosen similarly to the projection uΠ, see (12).

The strain energy function depends only on the gradient of the projection ∇uΠ|e, however, the construction of the inertia term requires the projected displacements uΠ, see e.g., [28]. Since linear ansatz functions are used, ∇p and ∇uΠ are constant at element domain Ωe and can be shifted out of the integral as

∇uΠ=1|Ωe|∫Ωe∇uhdΩ. (14)

Applying divergence theorem to (14) yields

∇uΠ|e=!1|Ωe|∫Γeuh⊗NdΓ, (15)

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 [14].

By employing the linear ansatz space, the left hand side of (15) takes the simple form

∇uΠ|e=[a12a13a14a22a23a24a32a33a34]. (16)

In the 2D case, the right hand side of (15) 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 [14]. However, in the 3D case, the element boundary consists of polygonal faces. Therefore the evaluation of the integral in (15) 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 [3], 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:

NhT=(ξ,η,1−ξ−η) (17)

uhT=NhTuTwhere uT=⋃uI∀I∈T, (18)

as outlined in [6]. Here uhT denotes the linear ansatz for the displacements at each triangle of the polygonal faces. uT is a list which contains the three nodal displacement vectors uI of the triangle T. ξ and η are the local dimensionless coordinates at the element level. The local nodes of T and the outward normal vector Ni are visualized in Fig. 2. Finally the right hand side of (15) can be computed. Using (18), the integral in (15) takes the form:

1|Ωe|∫Γeuh⊗NdΓ=1|Ωe|∑k=1nf∫ΓkuhT⊗NkdΓ=1|Ωe|∑k=1nf∑g=1ngwgNζuhgT⊗Ng (19)

Figure 2: Virtual element faces split into multiple triangles

Here nf is the number of element faces. ◻g denotes quantities which are evaluated at the Gauss point with the local coordinates ξ=1/3 and η=1/3. For an integration over triangles with linear shape functions (17) one point quadrature with ng = 1 Gauss point and wg = 1/2 Gauss weight is sufficient. The normal vector N and the Jacobian of the isoparametric mapping Nζ are evaluated as follows:

XT=NhTXI∀I∈T, (20)

gξ=∂XT∂ξ,gη=∂XT∂η,gζ=gξ×gη, (21)

Nζ=|gζ|,N=gζNζ. (22)

All quantities are related to the initial configuration. Comparing (16) and (19), the unknown virtual parameters aij|i∈(1,2,3)∧j∈(2,3,4) can be obtained by inspection, for further details see e.g., [14].

The constant part of the projection ai1|i∈(1,2,3) can be obtained from the first integral in (13). This L2-projection can be evaluated at the nodal points which yields for each virtual element Ωe

∑I=1nVuΠ(XI)=∑I=1nVuh(XI), (23)

where nV is the total number of boundary nodes and XI is the coordinate of nodal point I in the initial configuration. By substituting (12) in (23) the missing three parameters can be expressed in terms of the nodal displacements and the already known projection gradient ∇uΠ

(a11,a21,a31)=1nV∑I=1nV(uI−∇uΠXI). (24)

Finally with Eqs. (15) and (24) the ansatz function uΠ of the virtual element is completely defined in terms of the element unknowns, i.e., nodal displacements ue={u1,u2,…,unV}. The parameters aij in (16) and (24) can be related by a linear mapping to the nodal displacements and (12) can be rewritten

a=Π~∇ue⟶uΠ=H(X)Π~∇ue (25)

where H(X) is the matrix representation of the ansatz functions NΠ, see Eq. (12). It is defined in the three-dimensional case as

H(X)=1⋃1X⋃1Y⋃1Z=[100X00Y00Z000100X00Y00Z000100X00Y00Z]. (26)

3.2 Time Discretization

For the time integration scheme, we use the implicit Newmark method outlined in [23,24]. The equations for the velocity u˙=v and acceleration u¨ at time step tn+1 are given as

u˙(u)=u˙n+1(un+1)=γζΔt(un+1−un)−(γζ−1)u˙n−(γζ−1)Δtu¨n (27)

u¨(u)=u¨n+1(un+1)=1ζΔt2(un+1−un)−1ζΔtu˙n−(12ζ−1)u¨n, (28)

where γ and ζ are the so called Newmark parameters. ◻n are the known quantities from the previous time step tn. Δt=tn+1−tn is the time step. The Newmark parameters are chosen as ζ=1/4 and γ=1/2 which yields a solution without numerical damping, see e.g., [24].

3.3 Construction of the Virtual Element

As introduced in Section 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. The nodal degrees of freedom of an element are in each element projected to a polynomial projection function. Further each displacement component is approximated with the same interpolation function, but having unique set of ai parameters each. 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 [22]. For the construction of the virtual element method we start from the potential function (10). After summing up all element contributions for the ne virtual elements we obtain the following expression:

3.3.1 Consistency Part

For the consistency part, the projection uΠ as introduced in Section 3.1 is used in the total potential (10), thus the first part of Eq. (29) is given by

Uc(uΠ)=Ustat(uΠ)+Udyn(uΠ)=∫Ωe[Ψ(uΠ)−f¯⋅uΠ]Ω−∫ΓeNt¯⋅uΠΓ+∫Ωeρu¨Π⋅uΠdΩ (30)

The gradient of the projection ∇uΠ is constant on the entire domain Ωe thus all kinematic quantities, that stem from it, e.g., F=1+∇uΠ are constant as well. Hence the integration of the strain energy function can be simplified as

∫ΩeΨ(C)Ω=Ψ(C)|Ωe|, (31)

which is still nonlinear with respect to the unknown nodal degrees of freedom.

The acceleration can be evaluated from (28) as u¨Π=u¨(uΠ). Since the projected displacement uΠ is linear, the pseudo potential Udyn(uΠ) (30) is a quadratic function and can be computed in various ways as demonstrated next:

1.    First possibility is to evaluate the integral at the centroid Xc of the polygon in 2D and of the polyhedra in 3D. The displacements and the accelerations are then evaluated at the centroid and multiplied by the area (2D) or the volume (3D) of the element, which is an approximation (due to under-integration)

∫Ωeρu¨Π⋅uΠdΩ=ρu¨Π(Xc)⋅uΠ(Xc) |Ωe| . (32)

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 integration points Xg|T of each triangle (2D) or tetrahedron (3D)

∫Ωeρu¨Π⋅uΠdΩ=ρ∑TnT∑gngwg|Tdet(Xg,Ξ|T)u¨Π(Xg|T)⋅uΠ(Xg|T) (33)

Since the integral contains quadratic terms of X, Y and in 3D additionally of Z, the integration above approximates the integral.

3.    As a third option, the integral can be exactly computed using the nodal coordinates at the boundary via divergence theorem, see [29,30]

∫Ωef(X)dV =1ndim∫Γe[ ∫f(X)dX ∫f(X)dY∫f(X)dZ  ]︷fint⋅ NdΓ      =3D1ndim∑k=1nT∑g=1ngwgNζfint⋅Ng integration over nT triangles(34)

=2D1ndim∑k=1nV∑g=1ngwg|Xg, ξ|fint⋅ Ng integration over nV edges(35)

3.3.2 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 U^, with the help of which the stabilization term is redefined, as introduced in [14]

Ustab(uh−uΠ)=U^(uh)−U^(uΠ) (36)

We further define a stabilization parameter β∈(0,1] for the definition of the positive definite energy as U^=βUc, additionaly one can introduce a different β factor for each part of potential as

U^=βstatUstat+βdynUdyn (37)

Here, βstat and βdyn are the stabilization parameters for the static and dynamic part, where βstat,βdyn∈(0,1]. Thus the final stabilization term is given by

Ustab(uh−uΠ)=βstat[Ustat(uh)−Ustat(uΠ)]+βdyn[Udyn(uh)−Udyn(uΠ)] (38)

Applying Eq. (38) in Eq. (29) yields the final form of the total potential energy function

  +(1−βdyn)Udyn(uΠ)+βdynUdyn(uh) ] (40)

For the case where βstat=βdyn=β, expression (40) simplifies to

The consistency part in Eq. (40) can be computed in a similar way as described in Section 3.3.1. The stabilization part needs an approximation, see e.g., [14]. The displacement field is approximated by introducing an internal submesh of 3 noded triangles in 2D or 4 noded tetrahedra in 3D with linear ansatz functions NhT

NhT=(1−ξ−η−ζ,ξ,η,ζ) (42)

uhT=NhTuTwhere uT=⋃uI∀I∈T, (43)

where uhT denotes the ansatz for the displacements at each tetrahedron of the virtual element. uT is a matrix which contains the four nodal displacement vectors uI of the tetrahedron T. ξ, η and ζ are the local dimensionless coordinates at the element level. The nodes of the generated submesh belong to the set of nodes defining the virtual element (uT∈ue), hence no additional nodes have to be introduced. For the dynamic part, the integral can also be evaluated at the integration points of each triangle as shown for elements VEM H2S-I and H2S-II in Section 5. For this case, we would apply (33) but then replace the projected quantities with the displacement uhT and the acceleration u¨hT, see (45). The accelerations can be simply computed based on (27) and (28)

Ustat(uh)=∑T=1nT∑g=1ngwg|Tdet(Xg,Ξ|T)Ψ(uhT(Xg|T)) (44)

Udyn(uh)=∑T=1nT∑g=1ngwg|Tdet(Xg,Ξ|T)ρu¨hT(Xg|T)⋅uhT(Xg|T) (45)

The stabilization parameter 0<β≤1 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. For β=0 a rank deficient stiffness and mass matrix would be generated. The choice for the stabilization parameter β was analyzed in [6,16] and it has been shown that the optimal value is in the range β∈[0.2,0.6]. However, the split of the stabilization parameter β into an dynamic and static part, allows to set the stabilization for both parts individually. For our investigations we choose βstat=0.4 and for βdyn different values, which will be discussed in Chapter 5. Since the ansatz (43) is related to a finite element approximation, we call this stabilization procedure mixed VEM-FEM-Stabilization.

4  Solution Scheme and Linearization

A global Newton type iteration results to the following linearized system of equations:

which allows to determine at given global primary field of unknowns u (here displacements), by calculating their linear increment Δu in a classical Newton type iterative solution scheme.

To obtain the element residual vector Re and effective element tangent matrix Ke, the first and second derivative of the total energy U(uh) (40) have to be computed with respect to the element unknowns ue. In the first step the residual is obtained

Re=(1−βstat)∂Ustat(uΠ)∂ue+(1−βdyn)∂Udyn(uΠ)∂ue|u¨e=const.+βstat∂Ustat(uh)∂ue+βdyn∂Udyn(uh)∂ue|u¨e=const.=(1−βstat)RΠ,estat+(1−βdyn)RΠ,edyn+βstatRh,estat+βdynRh,edyn (47)

Note that ue has to be kept constant when evaluating the dynamic pseudo potential.

With (25) we can write for the projected displacement, velocity and acceleration

uΠ=HΠ~∇ue,u˙Π=HΠ~∇u˙e,andu¨Π=HΠ~∇u¨e (48)

and equivalently for the ansatz uhT for the subtriangularization in (43)

uhT=NhTuT,u˙hT=NhTu˙T,andu¨hT=NhTu¨T (49)

Inserting Eq. (48) into Udyn in (30) and using the procedure in (47) yieds the explicit form of the residual for the Newmark time integration with Eqs. (27) and (28) for the projected part

RΠ,edyn=(Π~∇)T∫ΩρHTHdΩΠ~∇[1ζΔt2ue,n+1−u^¨e,n] (50)

with u^¨e,n=1ζΔt2ue,n+1ζΔtu˙e,n+(12ζ−1)u¨e,n. In the same way the residual is obtained for the discretization in (49)

with u^¨nT=1ζΔt2unT+1ζΔtu˙nT+(12ζ−1)u¨nT.

The projection and stabilization part of the total element tangent for the dynamic part takes the form, note that differentiation is only performed with respect to nodal displacements at time tn+1

The total element tangent includes static and dynamics parts of the projection and stabilization

Ke=(1−βstat)KΠ,estat+(1−βdyn)KΠ,edyn+βstatKh,estat+βdynKh,edyn. (53)

All differentiations leading to the residual and tangent of the elastodynamic virtual element were obtained with the software tool AceGen, see [26]. It provides the most efficient element routines when a potential formulation is used. Let us note that the exact weak form (6) follows from Ustat and from the pseudo potential Udyn for fixed accelerations and thus the derivation above is equivalent to using the weak form directly. With expressions (27), (28), (47) and (53) at hand the global tangent matrix K and residual vector R can be assembled (46).

It is sufficient to use the consistency term alone (i.e., βdyn=1) for the construction of the dynamic part, without any stabilization, if the problem is not reaction dominated, as shown in [28,31].

The presented tangent matrix Ke (53) includes the mass matrix implicitly through the Newmark algorithm: Ke=Kestat+1ζΔt2Me. Thus the rank deficiency of mass is not a mayor factor in the simulation as will be shown in examples. Its calculation is not needed for transient boundary and initial value problems, it is however needed for the eigenvalue analysis, which is also shown in the examples. Different ways of how to integrate the mass matrix MΠ,e in (52) are explained in Section 3.3.1.

5  Numerical Examples

In this section we demonstrate the performance of the derived 2D and 3D virtual element formulation for dynamic problems at finite deformations. For comparison purposes results of the standard finite element method (FEM) are included. The material parameters used in this work are the same for all examples and provided in Table 1.

In this contribution, the following mesh types with first order virtual element discretizations are introduced. Different ways to evaluate the mass matrix MΠ,e in (52) are employed in the following element types:

•   VEM ◻: 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 (32). The ◻ denotes the element topology with:

•   VEM ◻ Stab: Elements which are additionally denoted with Stab have a stabilized mass matrix with βdyn=0.4, which means that the mass matrix is stabilized, as introduced in 3.3.2.

•   VEM ◻ BI: Elements which are additionally denoted with BI are using Eq. (35) to evaluate the dynamic part exactly on the boundary.

•   VEM ◻-I: Elements which are additionally denoted with I are using Eq. (33) to evaluate the dynamic part on the internal sub mesh.

•   VEM ◻-II: Elements which are additionally denoted with II are computed with βdyn=1, which result in a pure finite element mass matrix.

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 Q1: A regular 2D mesh with 4 noded quadliteral first order finite elements.

•   FEM Q2: A regular 2D mesh with 9 noded quadliteral second 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 all the simulations the stabilization parameter of the static part was chosen as βstat=0.4.

5.1 Eigenvalue Analysis

In this section, the eigenfrequencies of a single element and of a structural system are analyzed to check the correctness of the eigenmodes and the eigenfrequencies when compared to classical finite elements.

5.1.1 Single Element Analysis

The eigenfrequencies of a single quadliteral element which has a free-free boundary condition are shown in Table 2. Here eigenfrequencies related to rigid body motions are excluded. To investigate the effect of the stabilization parameters βstat and βdyn on the eigenfrequencies both stabilization parameters have been varied. It can be observed, that both stabilization parameters influence the eigenfrequencies of the element. In previous publications, an optimal the stabilization parameter βstat for the static part was found to be in the range 0.2–0.4, see [14]. Thus for transient analysis, βstat=0.4 is chosen as well which lead to results are close to the one obtained with a linear finite element FEM Q1. The results obtained with VEM Q1 Stab (evaluating the mass matrix at the cenroid) are slightly stiffer when compared with FEM Q1 and VEM Q1 BI Stab.

For βstat=βdyn=0.4, the first eight eigenfrequencies and eigenmodes are depicted in Fig. 3. The elements used in this analysis are:

The standard finite element FEM Q1 in Fig. 3a.

—   The virtual element VEM Q1 Stab in Fig. 3b, evaluating the dyamic part with Eq. (32).

—   The virtual element VEM Q1 BI Stab in Fig. 3c, evaluating the dynamic part with Eq. (35).

—   The virtual element VEM Q1-I Stab in Fig. 3d, evaluating the dynamic part with Eq. (33), using the submesh.

—   Note that all these virtual elements have a stabilized mass matrix, since the computation of eigenfrequencies with rank deficient mass matrix is not possible for a single element.

Figure 3: Eigenfrequencies ωn and modeshapes for single element with free-free boundary condition and stabilization parameters βstat=βdyn=0.4 (a) Eigenfrequencies ωn and modeshapes for FEM Q1 (b) Eigenfrequencies ωn and modeshapes for VEM Q1 Stab (c) Eigenfrequencies ωn and modeshapes for VEM Q1 BI Stab (d) Eigenfrequencies ωn and modeshapes for VEM Q1-I Stab

5.1.2 Structural Analysis

An eigenvalue analysis of the mass and stiffness matrix for the computation of specific initial boundary value problems is performed. A cantilever beam, which is clamped at one side and free at the other side (C-F) is considered. The material parameters are taken from Table 1. The beam has a length l=30 mm and a height h=0.3 mm. In this study the stabilization parameters for the static and the dynamic part are set to βdyn=βstat=0.4.

Fig. 4 shows the eigenvalues with respect to the mode numbers, obtained with FEM Q1 and VEM Q1 Stab. In all figures no distinction is made between the mode types. The eigenvalues are computed with a discretization that has 360 unknowns. The eigenvalues of the stiffness and mass matrix computed with VEM are very close to the eigenvalues obtained by FEM. Note that a mass matrix which is based purely on the projection part (βdyn=0) will yield in a rank deficient mass matrix und thus the eigenvalues are not computable. Nevertheless it will be shown later that mass matrices with βdyn=0 and βdyn>0) yield good results when applying them to transient initial boundary value problems.

Figure 4: Eigenvalues of stiffness and mass matrix for 2D beam (C-F) (a) Eigenvalues of stiffness matrix K (b) Eigenvalues of massmatrix M

Table 3 depicts the eigenfrequencies which are corresponding to the first six longitudinal (L) and transversal (T) modes for two different mesh densities. For the longitudinal modes, the eigenfrequencies computed with VEM Q1 Stab are nearly the same when compared to the eigenfrequencies obtained with FEM Q1 (see also Fig. 5). For the bending modes, the eigenfrequencies have some shift, but they are in a good agreement. However, when increasing the number of nodes, changing the element topology from Q1 to Q2S topology (from 4 to 8 nodes per virtual element), the quality of the computed eigenfrequencies increases for the virtual element. We note that the eigenvalues are converging to the analytical solution for all element types for refined meshes.

Figure 5: Eigenfrequencies for 2D beam (C-F) (a) Eigenfrequencies of transversal modes (b) Eigenfrequencies of longitudinal modes

5.2 2D Boundary Value Problems

5.2.1 Wave Propagation in Longitudinal Beams

In this example wave propagation in longitudinal beams is analyzed. The geometric setup and the loading conditions of the specimen are depicted in Fig. 6. Table 1 provides the material parameters. 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. Due to the high frequencies which appear in this specific example, the time increment is set to Δt=0.001 μs. The initial velocity of all nodes is set to v0,x=20 m/s. The virtual element method will be compared with the finite element method and the analytical solution. The analytical solution in (55) is obtained by solving the wave Eq. (54). Figs. 7a and 7b illustrate the wave propagation through the elastic body. The displacement field over time for different VEM and FEM formulations are compared with analytical results. The FEM results are computed for 4×200 elements, where the virtual element results are computed for 4×100 elements. We observe a good agreement of VEM compared with FEM solution and the analytical solution. In terms of the period and the amplitude of the wave, the virtual elements shows results that are close to the analytical solution. Furthermore, the time history of the displacements are nearly the same for both elements VEM Q2S Stab and VEM Q2S with stabilized and non stabilized mass matrix. This shows, that even a rank deficient mass matrix leads to sufficient accurate results.

∂2u∂t2=c2∂2u∂x2where c=Eρ (54)

u(x,t)=∑n=0∞2v0cℓωn2sin(wnxc)sin(wnt)with wn=(2n+1)πc2ℓ. (55)

For both virtual elements VEM Q2S and VEM Q2S Stab the integral for the dynamic part in Eq. (30) is evaluated at the centroid of the element, hence this simple and efficient scheme seems to be sufficient.

Figure 6: 2D example–-Wave propagation in longitudinal beams (boundary value problem)

Figure 7: Displacement over time response for 2D example–-Wave propagation in longitudinal beams (a) Response at X=ℓ/2 (b) Response at X=ℓ

5.2.2 Transversal Beam Vibration

The second benchmark test is concerned with the analysis of transversal vibrations in beams. The geometric setup and the loading conditions of the cantilever beam are depicted in Fig. 8a. For the material parameters see Table 1. The length of the bar is set to ℓ=30 mm and the height is h=5 mm. The force is applied transversal as a point load at the upper corner at X=ℓ as shown in Fig. 8b. The temporal course of the force is given by a half sine, where the maximum of the force is set to Pmax=100 kN. The period T of the applied force is adjusted to the bending stiffness of the beam and defined as

T=3.51562πℓ212ρEbh3 (56)

Figure 8: 2D example–-Transversal beam vibration. Boundary value problem in (a) and applied force in (b)

In order to analyze the position effect of the element centroid on evaluating the integral of the dynamic part, we used different type of meshes which can be seen in Fig. 9. The “animal” mesh (Fig. 9b) includes concave elements. To see the effect of using concave 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 (Fig. 9c).

Figure 9: 2D example–-Transversal beam vibration. VEM Q2S Mesh (a), VEM Animal-Mesh (b) and C-Mesh (c)

Figs. 10 and 11 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 VO and VO BI shows that it makes no difference if the integral of the dynamic part in Eq. (30) 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. Furthermore the elements VEM Q2S and VEM Q2S Stab are reproducing nearly the same response. Thus almost identical results can be reproduced with a non stabilized mass matrix.

Figure 10: Displacement over time response at X=ℓ for 2D example–-Transversal beam vibration (1)

Figure 11: Displacement over time response at X=ℓ for 2D example–-Transversal beam vibration (2)

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 FEM Q2 results and the virtual element VEM 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.

5.2.3 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 Fig. 12b. In this test a force driven scenario is applied at the right edge as a line load as depicted in Fig. 12b. The force is applied as shown in Fig. 8b with Pmax=10000 kN/mm. VEM VO mesh and regular VEM Q2S mesh are also plotted in Figs. 12a and 12c, respectively. The material data are provided in Table 1. The contour plots of the von Mises stress distribution for different elements at the time t=0.035 ms are sketched in Fig. 13. Both elements VEM Q2S and VEM Q2S Stab, which use the stabilized and non stabilized mass matrix, result in nearly the same von Mises stress distibution. The nonlinear behavior is clearly observed in the deformation process due to the dynamic effects at finite strains.

Figure 12: 2D example–-Cook’s membrane. VEM Voronoi Mesh (a), boundary value problem (b) and VEM Q2S Mesh (c)

Figure 13: 2D example–-Cook’s membrane. Von Mises stress distribution at time t=0.035 ms for different elements at same scale (a) FEM Q2 (b) VEM Q2S (c) VEM VO (d) VEM Q2S Stab (e) VEM VO Stab

Fig. 14 shows a mesh refinement study with the element division of 2N 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.

Figure 14: Convergence study–-Displacement over time response for 2D example–-Cook’s membrane. Element division 2N, where N increases from (a) to (d) (a) N = 1 (b) N = 2 (c) N = 3 (d) N = 4

This study shows that again, that the evaluation of the integral of the dynamic part in (30) at the element centroid is absolutely sufficient to compute the mass matrix and that an element with a rank deficient mass matrix reproduces almost identical responses.

5.3 3D Boundary Value Problems

5.3.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 ℓ=30 mm and the height is equal to the width h=b=5 mm. We apply an initial velocity of v0,x=20 m/s to all nodes in longitudinal direction. The material parameters can be taken from Table 1 and the time increment is set to Δt=0.01 μm. 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 VEM H2S, VEM H2S Stab, VEM H2S-I and VEM H2S-II with the finite element FEM H1 and the analytical solution which was obtained for the 1D case in Eq. (55). As already introduced before, the variable βdyn indicates how the dynamic part is going to be evaluated. For βdyn=0, the dynamic part is calculated using only the projection part. Whereas for βdyn=1 the computation of the dynamic part is carried out using the stabilization part. Fig. 15 depicts the displacement over time response in longitudinal direction at X=ℓ and X=ℓ/2. The computation of the dynamic part using VEM H2S-I and VEM H2S-II results in a very similar response. Further the computation using the projected part and evaluating the integral of the dynamic part in (30) at the element centroid (i.e., VEM H2S) produces nearly the same results as the finite element H1 and the analytical solution.

Figure 15: Displacement over time response for 3D example-Wave propagation in a bar that having an initial velocity of v0,x=20 m/s (a) Response at X=ℓ (b) Response at X=ℓ/2 (c) Response at X=ℓ/2

5.3.2 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 Fig. 16. Here a line load is applied along the upper edge at the end of the beam with Pmax=6 kN/mm. The temporal course of the force is again given by a half sine, as shown in Fig. 8b. In this example we used the same material parameters, see Table 1. 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 elements VEM H1, VEM H2S and VEM VO with non-stabilized mass matrix and VEM H1 Stab, VEM H2S Stab and VEM VO Stab with stabilized mass matrix with the finite elements FEM H1 and FEM H2. For this purpose a mesh refinement is employed from 8, 32, 128 to 1024 elements (N = 1, 2, 3, 4). The FEM H2 solution is computed with 3200 elements and can be seen as a reference solution. The maximum deformation state is sketched in Fig. 17b, representing the deflection w. Here the nonlinear behavior is clearly observed due to the dynamic effects at finite strains. Fig. 19 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.

Figure 16: 3D example–-Transversal vibration of a thick beam (boundary value problem)

Additionaly, we employ the virtual element VEM H2S computed with 256 elements in Fig. 17a 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.

Figure 17: 3D example–-Transversal vibration of a thick beam. (a) Displacement over time response at X=ℓ and (b) undeformed and maximal deformed mesh

This test also confirms that evaluating the integral of the dynamic part in (30) for the computation of the mass matrix only at the centroid of the polygon/polyhedra is absolutely enough to get satisfying results. Furthermore, the virtual elements with stabilized and non-stabilized mass matrix lead to similar results. This again shows, that a rank deficient mass matrix can also be employed to use this virtual elements for elastodynamic problems. However, for higher frequencies, further investigations need to be done.

Since the resulting stresses play an important role in engineering applications, Fig. 18 shows the von Mises stress distribution at time t=0.1 ms. The distribution of the von Mises stress shows a good agreement between all elements. Due to the unhomogenous distribution of the voronoi elements VEM VO and VEM VO Stab, the stresses are slightly lower but show a qualitatively similar distribution, compared to all other virtual and finite elements.

Figure 18: 3D example–-Transversal vibration of a thick beam. Von Mises stress distribution at time t = 0.1 ms for different elements. (a) FEM H1 (b) VEM H1 (c) VEM VO (d) FEM H2 (e) VEM H1 Stab (f) VEM VO Stab (g) VEM H2S (h) VEM H2S Stab

Figure 19: Convergence study–-Displacement over time response for 3D example–-Transversal vibration of a thick beam. Element division 2N, where N increases from (a) to (d). (a) N = 1 (b) N = 2 (c) N = 3 (d) N = 4

5.3.3 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 Fig. 20. The material parameters are the same as in the previous examples, see Table 1. We further set an initial condition such that the initial velocity of all nodes is set to v0,z=200 m/s, see Fig. 20. Fig. 22 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 Fig. 21 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 results. 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/H1 Stab, H2S/H2S Stab, VO/VO Stab and the finite elements H1 and H2. 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 of the dynamic part in (30) at the centroid of the element yields sufficiently accurate results. Further, as already observed in previous examples, the virtual elements with stabilized and non-stabilized mass matrix result in similar responses.

Figure 20: 3D example–-Vibration of a thick plate (boundary value problem)

Figure 21: 3D example–-Vibration of a thick plate. Vertical displacement over time response at the center of the plate

Figure 22: 3D example–-Vibration of a thick plate. Evolution of the displacement in the z-direction for different deformation states. (a) t = 0 s (b) t = 0.0000014798 (c) t = 0.00000290479 (d) t = 0.00000535609 (e) t = 0.000007597 (f) t = 0.00000989465

6  Summary and Conclusions

In this work, an efficient low order virtual element formulation for nonlinear elastodynamic was developed. The presented contribution does not consider the effect of damping, which can be included in further works. A formulation that derives single tangent matrix of dynamics problem is derived. The Newmark time discretization is performed within element itself. Thus the only unknowns of our problem are nodal displacements whereas the mass matrix is not required for solving simulations. However, the mass matrix can simply be exported, for eigenvalue analysis if required. Further solving the elastodynamics problem, employing the Newmark time integration on a global level will lead to similar results. Various schemes to integrate the dynamic part were shown, with and without stabilization. It was shown that the dynamic part does not need to be stabilized for the correctness and convergence of the procedure, unless eigenvalue analysis is needed. The virtual element results show a very good match with finite elements and analytical results for boundary and initial value problems. Arbitrary shaped elements with a various number of nodes could be used successfully for the simulations.

It shows that within this framework, the stabilization of the mass matrix is not needed. This is valid only for problems, where the equations are not reaction dominated [28,31]. To compute the integral of the dynamic part in Eq. (30), the argument can be evaluated at the element centroid. This is sufficiently 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 (35) 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. Further we want to emphasize, that further investigations to high frequency applications, such as impact problems would also be interesting.

1Starting from a potential form is advantageous when using automatic differentiation, like the software tool AceGen. In this way the most efficient code for the residuals and consistent tangent operators of the virtual element will be generated for elastodynamic problems, see [26].

Funding Statement: The authors gratefully acknowledges support for this research by the “German Research Foundation” (DFG) in (i) the Collaborative Research Center CRC 1153 and (ii) the Priority Program SPP 2020.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.