A Stress/Displacement Virtual Element Method for Plane Elasticity Problems

TL;DR

This paper introduces a low-order Virtual Element Method for 2D elasticity problems that ensures symmetric stresses, backed by stability analysis and numerical tests, improving computational approaches in plane elasticity.

Contribution

The paper presents a novel low-order VEM with symmetric stress approximation for 2D elasticity, including stability proof and convergence analysis.

Findings

Method achieves stable and convergent solutions.

Numerical tests confirm effectiveness and accuracy.

Ensures symmetric stress approximation in VEM.

Abstract

The numerical approximation of 2D elasticity problems is considered, in the framework of the small strain theory and in connection with the mixed Hellinger-Reissner variational formulation. A low-order Virtual Element Method (VEM) with a-priori symmetric stresses is proposed. Several numerical tests are provided, along with a rigorous stability and convergence analysis.

Equations

$$\left\{{\begin{aligned} &\mbox{Find }({\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt},\mathbf{u})~{}\mbox{such that}\\ &-\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}=\mathbf{f}\qquad\mbox{in $\Omega$}\\ &{\kern 0.20004pt\hbox{$\sigma$}\kern-5.71413pt\kern-0.20004pt\hbox{$\sigma$}\kern-5.71413pt\raise 0.29999pt\hbox{$\sigma$}\kern 0.20004pt}=\mathbb{C}{{\kern 0.20004pt\hbox{$\varepsilon$}\kern-4.66318pt\kern-0.20004pt\hbox{$\varepsilon$}\kern-4.66318pt\raise 0.29999pt\hbox{$\varepsilon$}\kern 0.20004pt}}(\mathbf{u})\qquad\mbox{in $\Omega$}\\ &\mathbf{u}_{|\partial\Omega}={\bf 0}\end{aligned}}\right.$$

$$\{\begin{array}{rl}& \text{mbox}Find(\sigma \sigma \text{raise}0.29999pt\sigma ,\mathbf{u})\in \Sigma \times U\text{mbox}suchthat\\ & a(\sigma \sigma \text{raise}0.29999pt\sigma ,\tau \tau \text{raise}0.29999pt\tau )+(\mathrm{div}\tau \tau \text{raise}0.29999pt\tau ,\mathbf{u})=0\forall \tau \tau \text{raise}0.29999pt\tau \in \Sigma \\ & (\mathrm{div}\sigma \sigma \text{raise}0.29999pt\sigma ,\mathbf{v})=-(\mathbf{f},\mathbf{v})\forall \mathbf{v}\in U\end{array}$$

$$||\sigma \sigma \text{raise}0.29999pt\sigma |{|}_{\Sigma }+||\mathbf{u}|{|}_U\le C||\mathbf{f}|{|}_0,$$

$$a(\sigma \sigma \text{raise}0.29999pt\sigma ,\tau \tau \text{raise}0.29999pt\tau )=\sum _{E\in {\mathcal{T}}_h}{a}_E(\sigma \sigma \text{raise}0.29999pt\sigma ,\tau \tau \text{raise}0.29999pt\tau )\text{ with }{a}_E(\sigma \sigma \text{raise}0.29999pt\sigma ,\tau \tau \text{raise}0.29999pt\tau ):={\int }_E\mathbb{D}\sigma \sigma \text{raise}0.29999pt\sigma :\tau \tau \text{raise}0.29999pt\tau$$

$$(\mathrm{div}\tau \tau \text{raise}0.29999pt\tau ,\mathbf{v})=\sum _{E\in {\mathcal{T}}_h}(\mathrm{div}\tau \tau \text{raise}0.29999pt\tau ,\mathbf{v}{)}_E\text{ with }(\mathrm{div}\tau \tau \text{raise}0.29999pt\tau ,\mathbf{v}{)}_E:={\int }_E\mathrm{div}\tau \tau \text{raise}0.29999pt\tau \cdot \mathbf{v},$$

$$\mathbb{C}\epsilon \epsilon \text{raise}0.29999pt\epsilon =2\mu \epsilon \epsilon \text{raise}0.29999pt\epsilon +\lambda \mathrm{tr}(\epsilon \epsilon \text{raise}0.29999pt\epsilon )Id\forall \text{mbox}\mathrm{symmetrictensor}\epsilon \epsilon \text{raise}0.29999pt\epsilon ,$$

$$K=\left\{\tau \tau \text{raise}0.29999pt\tau \in \Sigma :(\mathrm{div}\tau \tau \text{raise}0.29999pt\tau ,\mathbf{v})=0\forall \mathbf{v}\in U\right\}.$$

$$a(\tau \tau \text{raise}0.29999pt\tau ,\tau \tau \text{raise}0.29999pt\tau )\ge \alpha ||\tau \tau \text{raise}0.29999pt\tau |{|}_{\Sigma }^2\forall \tau \tau \text{raise}0.29999pt\tau \in K,$$

$$h_E:=\mathrm{diameter}(E),h:=\underset{E\in {\mathcal{T}}_h}{\sup }{h}_E.$$

$$RM(E)=\left\{\mathbf{r}(\mathbf{x})=\mathbf{a}+b(\mathbf{x}-{\mathbf{x}}_C{)}^{\perp }\mathbf{a}\in {\mathbb{R}}^2,b\in \mathbb{R}\right\}.$$

$$R(e)=\left\{\mathbf{t}(s)=\mathbf{c}+ds\mathbf{n}\mathbf{c}\in {\mathbb{R}}^2,d\in \mathbb{R},s\in [-1/2,1/2]\right\}.$$

$$\displaystyle\Sigma_{h}(E)=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in$$

$$\displaystyle({\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\,\mathbf{n})_{|e}\in R(e)\quad\forall e\in\partial E;\quad\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in RM(E)\Big{\}}.$$

$$\displaystyle\Sigma_{h}(E)=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in$$

$$\displaystyle({\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\,\mathbf{n})_{|e}\in R(e)\quad\forall e\in\partial E;\quad\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in RM(E)\Big{\}}.$$

$${\int }_E\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h\cdot \mathbf{r}={\int }_{\partial E}\phi \phi \text{raise}0.29999pt\phi \cdot \mathbf{r}\forall \mathbf{r}\in RM(E),$$

$$\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h={\alpha \alpha \text{raise}0.29999pt\alpha }_E+{\beta }_E(\mathbf{x}-{\mathbf{x}}_C{)}^{\perp },$$

$$\{\begin{array}{rl}& {\alpha \alpha \text{raise}0.29999pt\alpha }_E=\frac{1}{|E|}{\int }_{\partial E}\phi \phi \text{raise}0.29999pt\phi =\frac{1}{|E|}\sum _{e\in \partial E}{\int }_e{\mathbf{c}}_e\\ & {\beta }_E=\frac{1}{{\int }_E|\mathbf{x}-{\mathbf{x}}_C{|}^2}{\int }_{\partial E}\phi \phi \text{raise}0.29999pt\phi \cdot (\mathbf{x}-{\mathbf{x}}_C{)}^{\perp }=\frac{1}{{\int }_E|\mathbf{x}-{\mathbf{x}}_C{|}^2}\sum _{e\in \partial E}{\int }_e({\mathbf{c}}_e+d_{e}s\mathbf{n})\cdot (\mathbf{x}-{\mathbf{x}}_C{)}^{\perp }.\end{array}$$

$$U_{h}(E)=\Big{\{}\mathbf{v}_{h}\in L^{2}(E)^{2}\ :\ \mathbf{v}_{h}\in RM(E)\Big{\}}.$$

$${\int }_E\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h\cdot {\mathbf{v}}_h$$

$$a_E({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h)={\int }_E\mathbb{D}{\sigma \sigma \text{raise}0.29999pt\sigma }_h:{\tau \tau \text{raise}0.29999pt\tau }_h$$

$$\{\begin{array}{rl}& {\Pi }_E:{\Sigma }_h(E)\to {\mathcal{P}}_0(E{)}_s^{2\times 2}\\ & {\tau \tau \text{raise}0.29999pt\tau }_h\mapsto {\Pi }_E{\tau \tau \text{raise}0.29999pt\tau }_h\\ & a_E({\Pi }_E{\tau \tau \text{raise}0.29999pt\tau }_h,{\pi \pi \text{raise}0.29999pt\pi }_0)=a_E({\tau \tau \text{raise}0.29999pt\tau }_h,{\pi \pi \text{raise}0.29999pt\pi }_0)\forall {\pi \pi \text{raise}0.29999pt\pi }_0\in {\mathcal{P}}_0(E{)}_s^{2\times 2}\end{array}$$

$$a_E^h({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h)$$

$$={\int }_E\mathbb{D}({\Pi }_E{\sigma \sigma \text{raise}0.29999pt\sigma }_h):({\Pi }_E{\tau \tau \text{raise}0.29999pt\tau }_h)+s_E\left((Id-{\Pi }_E){\sigma \sigma \text{raise}0.29999pt\sigma }_h,(Id-{\Pi }_E){\tau \tau \text{raise}0.29999pt\tau }_h\right),$$

$$s_E({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h):={\kappa }_E{h}_E{\int }_{\partial E}{\sigma \sigma \text{raise}0.29999pt\sigma }_h\mathbf{n}\cdot {\tau \tau \text{raise}0.29999pt\tau }_h\mathbf{n},$$

$$s_E({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h):={\kappa }_E\sum _{e\in \partial E}{h}_e{\int }_e{\sigma \sigma \text{raise}0.29999pt\sigma }_h\mathbf{n}\cdot {\tau \tau \text{raise}0.29999pt\tau }_h\mathbf{n}.$$

$$(\mathbf{f},{\mathbf{v}}_h)={\int }_{\Omega }\mathbf{f}\cdot {\mathbf{v}}_h=\sum _{E\in {\mathcal{T}}_h}{\int }_E\mathbf{f}\cdot {\mathbf{v}}_h.$$

$$\Sigma_{h}=\Big{\{}{\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h}\in H(\mathop{\bf div}\nolimits;\Omega)\ :\ {\kern 0.20004pt\hbox{$\tau$}\kern-4.37154pt\kern-0.20004pt\hbox{$\tau$}\kern-4.37154pt\raise 0.29999pt\hbox{$\tau$}\kern 0.20004pt}_{h|E}\in\Sigma_{h}(E)\quad\forall E\in{\mathcal{T}}_{h}\Big{\}}.$$

$$U_{h}=\Big{\{}\mathbf{v}_{h}\in L^{2}(\Omega)^{2}\ :\ \mathbf{v}_{h|E}\in U_{h}(E)\quad\forall E\in{\mathcal{T}}_{h}\Big{\}}.$$

$$a_h({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h):=\sum _{E\in {\mathcal{T}}_h}{a}_E^h({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h).$$

$$\{\begin{array}{rl}& \text{mbox}Find({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\mathbf{u}}_h)\in {\Sigma }_h\times U_h\text{mbox}suchthat\\ & a_h({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h)+(\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h,{\mathbf{u}}_h)=0\forall {\tau \tau \text{raise}0.29999pt\tau }_h\in {\Sigma }_h\\ & (\mathrm{div}{\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\mathbf{v}}_h)=-(\mathbf{f},{\mathbf{v}}_h)\forall {\mathbf{v}}_h\in U_h.\end{array}$$

$${\mathcal{A}}_h({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\mathbf{u}}_h;{\tau \tau \text{raise}0.29999pt\tau }_h,{\mathbf{v}}_h):=a_h({\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\tau \tau \text{raise}0.29999pt\tau }_h)+(\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h,{\mathbf{u}}_h)+(\mathrm{div}{\sigma \sigma \text{raise}0.29999pt\sigma }_h,{\mathbf{v}}_h),$$

Full text

1 Introduction

The Virtual Element Method (VEM) is a new technology for the approximation of partial differential equation problems. VEM was born in 2012, see [7], as an evolution of modern mimetic schemes (see for instance [18, 5, 15]), which shares the same variational background of the Finite Element Method (FEM). The initial motivation of VEM is the need to construct an accurate conforming Galerkin scheme with the capability to deal with highly general polygonal/polyhedral meshes, including “hanging vertexes” and non-convex shapes. The virtual element method reaches this goal by abandoning the local polynomial approximation concept, and uses, instead, approximating functions which are solution to suitable local partial differential equations (of course, connected with the original problem to solve). Therefore, in general, the discrete functions are not known pointwise, but a limited information of them is at disposal. The key point is that the available information are indeed sufficient to implement the stiffness matrix and the right-hand side. We remark that VEM is not the only available technology for dealing with polytopal meshes: a brief representative sample of the increasing list of technologies that make use of polygonal/polyhedral meshes can be found in [17, 8, 15, 9, 11, 13, 28, 30, 32, 35, 39, 37, 40, 41, 42, 23, 34, 43, 20]. We here recall, in particular, the polygonal finite elements and the mimetic discretisation schemes. However, VEM is experiencing a growing interest towards Structural Mechanics problems, also in the engineering community. We here cite the recent works [24, 22, 2, 3, 44, 21, 1] and [4, 19], for instance.

In the present paper we apply the VEM concept to two-dimensional elasticity problems in the framework of small displacements and small deformations. More precisely, we consider the (mixed) Hellinger-Reissner functional (see, for instance, [12, 14]) as the starting point of the discretization procedure. Thus, the numerical scheme approximates both the stress and the displacement fields.

It is well-known that in the Finite Element practice, designing a stable and accurate element for the Hellinger-Reissner functional, is not at all a trivial task. Essentially, one is led either to consider quite cumbersome schemes, or to relax the symmetry of the Cauchy stress field, or to employ composite elements (a discussion about this issue can be found in [12], for instance). We here exploit the flexibility of the VEM approach to propose and study a low-order scheme, with a-priori symmetric Cauchy stresses, that can be used for general polygons, from triangular shapes on. Furthermore, the method is robust with respect to the compressibility parameter, and therefore can be used for nearly incompressible situations. Our scheme approximates the stress field by using traction degrees of freedom (three per each edge), while the displacement field inside each polygon is essentially a rigid body motion. The VEM concept is then applied essentially for the stress field. We also remark that the construction of the discrete stress field is somehow similar to the construction of the discrete velocity field used for the Stokes problem in [10]. Instead, the displacement field is modelled with polynomial functions, in accordance with the classical Finite Element procedure.

An outline of the paper is as follows. In Section 2 we briefly introduce the Hellinger-Reissner variational formulation of the elasticity problem. Section 3 concerns with the discrete problem: all the bilinear and linear forms are introduced and detailed. Numerical experiments are reported in Section 4, where suitable error measures are considered. These numerical tests are supported by the stability and convergence analysis developed in Section 5. Finally, Section 6 draws some conclusions, including possible future extensions of the present study.

Throughout the paper, given two quantities $a$ and $b$, we use the notation $a\lesssim b$ to mean: there exists a constant $C$, independent of the mesh-size, such that $a\leq C\,b$. Moreover, we use standard notations for Sobolev spaces, norms and semi-norms (cf. [27], for example).

2 The elasticity problem in mixed form

In this section we briefly present the elasticity problem as it stems from the Hellinger-Reissner principle. More details can be found in [12, 14].

[TABLE]

Defining $(\cdot,\cdot)$ as the scalar product in $L^{2}$, and a({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}):=(\mathbb{D}{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}), a mixed variational formulation of the problem reads:

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ is a polygonal domain, $\Sigma=H(\mathop{\bf div}\nolimits;\Omega)$, $U=\times L^{2}(\Omega)^{2}$, and the loading $\mathbf{f}\in L^{2}(\Omega)^{2}$. We recall that $\mathop{\bf div}\nolimits$ is the vector-valued divergence operator, acting on a second order tensor field. Thus, \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt} is, in Cartesian components: $\frac{\partial\tau_{ij}}{\partial x_{j}}$ (Einstein’s summation convention is here adopted). The elasticity fourth-order symmetric tensor $\mathbb{D}:=\mathbb{C}^{-1}$ is assumed to be uniformly bounded and positive-definite. It is well known that problem (2) is well-posed (see [12], for instance). in particular, it holds:

[TABLE]

where $C$ is a constant depending on $\Omega$ and on the material tensor $\mathbb{C}$.

Note also that the bilinear form $a(\cdot,\cdot)$ in (2) can obviously be split as

[TABLE]

for all {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma. Above, ${\mathcal{T}}_{h}$ is a polygonal mesh of meshsize $h$.

Similarly, it holds

[TABLE]

for all ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt},\mathbf{v})\in\Sigma\times U.

Remark 1**.**

As discussed in [12], estimate (3) does not break down for nearly incompressible materials. More precisely, considering the constitutive law:

[TABLE]

with $\lambda,\mu>0$ the Lame’s parameters and ${\rm tr}(\cdot)$ the trace operator, the constant $C$ in (3) can be chosen independent of $\lambda$. The key point is that it is sufficient to check the $\Sigma$-coercivity of the bilinear form $a(\cdot,\cdot)$ in (2) for the subspace:

[TABLE]

In fact, there exists a positive constant $\alpha$ such that (see [12]):

[TABLE]

with $\alpha$ independent of $\lambda$.

3 The Virtual Element Method

We outline the Virtual Element discretization of problem (2). Let $\{\mathcal{T}_{h}\}_{h}$ be a sequence of decompositions of $\Omega$ into general polygonal elements $E$ with

[TABLE]

In what follows, $|E|$ and $|e|=h_{e}$ will denote the area of $E$ and the length of the side $e\in\partial E$, respectively.

We suppose that for all $h$, each element $E$ in $\mathcal{T}_{h}$ fulfils the following assumptions:

• •

$\mathbf{(A1)}$ $E$ is star-shaped with respect to a ball of radius $\geq\,\gamma\,h_{E}$,

• •

$\mathbf{(A2)}$ the distance between any two vertexes of $E$ is $\geq c\,h_{E}$,

where $\gamma$ and $c$ are positive constants. We remark that the hypotheses above, though not too restrictive in many practical cases, can be further relaxed, as noted in  [7].

In addition, we suppose that the tensor $\mathbb{D}$ is piecewise constant with respect to the underlying mesh ${\mathcal{T}}_{h}$.

3.1 The local spaces

Given a polygon $E\in{\mathcal{T}}_{h}$ with $n_{E}$ edges, we first introduce the space of local infinitesimal rigid body motions:

[TABLE]

Here above, given $\mathbf{c}=(c_{1},c_{2})^{T}\in\mathbb{R}^{2}$, $\mathbf{c}^{\perp}$ is the counterclock-wise rotated vector $\mathbf{c}^{\perp}=(c_{2},-c_{1})^{T}$, and $\mathbf{x}_{C}$ is the baricenter of $E$. For each edge $e$ of $\partial E$, we introduce the space

[TABLE]

Here above, $s$ is a local linear coordinate on $e$, such that $s=0$ corresponds to the edge midpoint. Furthermore, $\mathbf{n}$ is the outward normal to the edge $e$. Hence, $R(e)$ consists of vectorial functions which have the edge tangential component constant, and the edge normal component linear along the edge. Our local approximation space for the stress field is then defined by

[TABLE]

Remark 2**.**

Alternatively, the space (11) can be defined as follows.

[TABLE]

Here above, the equation \mathop{\rm curl}\nolimits\mathop{\bf curl}\nolimits(\mathbb{D}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h})=0 is to be intended in the distribution sense.

We remark that, once ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\,\mathbf{n})_{|e}=\mathbf{c}_{e}+d_{e}s\,\mathbf{n} is given for all $e\in\partial E$, cf. (10), the quantity \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in RM(E) is determined. Indeed, denoting with {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}:\partial E\to\mathbb{R}^{2} the function such that {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}_{|e}:=\mathbf{c}_{e}+d_{e}s\,\mathbf{n}, the obvious compatibility condition

[TABLE]

allows to compute \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h} using the $\mathbf{c}_{e}$’s and the $d_{e}$’s. More precisely, setting (cf (9))

[TABLE]

from (13) we infer

[TABLE]

The local approximation space for the displacement field is simply defined by, see (9):

[TABLE]

We notice that $\dim(\Sigma_{h}(E))=3\,n_{E}$, while $\dim(U_{h}(E))=3$.

3.2 The local bilinear forms

Given $E\in{\mathcal{T}}_{h}$, we first notice that, for every {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in\Sigma_{h}(E) and $\mathbf{v}_{h}\in U_{h}(E)$, the term

[TABLE]

is computable from the knowledge of the degrees of freedom. Therefore, there is no need to introduce any approximation in the structure of the terms (\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt},\mathbf{u}) and (\mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{v}) in problem (2). Instead, the term

[TABLE]

is not computable for a general couple ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h},{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h})\in\Sigma_{h}(E)\times\Sigma_{h}(E). As usual in the VEM approach (see [7], for instance), we then need to introduce a suitable approximation $a_{E}^{h}(\cdot,\cdot)$ of $a_{E}(\cdot,\cdot)$. To this end, we first define the projection operator

[TABLE]

Above and in the sequel, given a domain $\omega$ and an integer $k\geq 0$, the space ${\mathcal{P}}_{k}(\omega)$ denotes the polynomials up to degree $k$, defined on $\omega$. Furthermore, given a functional space $X$, $X^{2\times 2}_{s}$ denotes the $2\times 2$ symmetric tensors whose components belong to $X$. Therefore, the operator in (19) is a projection onto the piecewise constant symmetric tensors.

We then set

[TABLE]

where $s_{E}(\cdot,\cdot)$ is a suitable stabilization term. We propose the following choice:

[TABLE]

where $\kappa_{E}$ is a positive constant to be chosen. For instance, in the numerical examples of Section 4, $\kappa_{E}$ is set equal to $\frac{1}{2}{\rm tr}(\mathbb{D}_{|E})$; however, any norm of $\mathbb{D}_{|E}$ can be used. A possible variant of (21) is provided by

[TABLE]

3.3 The local loading terms

We need to consider the term, see (2):

[TABLE]

We remark that, since $\mathbf{v}_{h}\in RM(E)$, computing (23) is possible once a suitable quadrature rule is available for polygonal domains. For such an issue, see for instance [31, 36, 30].

3.4 The discrete scheme

We are now ready to introduce the discrete scheme. We introduce a global approximation space for the stress field, by glueing the local approximation spaces, see (11):

[TABLE]

For the global approximation of the displacement field, we take, see (16):

[TABLE]

Furthermore, given a local approximation of $a_{E}(\cdot,\cdot)$, see (20), we set

[TABLE]

The method we consider is then defined by

[TABLE]

Introducing the bilinear form ${\mathcal{A}}_{h}:(\Sigma_{h}\times U_{H})\times(\Sigma_{h}\times U_{h})\to\mathbb{R}$ defined by

[TABLE]

problem (27) can be written as

[TABLE]

We will prove in Section 5 that our method is first order convergent with respect to the natural norms, see in particular Theorem 5.8. More precisely, the following error estimate holds true.

[TABLE]

where C=C(\Omega,{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u}) is independent of $h$ but depends on the domain $\Omega$ and on the Sobolev regularity of \sigma$$\sigma$$\sigma and $\mathbf{u}$.

4 Numerical results

The present section is devoted to the validation of the proposed methodology through the assessment of accuracy on a selected number of test problems. Applicability to structural analysis is then demonstrated through a classical benchmark.

4.1 Accuracy assessment

We consider two boundary value problems on the unit square domain $\Omega=[0,1]^{2}$, with known analytical solution, discussed in [22, 2]. The material obeys to a homogeneous isotropic constitutive law, see (6), with material parameters assigned in terms of the Lamé constants, here set as $\lambda=1$ and $\mu=1$. Plane strain regime is invoked throughout. The tests are defined by choosing a required solution and deriving the corresponding body load $\mathbf{f}$, as synthetically indicated in the following:

• $\bullet$

Test $a$

[TABLE]

• $\bullet$

Test $b$

[TABLE]

As it can be observed, Test $a$ is a problem with Dirichlet non-homogeneous boundary conditions, zero loading and a polynomial solution; whereas Test $b$ has homogeneous Dirichlet boundary conditions, trigonometric distributed loads with a trigonometric solution.

In order to test the robustness of the proposed procedure with respect to element topology and mesh distortion, six different meshes are considered, as can be inspected in Fig. 1. Three are structured meshes composed of triangles, quadrilaterals and a set of quads, pentagons and hexagons. In the following, such meshes are denoted by the letter "S". Three unstructured meshes are considered as well, comprising triangles, quadrilaterals and random polygons; these are denoted by the letter "U". In the numerical campaign the mesh size parameter is chosen to be the average edge length, denoted with $\bar{h}_{e}$. We remark that, under mesh assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ and for a quasi-uniform family of mesh, $\bar{h}_{e}$ is indeed equivalent to both $h_{E}$ and $h$. The accuracy and the convergence rate assessment is carried out using the following error norms:

• $\bullet$

Discrete error norms for the stress field:

[TABLE]

where $\kappa_{e}=\kappa=\frac{1}{2}{\rm tr}(\mathbb{D})$ (the material is here homogeneous). We remark that the quantity above scales like the internal elastic energy, with respect to the size of the domain and of the elastic coefficients.

We make also use of the $L^{2}$ error on the divergence:

[TABLE]

• $\bullet$

$L^{2}$ error norm for the displacement field:

[TABLE]

Figure 2 reports the $\bar{h}_{e}-$convergence of the proposed method for Test $a$. As expected, the asymptotic convergence rate is approximately equal to $1$ for all the considered error norms and meshes. It is noted that, in this case, the E_{{\kern 0.14003pt\hbox{\sigma}\kern-3.9999pt\kern-0.14003pt\hbox{\sigma}\kern-3.9999pt\raise 0.20999pt\hbox{\sigma}\kern 0.14003pt},\mathop{\bf div}\nolimits} plots are not reported because such a quantity is captured up to machine precision for all the considered computational grids.

Figure 3 reports $\bar{h}_{e}-$convergence for Test $b$. Asymptotic converge rate is approximately equal to $1$ for all investigated mesh types and error measures, including E_{{\kern 0.14003pt\hbox{\sigma}\kern-3.9999pt\kern-0.14003pt\hbox{\sigma}\kern-3.9999pt\raise 0.20999pt\hbox{\sigma}\kern 0.14003pt},\mathop{\bf div}\nolimits}. These results highlight the expected optimal performance of the proposed VEM approach and its robustness with respect to the adopted computational grid.

4.1.1 Nearly incompressibility regime

A problem on the unit square domain $\Omega=[0,1]^{2}$, with known analytical solution, is considered. A nearly incompressible material is chosen by selecting Lamé constants as $\lambda=10^{5}$, $\mu=0.5$. The test is designed by choosing a required solution for the displacement field and deriving the load $\mathbf{f}$ accordingly. The displacement solution is as follows:

[TABLE]

Figure 4 reports the results obtained for both structured and unstructured meshes. In can be clearly seen that the proposed method shows the expected asymptotic rate of convergence also in this case.

4.2 Structural analysis benchmark: Cook’s membrane

The present section deals with the classical Cook’s membrane 2D problem [45]. The geometry of the domain $\Omega$ is presented in Fig. 5 with length data $H_{1}=44$, $H_{2}=16$, $L=48$. The loading is given by a constant tangential traction $q=6.25$ on the right edge of the domain. The Young modulus, $E$, is set equal to $70$ and two Poisson ratios are considered, one corresponding to $\nu=1/3$ and one corresponding to a nearly incompressible case characterized by $\nu=0.499995$.

The problem is solved using three types of meshes: an evenly distributed quadrilateral mesh denoted as Quad, a centroid based Voronoi tessellation, denoted as CVor, and a random based Voronoi tessellation indicated as RVor. An overview of the adopted meshes is reported in Fig. 6.

Convergence results are reported in terms of mesh refinement monitoring $v_{A}$, the vertical displacement of point A (see Fig. 5), approximated as the vertical displacement at the centroid of the closest polygon. In particular, Fig. 7(a) corresponds to the case in which $\nu=1/3$ while Fig. 7(b) reports the results obtained for the nearly incompressible case. The reference solution is indicated with a dotted red line corresponding to an overkilling accurate solution obtained with the hybrid-mixed CPE4I element [38]. In accordance with the results of Section 4.1.1, it can be clearly observed that the proposed formulation is robust with respect to the compressibility parameter, as the convergence behaviour of both cases (a) and (b) is almost the same.

Finally, contours representing the von Mises equivalent stress distributions are reported in Fig. 7. We remark that, inside the polygons, the stress distribution {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h} is not known, but its projection \Pi_{E}{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h} onto the constant tensors is (cf. (19)). Thus, we have used this latter quantity to compute the von Mises equivalent stress displayed in Fig. 7. Finally, the results refer to the case $\nu=1/3$, being the nearly incompressible case extremely similar.

5 Stability and convergence analysis

In this section, we provide a rigorous analysis of the proposed VEM method. For all $E\in{\mathcal{T}}_{h}$, we first introduce the space:

[TABLE]

The global space $\widetilde{\Sigma}$ is defined as

[TABLE]

In the sequel, given a measurable subset $\omega\subseteq\Omega$ and $r>2$, we will use the space

[TABLE]

equipped with the obvious norm. Under our assumptions on the mesh, we recall the following version of the Korn’s inequality:

[TABLE]

Given $\mathbf{v}\in H^{1}(E)^{2}$, the above inequality can be derived by classical results (see [33], for instance), and by choosing $\mathbf{r}_{\mathbf{v}}\in RM(E)$ such that $\int_{E}(\mathbf{v}-\mathbf{r}_{\mathbf{v}})={\bf 0}$.

We will also use the following result.

Lemma 5.1**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. Given $E\in{\mathcal{T}}_{h}$, let $\mathbf{w}\in H^{1}(E)^{2}$ be a solution of the problem:

[TABLE]

where $\mathbf{g}\in L^{2}(E)^{2}$ and $\mathbf{h}\in L^{2}(\partial E)^{2}$ satisfy the compatibility condition

[TABLE]

Then it holds:

[TABLE]

Proof.

For every $\mathbf{r}\in RM(E)$, we have

[TABLE]

by which we get

[TABLE]

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, the Agmon’s inequality then gives

[TABLE]

Estimate (50) now follows from (47).

∎

5.1 An interpolation operator for stresses

We now introduce the local interpolation operator ${\mathcal{I}}_{E}:W^{r}(E)\to\Sigma_{h}(E)$, defined as follows. Given {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in W^{r}(E), {\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma_{h}(E) is determined by:

[TABLE]

where

[TABLE]

If \tau$$\tau$$\tau is not sufficiently regular, the integral in the right-hsnd side of (54) is intended as a duality between $W^{-\frac{1}{r},r}(\partial E)^{2}$ and $W^{\frac{1}{r},r^{\prime}}(\partial E)^{2}$. If \tau$$\tau$$\tau is a regular function, the above condition is equivalent to require:

[TABLE]

The following result shows, in particular, that {\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma_{h}(E) is well-defined by conditions (54).

Lemma 5.2**.**

If {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in\Sigma_{h}(E), then

[TABLE]

imply {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}={\bf 0}.

Proof.

First, recall that for {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in\Sigma_{h}(E) it holds ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\mathbf{n})_{|e}=\mathbf{c}_{e}+d_{e}s\,\mathbf{n} for each edge $e\in\partial E$, cf. (10) and (11). By (57), choosing {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}_{\ast} such that {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}_{\ast|e}={\kern 0.20004pt\hbox{\gamma}\kern-5.1773pt\kern-0.20004pt\hbox{\gamma}\kern-5.1773pt\raise 0.29999pt\hbox{\gamma}\kern 0.20004pt}_{e} for each $e\in\partial E$, it follows that ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\mathbf{n})_{|e}=d_{e}s\,\mathbf{n}. Choosing now {\kern 0.20004pt\hbox{\varphi}\kern-6.54167pt\kern-0.20004pt\hbox{\varphi}\kern-6.54167pt\raise 0.29999pt\hbox{\varphi}\kern 0.20004pt}_{\ast|e}=\delta_{e}(\mathbf{x}-\mathbf{x}_{C})^{\perp}, conditions (57) then give

[TABLE]

A direct computation (for instance by using the Cavalieri-Simpson rule) shows that (58) is equivalent to

[TABLE]

Above, $\mathbf{p}_{e}$ and $\mathbf{q}_{e}$ denote the endpoints of $e$. From (59) we infer $d_{e}=0$ for each $e\in\partial E$, which concludes the proof. ∎

The global interpolation operator ${\mathcal{I}}_{h}:W^{r}(\Omega)\to\Sigma_{h}$ is then defined by simply glueing the local contributions provided by ${\mathcal{I}}_{E}$. More precisely, we set ({\mathcal{I}}_{h}\tau)_{|E}:={\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{|E} for every $E\in{\mathcal{T}}_{h}$ and {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in W^{r}(\Omega).

5.2 Approximation estimates

Proposition 5.3**.**

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, for the interpolation operator ${\mathcal{I}}_{E}$ defined in (56), the following estimates hold:

[TABLE]

[TABLE]

Proof.

Let {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\widetilde{\Sigma}(E)\cap H^{1}(E)^{2\times 2}_{s}, and let $\mathbf{w}\in H^{1}(E)^{2}$ be such that {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}=\mathbb{C}{{\kern 0.20004pt\hbox{\varepsilon}\kern-4.66318pt\kern-0.20004pt\hbox{\varepsilon}\kern-4.66318pt\raise 0.29999pt\hbox{\varepsilon}\kern 0.20004pt}}(\mathbf{w}), see (44). Furthermore, consider {\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}\in\Sigma_{h}(E) and $\mathbf{w}^{\ast}\in H^{1}(E)^{2}$ such that {\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}=\mathbb{C}{{\kern 0.20004pt\hbox{\varepsilon}\kern-4.66318pt\kern-0.20004pt\hbox{\varepsilon}\kern-4.66318pt\raise 0.29999pt\hbox{\varepsilon}\kern 0.20004pt}}(\mathbf{w}^{\ast}), see (11). Hence, setting {\kern 0.20004pt\hbox{\delta}\kern-4.44444pt\kern-0.20004pt\hbox{\delta}\kern-4.44444pt\raise 0.29999pt\hbox{\delta}\kern 0.20004pt}:=(\mathbf{w}-\mathbf{w}^{\ast})\in H^{1}(E)^{2}, it holds:

[TABLE]

Furthermore, using (56), (14) and (15), we infer that {\kern 0.20004pt\hbox{\delta}\kern-4.44444pt\kern-0.20004pt\hbox{\delta}\kern-4.44444pt\raise 0.29999pt\hbox{\delta}\kern 0.20004pt}\in H^{1}(E)^{2} satisfies:

[TABLE]

where $\chi_{e}$ denotes the characteristic function of the edge $e$. Applying Lemma 5.1 with:

[TABLE]

we get

[TABLE]

We now estimate $\mathbf{g}$ and $\mathbf{h}$. We denote respectively with $\Pi_{0,E}$, $\Pi_{RM,E}$ and $\Pi_{0,\partial E}$ the $L^{2}$-projection operators onto the constant functions on $E$, onto the space $RM(E)$ (see (9)), and on the piecewise constant functions on $\partial E$ (with respect to the edge subdivision of $\partial E$).

The divergence theorem and a direct computation show that:

[TABLE]

Therefore, from the first equation of (64), we have

[TABLE]

Noting that ${\mathcal{P}}_{0}(E)^{2}\subset RM(E)$, from the properties of the $L^{2}$ projection operator, we then get

[TABLE]

and

[TABLE]

For the second equation of (64), we remark that:

[TABLE]

Hence, using a standard approximation estimate and a trace inequality, we get

[TABLE]

Taking into account (68) and (71), from (65) we obtain estimate (60).

We now notice that from (62), (63) and (64), we have:

[TABLE]

Then, using (69), we immediately get (61).

∎

5.3 Proving the ellipticity-on-the-kernel condition

We first notice that by (19), (20) and (21), using the techniques of [7, 16], one has:

[TABLE]

We also notice that (see (24), (11) and (25), (16)):

[TABLE]

As a consequence, introducing the discrete kernel $K_{h}\subseteq\Sigma_{h}$:

[TABLE]

we infer that {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in K_{h} implies \mathop{\bf div}\nolimits{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}={\bf 0}. Hence, it holds:

[TABLE]

We are now ready to prove the following ellipticity-on-the-kernel condition.

Proposition 5.4**.**

For the method described in Section 3, there exists a constant $\alpha>0$ such that

[TABLE]

Proof.

By recalling (26), from (73) we get the existence of $\alpha>0$ such that

[TABLE]

Estimate (77) now follows by recalling (76). ∎

Remark 3**.**

Notice that for our method it holds $K_{h}\subset K$, where $K$ is defined by (7). Considering an isotropic material, see (6), from Remark 1 we infer that the coercivity constant $\alpha$ can be chosen independent of $\lambda$. Therefore, our numerical method does not suffer from volumetric locking (see [26], for instance) and can be used also for nearly incompressible materials. This feature is confirmed by the numerical tests presented in Section 4.

5.4 Proving the inf-sup condition

We start by stating the following proposition, which can be derived by regularity results for the elasticity problem on Lipschitz domains (see [25], for example).

Proposition 5.5**.**

Given the polygonal domain $\Omega$, there exist $s>2$ and $\beta^{\ast}>0$ such that

[TABLE]

where $W^{s}(\Omega)$ is the Banach space defined by (46).

We are now ready to prove the discrete inf-sup condition for our choice of the approximation spaces.

Proposition 5.6**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. There exists $\beta>0$ such that

[TABLE]

Proof.

We will apply Fortin’s criterion (see [12]), using the operator ${\mathcal{I}}_{h}:W^{s}(\Omega)\to\Sigma_{h}$, see (56) for the definition of the local contributions. More precisely, we will show that it holds:

[TABLE]

Together with (79), conditions (81) imply (80), see [12].

To prove the first condition in (81), recalling that $\mathbf{v}_{h|E}\in RM(E)$, it is sufficient to show that:

[TABLE]

The above equation directly follows from the divergence theorem and definition (54).

We now prove the continuity estimate (i.e. the second equation in (81)). We will exploit again Lemma 5.1. More precisely, we take $\mathbf{w}^{\ast}\in H^{1}(E)^{2}$ such that {\mathcal{I}}_{E}{\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}=\mathbb{C}{{\kern 0.20004pt\hbox{\varepsilon}\kern-4.66318pt\kern-0.20004pt\hbox{\varepsilon}\kern-4.66318pt\raise 0.29999pt\hbox{\varepsilon}\kern 0.20004pt}}(\mathbf{w}^{\ast}). It follows that $\mathbf{w}^{\ast}$ solves, cf. (66):

[TABLE]

where the $\mathbf{c}_{e}$’s are given by the dualities for the couple $<W^{-\frac{1}{s},s}(\partial E),W^{\frac{1}{s},s^{\prime}}(\partial E)>$:

[TABLE]

From (83) we obviously deduce

[TABLE]

We now apply Lemma 5.1 with:

[TABLE]

and estimate $||\mathbf{h}||_{0,\partial E}$. We start by noting that:

[TABLE]

A duality estimate and a trace bound shows that

[TABLE]

Similarly, it holds:

[TABLE]

From (84), (88) and (89) we get

[TABLE]

by which we deduce, see (87):

[TABLE]

Lemma 5.1 thus gives

[TABLE]

The continuity estimate in (81) now follows by collecting all the local estimates (92).

∎

5.5 Error estimates

We denote with ${\mathcal{P}}_{0}({\mathcal{T}}_{h})$ the space of piecewise constant functions with respect to the given mesh ${\mathcal{T}}_{h}$. We can prove the Proposition:

Proposition 5.7**.**

Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. For every ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I},\mathbf{u}_{I})\in\Sigma_{h}\times U_{h} and every {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{\pi}\in{\mathcal{P}}_{0}({\mathcal{T}}_{h})^{2\times 2}_{s}, the following error equation holds:

[TABLE]

Proof.

Given ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I},\mathbf{u}_{I})\in\Sigma_{h}\times U_{h}, we form ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h}-{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I},\mathbf{u}_{h}-\mathbf{u}_{I})\in\Sigma_{h}\times U_{h}. Then, using the ellipticity-on-the-kernel condition of Proposition 5.4 and the inf-sup condition of Proposition 5.6, there exists ({\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h},\mathbf{v}_{h})\in\Sigma_{h}\times U_{h} such that (see [12] and [14], for instance):

[TABLE]

and

[TABLE]

We have

[TABLE]

Concerning $T_{1}$, it holds:

[TABLE]

We have, using the continuity of $a_{E}(\cdot,\cdot)$ and of $\Pi_{E}$:

[TABLE]

Furthermore, it holds:

[TABLE]

Under assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$, we notice that, given {\kern 0.20004pt\hbox{\tau}\kern-4.37154pt\kern-0.20004pt\hbox{\tau}\kern-4.37154pt\raise 0.29999pt\hbox{\tau}\kern 0.20004pt}_{h}\in\Sigma_{h}(E), we have the 1D inverse estimate on $\partial E$:

[TABLE]

Using the techniques developed in [6], we deduce the scaled trace estimate:

[TABLE]

Hence, we get:

[TABLE]

From (99) and (102) we then deduce

[TABLE]

Since it holds, using also the $L^{2}$ continuity of $\Pi_{E}$:

[TABLE]

Therefore, we get:

[TABLE]

Combining (97), (98) and (105), we infer

[TABLE]

Regarding $T_{2}$, $T_{3}$ and $T_{4}$, one obviously have:

[TABLE]

From (95), (96), (106) and (107), we get:

[TABLE]

Estimate (93) follows from the triangle inequality, estimate (108) and bound (94). ∎

We are now ready to state and prove our main convergence result.

Theorem 5.8**.**

Let ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u})\in\Sigma\times U be the solution of Problem (2), and let ({\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{h},bbu_{h})\in\Sigma_{h}\times U_{h} be the solution of the discrete problem (27). Suppose that assumptions $\mathbf{(A1)}$ and $\mathbf{(A2)}$ are fulfilled. Assuming {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{|E}\in H^{1}(E)^{2\times 2}_{s} and (\mathop{\bf div}\nolimits\,{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt})_{|E}\in H^{1}(E)^{2}, the following estimate holds true:

[TABLE]

where C(\Omega,{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt},\mathbf{u}) is independent of $h$ but depends on the domain $\Omega$ and on the Sobolev regularity of \sigma$$\sigma$$\sigma and $\mathbf{u}$.

Proof.

In Proposition 5.7 let us choose {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{I}={\mathcal{I}}_{h}{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}\in\Sigma_{h} as detailed in Section 5.1, $\mathbf{u}_{I}=P_{0}\mathbf{u}\in U_{h}$ and {\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}_{\pi}=P_{0}{\kern 0.20004pt\hbox{\sigma}\kern-5.71413pt\kern-0.20004pt\hbox{\sigma}\kern-5.71413pt\raise 0.29999pt\hbox{\sigma}\kern 0.20004pt}\in\Sigma_{h}. Estimate (109) easily follows from Proposition 5.3 and standard approximation results. ∎

Remark 4**.**

An alternative way to develop the stability and error analysis might be the use of suitable mesh-dependent norms, as detailed in [29] for the Poisson problem in mixed form.

6 Conclusions

We have proposed, numerically tested and analysed a new Virtual Element Method for the Hellinger-Reissner formulation of two-dimensional elasticity problems. Our scheme is low-order, it has a-priori symmetric stresses and it optimally converges. Possible future developments of the present study include the design of higher-order schemes in the framework of the same variational principle. In addition, accurate post-processed displacements might be considered and used for mesh adaptive strategies, based on suitable a-posteriori error estimators.

Aknowledgements

EA gratefully acknowledges the partial financial support of the Italian Minister of University and Research, MIUR (Program: Consolidate the Foundations 2015; Project: BIOART; Grant number (CUP): E82F16000850005).