A virtual element method for frictional contact including large deformations

Peter Wriggers (Department of Mechanical Engineering, Leibniz University Hanover, Hannover, Germany)
Wilhelm T. Rust (Department of Mechanical Engineering, Leibniz University Hanover, Hannover, Germany)

Engineering Computations

ISSN: 0264-4401

Article publication date: 20 June 2019

Issue publication date: 12 September 2019



This paper aims to describe the application of the virtual element method (VEM) to contact problems between elastic bodies.


Polygonal elements with arbitrary shape allow a stable node-to-node contact enforcement. By adaptively adjusting the polygonal mesh, this methodology is extended to problems undergoing large frictional sliding.


The virtual element is well suited for large deformation contact problems. The issue of element stability for this specific application is discussed, and the capability of the method is demonstrated by means of numerical examples.


This work is completely new as this is the first time, as per the authors’ knowledge, the VEM is applied to large deformation contact.



Wriggers, P. and Rust, W.T. (2019), "A virtual element method for frictional contact including large deformations", Engineering Computations, Vol. 36 No. 7, pp. 2133-2161. https://doi.org/10.1108/EC-02-2019-0043



Emerald Publishing Limited

Copyright © 2019, Emerald Publishing Limited

1. Introduction

Almost all engineering applications include parts acting on other elements through an area of contact. In many cases the behaviour of the contact has an influence on the performance of machines, looking for example at friction and wear. Therefore, it is necessary to capture these effects also in the predictive numerical simulations during development phase. In the field of finite element methods (FEM), the computation of contact problems has a long history. Different approaches and discretization techniques where developed over the years by many scientists. An overview can be found in the textbooks of Laursen (2002) and Wriggers(2006).

Generally contact falls in the category of unilateral problems meaning that no penetration can take place and the deformation is restricted. These problems have to be treated as constraint optimisation problems as the contact area is a priori unknown. Within this methodology, one has to enforce contact constraints. This means that the surfaces in contact have to be coupled and contact tractions need to be transmitted at the interface. This can be achieved by the method of Lagrange multipliers, the penalty method or other techniques like the augmented Lagrangian formulation or barrier techniques.

Within a numerical simulation method another essential topic is the discretization of the contact area. Classical works started from establishing contact pointwise by evaluating contact along contacting nodes or segments Hallquist, 1979; Wriggers and Simo, 1985), and (Hallquist et al., 1985). To increase the robustness of solution schemes, the interpolation was enhanced in (Padmanabhan and Laursen, 2001) or (Wriggers et al., 2001) to reach C1-continuity which leads to a continuous and smooth field of the normal vector. With the method of isogeometric formulation, the surface quality was increased even more, this approach is discussed by several authors (Temizer et al., 2011) and (de Lorenzis et al., 2012). To have a continuous discretization, the so-called mortar methods were introduced in recent years. Laursen et al. (2012) provide a good overview. Here, by integrating on a virtual mortar domain, the contact discretisation represents both non-matching sides exactly. However, the additional discretisation of non-matching surface patches adds to the complexity of the problem. Another approach to cope with non-matching interfaces was the use of adaptive meshing. By simply moving the underlying mesh to match both meshes at the interface, as shown for cracks in Kaczmarczyk et al. (2014) and contact in Xu and Hjelmstad (2008), leads to a poor mesh quality. This in turn requires additional effort on mesh smoothing. This problem can be overcome by using polygonal elements with less requirements on the element topology. Polygonal bases were developed in Sukumar (2004) and Sukumar and Malsch (2006) and also applied to contact problems and extended to finite deformations, see Biabanaki and Khoei (2012) and Biabanaki et al. (2014).

An even more flexible polygonal method, the virtual element method (VEM), was developed in the past years. Initiated in the work by Beirão Da Veiga et al. (2013a) and Beirão Da Veiga et al. (2014) for Laplace equations further applications followed in the area of elasticity in Beirão Da Veiga et al. (2013b) and Gain et al. (2014). The method is based on a projection of polynomial ansatz functions on the physical space. This projection does not depend on the number of vertices of one element which makes the method flexible in terms of the element topology. The number of nodes can be arbitrary large permitting aligned nodes along edges or non-convex elements. Even though the VEM is in general flexible with respect to the polygonal and polyhedral shapes it requires stabilisation. Starting from a simple nodal based error in Beirão Da Veiga et al. (2013b), it occurred that for non-homogeneous meshes special care has to be taken in the stabilisation. A first, discussion was carried out by Cangiani et al. (2015) followed by a study in Beirão Da Veiga et al. (2017). Later the method was extended to finite deformations by adjusting the node based stabilisation to non-linear problems, see Beirão Da Veiga et al. (2015) and Chi et al. (2017). A different energy-based approach that is also used in this work was derived in Wriggers et al. (2017) for finite deformations, see also Nadler and Rubin (2003), Mueller et al. (2009) and Krysl (2015a, 2015b).

In this paper, we extend the contact algorithm described in Wriggers et al. (2016) to finite deformation contact problems including friction. The methodology is based on adjusting non-matching meshes by simply inserting nodes where a nodal pair along the contact interface is needed. In this way, polygonal elements with new nodes are created along the contact surface. As the elements are only changed locally the rest of the mesh is not affected. A special variant that updates nodal position and nodes within a non-linear solution procedure allows and application to problems with large deformations and large tangential sliding. The possibility to adaptively change the inserted nodes provides the necessary flexibility for large relative tangential displacements along the interface

In Section 2, the governing equations for elastic problems undergoing finite deformations including contact are summarised.

This is followed in Section 3 by the derivation of the VEM for finite strains as developed in Wriggers et al. (2017). This section also focusses on the necessary element stabilisation. Here three different stabilisation methods are discussed that are based on the element displacement error, edge based displacement error and an energy based error. These formulations are discussed in detail with regard to the application to contact problems. In case of general sliding the distribution of nodes is not always optimal, e.g. nodes can be placed very close to each other on aligned edges. Even though the VEM, in general, is shape independent these element topologies still require special focus with respect to the stabilisation. The differences are demonstrated with a test setup.

With the basic method developed, Section 4 describes the derivation of the stiffness and residual terms for contact using a nodal based penalty method. Also the handling of the frictional problem is described. The combination of a contact for small changes during one load step with an adaptive mesh adjustment leads to a method that can represent frictional problems with large tangential movement.

Finally the performance of the method is demonstrated in Section 5 with numerical examples. First, the contact formulation is validated using the classical Hertzian contact problem. Here, a restriction to small deformations allows to compare the results to an analytical solution. This is followed by the same setup, but now for large deformations. The last example is an ironing problem of a stiff indenter and a soft block. Contact tractions are compared for the frictional and frictionless case as well as two different meshes.

2. Governing equations for finite elasticity and contact

The problem description starts from two elastic bodies with the domains Ωi (i =1,2) ⊂ ℝ. Every material point has an initial position X and is mapped to the position x at time t by:

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

Here, u(X, t) denotes the displacement field. Once the deformation is defined the deformation gradient F follows from:

(2) F=Grad(φ)=1+Grad(u(X,t))
where the gradient Grad is computed with respect to X. This gradient includes stretching and rotating of the body. As only the stretching parts and not the rotations lead to internal stresses the right Cauchy–Green tensor can be defined as a strain measure and then considered in the internal energy. The Cauchy–Green tensor C is defined as:
(3) C=FTF.

To set-up a boundary value problem boundary conditions, balance and constitutive equations have to be formulated.

Each of the bodies in Figure 1 has a boundary Γi that can be split into the parts ΓiD, ΓiN and ΓiC. On the Dirichlet boundaries ΓId, the displacements are prescribed:

(4) u=u¯onΓiD.

The Neumann boundaries are related to surface tractions t that can be, together with the Cauchy theorem t = P N, formulated as:

(5) PN=t¯
where N the surface normal vector, P the first Piola–Kirchhoff stress tensor and t̄ are the prescribed surface tractions.

Finally, the boundary ΓC is related to the contact part which is an a priori unknown area where the two bodies come in contact (Γ1C = Γ2C), see Figure 1.Both bodies satisfy on Ω1 and Ω2 the equilibrium:

(6) DivP=f
with the body force f.

The stresses in Ω1 and Ω2 follow from the strain energy function Ψ(u) and the deformation gradient:

(7) P(u)=Ψ(u)F.

Here, a homogeneous compressible isotropic hyperelastic material is assumed and a neo-Hookean strain energy function is used:

(8) Ψ(u)=λ4(J(u)212ln(J(u)))+μ2(tr(C(u))22ln(J(u))),
where λ and μ are the Lamé constants and the Jacobian J = det F represents the change in volume.An equivalent formulation to equation (6) is the potential energy U(u) that yields as Euler equations equilibrium conditions. It is often used as a basis for the derivation of a discretization scheme and will here also be applied for the VEM. For the hyperelastic case, using equation (8), the potential energy function is given by:
(9) U(u)=Ω[Ψ(u)f·u]dΩΓNt·udΓ+Uc.

The energy contribution for the contact is only present for segments actually in contact. These can be defined from the gap between the two surfaces:

(10) gn=(u2u1)·n1+g0
where g0 is the initial gap. With the contact normal pressure σn = σnn · n, which has to be computed in finite deformation with variables related to the deformed configuration, the so-called Kuhn–Tucker–Karush conditions can be formulated:
(11) gn0,σn0andgnσn=0.

These conditions state that the bodies do not interpenetrate each other, the contact pressure is always negative and that the gap and normal pressure are not non-zero at the same time. For a penalty contact formulation the energy Uc contributing to (9) is:

(12) Uc=12ΓC[ϵn(gn)2+ϵt(gt)2]dΓ
with the contact penalty parameters ϵn for the normal and ϵt for the tangential direction. The tangential gap gt will be defined later. It is here formulated for the stick case of no sliding at the contact interface.

3. Virtual element method for finite strains

The VEM offers an element formulation for a mesh of polygons with a variable number of nodes. This is possible by choosing the ansatz space directly in the physical domain. In this case the displacements are directly approximated using a simple split of the ansatz space: The basis is formed by a polynomial projection function and a remainder term:

(13) uh=Πuh+(uhΠuh).

To obtain a low order approximation the projection Π is selected as a linear function:

(14) Πuh=Ha=[10X0Y0010X0Y][a1a2a6]
for every 2D element with the set of six unknown parameters a. The remainder term guarantees consistency and serves as a stabilisation part. The goal is now to compute the set of parameters a as a map of the displacement values and thus the nodal degrees of freedom of a virtual element. Here we can remark that the ansatz, while looking like the ansatz for a three-noded triangular element, is used for a virtual element with arbitrary number of nodes and thus is different from a finite element ansatz space for triangles. However, in case that the virtual element only has three nodes, it coincides with the classical linear triangular element.

3.1 Definition of the virtual element method ansatz function

The computation of this polynomial function is based on the requirement, that the gradient of the remainder is orthogonal to the gradient of a polynomial up with the same ansatz order as Πuh. This yields:

(15) Ωe[up·(uhΠuh)]dΩe=0.

As ∇up= const. and ∇(Πuh) = const. Equation (15) reduces to the condition that the gradients computed from the projected and real displacements must be equal:

(16) Ωe(Πuh)dΩe=!ΩeuhdΩe.

While the left-hand-side is easily computable, the integral on the right-hand-side cannot be computed because uh is not known within the element Ωe. Therefore, the integral on the right-hand side is transformed to a boundary integral:

(17) Πuh|e=!1ΩeΩeuhNdΓ
where N is the outward normal to the element surface. For the two-dimensional case this results to:
(18) [a3a5a4a6]=!1ΩekNVΓk[ux(Xk)Nxux(Xk)Nyuy(Xk)Nxuy(Xk)Ny]LkdΓ,
where the left hand side is related to the gradient of Πuh and the right hand side represents the integration of a general polygonal virtual element with NV vertices. Nx and Ny denote the components of the normal vector N. Here, ux(Xk) and uy(Xk) are the components of uh and linear functions defined at the boundary of the virtual element.

From this equation, the parameters a3 to a6 can be found by simple comparison; see Wriggers et al. (2017). The parameters a3 to a6 represent the strain modes. For completeness of the polynomial in equation (14), the two parameters a1 and a2 have to be determined. They follow from the condition that the average displacement when computed using the projection Πuh has to be equal to the average displacement of the ansatz function uh:

(19) 1NVi=1NVΠuh(Xi)=!1NVi=1NVuh(Xi).

Now all kinematic quantities can be expressed in terms of the projection function. The deformation gradient is therefore:

(20) F=I+Πuh|e.

For the linear ansatz in equation (14), the deformation gradient is constant; see equation (18).

3.2 Construction of the virtual element

Knowing the potential equation (12) and the kinematical quantities equation (20), a virtual element formulation can be derived as described in Wriggers et al. (2017). With the split in equation (13), the energy can be also split in a part governed by the deformation gradient of the projection and a stabilisation part:

(21) U(u)=Ae=1Ne[Up(Πuh|e)+Ustab(uh|eΠuh|e)].
where A is the assembly operator. The energy contribution of each virtual element resulting from the deformation is found by inserting the basis from Section 3.1 into the potential equation (9):
(22) Up(Πuh|e)=Ωe[Ψ(Πuh|e)f·Πuh|e]dΩΩet·Πuh|edΓ.

For a linear projection function, the deformation gradient equation (20) is constant in Ωe. However, the parameters from equation (14) still depend on the displacement and thus a non-linear strain energy function leads to a non-linear equations system for the unknowns Πuh|e that depend on uh through equations (18) and (19). Now the residual vector Re and the stiffness matrix Ke can be computed at element level using Up. This yields:

(23) Rep=Up(Πuh|e)uh|eandKep=Rep(uh|e)uh|e
with respect to the unknown displacements uh. All derivations are executed using the automatic differentiation and code generation package AceGen as in Wriggers et al. (2017). More details on AceGen and its applications can be found in Korelc and Wriggers (2016).

3.3 Stabilisation of the method

A stabilisation is needed wherever not all possible deformation modes can be represented by the basis. This is the case for the virtual element where the projection Πuh |e can only cover the basic modes. Hence, a stabilisation of these elements is necessary. Many stabilisation techniques can be used. Here, we will restrict ourselves to three different approaches to formulate the stabilisation part Ustab in equation (21).

3.3.1 Element stabilisation.

The easiest approach which was described in Beirão Da Veiga et al. (2013a) is to introduce a pointwise error measure between the approximation function and the real displacement values ûI:

(24) du(XI)=[u^IΠuh|e(XI)]
where ûI are just the nodal values at the vertices of the virtual element. When summed up over all element vertices NV, this leads to a stabilisation term:
(25) Ustab=γ2I=1NV[u^IΠuh|e(XI)]T[u^IΠuh|e(XI)]
which penalises the error. Here, a stabilisation parameter γ is introduced that for small deformational problems can be computed by:
(26) γ=γ0hD2tr(D)tr(HT(XI)H(XI))
which includes the element diameter hD as well as a part of the constitutive tensor 𝔻. Here only the parameter γ0 has to be chosen. Unfortunately, independent of the factor γ, this element stabilisation leads to large errors in the vicinity of nodes, where vertices have a small distance to each other. Especially in adaptive schemes and at contact surfaces, as shown in Figure 2, this situation can appear quite frequently leading to erroneous results. A demonstration of this including a comparison to other stabilisations can be found at the end of this section.

3.3.2 Edge stabilisation.

Another type of stabilisation is taylored to discretisations with collapsing nodes. Here again the nodal-wise error from equation (24) is computed, but then the difference in errors of two neighbouring nodes is used within the penalty term. This formulation takes into account the length of the segment (Figure 3). In the term:

(27) Ustab=γI=1NV1|XI+1XI|×[du(XI+1)du(XI)]T[du(XI+1)du(XI)]
small edges have a larger impact on the stabilisation. Here, again the stabilisation factor γ0 from equation (26) needs to be determined. A problem dependent parameter optimisation is discussed in Wriggers et al. (2016) for the small deformational contact comparing the normal pressure and tangential gap. In summary, it was observed that with an increasing factor γ0 oscillations in the tangential gap are suppressed and the contact pressure is represented more accurately. However, by exceeding an optimal value for γ0, the results show a stiffening behaviour. Thus, penalty parameters that yield optimal results lie in a small interval which is in practice not sufficient, in particular when working with unknown problems. Hence, the penalty factor appears to be problem dependent and it is not physically motivated. Even though the results are better than for the stabilisation in Section 3.3.1 they are not robust, especially for finite deformation problems.

In non-linear problems the stiffness terms depend on the deformation. Thus, the stabilisation has to consider such dependencies and hence, different approaches for non-linear problems were developed taking into account the behaviour of the constitutive tensor. In Beirão Da Veiga et al. (2015), the norm ||PF(Fe)|| and in Chi et al. (2015), the trace 14tr[PF(Fe)] is used to automatically adjust the parameter to the non-linear material behaviour. However, this factor is not easily linearisable, so it is mostly computed from the previous load step within a Newton iteration.

3.3.3 Energy stabilisation.

To create a stabilisation that will lead to a consistent linearisation for non-linear problems, a third approach was developed. For the stabilisation a rank deficient formulation (Krysl, 2015a, 2015b) introduced an energy sampling stabilisation that was earlier used in Nadler and Rubin (2003) for Cosserat point elements. This methodology was adopted for the VEM in Wriggers et al. (2017). In contrast to the other stabilisation methods, this approach is based on the error in the resulting deformation energy rather than the error in just the displacement values. The strain energy Up, see equation (21), is extended by a supplementary energy Û:

(28) U=Up+U^(uh)U^(Πuh)Ustab
so that for a convergent solution with very small elements, Û will cancel out, guaranteeing a consistent formulation. For coarser meshes, the supplementary energy will work as a non-linear penalty term. In this work, the additional energy is computed first from the approximation function and second from the real displacement values, see equation (28) where the energy is computed for all elements Ne with:
(29) U^(uh)=Ae=1NeΩeΨ^(uh|e)dΩ
and for Ûuh |e) respectively. The energy term related to the projected field Πuh can be evaluated as in Section 3.2. Again, as described in Section 3.1, the contribution of uh cannot be directly computed. Hence, a triangulation as shown in Figure 4 is performed and a separate linear ansatz is used for every triangle. Note that this triangulation only uses the nodes that define the polygonal element. Hence, no extra degrees of freedom are introduced. Based on this ansatz, the total energy in Ωe is computed as the sum over the Nint internal triangles:
(30) U^(uh|e)=m=1NintΩmψ^(um)dΩ.

The strain energy Ψ^ can be freely selected. Here, we use a Neo-Hookean material model, see also (Wriggers et al., 2017). It leads to a simpler version of equation (8):

(31) Ψ^(u)=λ^2(J(u)1)2+μ^2(tr(C(u))22ln(J(u))).

For the internal triangles, an ansatz function similar to equation (14) can be selected:

(32) um=Hmdm.

As the linear ansatz yields a constant gradient within each sub-triangle m the deformation gradient Fmd = 1 + Grad(um) can be computed directly. With the Jacobi determinant Jm = det (Fm) and the Cauchy–Green tensor Cm=FmTFm the energy equation (31) can be evaluated for each triangle in Figure 4 and then summed up with equation (30). Note that no numerical integration is necessary as all kinematical quantities in the sub-triangles are constant.

The derivation of the energy contribution to result in the residual vector and the stiffness matrix is again carried out as in equation (23) by the automatic differentiation tool AceGen from Korelc and Wriggers (2016):

(33) Res=U^(uh|e)U^(Πuh|e)uh|eandKes=R(uh|e)uh|e.

Together with the continuum part, the total residual and stiffness therefore add up to:

(34) Re=Rep+ResandKe=Kep+Kes.

The last part in this stabilisation procedure is the determination of the material parameters λ^ and μ^. The stabilisation can be optimised by adjusting the Lamé parameters in the stabilisation energy from equation (31). As described by Krysl (2016) and applied in Wriggers et al. (2017) for virtual elements, the correction can be done to favour locking free and good bending behaviour. The Poisson ratio is set to:

(35) ν^=0.3
for all elements to avoid locking. To guarantee a good bending behaviour, Krysl (2016) compares the bending energy of a hexahedron element with a beam. This implies a reduced Young’s modulus E^ for slender elements:
(36) E^=Eβ1+β
where the factor β takes into account the height to length ratio of the rectangular element. For arbitrary shaped virtual elements, this was adopted in Wriggers et al. (2017) by taking into account the ratio between an inner radius Ri and outer radius Ro of the elements convex hull, shown in Figure 5. The height hz = 2 Ri and length hx=2Ro2Ri2 lead to a ratio of one for square shaped elements. In total, this leads to the factor β as:
(37) β=2(1+ν)2Ri2Ro2Ri2.

With these values, the Lamé parameters are defined by:

(38) λ^=μ^=0,385E^.

3.3.4 Stabilisation test.

For the qualitative assessment of different stabilisation methods, Beirão Da Veiga et al. (2017) developed the test depicted in Figure 6. A regular mesh is connected to a Voronoi mesh. The displacements of both meshes are fixed at the sides and the bottom and then loaded at the top. This discretisation creates irregular spaced nodes at the interface. In Figure 7, the displacements in x-direction are plotted which show the quality of the used stabilisations.

Starting with the element stabilization, it can well be seen that the displacements at the interface of both meshes oscillate heavily, while the overall result is correct. This can also not be avoided completely by increasing the order of the polynomial basis. Therefore, the stabilisation is not sufficient as the oscillations will lead to problems when dealing with contact. Incorrect surface displacements lead to wrong gaps and incorrect contact stresses which yield erroneous contact or sliding states. The edge stabilisation in the second graph exibits a much smoother result for the surface displacements. Because of this, it is well suited for the virtual element formulation in small deformational contact computations. However, as was discussed in this section, the method still requires a parameter fitting to perform well, making it unusable for non-linear problems. This is overcome by the third method depicted by the results in the last graph in Figure 7 which shows a smooth surface displacement already for the linear basis.

4. Contact for large deformations including friction

The virtual element discretisation allows a very simple contact formulation, even in the case of large strain. The basic idea is to define additional nodes of an element at the contact interface when contact is detected. As virtual elements can have an unlimited number of nodes the addition of nodes does not change the general discretisation. It is only local. This idea which was presented first in Wriggers et al. (2016) for geometrically linear problems is now adopted for contact problems undergoing finite deformations. Surprisingly enough, no locking occurs when the number of nodes are increased in virtual elements. This was shown in Wriggers et al. (2017) for incompressible material response including finite deformations and will work also for finite deformation contact as the non-to-node contact interface does not add additional constraints.

4.1 Node insertion algorithm

One of the main advantages of the VEM contact is, that the continuum discretisation does not change, only new nodes have to be generated to formulate the contact. For this task a node insertion algorithm has to be defined. Here, the scheme described in Wriggers et al. (2016) is used. The starting point for the contact search are two contacting surfaces, Γ1 and Γ2, each divided into contact segments which lie between the nodes defining the possible contact interface. Each point x(ξ)i on the segment c can be described in a parametrised configuration with respect to the current coordinates:

(39) x1(ξ)i,0=(1ξ)xi11+ξxi+11.

Following a similar approach as in the node-to-surface description, an orthogonal projection of each surface node onto the opposite surface provides the correct contact position. For the projection of a node x02 onto a segment c defined by node xi11 to node xi+11, as shown in Figure 8, this leads to:

(40) (xi+11xi11)·[x02x1(ξ)i,0]=!0
(41) ξ=1c2(xi+11xi11)·[x02xi11]

where the length of segment c is c=||xi+11xi11|| which can be used to define the normalised tangential vector:

(42) ac=1/c(xi+11xi11).

Note that the reference for the normal and tangential vector during the projection is related to surface c.

In the general case of finite deformation contact there will be always node to segment contact. This is now adapted in VEM by the insertion of nodes. As the changes only affect one element a node can be inserted in the mesh at the found position. This creates what is generally known as hanging nodes in standard FEM, but here only the polygonal mesh is extended by an additional node. Differing from the standard contact, the projection carried out for both surfaces, see Figure 9, leading to a symmetric treatment of contact. As can be seen from Figure 9 the projection creates contact elements along each side of the current contact interface. These can be described by an adequate interpolation.

This node insertion algorithm is executed for the complete contact surface within an iteration step of the global Newton method before the actual detection step for active/in-active contact is carried out.

4.2 Contact discretization: frictionless

An approach that enforces the contact constraints gn = 0 without introducing new variable is the penalty method. It penalises the constraint equation – non-penetration condition – using the potential (12):

(43) Ucn=Γcϵn2gn2dΓ
with the penalty factor ϵn. The normal contact traction can be obtained as:
(44) tn=ϵngn.

For the case of active contact the segments from both sides coincide and therefore the interpolation is the same. The nodal gap can therefore be written in short by:

(45) gncA=[xA2xA1]·nc1.

The gap function in a contact segment can be interpolated by choosing ansatz functions for the current surface coordinates. For a segment c between two nodes the functions N1(ξc) = 1 – ξc and N2 (ξ ) = ξc interpolate the current coordinates at the surface Γ1 as:

(46) x1(ξc)=A=12NA(ξc)xA1
and for the surface Γ2 the similar functions NA (ηc) using ηc as surface coordinate lead to the interpolation within the segment c:
(47) x2(ηc)=A=12NA(ηc)xA2.

As the normal vector nc1 is constant within a segment c the gap function can be written as:

(48) gn=[x2(ηc)x1(ξc)]·nc1

Due to the node projection in each step the coordinates ξc and ηc coincide along a segment c the gap function can be written as:

(49) gn(ξc)=[A=12NA(ξc)xA2A=12NA(ξc)xA1]·nc1

Note that the current coordinates can be expressed in terms of the displacements as xAα=XAα+uAα for α ∈ 1, 2. Furthermore, the surface coordinate ξ depends on the deformation and with that on the nodal displacements, see equation (41). This is also true for the normal vector nc1=e3×ac=e3×1/c(xi+11xi11) which is related to the current configuration. Here e3 is the unit vector in the third direction. All these dependencies have to be considered when the residual and tangent matrix of the contact potential equation (12) is computed.

Again AceGen is used to compute the residual and tangent matrix:

(50) Rc=UcucandKc=Rcuc
where differentiation is performed with respect to the nodal displacements uc={u12,u22,u11,u21} that are related to a segment.

4.3 Contact discretization: friction

At the contact interface the sliding in tangential direction can be split into a stick and a sliding phase. This is related to friction which affects tangential movements and introduces further complexity in the solution of contact problems. No tangential relative movement at the contact interface is associated with stick and leads to a constraint as in normal contact. Thus to include frictional stick it is necessary to constrain the relative tangential movement. In this contribution, it is achieved by using the penalty method, see the second term in equation (12) where gt denotes the tangential relative displacement. The relative tangential displacement can be computed as:

(51) gt(ξ)=[x2(ξc)x1(ξc)]·ac
where ac is the tangent vector at the interface c; see equation (42). This discretization leads to a contribution to the potential energy:
(52) Uct=c=1nc01ϵtgt2(ξc)lcdξ
where the stick constraint gt = 0 is enforced by a penalty factor ϵt.

A Coulomb friction law provides a relation between the normal and tangential tractions depending on the friction coefficient μ̄. The classical algorithm for frictional problems is based on an operator split technique which leads to a return mapping scheme, like in elasto-plasticity, see e.g (Laursen, 2002). and (Wriggers, 2006). Here first a ‘trial’ step, assuming stick is computed. After this trial step, the contact state is checked by inserting the computed contact tractions in normal and tangential direction; see equation (44) and tt = ∈t gt,trial, into the slip function:

(53) fstr=|ϵtgt,trial||μ¯ϵngn|0.

If the slip function fstr is ≤ 0, then the contact state is related to stick in segment c and the computed variables can be used in the next iteration step. Once the slip function is violated, the contact segment c is in slip state. By assuming that relative gap function gt can be split into a stick (elastic) and sliding (plastic) part:

(54) gt=gt,stick+gt,slip
it is possible to compute within the return mapping algorithm the relative tangential sliding gt,slip. This yields (Wriggers, 2006):
(55) gt,slip=fstrϵt=|ϵtgt,trial||μ¯ϵngn|ϵt,fs>0.

With the updated tangential gap, the contact is recomputed until convergence is reached.

Due to the relative tangential slip the projected nodes do not coincide anymore. This error can be neglected for small changes during one load step. However, after each converged load step the nodes in contact are newly projected onto the surface to keep the error low. Thus, the next load step is based again on a geometry in which again node-to-node contact is realised. After reprojecting the nodes, the history variable due to friction (the slip amount gt,slip) is reset to zero. This is permitted as the Coulomb friction law does not directly depend on the slip amount gt,slip.

Let us note, that for a sufficiently small step size (which may be necessary anyway to determine the correct physical behaviour of a frictional problem), good results can be achieved by fixing the normal and tangential vectors within each load step which simplifies the linearisation.

Again the additions to the stiffness matrix and residuum directly follow by using equation (50); see Section 4.2. As sliding has no potential formulation, it is necessary to compute the residual by keeping the relative tangential slip constant during the derivation:

(56) Rc=Ucuc|gt,slip=const.andKc=Rcuc.

This yields a non-symmetric tangent matrix in case of slip.

5. Numerical examples

In this section, the new virtual element formulation for large deformation contact problems is verified by means of several examples. These are the classical geometrically linear Hertz and Mindlin problem, its large deformation extension and an ironing problem.

5.1 Hertzian contact problem for small deformations

5.1.1 Frictionless contact.

The Hertzian contact problem is widely used to verify contact discretization, as an analytical solution exists for the small deformational case. In a two-dimensional setting the Hertz problem describes the contact of an elastic cylinder with an elastic half spaces. Both bodies are considered to have an infinite extension in longitudinal direction and so a plane strain condition is assumed. The feature of this test is, that the two surfaces are geometrically non-conforming which makes the contact area considerably smaller than the total surface.

The analytical solution can be found in Johnson (1987). For the two bodies, it considers a relative curvature from the two surface radii R1 and R2 and a relative Young’s modulus:

(57) 1R*=1R1+1R2and1E*=1ν12E1+1ν22E2.

With these definitions, the width a of the contact area can be computed as a function of the load F by:

(58) a2=4R*FπE*.

Furthermore, the distribution of the surface stress follows as:

(59) p(x)=2Fπa2a2x2.

The original setup in 2D is depicted in Figure 10 where a circular cross-section is loaded with a single point load.

As the problem is symmetric, only half of the original geometry is taken into account for the discretization where only one quarter of the circle is sufficient when the single point load is approximated as a distributed load. Figure 10 demonstrates the geometry, the boundary conditions and the loading. Both bodies are additionally fixed in tangential direction along the symmetry axis. The top body of radius R1 = 10 is loaded by a distributed load pn = 25 which leads back to a total original force of F =2 R.pn = 500 for the full circle. For the foundation a Young’s modulus of E2 = 7 × 103 is chosen and a Poisson’s ratio of ν = 0.3. The indenter is stiffer with a Young’s modulus of E1 = 7 × 104 and a Poisson’s ratio of ν = 0.3.

The given constitutive and geometrical parameters lead with equation (59) to an analytical value of the contact pressure: pn = 333.

For this test which is used to check the accuracy of the node insertion algorithm in the VEM the geometry is discretised using 4 noded virtual elements. Starting from the mesh in Figure 10 where the number of elements Ne for each section is marked, the mesh is refined by homogeneously dividing the elements in N parts along each direction. The maximum pressure, pn = 333, can now be compared with the maximum pressure resulting from the mesh sequence (Figure 11). After five homogeneous refinements convergence is achieved. Looking closer at the contact pressure of this refined mesh, Figure 12 shows the contour of the surface normal pressure in comparison with the analytical solution. The result for the linear case is already in good agreement with the analytical solution. The stress contour plots in Figure 13 show the typical Hertzian stress distribution. A comparison with a classical node-to-segment contact shows the superiority of the new approach. The two-sided node-to-node contact used in the VEM contact procedure easily captures interfaces with non-matching meshes.

5.1.2 Frictional contact.

The node-to-node discretization using virtual elements is verified using the same example and the converged mesh of the frictionless case. Frictional response is triggered by an additional tangential loading Q at the top surface. For the Coulomb friction, a friction coefficient of μ̄ = 0.2 is chosen. The analytical solution distinguishes between three cases: full sticking, full sliding and a mixed stick/slip state. In this example, the tangential tractions are chosen to cause a mixed state where the outer part of the surface reaches a sliding state, but the centre still sticks with the effect that the bodies do not move relative to each other.

The analytical solution includes the sticking area –c  xc as a part of the full contact area a, see equation (58):

(60) c=a(1Qμ¯F).

Detailed derivations can be found in Johnson (1987). The resulting tangential traction tt within the sticking area is then given by:

(61) tt(x)=μ¯2Fπa2(a2x2c2x2),cxc
while the traction outside the sticking area simply reduces to tt(x) = μ̄ p(x). This implies that the tangential traction has no influence on the normal contact pressure. However, Bufler (1959) pointed out that this is only true for examples where the material properties of both bodies are equal. Otherwise the difference in Young’s modulus will have an influence on the result. To rule out this effect, the Young’s moduli are chosen equal E1 = E2 = 7 × 103 for this verification. In Figure 14, the resulting tractions are compared with analytically computed tractions. Overall the values are in good agreement. In the transition zone from the sticking area to the sliding part of the surface the analytical solution depicts a kink while the numerical result is smooth. The reason for that is that no contact node lies exactly at the stick-slip border. As the discretization with linear ansatz functions cannot produce a kink within the contact segment, the sharp edge is smoothed out. The same is true for the increasing small oscillations in normal pressure towards the edge of the contact area. On the edge of the contact area, there are always elements only partly in contact where the contact state is then reproduced not correctly.

5.2 Large deformational contact

In the second example, using the node adjustment algorithm, the Hertzian contact problem sketched in Figure 10 is computed allowing large deformations. The top surface is loaded in ten load steps with a distributed pressure of pn = 4000. Friction is considered with a Coulomb friction coefficient of μ̄ = 0.3.

The deformed configuration is shown in Figure 15 for three load steps starting from Step 1 over Step 5 up to the final configuration in Step 10. From the contour plot, it can be observed that the stresses σyy are transmitted correctly.

5.3 Ironing problem

The last example shows the behaviour in a case with large tangential sliding. This so-called ironing problem is similar to those in Yang et al. (2005) and Tur et al. (2009) for mortar contact. A deformable indenter is pressed into a deformable block.

The Neo-Hooke material in (8) is used to model the material response of the system. The block has the Lamé constants of λ2 = μ2 = 2,692.103 while the indenter is ten times stiffer with the Lamé constants λ1 = μ1 = 2,692 · 104.

Contact is established using the penalty formulation with re-projection as described in Section 4.3. The penalty factors are chosen as n = 7 · 109 and t = 7 · 104. For the Coulomb friction law, a coefficient μ̄ = 0.3 is selected. The computations are performed and compared as well for frictionless as for frictional contact. The frictional case is computed for a quadrilateral mesh and a Voronoi mesh using the VEM. The size of the block’s cross-section is 10 × 7. The indenter has a top surface of length 1 and a radius of r=2. Figure 16 shows the setup of the problem with the block being meshed by Voronoi cells.

First the indenter is pushed down into the softer block at surface position x =3 by means of a prescribed displacement uy = –1 at the top surface which is applied in 20 load steps. The vertical prescribed displacement is then fixed and the indentor is dragged along the surface by applying at the top surface the horizontal ux = 8 in 1,000 load steps.

The contour plots on the left of Figure 17 depict the deformed configuration for a displacement ux while the graphs on the right-hand side show the normal and tangential reaction forces recorded at the top surface of the indenter. For all cases, an increase in normal reaction during the indentation, see point A, can be observed. The small tangential reaction in the frictionless case is caused by the contact forces acting on the curved indenter surface. It is only zero once the indenter is in the symmetric position in the middle of the block.

For the frictional case a considerable increase in the tangential force can be noted up to the point marked with B where the contact changes its status from sticking to sliding. For an infinitely large surface the tangential force during a steady sliding motion is expected to remain constant after this point. However, for the case shown here the normal reaction and thus also the tangential reaction slowly decreases. The reason for this decrease can be found in the detailed picture of the deformed configuration in Figure 17 in the graph for the frictionless contact reactions. As the indenter approaches and then leaves the blocks edge the loading direction and area of the contact surface is changed. Resulting from this specific deformation, the reaction of the elastically deformed block acts increasingly in the direction of movement.

A detailed study of the reaction plot for the quadrilateral mesh shows a smooth normal reaction and very small oscillations in the tangential direction. These oscillations where also noted in Puso and Laursen (2004) and de Lorenzis et al. (2011) where this issue was resolved with a segment-to-segment mortar approach or a smooth NURBS discretisation. Compared to the results stated in these papers the results obtained in this work using the re-projection node-to-node contact are comparable and thus the new approach can be applied for finite tangential sliding problems. The reactions computed with the Voronoi mesh show slightly larger oscillations for the normal reaction and considerable oscillations for the frictional response, see last graph in Figure 17. This can be explained with the coarser mesh and a Voronoi discretization that does not produce (by purpose) an evenly spaced surface mesh. As discussed in (Wriggers et al., 2017) Voronoi meshes with a large difference in cell sizes depict a lower convergence rate than homogeneous quadrilateral or polygonal mesh.

Even for cases with size differences in the two contacting meshes, the new re-projection contact algorithm performs stable.

6. Conclusion

The work presented here shows the application of the VEM for large deformations to contact problems. Within the formulation of the virtual element the importance of the stabilisation term was demonstrated, especially for the treatment of contact. By using a node insertion algorithm it was possible to create a stable node-to-node contact for friction. By adjusting the mesh adaptively, using re-projection, large tangential movements are allowed. In the contact discretization the advantage of the VEM was used regarding to geometrical flexibility in shape and number of vertices. The results were verified first with an analytical solution in the small deformation regime. Then the application to a finite deformation contact problem was shown. It was observed that for large tangential sliding the node-to-node contact leads to smooth traction results and due to its simplicity is advantageous compared to classical segment based contact enforcements.

Furthermore, due to their construction virtual elements can easily be coupled with finite elements. In case of contact, it is possible to use a finite element discretization for the solids and the virtual elements with the new re-projection algorithm only at the contact interface.

In the future recent developments, in the VEM can be included in the method. With the introduction of curved surfaces for VEM, started by Beirão Da Veiga et al. (2017), the discretisation quality could be increased. Also based on the work of Hudobivnik et al. (2019), the contact formulation using node insertion can be extended to three dimensional problems.


Two continuum bodies Ω(1) and Ω(2) in contact

Figure 1.

Two continuum bodies Ω(1) and Ω(2) in contact

Small edges and collapsing nodes in adaptive mesh methods or arbitrary mesh geometries are creating stability problems

Figure 2.

Small edges and collapsing nodes in adaptive mesh methods or arbitrary mesh geometries are creating stability problems

Sketch of the terms for the linear stabilisation approaches.

Figure 3.

Sketch of the terms for the linear stabilisation approaches.

Supplementary triangulation of the element area as a basis for the integration

Figure 4.

Supplementary triangulation of the element area as a basis for the integration

The ratio between inner radius Ri and outer radius Ro of the convex hull is used to introduce the influence of the element shape into the stabilisation

Figure 5.

The ratio between inner radius Ri and outer radius Ro of the convex hull is used to introduce the influence of the element shape into the stabilisation

Setup for the testing of the different stabilisations with a coupled quadrilateral and voronoi mesh

Figure 6.

Setup for the testing of the different stabilisations with a coupled quadrilateral and voronoi mesh

Plot of x-displacement along the interface of both meshes in the stability test

Figure 7.

Plot of x-displacement along the interface of both meshes in the stability test

Projection of a node 

x02 contacting a segment

Figure 8.

Projection of a node x02 contacting a segment

With the nodal projection and insertion algorithm nodes are projected and inserted symmetrically onto both contact surfaces

Figure 9.

With the nodal projection and insertion algorithm nodes are projected and inserted symmetrically onto both contact surfaces

Hertz contact problem with a two-dimensional contact of a cylindrical body and an elastic halfspace

Figure 10.

Hertz contact problem with a two-dimensional contact of a cylindrical body and an elastic halfspace

Convergence to the maximum contact pressure pn, max for an increasing number of element divisions going from a coarse to a fine mesh

Figure 11.

Convergence to the maximum contact pressure pn, max for an increasing number of element divisions going from a coarse to a fine mesh

Comparison of the normal contact pressure computed with the VEM contact with the analytical solution for the frictionless Hertzian contact problem

Figure 12.

Comparison of the normal contact pressure computed with the VEM contact with the analytical solution for the frictionless Hertzian contact problem

Stress contours for Hertzian contact problem with a detailed view of the contact zone for different mesh sizes

Figure 13.

Stress contours for Hertzian contact problem with a detailed view of the contact zone for different mesh sizes

Comparison of normal and tangential contact tractions with the analytical solution for the frictional Hertzian contact problem

Figure 14.

Comparison of normal and tangential contact tractions with the analytical solution for the frictional Hertzian contact problem

Sequence of loading steps from 1 over 5 to 10 for the finite deformational Hertzian contact problem

Figure 15.

Sequence of loading steps from 1 over 5 to 10 for the finite deformational Hertzian contact problem

Ironing problem of a stiff indenter and a larger and softer base

Figure 16.

Ironing problem of a stiff indenter and a larger and softer base

Stress contours (left) and contact tractions (right) for the ironing problem with different meshes and friction states

Figure 17.

Stress contours (left) and contact tractions (right) for the ironing problem with different meshes and friction states


Beirão Da Veiga, L., Brezzi, F. and Marini, L. (2013b), “Virtual elements for linear elasticity problems”, SIAM Journal on Numerical Analysis, Vol. 51 No. 2, pp. 794-812.

Beirão Da Veiga, L., Lovadina, C. and Mora, D. (2015), “A virtual element method for elastic and inelastic problems on polytope meshes”, Computer Methods in Applied Mechanics and Engineering, Vol. 295, pp. 327-346.

Beirão Da Veiga, L., Lovadina, C. and Russo, A., (2017), “Stability analysis for the virtual element method”, Mathematical Models and Methods in Applied Sciences, Vol. 27 No. 13, pp. 2557-2594.

Beirão Da Veiga, L., Russo, A. and Vacca, G. (2017), “The virtual element method with curved edges”, arXiv preprint arXiv:1711.04306.

Beirão Da Veiga, L., Brezzi, F., Marini, L. and Russo, A. (2014), “The hitchhiker’s guide to the virtual element method”, Mathematical Models and Methods in Applied Sciences, Vol. 24 No. 8, pp. 1541-1573.

Beirão Da Veiga, L., Brezzi, F., Cangiani, A., Manzini, G., Marini, L. and Russo, A. (2013a), “Basic principles of virtual element methods”, Mathematical Models and Methods in Applied Sciences, Vol. 23 No. 1, pp. 199-214.

Biabanaki, S. and Khoei, A. (2012), “A polygonal finite element method for modeling arbitrary interfaces in large deformation problems”, Computational Mechanics, Vol. 50 No. 1, pp. 19-33.

Biabanaki, S., Khoei, A. and Wriggers, P. (2014), “Polygonal finite element methods for contact-impact problems on non-conformal meshes”, Computer Methods in Applied Mechanics and Engineering, Vol. 269, pp. 198-221.

Bufler, H. (1959), “Zur theorie der rollenden reibung”, Ingenieur-Archiv, Vol. 27 No. 3, pp. 137-152.

Cangiani, A., Manzini, G., Russo, A. and Sukumar, N. (2015), “Hourglass stabilization and the virtual element method”, International Journal for Numerical Methods in Engineering, Vol. 102 Nos 3/4, pp. 404-436.

Chi, H., Beirãoda Veiga, L. and Paulino, G. (2017), “Some basic formulations of the virtual element method (vem) for finite deformations”, Computer Methods in Applied Mechanics and Engineering, Vol. 318, pp. 148-192.

Chi, H., Talischi, C., Lopez-Pamies, O. and Paulino, G.H. (2015), “Polygonal finite elements for finite elasticity”, International Journal for Numerical Methods in Engineering, Vol. 101 No. 4, pp. 305-328.

de Lorenzis, L., Wriggers, P. and Zavarise, G. (2012), “Isogeometric analysis of 3D large deformation contact problems with the augmented lagrangian formulation”, Computational Mechanics, Vol. 49 No. 1, pp. 1-20.

de Lorenzis, L., Temizer, I., Wriggers, P. and Zavarise, G. (2011), “A large deformation frictional contact formulation using NURBS-based isogeometric analysis”, International Journal for Numerical Methods in Engineering, Vol. 87 No. 13, pp. 1278-1300.

Gain, A.L., Talischi, C. and Paulino, G.H. (2014), “On the virtual element method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes”, Computer Methods in Applied Mechanics and Engineering, Vol. 282, pp. 132-160.

Hallquist, J.O. (1979), NIKE2d: An Implicit, Finite-Deformation, Finite Element Code for Analysing the Static and Dynamic Response of Two-dimensional Solids,, Tech. Rep. UCRL–52678, University of California, Lawrence Livermore National Laboratory.

Hallquist, J.O., Goudreau, G.L. and Benson, D.J. (1985), “Sliding interfaces with contact-impact in large-scale lagrange computations”, Computer Methods in Applied Mechanics and Engineering, Vol. 51 Nos 1/3, pp. 107-137.

Hudobivnik, B., Aldakheel, F. and Wriggers, P. (2019), “Low order 3d virtual element formulation for finite elasto-plastic deformations”, Computational Mechanics, Vol. 63 No. 2, pp. 253-269. available at: https://doi.org/10.1007/s00466-018-1593-6

Johnson, K.L. (1987), Contact Mechanics, Cambridge university press.

Kaczmarczyk, Ł., Nezhad, M.M. and Pearce, C. (2014), “Three-dimensional brittle fracture: configurational-force-driven crack propagation”, International Journal for Numerical Methods in Engineering, Vol. 97 No. 7, pp. 531-550.

Korelc, J. and Wriggers, P. (2016), Automation of Finite Element Methods, Springer, Berlin.

Krysl, P. (2015a), “Mean-strain eight-node hexahedron with optimized energy-sampling stabilization for large-strain deformation”, International Journal for Numerical Methods in Engineering, Vol. 103 No. 9, pp. 650-670.

Krysl, P. (2015b), “Mean-strain eight-node hexahedron with stabilization by energy sampling”, International Journal for Numerical Methods in Engineering, Vol. 102 Nos 3/4, pp. 437-449.

Krysl, P. (2016), “Mean-strain 8-node hexahedron with optimized energy-sampling stabilization”, Finite Elements in Analysis and Design, Vol. 108, pp. 41-53.

Laursen, T.A. (2002), Computational Contact and Impact Mechanics, Springer, Berlin, New York, NY, Heidelberg.

Laursen, T.A., Puso, M.A. and Sanders, J. (2012), “Mortar contact formulations for deformable–deformable contact: past contributions and new extensions for enriched and embedded interface formulations”, Computer Methods in Applied Mechanics and Engineering, Vol. 205, pp. 3-15.

Mueller-Hoeppe, D.S., Loehnert, S. and Wriggers, P. (2009), “A finite deformation brick element with inhomogeneous mode enhancement”, International Journal for Numerical Methods in Engineering, Vol. 78, pp. 1164-1187.

Nadler, B. and Rubin, M. (2003), “A new 3-d finite element for nonlinear elasticity using the theory of a cosserat point”, Int. J. of Solids and Structures, Vol. 40 No. 17, pp. 4585-4614.

Padmanabhan, V. and Laursen, T. (2001), “A framework for development of surface smoothing procedures in large deformation frictional contact analysis”, Finite Elements in Analysis and Design, Vol. 37 No. 3, pp. 173-198.

Puso, M.A. and Laursen, T.A. (2004), “A mortar segment-to-segment contact method for large deformation solid mechanics”, Computer Methods in Applied Mechanics and Engineering, Vol. 193 Nos 6/8, pp. 601-629.

Sukumar, N. (2004), “Construction of polygonal interpolants: a maximum entropy approach”, International Journal for Numerical Methods in Engineering, Vol. 61 No. 12, pp. 2159-2181.

Sukumar, N. and Malsch, E.A. (2006), “Recent advances in the construction of polygonal finite element interpolants”, Archives of Computational Methods in Engineering, Vol. 13 No. 1, pp. 129-163.

Temizer, I., Wriggers, P. and Hughes, T.J.R. (2011), “Contact treatment in isogeometric analysis with NURBS”, Computer Methods in Applied Mechanics and Engineering, Vol. 200 Nos 9/12, pp. 1100-1112.

Tur, M., Fuenmayor, F. and Wriggers, P. (2009), “A mortar-based frictional contact formulation for large deformations using lagrange multipliers”, Computer Methods in Applied Mechanics and Engineering, Vol. 198 Nos 37/40, pp. 2860-2873.

Wriggers, P. (2006), Computational Contact Mechanics, 2nd ed., Springer, Berlin, Heidelberg, New York, NY.

Wriggers, P. and Simo, J. (1985), “A note on tangent stiffnesses for fully nonlinear contact problems”, Communications in Applied Numerical Methods, Vol. 1 No. 5, pp. 199-203.

Wriggers, P., Krstulovic-Opara, L. and Korelc, J. (2001), “Smooth C1- interpolations for two-dimensional frictional contact problems”, International Journal for Numerical Methods in Engineering, Vol. 51 No. 12, pp. 1469-1495.

Wriggers, P., Reddy, B., Rust, W. and Hudobivnik, B. (2017), “Efficient virtual element formulations for compressible and incompressible finite deformations”, Computational Mechanics, Vol. 60 No. 2, pp. 253-268.

Wriggers, P., Rust, W. and Reddy, B. (2016), “A virtual element method for contact”, Computational Mechanics, Vol. 58 No. 6, pp. 1039-1050.

Xu, D. and Hjelmstad, K.D. (2008), A New Node-node to-node Approach to Contact/Impact Problems for Two-dimensional Elastic Solids Subject to Finite Deformation,, Technical report, Newmark Structural Engineering, Laboratory. University of Illinois at Urbana-Champaign.

Yang, B., Laursen, T.A. and Meng, X. (2005), “Two dimensional mortar contact methods for large deformation frictional sliding”, International Journal for Numerical Methods in Engineering, Vol. 62 No. 9, pp. 1183-1225.


The first author would like to acknowledge the support of the “Deutsche Forschungsgemeinschaft” under contract “Novel finite element technologies for anisotropic media” via the priority program “Reliable Simulation Techniques in Solid Mechanics, Development of Non-standard Discretization Methods, Mechanical and Mathematical Analysis” (SPP 1748).

Corresponding author

Peter Wriggers can be contacted at: wriggers@ikm.uni-hannover.de