A priori error analysis of virtual element method for contact problem

As an extension of the finite element method, the virtual element method (VEM) can handle very general polygonal meshes, making it very suitable for non-matching meshes. In (Wriggers et al. in Comput. Mech. 58:1039–1050, 2016), the lowest-order virtual element method was applied to solve the contact problem of two elastic bodies on non-matching meshes. The numerical experiments showed the robustness and accuracy of the virtual element scheme. In this paper, we establish a priori error estimate of the virtual element method for the contact problem and prove that the lowest-order VEM achieves linear convergence order, which is optimal.


Introduction
Contact phenomena among deformable bodies are abundant in structural and mechanical systems, so a considerable effort has been made in modeling and numerical analysis of the contact processes. For the formulation of contact problems, variational inequality is a powerful mathematical tool. The well-known Signorini problem, one variational inequality of the first kind, is an elastostatic problem describing the contact of a deformable body with a rigid frictionless foundation [19]. When friction effects are considered, we need other formulations to describe the frictional contact problems, which are variational inequalities of the second kind, featured by the presence of non-differentiable terms in the formulation [15]. No penetration into the foundation is allowed for a rigid foundation, so a normal compliance contact condition was proposed for studying the interpenetration of the body's surface into the foundation [1,17,20,22,24,28]. This contact problem can be described by a quasi-variational inequality, and its well-posedness is proved by the fixed-point argument [18,30]. Actually, the fixed-point argument can be applied to study the existence and uniqueness of many contact problems, see, e.g., [6,21,25,29] and the references therein.
Because the finite element method (FEM) is based on variational formulation, it is a natural numerical discretization method for variational inequalities [10,15]. However, the classical FEM works on the elements with simple geometries, like triangles and rectan- gles, so it is very hard to handle non-matching meshes to discretize contact problems. Hence, some techniques have been developed to discretize contact constraints, for example, mortar methods [13,27,39,40] impose the contact constraints in a weak sense. Due to the flexibility in constructing the local function space, the discontinuous Galerkin (DG) method [9,16,[31][32][33][34] can handle very general meshes with hanging nodes, which make them very suitable for hp adaptivity and non-matching grids. However, the DG method needs a large number of degrees of freedom.
Recently, borrowing the ideas from the mimetic finite difference method, the virtual element method (VEM), which can be regarded as a generalization of the classic FEM to general polygonal meshes, has been developed for solving a variety of partial differential equations, cf., e.g., [2,3,5,7,14,26,43]. The virtual element method can handle very general polygonal elements with an arbitrary number of edges. In addition, it allows geometrical hanging nodes in the meshes. Actually, they are treated as vertices of the polygonal elements in practice, so it is easy to implement the h adaptive strategy for the VEM. In [42], the lowest-order virtual element scheme was applied to solve the contact problem of two elastic bodies on initially non-matching meshes, but the error analysis was not given. Even the initial meshes for the two subdomains are non-matching, due to the feature of VEM, new nodes can be inserted into the contact interface without difficulty, then the non-matching meshes are transformed into matching meshes. Recently, VEMs were studied to solve various variational inequalities [8,11,35,36,38] and hemivariational inequalities [12,23,37].
In this paper, following the setup in [42], we give a priori error estimate of the virtual element method for the contact problem. Furthermore, we also consider the contact problem with the Tresca friction law. The rest of the paper is organized as follows: In Sect. 2, we describe the contact problems in variational inequality formulations. In Sect. 3, we introduce the abstract framework of the virtual element method. Section 4 establishes a priori error analysis, which shows that the lowest-order virtual element achieves optimal convergence order.

Contact problem
The contact problem of two elastic bodies without friction is an elastostatics problem describing the contact between two deformable bodies. Let i ⊂ R d (i = 1, 2 and d = 2, 3) be an open bounded connected domain with a Lipschitz boundary i that is divided into three parts iD , iN and iC with iD , iN and iC relatively open and mutually disjoint such that meas( iD ) > 0. Note that 1C coincides with 2C , so we use C to represent them. The displacement and stress tensors are second-order symmetric tensors, which take values in S d , the space of second-order symmetric tensors on R d with the inner product σ : τ = σ ij τ ij . Let ν i be the unit outward normal to i . For a vector v, denote its normal component and tangential component by v ν = v · ν and v τ = vv ν ν on the boundary. Similarly, for a tensor-valued function σ : → S d , we define its normal component σ ν = (σ ν) · ν and tangential component σ τ = σ νσ ν ν. We have the decomposition formula For a tensor valued function σ , define its divergence by Then, for any symmetric tensor σ and any vector field v, both being continuously differentiable over D, we have the following integration by parts formula: the frictionless contact problem of two elastic bodies is to find displacement fields u i : i → R d and stress fields σ i : i → S d such that [41] In the above problem, (2.2) follows from the constitutive relation of the elastic material, (2.3) is the equilibrium equation, in which volume forces of density f i acts in i . Boundary condition (2.4) means that the body is clamped on iD , and so the displacement field vanishes there. Surface tractions of density g i act on iN in (2.5). In the contact boundary condition (2.6), σ iτ = 0 implies that it is a frictionless contact. Here, [u] ν = u 1 · ν 1 + u 2 · ν 2 , and h 0 ∈ H 1 ( C ) is the initial gap. The condition (2.7) means no interpenetration between two bodies. Note that σ 1ν 1 = σ 2ν 2 on C , so we use σ ν instead. The fourth-order elasticity tensor C : i × S d → S d is assumed to be bounded, symmetric and positive definite in i , i.e., If the elastic behavior of the material is homogeneous and isotropic, then the elasticity tensor is given by where λ > 0 and μ > 0 are the Lamé coefficients. To give the weak formulation of the contact problem, we define The admissible set K is non-empty, closed, and convex. Following a standard argument [4,41], the variational formulation of the frictionless contact problem of two elastic bodies (2.2)-(2.7) is: Hence, the bilinear form is coercive and bounded for any u, v ∈ V i . For considering the effect of friction, we can replace the frictionless condition (2.6) by the Tresca friction law Here, [u] 1τ = u 1 -u 2 -[u] ν ν 1 = -[u] 2τ , and note that σ 1τ = -σ 2τ . Then the variational formulation of the frictional contact problem of two elastic bodies is: where Note that the contact problem with Coulomb's law can be approximated by a succession of states of the Tresca friction in which the friction threshold is fixed at each fixed-point iteration.

Abstract framework of VEM
In this section, following the ideas in [5,14,42], we present an abstract framework on the virtual element methods for solving the variational inequality problems (2.12) and (2.18).
With i = 1, 2, for a polygonal domain i , let T h i be a subdivision of i into elements denoted by T, set h T = diam(T), h i = max{h T : T ∈ T h i }, and h = max{h 1 , h 2 }. All the subdivisions are compatible with the boundary splitting: i = iD ∪ iN ∪ C . Note that the two meshes T h 1 and T h 2 may not be matching on the contact boundary C . However, in the light of VEM's capability of adding arbitrarily new nodes to existing discretization, some new nodes on the contact interface can be inserted into elements without difficulty so that the non-matching meshes are transformed into matching meshes, see details in [42]. Without ambiguity, we still use T h i , i = 1, 2 to denote the new matching meshes.
Following [5], we make an assumption on the decompositions {T h i } h as follows: A1. For each T h i , and every T ∈ T h i , there exists a positive constant γ such that • the length of the shortest edge is greater or equal to γ h T , • T is star-shaped with respect to a ball of radius ≥ γ h T .
Note that the bilinear form can be split as Following [5,14], we assume that there exists a finite-dimensional space V h i ⊂ V i and bilinear form a h i : V h i × V h i → R satisfying the following assumption: A2. For each T h i , we assume that there exists a symmetric bilinear form a h i (·, ·) : V h i × V h i → R, which can be split as • Stability: there exist two positive constants α * and α * , independent of h i and T, such The symmetry of a h (·, ·), stability (3.2), and the continuity (2.16) of a(·, ·) easily imply

A3. For each h, we assume that there exists an element f
Under the assumption (A1), we have the following approximation properties. For 2. For the specific constructions of the bilinear form a h i (·, ·), and f h i ∈ (V h i ) such that the assumptions (A2)-(A3) are satisfied, we refer the reader to [5,37,42].
The VE scheme for solving Problem (P 1 ) is: and Remark 3.2 Note that the functions of the lowest-order VEM are linear on the element boundary, so the inequality condition in This contact constraint can be imposed by the Lagrangian multiplier approach or penalty formulation [41,42].
The VE scheme for solving Problem (P 2 ) is:

Error estimates
In this section, we establish a priori error analysis of the virtual element method for solving the two contact problems (2.12) and (2.18).

Theorem 4.1
Let u i ∈ H 2 ( i ) and u h i be the solutions of (2.12) and (3.8), respectively. Furthermore, assume that u i | C ∈ H 2 ( C ) and h 0 ∈ H 2 ( C ), then i=1,2 where the constant C depends on C s , Proof First, we split the error as two parts From (2.15), the properties (3.1)-(3.2) of the bilinear form a h , and the discrete scheme where 3) Next, we estimate the above inequality term by term. By the boundedness (2.16), (3.4) of the bilinear forms, and the approximation properties (3.6), (3.7), we get