A three-dimensional Hellinger-Reissner Virtual Element Method for linear elasticity problems

TL;DR

This paper introduces a 3D Virtual Element Method for linear elasticity problems using the Hellinger-Reissner principle, ensuring symmetric stresses and continuous tractions, with proven convergence and stability.

Contribution

It develops a low-order 3D Virtual Element Method based on the Hellinger-Reissner principle, with theoretical analysis and numerical validation.

Findings

Method achieves convergence and stability.

Numerical tests confirm theoretical predictions.

Ensures symmetric stresses and continuous tractions.

Abstract

We present a Virtual Element Method for the 3D linear elasticity problems, based on Hellinger-Reissner variational principle. In the framework of the small strain theory, we propose a low-order scheme with a-priori symmetric stresses and continuous tractions across element interfaces. A convergence and stability analysis is developed and we confirm the theoretical predictions via some numerical tests.

Figures

Tables

Table 1. Table 1: Example 3 (unloaded body): h ℎ h -convergence results for error E σ σ σ , 𝐝𝐢𝐯 subscript 𝐸 σ σ σ 𝐝𝐢𝐯 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} .

Step	Cube	Tetra	CVT	Random
1	2.1904e-14	4.7821e-14	2.2369e-14	2.9712e-14
2	4.7351e-14	1.2589e-13	4.5054e-14	6.4063e-14
3	1.0024e-13	2.0757e-13	8.7751e-14	1.3174e-13
4	1.1793e-13	2.7652e-13	1.0922e-13	1.6482e-13

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})\leavevmode\nobreak\ \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}\ &\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})\ &\mbox{in $\Omega$}\\ &\mathbf{u}={\bf 0}\ &\mbox{in $\partial\Omega$}\end{aligned}}\right.$$

$$\{\begin{array}{rlr}& \text{mbox}Find(\sigma \sigma \text{raise}0.29999pt\sigma ,\mathbf{u})\in \Sigma \times U\text{leavevmode}\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}$$

$$U={\left[L^2(\Omega )\right]}^3$$

$$\Sigma =\left\{\tau \tau \text{raise}0.29999pt\tau \in H(\mathrm{div};\Omega ):\tau \tau \text{raise}0.29999pt\tau \text{mbox}issymmetric\right\}.$$

$$||\sigma \sigma \text{raise}0.29999pt\sigma |{|}_{\Sigma }^2:={\int }_{\Omega }|\sigma \sigma \text{raise}0.29999pt\sigma {|}^2\text{leavevmode}\text{d}\Omega +{\int }_{\Omega }|\mathrm{div}(\sigma \sigma \text{raise}0.29999pt\sigma ){|}^2\text{leavevmode}\text{d}\Omega ,||\mathbf{u}|{|}_U^2:={\int }_{\Omega }|\mathbf{u}{|}^2\text{leavevmode}\text{d}\Omega ,$$

$$||\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 \text{leavevmode}\text{d}E$$

$$(\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}\text{leavevmode}\text{d}E,$$

$$h:=\underset{E\in {\mathcal{T}}_h}{\sup }{h}_E.$$

$$RM(E):=\left\{\mathbf{r}(\mathbf{x})={\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}+{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\ \ \text{s.t.}\ {\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt},\ {\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\in\mathbb{R}^{3}\right\},$$

$$T_{h}(f):=\left\{{\kern 0.20004pt\hbox{$\psi$}\kern-6.51392pt\kern-0.20004pt\hbox{$\psi$}\kern-6.51392pt\raise 0.29999pt\hbox{$\psi$}\kern 0.20004pt}({\tilde{\mathbf{x}}})=\,\mathbf{t}_{f}+a\big{[}\mathbf{n}_{f}\wedge(\mathbf{x}({\tilde{\mathbf{x}}})-\mathbf{x}_{f})\big{]}+p_{1}({\tilde{\mathbf{x}}})\mathbf{n}_{f},\ \text{s.t.}\ a\in\mathbb{R},\ p_{1}({\tilde{\mathbf{x}}})\in{\mathcal{P}}_{1}(f)\right\},$$

$$\mathbf{t}=b_1{\mathbf{t}}_1+b_2{\mathbf{t}}_2.$$

$$p_1(\tilde{\mathbf{x}})=c_1+c_2(\tilde{x}-{\tilde{x}}_f)+c_3(\tilde{y}-{\tilde{y}}_f).$$

$$\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})_{|f}\in T_{h}(f)\quad\forall f\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{\}}.$$

$$\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}={\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}_{E}+{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)},$$

$${\alpha \alpha \text{raise}0.29999pt\alpha }_E=\frac{1}{|E|}\left(\sum _{f\in \partial E}{\int }_f({\tau \tau \text{raise}0.29999pt\tau }_h\mathbf{n})\text{leavevmode}\text{d}f\right),$$

$$\int_{E}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\leavevmode\nobreak\ \text{d}E=\sum_{f\in\partial E}\int_{f}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge({\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})\leavevmode\nobreak\ \text{d}f.$$

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

$$\int_{E}{\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}_{E}\cdot{\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}\leavevmode\nobreak\ \text{d}E=\int_{\partial E}{\kern 0.20004pt\hbox{$\varphi$}\kern-6.54167pt\kern-0.20004pt\hbox{$\varphi$}\kern-6.54167pt\raise 0.29999pt\hbox{$\varphi$}\kern 0.20004pt}_{|f}\cdot{\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}\leavevmode\nobreak\ \text{d}f=\Big{(}\sum_{f\in\partial E}\int_{f}({\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})\leavevmode\nobreak\ \text{d}f\Big{)}\cdot{\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}\qquad\forall{\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}\in\mathbb{R}^{3}.$$

$${\kern 0.20004pt\hbox{$\alpha$}\kern-6.39702pt\kern-0.20004pt\hbox{$\alpha$}\kern-6.39702pt\raise 0.29999pt\hbox{$\alpha$}\kern 0.20004pt}_{E}=\frac{1}{|E|}\Big{(}\sum_{f\in\partial E}\int_{f}({\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})\leavevmode\nobreak\ \text{d}f\Big{)}.$$

$$\int_{E}\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\cdot\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\leavevmode\nobreak\ \text{d}E=\sum_{f\in\partial E}\int_{f}{\kern 0.20004pt\hbox{$\varphi$}\kern-6.54167pt\kern-0.20004pt\hbox{$\varphi$}\kern-6.54167pt\raise 0.29999pt\hbox{$\varphi$}\kern 0.20004pt}_{|f}\cdot\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\leavevmode\nobreak\ \text{d}f\quad\forall{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\in\mathbb{R}^{3},$$

$${\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\cdot\left(\int_{E}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\leavevmode\nobreak\ \text{d}E\right)={\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\cdot\left(\sum_{f\in\partial E}\int_{f}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge{\kern 0.20004pt\hbox{$\varphi$}\kern-6.54167pt\kern-0.20004pt\hbox{$\varphi$}\kern-6.54167pt\raise 0.29999pt\hbox{$\varphi$}\kern 0.20004pt}_{|f}\leavevmode\nobreak\ \text{d}f\right)\quad\forall{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}\in\mathbb{R}^{3}.$$

$$\int_{E}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\leavevmode\nobreak\ \text{d}E=\sum_{f\in\partial E}\int_{f}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge({\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})\leavevmode\nobreak\ \text{d}f,$$

$${\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\longmapsto\int_{E}\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\wedge\left[{\kern 0.20004pt\hbox{$\omega$}\kern-6.22453pt\kern-0.20004pt\hbox{$\omega$}\kern-6.22453pt\raise 0.29999pt\hbox{$\omega$}\kern 0.20004pt}_{E}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}\right]\text{d}E$$

$$\dim ({\Sigma }_h(E))=6n_f^E.$$

$$U_h(E)=\left\{{\mathbf{v}}_h\in {\left[L^2(E)\right]}^3:{\mathbf{v}}_h\in RM(E)\right\},$$

$$\dim (U_h(E))=6.$$

$${\left(\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h,{\mathbf{v}}_h\right)}_E={\int }_E\mathrm{div}{\tau \tau \text{raise}0.29999pt\tau }_h\cdot {\mathbf{v}}_h\text{leavevmode}\text{d}E$$

$$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\text{leavevmode}\text{d}E$$

$${\int }_E{\Pi }_E{\tau \tau \text{raise}0.29999pt\tau }_h:{\pi \pi \text{raise}0.29999pt\pi }_0={\int }_E{\tau \tau \text{raise}0.29999pt\tau }_h:{\pi \pi \text{raise}0.29999pt\pi }_0\text{leavevmode}\text{d}E\forall {\pi \pi \text{raise}0.29999pt\pi }_0\in {\left[{\mathcal{P}}_0(E)\right]}_s^{3\times 3}.$$

Full text

1 Introduction

The Virtual Element Method (VEM), introduced in [11, 13], is a recent technology for the approximation of partial differential equation problems. This method is a generalization of the Finite Element Method (FEM) which allows to deal with arbitrary polygonal/polyhedral meshes, also including non convex and distored elements. To garantee this flexibility, the virtual element method abandons the idea of the local polynomial approximation, typical of FEM, to use approximating functions which are solution of suitable local PDE. In general, these non-polynomial functions are not explicitly known. Therefore, the main idea of this method is to exploit the available information (the degrees of freedom) to compute the stiffness matrix and the right-hand side of the discretized problem.

During these years, VEM have been employed with success both in mathematical and engineering communities. Here we mention, as a rapresentative non-exaustive sample, a brief list of papers [12, 14, 16, 17, 18, 19, 23, 24, 27, 29, 36, 37]. In the framework of structural mechanics problems, we cite the recent works [1, 5, 6, 7, 8, 9, 41, 42, 43] and [10, 15, 26, 33], for instance. However, we remark that VEM is not the only technology which can make use the polytopal meshes. Considering elasticity problems, we mention [21, 28, 30, 32, 34] as representative examples.

In the present paper we extend the study presented in [7] to the three dimensional case. More precisely, we design and analyze a low-order virtual element method for linear elasticity problems. Within the framework of small displacements and small deformations, we consider the Hellinger-Reissner variational principle as the basis of our discretization procedure. This mixed formulation describes the problem by means of both the displacement and the stress fields. In the Finite Element practice it is a difficult task designing a stable and accurate scheme that preserves both the symmetry of the stress tensor and the continuity of the tractions at the inter-elements, see for instance [20] and [4]. The fundamental reason behind this difficulty lies in the local polynomial approximation, which forces the introduction of nodal degrees of freedom for the stress unknown. The resulting finite element schemes are typically quite cumbersome, especially in the three dimensional case, see [2]. Furthermore, the presence of nodal degrees of freedom introduces an additional complication if one aims at using the hybridization procedure to solve the discrete linear system, see [3]. We exploit the flexibility of virtual element methods to avoid these drawbacks and to develop an optimal scheme which is reasonably cheap with respect to the delivered accuracy. Since the approximated stresses does not have nodal degrees of freedom (on the contrary, the degress of freedom are entirely local to each polyhedron face), the hybridization procedure could be easily applied to our VEM scheme. This aspects show that, even for tetrahedral or hexahedral meshes, the proposed VEM method is a valid alternatives to FEM schemes.

The paper is organized as follows. In Section 2 we briefly introduce the classical Hellinger-Reissner formulation of the 3D elasticity problem. Section 3 describes the Virtual Element approximation we propose, while Section 4 is about the convergence analysis of the method. The numerical experiments, which confirm the theoretical predictions, are detailed in Section 5. Finally, we draw some conclusions.

Space notation.

Throughout the paper, we will make use of standard notations regarding Sobolev spaces, norms and seminorms (cf. [35] for example). In addition, given two quantities $a$ and $b$, we write $a\lesssim b$ when there exists a constant $C$, indipendent of the meshsize, such that $a\leq C\,b$. Finally, given any subset $A\subset\mathbb{R}^{n}$ and an integer $k\geq 0$, ${\mathcal{P}}_{k}(A)$ denotes the space of polynomials up to degree $k$, defined on $A$; whereas, given a functional space $X$, we denote by $\left[X\right]^{3\times 3}_{s}$ the $3\times 3$ symmetric tensors whose components belong to the space $X$.

Mesh notation.

Given a polyhedron $E$ with $n_{f}^{E}$ faces we denote its volume, diameter and barycenter by $|E|$, $h_{E}$ and $\mathbf{x}_{E}$, respectively. In a similar way we refer to the area, diameter and barycenter of a face $f$, while $|e|$ denotes the length of the edge $e$. Given a polygonal face $f$, we use $\mathbf{x}$ and ${\tilde{\mathbf{x}}}$ to indicate the global and local coordinates of a generic point of $f$, respectively.

2 The elasticity problem in mixed form

In this section we briefly present the 3D elasticity problem based on Hellinger-Reissner principle [20, 22]:

[TABLE]

where $\Omega\subseteq\mathbb{R}^{3}$ is a polyhedral domain, \sigma$$\sigma$$\sigma and $\mathbf{u}$ represent the stress and the displacement fields, respectively. We define the bilinear form 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}), where $(\cdot,\cdot)$ is the scalar product in $L^{2}$, and $\mathbb{D}:=\mathbb{C}^{-1}$. Then, the mixed variational formulation reads:

[TABLE]

where $\mathbf{f}\in\left[L^{2}(\Omega)\right]^{3}$ is the loading term, while the spaces $U$ and $\Sigma$ are

[TABLE]

and

[TABLE]

As usual, $H(\mathop{\bf div}\nolimits;\Omega)$ is the space of tensor in $\left[L^{2}(\Omega)\right]^{3\times 3}$ whose divergence is the vector-valued operator in $\left[L^{2}(\Omega)\right]^{3}$. The elasticity fourth-order symmetric tensor $\mathbb{D}$ is assumed to be uniformly bounded and positive-definite.

It is well known that Problem (1) is well posed, cf. [20] for example. In particular, considering the natural norms

[TABLE]

it holds:

[TABLE]

where C is a constant depending on $\Omega$ and on the material tensor $\mathbb{D}$, which, however, does not degenerate in the incompressible limit. We also remark that a possible interesting variant of the variational formulation (1) has been recently proposed and studied in [38] and [39].

When considering a polyhedral mesh ${\mathcal{T}}_{h}$ of the domain $\Omega$, the bilinear form $a(\cdot,\cdot)$ is 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. 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.

3 The Virtual Element Method

In this section we define our Virtual Element discretization of Problem (1). Let $\{\mathcal{T}_{h}\}_{h}$ be a sequence of decompositions of $\Omega$ into general polyhedral elements $E$ with

[TABLE]

We suppose that for all $h$, each element $E$ in $\mathcal{T}_{h}$ fulfils the following standard assumptions (cf. [11]):

• •

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

• •

$\mathbf{(A2)}$ for every face $f\in\partial E$ we have $h_{f}\geq\gamma h_{E}$ and $f$ is star-shaped with respect to a disk of radius $\geq\,\gamma\,h_{f}$,

• •

$\mathbf{(A3)}$ for every edge $e\in\partial f$, we have $|e|\geq\gamma\,h_{f}\geq\gamma^{2}\,h_{E}$,

where $\gamma$ is a suitable positive constant. We also assume that the elasticity tensor $\mathbb{D}$ is piecewise constant with respect to the decomposition ${\mathcal{T}}_{h}$, i.e. $\mathbb{D}$ is constant on each polyhedron of the mesh [7].

3.1 The local spaces

To describe the local spaces employed in the VEM proposed scheme, we introduce two spaces: $RM(E)$ and $T_{h}(f)$.

Space $RM(E)$.

It is the space of local infinitesimal rigid body motions:

[TABLE]

whose dimension is $6$.

Space $T_{h}(f)$.

For each face $f\in\partial E$, we introduce

[TABLE]

where $\mathbf{n}_{f}$ the outward normal to the face $f$, and $\mathbf{t}_{f}$ is an arbitrary vector tangent to the face $f$. Above, $\mathbf{x}({\tilde{\mathbf{x}}})$ is the three dimensional position vector of a point on $f$, determined by the two local coordinates ${\tilde{\mathbf{x}}}$. The dimension of such a space is $6$:

• •

The three dimensional tangent vector $\mathbf{t}$ is determined by a linear combination of two given linearly independent tangential vectors $\mathbf{t}_{1}$ and $\mathbf{t}_{2}$, i.e.

[TABLE]

• •

The rotational term a\big{[}\mathbf{n}_{f}\wedge(\mathbf{x}({\tilde{\mathbf{x}}})-\mathbf{x}_{f})\big{]} is determined by a single scalar value $a\in\mathbb{R}$.

• •

The polynomial $p_{1}({\tilde{\mathbf{x}}})\in{\mathcal{P}}_{1}(f)$ is a two variable polynomial with respect to the local face coordinate system so it is determined by three parameters, for instance:

[TABLE]

Therefore, this space consists of vector functions whose tangential component is a 2D face rigid body motion (the first two terms of (3)), while the normal component is a linear two-variable polynomial (the last term of (3)).

Stress space.

We are now ready to introduce our local approximation space for the stress field:

[TABLE]

We have the following Proposition.

Proposition 3.1**.**

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}_{h}\in\Sigma_{h}(E), then \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} is completely determined by ({\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})_{|f}, with $f\in\partial E$ face of $E$. More precisely, setting (cf (2))

[TABLE]

it holds

[TABLE]

and {\kern 0.20004pt\hbox{\omega}\kern-6.22453pt\kern-0.20004pt\hbox{\omega}\kern-6.22453pt\raise 0.29999pt\hbox{\omega}\kern 0.20004pt}_{E} is the unique solution of the $3\times 3$ linear system

[TABLE]

Proof.

Since {\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 \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) 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}_{h}\,\mathbf{n})_{|f}={\kern 0.20004pt\hbox{\psi}\kern-6.51392pt\kern-0.20004pt\hbox{\psi}\kern-6.51392pt\raise 0.29999pt\hbox{\psi}\kern 0.20004pt}({\tilde{\mathbf{x}}})_{|f}\in T_{h}(f), for $f\in\partial E$. Then, 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}^{3} 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}_{|f}:={\kern 0.20004pt\hbox{\psi}\kern-6.51392pt\kern-0.20004pt\hbox{\psi}\kern-6.51392pt\raise 0.29999pt\hbox{\psi}\kern 0.20004pt}({\tilde{\mathbf{x}}})_{|f}, the integration by parts:

[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 degrees of freedom of the space $T_{h}(f)$.

Testing (7) with constant functions \mathbf{r}(\mathbf{x})={\kern 0.20004pt\hbox{\alpha}\kern-6.39702pt\kern-0.20004pt\hbox{\alpha}\kern-6.39702pt\raise 0.29999pt\hbox{\alpha}\kern 0.20004pt} and recalling (5), we have

[TABLE]

Hence, we obtain

[TABLE]

Now, we test (7) selecting \mathbf{r}(\mathbf{x})={\kern 0.20004pt\hbox{\omega}\kern-6.22453pt\kern-0.20004pt\hbox{\omega}\kern-6.22453pt\raise 0.29999pt\hbox{\omega}\kern 0.20004pt}\wedge\big{(}\mathbf{x}-\mathbf{x}_{E}\big{)}. We have

[TABLE]

i.e.

[TABLE]

We then infer that {\kern 0.20004pt\hbox{\omega}\kern-6.22453pt\kern-0.20004pt\hbox{\omega}\kern-6.22453pt\raise 0.29999pt\hbox{\omega}\kern 0.20004pt}_{E} satisfies

[TABLE]

i.e. {\kern 0.20004pt\hbox{\omega}\kern-6.22453pt\kern-0.20004pt\hbox{\omega}\kern-6.22453pt\raise 0.29999pt\hbox{\omega}\kern 0.20004pt}_{E} solves system (6). We also notice that from (8) we deduce that the linear operator

[TABLE]

is symmetric and positive definite. ∎

From Proposition 3.1 and (3) we infer that the dimension of the space (4) is

[TABLE]

Displacement space.

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

[TABLE]

and it follows that

[TABLE]

3.2 The local forms

In this section we introduce the VEM counterparts of the local forms associated with the continuous problem.

The local mixed term.

Given $E\in{\mathcal{T}}_{h}$, we begin by noticing that the term

[TABLE]

is computable 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)$ via degrees of freedom. For this reason, we do not need to introduce any approximation 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 (1).

The local bilinear form $a_{E}(\cdot,\cdot)$.

The local bilinear form

[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 it is standard in the VEM procedure (cf. [11]), we build a computable approximation of the bilinear form by defining a suitable projection operator onto local polynomial functions. In our case, we introduce $\Pi_{E}:\Sigma_{h}(E)\to\left[{\mathcal{P}}_{0}(E)\right]^{3\times 3}_{s}$, by requiring

[TABLE]

This is a projection operator onto the constant symmetric tensor functions and it is computable. Indeed, we notice that each {\kern 0.20004pt\hbox{\pi}\kern-5.70027pt\kern-0.20004pt\hbox{\pi}\kern-5.70027pt\raise 0.29999pt\hbox{\pi}\kern 0.20004pt}_{0}\in\left[{\mathcal{P}}_{0}(E)\right]^{3\times 3}_{s} can be written as the symmetric gradient of a linear vectorial function, i.e. {\kern 0.20004pt\hbox{\pi}\kern-5.70027pt\kern-0.20004pt\hbox{\pi}\kern-5.70027pt\raise 0.29999pt\hbox{\pi}\kern 0.20004pt}_{0}={{\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{p}_{1}), with $\mathbf{p}_{1}\in\left[{\mathcal{P}}_{1}(E)\right]^{3}$. Hence, using the divergence theorem, the right-hand side of (10) becomes

[TABLE]

which is clearly computable (see also Proposition 3.1). Then, the approximation of $a_{E}(\cdot,\cdot)$ reads:

[TABLE]

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

[TABLE]

Above, $\kappa_{E}$ is a positive constant to be chosen according to $\mathbb{D}$. For instance, in the numerical examples of Section 5, $\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 (12) is

[TABLE]

The local loading term.

We split the load term on each element and we have

[TABLE]

Since $\mathbf{v}_{h}\in RM(E)$, the right-hand side is computable via quadrature rules for polyhedral domains.

Remark 1**.**

Since this integral involves a sufficiently regular function $f$, to get a “good” approximation we exploit a quadrature rule of high degree. For the numerical examples in Section 5, we use a quadrature rule of degree 4.

3.3 The discrete scheme

Starting from the local spaces and local terms introduced in the previous sections, we can set the global problem. More specifically, we introduce a global approximation space for the stress field, by glueing the local approximation spaces, see (4):

[TABLE]

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

[TABLE]

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

[TABLE]

The method we consider is then defined by

[TABLE]

4 Stability and convergence analysis

Since some results of the analysis follows the guidelines of the theory developed in [7], in this section we do not provide full details of all the proofs. First, for all $E\in{\mathcal{T}}_{h}$, we introduce the space:

[TABLE]

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

[TABLE]

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

[TABLE]

equipped with the obvious norm.

4.1 An interpolation operator for stresses

We introduce the local interpolation operator ${\mathcal{I}}_{E}:W^{r}(E)\to\Sigma_{h}(E)$, defined by:

[TABLE]

where:

[TABLE]

We remark that if \tau$$\tau$$\tau is not sufficiently regular, the integral in the right-hand side of (17) is intended as a duality between $\left[W^{-\frac{1}{r},r}(\partial E)\right]^{3}$ and $\left[W^{\frac{1}{r},r^{\prime}}(\partial E)\right]^{3}$. Instead, if \tau$$\tau$$\tau is a regular function, the above condition is equivalent to require:

[TABLE]

for each face $f\in\partial E$.

We now show 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 (17). Indeed this is an immediate consequence of the following Lemma.

Lemma 4.1**.**

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) is such that

[TABLE]

then {\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.

Recalling (4), Proposition 3.1 and (18), it is sufficient to prove that, given a face $f\in\partial E$, conditions

[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}\,\mathbf{n})_{|f}={\bf 0}. To this end, we first set (cf. (3))

[TABLE]

with $\mathbf{t}_{f}$ an arbitrary constant vector tangent to the face $f$, $a\in\mathbb{R}$, $\mathbf{n}_{f}$ the outward normal vector and

[TABLE]

Choosing {\kern 0.20004pt\hbox{\delta}\kern-4.44444pt\kern-0.20004pt\hbox{\delta}\kern-4.44444pt\raise 0.29999pt\hbox{\delta}\kern 0.20004pt}_{f}={\bf 0} and {\kern 0.20004pt\hbox{\gamma}\kern-5.1773pt\kern-0.20004pt\hbox{\gamma}\kern-5.1773pt\raise 0.29999pt\hbox{\gamma}\kern 0.20004pt}_{f} arbitrary in (20), and considering (21), we infer $\mathbf{t}_{f}={\bf 0}$ and $c_{1}=0$. Hence, it holds:

[TABLE]

We now select {\kern 0.20004pt\hbox{\gamma}\kern-5.1773pt\kern-0.20004pt\hbox{\gamma}\kern-5.1773pt\raise 0.29999pt\hbox{\gamma}\kern 0.20004pt}_{f}={\bf 0} and {\kern 0.20004pt\hbox{\delta}\kern-4.44444pt\kern-0.20004pt\hbox{\delta}\kern-4.44444pt\raise 0.29999pt\hbox{\delta}\kern 0.20004pt}_{f}=\mathbf{n}_{f} in (20). From (22) we get

[TABLE]

We have:

[TABLE]

Above, we have used that $\mathbf{n}_{f}\wedge(\mathbf{x}_{f}-\mathbf{x}_{E})$ is a constant vector and that $\mathbf{n}_{f}\wedge(\mathbf{x}({\tilde{\mathbf{x}}})-\mathbf{x}_{f})$ has zero mean value over the face $f$, to infer that:

[TABLE]

From (23) we deduce $a=0$, and therefore we get:

[TABLE]

We finally select {\kern 0.20004pt\hbox{\gamma}\kern-5.1773pt\kern-0.20004pt\hbox{\gamma}\kern-5.1773pt\raise 0.29999pt\hbox{\gamma}\kern 0.20004pt}_{f}={\bf 0} and {\kern 0.20004pt\hbox{\delta}\kern-4.44444pt\kern-0.20004pt\hbox{\delta}\kern-4.44444pt\raise 0.29999pt\hbox{\delta}\kern 0.20004pt}_{f}=\mathbf{t}_{f} in (20), with $\mathbf{t}_{f}$ an arbitrary vector tangential to the face $f$. Thus, we obtain:

[TABLE]

Using again $\mathbf{x}({\tilde{\mathbf{x}}})-\mathbf{x}_{E}=(\mathbf{x}({\tilde{\mathbf{x}}})-\mathbf{x}_{f})+(\mathbf{x}_{f}-\mathbf{x}_{E})$, we infer

[TABLE]

We now choose $\mathbf{t}_{f}$ such that the (tangential to the face) vector $\mathbf{n}_{f}\wedge\mathbf{t}_{f}$ has components $(c_{2},c_{3})$ with respect to the local coordinate system $({\tilde{x}},{\tilde{y}})$. Then, from (25) we get

[TABLE]

which implies $c_{2}=c_{3}=0$. Therefore, ({\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})_{|f}={\bf 0}, see (24), and the proof is complete.

∎

The global interpolation operator ${\mathcal{I}}_{h}:W^{r}(\Omega)\to\Sigma_{h}$ is then defined by glueing the local contributions provided by ${\mathcal{I}}_{E}$. More precisely, 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}\in W^{r}(\Omega) and $E\in{\mathcal{T}}_{h}$, we set

[TABLE]

Using the same steps detailed in [7], we have the following error estimates for the interpolation operator ${\mathcal{I}}_{h}$.

Proposition 4.2**.**

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

[TABLE]

and

[TABLE]

4.2 The ellipticity-on-the-kernel and inf-sup conditions

The proposed approach satisfies the compatibility conditions. First, we notice that (see (13), (4) and (14), (9)):

[TABLE]

Then, we introduce the discrete kernel $K_{h}\subseteq\Sigma_{h}$:

[TABLE]

and we infer from (26) 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]

This is essentially the property that leads to the following ellipticity-on-the-kernel condition.

Proposition 4.3**.**

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

[TABLE]

Proof.

Fix $E\in{\mathcal{T}}_{h}$. By (10), (11) and (12), using the techniques of [11, 25], one has:

[TABLE]

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

[TABLE]

Estimate (28) now follows by recalling (27). ∎

For the discrete inf-sup condition, we need the following commuting diagram property.

Proposition 4.4**.**

For the operator ${\mathcal{I}}_{h}:W^{r}(\Omega)\rightarrow\Sigma_{h}$ it holds:

[TABLE]

where $\Pi_{RM}$ denotes the $L^{2}$-projection operator onto the space of the rigid body motions, see (2).

Proof.

It is sufficient to prove property (29) locally, in each element $E\in{\mathcal{T}}_{h}$. Fix now $\mathbf{r}\in RM(E)$ 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}(E). We have:

[TABLE]

From (30) and the definition of $L^{2}$-projection operator, we get \mathop{\bf div}\nolimits({\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})=\Pi_{RM}(\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}) on $E$. ∎

Using Proposition 4.4 the following discrete inf-sup condition follows from the theory developed in [7].

Proposition 4.5**.**

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

[TABLE]

4.3 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}$. Using the techniques developed in [7], one can prove the following result.

Proposition 4.6**.**

Suppose that assumptions $\mathbf{(A1)}$, $\mathbf{(A2)}$ and $\mathbf{(A3)}$ 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\left[{\mathcal{P}}_{0}({\mathcal{T}}_{h})\right]^{3\times 3}_{s}, the following error equation holds:

[TABLE]

A suitable choice of $\mathbf{u}_{I}$, {\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}, 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} leads to the following error estimate, see [7].

Theorem 4.7**.**

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 (1), 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},\mathbf{u}_{h})\in\Sigma_{h}\times U_{h} be the solution of the discrete problem (16). Suppose that assumptions $\mathbf{(A1)}$, $\mathbf{(A2)}$ and $\mathbf{(A3)}$ 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\left[H^{1}(E)\right]^{3\times 3} 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\left[H^{1}(E)\right]^{3}, 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}$.

5 Numerical results

In this section we numerically assess the proposed VEM approach through the study of the method accuracy on a selected number of test problems. The numerical results confirm the proved theoretical results.

We consider the standard unit cube $\Omega=[0,1]^{3}$ as the domain of our problems and we take the following three types of mesh:

• •

Cube, a mesh composed by standard structured cubes;

• •

Tetra, a Delaunay tetrahedralization of the domain $\Omega$;

• •

CVT, a Voronoi tassellation obtained by the Lloyd algorithm [31];

• •

Random, a Voronoi tassellation achieved with random control points.

We remark that the meshes CVT and Random are very challenging. Indeed, they could have some elements with small faces and edges, and we remark that such case is not covered by the developed theory, i.e. assumptions $\mathbf{(A1)}$, $\mathbf{(A2)}$ and $\mathbf{(A3)}$. However, the numerical results show that the proposed methods are fairly robust with respect to this geometric situation. These two type of meshes are build via the voro++ library [40]. Moreover, the whole numerical scheme is developed inside the vem++ library, a c++ code realized at the Univeristy Milano - Bicocca during the CAVE project (https://sites.google.com/view/vembic/home).

In order to assess the convergence rate, for each type of mesh, we define the following mesh-size $h$:

[TABLE]

where we recall that $N_{E}$ is the number of elements in the mesh, and $h_{E}$ is the diameter of the polyhedron $E$. The accuracy and the convergence rate assessment is carried out using the following error norms:

• $\bullet$

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

[TABLE]

• $\bullet$

$L^{2}$ error on the divergence:

[TABLE]

• $\bullet$

$L^{2}$ error on the projection:

[TABLE]

• $\bullet$

Discrete error norms for the stress field:

[TABLE]

where ${\mathcal{F}}_{h}$ is the set of faces for ${\mathcal{T}}_{h}$ and $\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, i.e. $\sim h$.

Example 1 (compressible material).

We consider an elastic problem with a trigonometric solution and homogeneous Dirichlet boundary conditions. The material of this problems obeys to a homogeneous isotropic constituive law, with material parameters assigned in terms of the Lamé constants, here set as $\lambda=1$ and $\mu=1$. Applied loads are accordingly computed. More precisely, the test details are as follows:

[TABLE]

where $S(x,y,z)=\sin(\pi x)\sin(\pi y)\sin(\pi z)$.

Figure 2 reports $h$-convergence of the proposed methods for Example $1$. As expected, for the considered methods, the asymptotic convergence rate is approximately equal to $1$ for all error norms and meshes. In addition, the convergence graphs of each type of mesh are close to each others and this fact confirms the robustness of the proposed method with respect to the element shape.

Example 2 (nearly incompressible material).

We consider a problem with known analytical solution. 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]

In Figure 3 we report the convergence results for the proposed VEM approach. It can be clearly seen that our method shows the expected asymptotic rate of convergence for each kind of mesh. Moreover, also in this case the convergence lines are close to each other and this fact further confirms the robustness of the proposed scheme with respect to element shape.

Example 3 (unloaded body).

We consider a problem with polynomial solution, non-homogeneous Dirichlet boundary conditions and zero loading. We take a homogeneous and isotropic material with Lamé constants $\lambda=1$ and $\mu=1$ (compressible case). As in the previous examples, the test is defined by choosing a required solution and deriving the correspondig body load $\mathbf{f}$, as indicated in the following:

[TABLE]

We remark that this is a typical example where the displacement field is nontrivial, while stresses are divergence-free. As expected, this latter feature is numerically satisfied by our VEM scheme. Indeed, in Table 1 the errors 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} are close to the machine precision at each mesh refinement step. The other error behaviours are similar to the ones showed in the previous examples so we do not report such graphichs.

6 Conclusions

We have proposed a Virtual Element Method for the linear elasticity 3D problem, based on the mixed Hellinger-Reissner variational principle. The scheme takes advantage of low-order approximation spaces for both the stresses and the displacements. In addition, the stresses are a-priori symmetric and with continuous normal component across the element interfaces. The convergence and stability analysis has been confirmed by some numerical results. A possible future development of the present paper may concern the design of schemes with reduced (minimal) degrees of freedom, by exploiting different variational principles (e.g. suitable augmented lagrangian formulations).