A virtual element method for the vibration problem of Kirchhoff plates

TL;DR

This paper introduces a virtual element method for accurately approximating the vibration spectra of thin Kirchhoff plates on polygonal meshes, with proven optimal error estimates and confirmed through numerical experiments.

Contribution

It develops a simple, conforming VEM scheme for Kirchhoff plate vibrations that achieves optimal spectral approximation and error estimates, extending to non-convex polygons.

Findings

The scheme accurately approximates eigenvalues and eigenfunctions.

Numerical experiments confirm theoretical error estimates.

The method works on various mesh types, including non-convex polygons.

Abstract

The aim of this paper is to develop a virtual element method (VEM) for the vibration problem of thin plates on polygonal meshes. We consider a variational formulation relying only on the transverse displacement of the plate and propose an $H^2(\Omega)$ conforming discretization by means of the VEM which is simple in terms of degrees of freedom and coding aspects. Under standard assumptions on the computational domain, we establish that the resulting schemeprovides a correct approximation of the spectrum and prove optimal order error estimates for the eigenfunctions and a double order for the eigenvalues. The analysis restricts to simply connected polygonal clamped plates, not necessarily convex. Finally, we report several numerical experiments illustrating the behaviour of the proposed scheme and confirming our theoretical results on different families of meshes. Additional examples of cases not covered by our theory are also presented.

Tables

Table 1. Table 1: Lowest vibration frequencies of a simply supported square plate computed on different meshes with the method analyzed in this paper.

	Mesh	$N=32$	$N=64$	$N=128$	Order	Extrapolated	Exact
${\lambda }_1$		390.0184	389.7307	389.6599	2.02	389.6366	389.6364
${\lambda }_2$	${𝒯}_h^1$	2430.2171	2433.9024	2434.8914	1.90	2435.2523	2435.2273
${\lambda }_3$		2430.2171	2433.9024	2434.8914	1.90	2435.2523	2435.2273
${\lambda }_4$		6259.8318	6240.2949	6235.6906	2.09	6234.2872	6234.1818
${\lambda }_1$		389.0957	389.4908	389.5987	1.87	389.6395	389.63634
${\lambda }_2$	${𝒯}_h^2$	2412.1885	2429.0389	2433.6393	1.87	2435.3783	2435.2273
${\lambda }_3$		2433.8095	2434.8277	2435.1197	1.80	2435.2376	2435.2273
${\lambda }_4$		6199.2905	6224.8431	6231.7684	1.88	6234.3634	6234.1818

Table 2. Table 2: Lowest vibration frequencies of a clamped square plate computed on different meshes with the VEM method analyzed in this paper and the one in [ 35 ] .

	Mesh	$N=32$	$N=64$	$N=128$	Order	Extrapolated	[35]
${\lambda }_1$		1283.2286	1291.5607	1294.0225	1.76	1295.0526	1294.9369
${\lambda }_2$	${𝒯}_h^3$	5268.2854	5353.1383	5377.7322	1.79	5387.6973	5386.6675
${\lambda }_3$		5326.6504	5368.5269	5381.6191	1.68	5387.5416	5386.6675
${\lambda }_4$		11406.3068	11622.9583	11686.8981	1.76	11713.7035	11710.9076
${\lambda }_1$		1289.7221	1293.6088	1294.6010	1.97	1294.9410	1294.9369
${\lambda }_2$	${𝒯}_h^4$	5318.4039	5368.9773	5382.1939	1.94	5386.8279	5386.6675
${\lambda }_3$		5351.6510	5377.6743	5384.3950	1.95	5386.7517	5386.6675
${\lambda }_4$		11664.9586	11698.2652	11707.5973	1.84	11711.1942	11710.9076

Table 3. Table 3: Lowest vibration frequencies of an L-shaped clamped plate computed on uniform triangular meshes with the VEM method analyzed in this paper and the one in [ 35 ] .

	$N=32$	$N=64$	$N=128$	Order	Extrapolated	[35]
${\lambda }_1$	6827.5421	6753.6207	6725.1315	1.28	6707.4264	6704.2982
${\lambda }_2$	11128.5787	11073.4576	11059.3867	1.97	11054.5647	11055.5189
${\lambda }_3$	14989.9367	14926.5156	14910.6489	2.00	14905.3676	14907.0816
${\lambda }_4$	26325.7078	26195.9206	26163.4597	2.00	26152.6488	26157.9673

Table 4. Table 4: Lowest vibration frequencies of an L-shaped clamped-free plate computed on triangular meshes with the VEM method analyzed in this paper.

	$N=32$	$N=64$	$N=128$	Order	Extrapolated
${\lambda }_1$	1198.2579	1195.3003	1194.4606	1.82	1194.1302
${\lambda }_2$	4576.4950	4556.9217	4551.0233	1.73	4548.4764
${\lambda }_3$	6807.0921	6785.8226	6780.4745	2.00	6778.6710
${\lambda }_4$	15094.2896	15019.3352	14998.6457	1.86	14990.8077

Equations

$$\left\{\begin{array}[]{ll}\Delta^{2}w=\lambda w&\qquad\text{in }\Omega,\\ w=\partial_{n}w=0&\qquad\text{on }\Gamma,\end{array}\right.$$

$$a(w,v)=\lambda b(w,v)\forall v\in H_0^2(\Omega ),$$

$$a(w,v)$$

$$b(w,v)$$

$$T:H_0^2(\Omega )$$

$$f$$

$$a(u,v)=b(f,v)\forall v\in H_0^2(\Omega ).$$

$$a(v,v)\ge \alpha {\Vert v\Vert }_{2,\Omega }^2\forall v\in H_0^2(\Omega ).$$

$$a(Tf,g)=b(f,g)=b(g,f)=a(Tg,f)=a(f,Tg).$$

$$\Vert u{\Vert }_{2+s,\Omega }\le C\Vert f{\Vert }_{0,\Omega }.$$

$$V_h^K:=\{v_h\in H^2(K):{\Delta }^2{v}_h\in {\mathbb{P}}_2(K),v_h{|}_{\partial K}\in C^0(\partial K),v_h{|}_e\in {\mathbb{P}}_3(e)\forall e\in \partial K,$$

$$\nabla v_h{|}_{\partial K}\in C^0(\partial K{)}^2,{\partial }_n{v}_h{|}_e\in {\mathbb{P}}_1(e)\forall e\in \partial K\},$$

$$a(u,v)=\sum _{K\in {\mathcal{T}}_h}{a}_K(u,v),u,v\in H_0^2(\Omega ),$$

$$b(u,v)=\sum _{K\in {\mathcal{T}}_h}{b}_K(u,v),u,v\in H_0^2(\Omega ),$$

$$a_K(u,v):={\int }_K{D}^{2}u:D^{2}v,u,v\in H^2(K),$$

$$b_K(u,v):={\int }_{K}uv,u,v\in H^2(K).$$

$$\displaystyle a_{K}\big{(}\Pi_{K}^{\Delta}v,q\big{)}$$

$$(({\Pi }_K^{\Delta }v,q){)}_K$$

$$((u,v){)}_K=\sum _{i=1}^{N_K}u(P_i)v(P_i)\forall u,v\in C^0(\partial K),$$

$$W_h^K:=\left\{v_h\in V_h^K:{\int }_K({\Pi }_K^{\Delta }{v}_h)q={\int }_K{v}_{h}q\forall q\in {\mathbb{P}}_2(K)\right\}.$$

$${\int }_K({\Pi }_K^0{v}_h)q={\int }_K{v}_{h}q\forall q\in {\mathbb{P}}_2(K).$$

$${\int }_K({\Pi }_K^0{v}_h)q={\int }_K{v}_{h}q={\int }_K({\Pi }_K^{\Delta }{v}_h)q\forall q\in {\mathbb{P}}_2(K).$$

$$W_h:=\left\{v_h\in H_0^2(\Omega ):v_h{|}_K\in W_h^K\right\}.$$

$$c_0{a}_K(v_h,v_h)\le s_K(v_h,v_h)\le c_1{a}_K(v_h,v_h)$$

$$c_2{b}_K(v_h,v_h)\le s_K^0(v_h,v_h)\le c_3{b}_K(v_h,v_h)$$

$$a_h(u_h,v_h):=\sum _{K\in {\mathcal{T}}_h}{a}_{h,K}(u_h,v_h),u_h,v_h\in W_h,$$

$$b_h(u_h,v_h):=\sum _{K\in {\mathcal{T}}_h}{b}_{h,K}(u_h,v_h),u_h,v_h\in W_h,$$

$$\displaystyle a_{h,K}(u_{h},v_{h}):=a_{K}\big{(}\Pi_{K}^{\Delta}u_{h},\Pi_{K}^{\Delta}v_{h}\big{)}+s_{K}\big{(}u_{h}-\Pi_{K}^{\Delta}u_{h},v_{h}-\Pi_{K}^{\Delta}v_{h}\big{)},\qquad u_{h},v_{h}\in W^{K}_{h},$$

$$\displaystyle b_{h,K}(u_{h},v_{h}):=b_{K}\big{(}\Pi_{K}^{\Delta}u_{h},\Pi_{K}^{\Delta}v_{h}\big{)}+s_{K}^{0}\big{(}u_{h}-\Pi_{K}^{\Delta}u_{h},v_{h}-\Pi_{K}^{\Delta}v_{h}\big{)},\qquad u_{h},v_{h}\in W^{K}_{h}.$$

$$a_{h,K}(p,v_h)=a_K(p,v_h)\forall p\in {\mathbb{P}}_2(K),\forall v_h\in W_h^K,$$

$$b_{h,K}(p,v_h)=b_K(p,v_h)\forall p\in {\mathbb{P}}_2(K),\forall v_h\in W_h^K.$$

$${\alpha }_1{a}_K(v_h,v_h)\le a_{h,K}(v_h,v_h)\le {\alpha }_2{a}_K(v_h,v_h)\forall v_h\in W_h^K,$$

$${\alpha }_3{b}_K(v_h,v_h)\le b_{h,K}(v_h,v_h)\le {\alpha }_4{b}_K(v_h,v_h)\forall v_h\in W_h^K.$$

$$a_h(w_h,v_h)={\lambda }_h{b}_h(w_h,v_h)\forall v_h\in W_h.$$

$$a_h(v_h,v_h)\ge \beta {\Vert v_h\Vert }_{2,\Omega }^2\forall v_h\in W_h.$$

Full text

1 Introduction

The Virtual Element Method (VEM), introduced in [6, 7], is a recent generalization of the Finite Element Method which is characterized by the capability of dealing with very general polygonal/polyhedral meshes. The interest in numerical methods that can make use of general polytopal meshes has recently undergone a significant growth in the mathematical and engineering literature; among the large number of papers on this subject, we cite as a minimal sample [6, 20, 27, 28, 39, 40].

Indeed, polytopal meshes can be very useful for a wide range of reasons, including meshing of the domain (such as cracks) and data (such as inclusions) features, automatic use of hanging nodes, use of moving meshes, adaptivity. Moreover, the VEM presents the advantage to easily implement highly regular discrete spaces. Indeed, by avoiding the explicit construction of the local basis functions, the VEM can easily handle general polygons/polyhedrons without complex integrations on the element (see [7] for details on the coding aspects of the method). The Virtual Element Method has recently been applied successfully to a wide range of problems, see for instance [1, 3, 6, 7, 9, 12, 13, 18, 19, 21, 23, 29, 30, 34, 36, 37, 41, 42, 43].

The numerical approximation of eigenvalue problems for partial differential equations derived from engineering applications, is object of great interest from both, the practical and theoretical points of view. We refer to [14, 15] and the references therein for the state of the art in this subject area. In particular, this paper focus on the so called thin plate vibration problem, which involves the biharmonic operator. Among the existing techniques to solve this problem, various finite element methods have been introduced and analyzed. In particular, we mention nonconforming methods and different mixed formulations for the Kirchhoff model, see for instance [5, 22, 33, 35, 38]. More recently, in [16] a discontinuous Galerkin method has been proposed and analyzed for the vibration and buckling problems of thin plates. On the other hand, the construction of high regularity conforming finite elements for $H^{2}(\Omega)$ is difficult in general, since them usually involve a large number of degrees of freedom (see [24]).

Recently, thanks to the flexibility of VEM, it has been presented in [10, 18] that virtual elements can be used to build global discrete spaces of arbitrary regularity that are simple in terms of degrees of freedom and coding aspects (see also [4, 11]). Thus, in the present contribution, we follow a similar approach in order to solve an eigenvalue problem modelling the two-dimensional plate vibration problem considering a conforming $C^{1}$ discrete formulation.

The aim of this paper is to introduce and analyze a $C^{1}$-VEM which applies to general polygonal meshes, made by possibly non-convex elements, for the two-dimensional plate vibration problem. We begin with a variational formulation of the spectral problem relying only on the transverse displacement of the plate. Then, we exploit the capability of VEM to built highly regular discrete spaces and propose a conforming $H^{2}(\Omega)$ discrete formulation. The resulting discrete bilinear form is continuous and elliptic. This method makes use of a very simple set of degrees of freedom, namely 3 degrees of freedom per vertex of the mesh. By using the abstract spectral approximation theory (see [25, 26]), under rather mild assumptions on the polygonal meshes, we establish that the resulting scheme provides a correct approximation of the spectrum and prove optimal order error estimates for the eigenfunctions and a double order for the eigenvalues. We remark that the present method is new on triangular meshes, and in this case the computational cost is almost $3N_{v}$, where $N_{v}$ denotes the number of vertices, thus providing a very competitive alternative in comparison to other classical techniques based on finite elements.

The outline of this article is as follows: We introduce in Section 2 the variational formulation of the vibration eigenvalue problem, define a solution operator and establish its spectral characterization. In Section 3, we introduce the virtual element discrete formulation, describe the spectrum of a discrete solution operator and prove some auxiliary results. In Section 4, we prove that the numerical scheme provides a correct spectral approximation and establish optimal order error estimates for the eigenvalues and eigenfunctions. Several numerical tests that allow us to assess the convergence properties of the method, to confirm that it is not polluted with spurious modes and to check whether the experimental rates of convergence agree with the theoretical ones are reported in Section 5. Finally, we summarize some conclusions in Section 6.

Throughout the article we will use standard notations for Sobolev spaces, norms and seminorms. Moreover, we will denote by $C$ a generic constant independent of the mesh parameter $h$, which may take different values in different occurrences.

Finally, given a linear bounded operator $T:\,X\to X$, defined on a Hilbert space $X$, we denote its spectrum by $\mathop{\mathrm{sp}}\nolimits(T):=\left\{z\in\mathbb{C}:\,\left(zI-T\right)\mbox{ is not invertible}\right\}$ and by $\rho(T):=\mathbb{C}\setminus\mathop{\mathrm{sp}}\nolimits(T)$ the resolvent set of $T$. Moreover, for any $z\in\rho(T)$, $R_{z}(T):=\left(zI-T\right)^{-1}:\,X\to X$ denotes the resolvent operator of $T$ corresponding to $z$.

2 The spectral problem

Let $\Omega\subset\mathbb{R}^{2}$ be a polygonal bounded simply-connected domain occupied by the mean surface of a plate, clamped on its whole boundary $\Gamma$. The plate is assumed to be homogeneous, isotropic, linearly elastic, and sufficiently thin as to be modeled by Kirchhoff-Love equations. We denote by $w$ the transverse displacement of the mean surface of the plate.

The plate vibration problem reads as follows:

Find $(\lambda,w)\in\mathbb{R}\times{H^{2}(\Omega)}$, $w\neq 0$, such that

[TABLE]

where $\lambda=\omega^{2}$, with $\omega>0$ being the vibration frequency, and $\partial_{n}$ denotes the normal derivative. To simplify the notation we have taken the Young modulus and the density of the plate, both equal to $1$.

To obtain a weak formulation of the spectral problem (2.1), we multiply the corresponding equation by $v\in{H_{0}^{2}(\Omega)}$ and integrate twice by parts in $\Omega$. Thus, we obtain:

Find $(\lambda,w)\in\mathbb{R}\times{H_{0}^{2}(\Omega)}$, $w\neq 0$, such that

[TABLE]

in (2.2) the bilinear forms are defined for any $w,v\in{H_{0}^{2}(\Omega)}$ by

[TABLE]

with $":"$ denotes the usual scalar product of $2\times 2$-matrices, $D^{2}v:=(\partial_{ij}v)_{1\leq i,j\leq 2}$ denotes the Hessian matrix of $v$. We note that those are bounded bilinear symmetric forms. Moreover, it is immediate to prove that the eigenvalues of the problem above are real and positive.

Next, we define the solution operator associated with the variational eigenvalue problem (2.2):

[TABLE]

where $u\in{H_{0}^{2}(\Omega)}$ is the solution of the corresponding source problem:

[TABLE]

The following lemma allows us to establish the well-posedness of this source problem.

Lemma 2.1**.**

There exists a constant $\alpha>0$, depending on $\Omega$, such that

[TABLE]

Proof.

The result follows immediately from the fact that $\|D^{2}v\|_{0,\Omega}$ is a norm on ${H_{0}^{2}(\Omega)}$, equivalent with the usual norm. ∎

We deduce from Lemma 2.1 that the linear operator $T$ is well defined and bounded. Notice that $(\lambda,w)\in\mathbb{R}\times{H_{0}^{2}(\Omega)}$ solves problem (2.2) (and hence problem (2.1)) if and only if $Tw=\mu w$ with $\mu\neq 0$ and $w\neq 0$, in which case $\mu:=\frac{1}{\lambda}$. Moreover, it is easy to check that $T$ is self-adjoint with respect to the $a(\cdot,\cdot)$ inner product. Indeed, given $f,g\in{H_{0}^{2}(\Omega)}$,

[TABLE]

The following is an additional regularity result for the solution of problem (2.3) and consequently, for the eigenfunctions of $T$.

Lemma 2.2**.**

There exist $s\in(\frac{1}{2},1]$ and $C>0$ such that, for all $f\in L^{2}(\Omega)$, the solution $u$ of problem (2.3) satisfies $u\in{H^{2+s}(\Omega)}$ and

[TABLE]

Proof.

The proof follows from the classical regularity result for the biharmonic problem with its right-hand side in $L^{2}(\Omega)$ (cf. [32]). ∎

The constant $s$ in the lemma above is the Sobolev regularity for the biharmonic equation with homogeneous Dirichlet boundary conditions. This constant only depends on the domain $\Omega$. If $\Omega$ is convex, then $s=1$. Otherwise, the lemma holds for all $s<s_{0}$, where $s_{0}\in(\frac{1}{2},1)$ depends on the largest reentrant angle of $\Omega$ (see [32] for the precise equation determining $s_{0}$). Hence, because of the compact inclusion ${H^{2+s}(\Omega)}\hookrightarrow{H_{0}^{2}(\Omega)}$, $T$ is a compact operator. Therefore, we have the following spectral characterization result.

Lemma 2.3**.**

The spectrum of $T$ satisfies $\mathop{\mathrm{sp}}\nolimits(T)=\{0\}\cup\left\{\mu_{k}\right\}_{k\in\mathbb{N}}$, where $\left\{\mu_{k}\right\}_{k\in\mathbb{N}}$ is a sequence of real positive eigenvalues which converges to [math]. The multiplicity of each eigenvalue is finite.

3 Spectral approximation

In this section, first we recall the mesh construction and the assumptions considered to introduce the discrete virtual element spaces. Then, we will introduce a virtual element discretization of the eigenvalue problem (2.2) and provide a spectral characterization of the resulting discrete eigenvalue problem.

Let $\left\{\mathcal{T}_{h}\right\}_{h}$ be a sequence of decompositions of $\Omega$ into polygons $K$. Let $h_{K}$ denote the diameter of the element $K$ and $h$ the maximum of the diameters of all the elements of the mesh, i.e., $h:=\max_{K\in\mathcal{T}_{h}}h_{K}$. In what follows, we denote by $N_{K}$ the number of vertices of $K$.

For the analysis, we will make the following assumptions as in [6, 11]: there exists a positive real number $C_{\mathcal{T}}$ such that, for every $h$ and every $K\in\mathcal{T}_{h}$,

• A1:

the ratio between the shortest edge and the diameter $h_{K}$ of $K$ is larger than $C_{\mathcal{T}}$;

• A2:

$K\in\mathcal{T}_{h}$ is star-shaped with respect to every point of a ball of radius $C_{\mathcal{T}}h_{K}$.

For any subset $S\subseteq\mathbb{R}^{2}$ and nonnegative integer $k$, we indicate by $\mathbb{P}_{k}(S)$ the space of polynomials of degree up to $k$ defined on $S$.

We consider now a simple polygon $K$ (meaning open simply connected sets whose boundary is a non-intersecting line made of a finite number of straight line segments) and we define the following finite-dimensional space

[TABLE]

where $\Delta^{2}$ represents the biharmonic operator and $\partial_{n}$ denotes the normal derivative. We observe that any $v_{h}\in V^{K}_{h}$ satisfy the following conditions:

the trace (and the trace of the gradient) on the boundary of $K$ is continuous;

• 2.

$\mathbb{P}_{2}(K)\subseteq V^{K}_{h}$.

We now introduce two sets ${\bf D_{1}}$ and ${\bf D_{2}}$ of linear operators from $V^{K}_{h}$ into $\mathbb{R}$. For all $v_{h}\in V^{K}_{h}$ they are defined as follows:

${\bf D_{1}}$ contains linear operators evaluating $v_{h}$ at the $N_{K}$ vertices of $K$;

• 2.

${\bf D_{2}}$ contains linear operators evaluating $\nabla v_{h}$ at the $N_{K}$ vertices of $K$.

Note that, as a consequence of definition of $V^{K}_{h}$, the output values of the two sets of operators ${\bf D_{1}}$ and ${\bf D_{2}}$ are sufficient to uniquely determine $v_{h}$ and $\nabla v_{h}$ on the boundary of $K$.

In order to construct the discrete scheme, we need some preliminary definitions. First, we split the bilinear forms $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ introduced in the previous section as follows:

[TABLE]

[TABLE]

with

[TABLE]

and

[TABLE]

Now, we define the projector $\Pi_{K}^{\Delta}:\ V^{K}_{h}\longrightarrow\mathbb{P}_{2}(K)\subseteq V^{K}_{h}$ for each $v\in V^{K}_{h}$ as the solution of

[TABLE]

where $((\cdot,\cdot))_{K}$ is defined as follows:

[TABLE]

where $P_{i},1\leq i\leq N_{K}$, are the vertices of $K$. We note that the bilinear form $a_{K}(\cdot,\cdot)$ has a non-trivial kernel, given by $\mathbb{P}_{1}(K)$. Hence, the role of condition (3.1b) is to select an element of the kernel of the operator.

Now, we introduce our local virtual space:

[TABLE]

It is easy to check that $W^{K}_{h}\subseteq V^{K}_{h}$. Therefore, the operator $\Pi_{K}^{\Delta}$ is well defined on $W^{K}_{h}$ and computable only on the basis of the output values of the operators in ${\bf D_{1}}$ and ${\bf D_{2}}$.

In [4, Lemma 2.1] has been established that the set of operators ${\bf D_{1}}$ and ${\bf D_{2}}$ constitutes a set of degrees of freedom for the space $W^{K}_{h}$. Moreover, it is easy to check that $\mathbb{P}_{2}\subseteq W^{K}_{h}$. This will guarantee the good approximation properties for the space.

Additionaly, we have that the $L^{2}(\Omega)$ projector operator $\Pi_{K}^{0}:W^{K}_{h}\to\mathbb{P}_{2}(K)$ is computable from the set of degrees freedom. In fact, for all $v_{h}\in W^{K}_{h}$, the function $\Pi_{K}^{0}v_{h}\in\mathbb{P}_{2}(K)$ is defined by:

[TABLE]

Now, due to the particular property appearing in definition of the space $W^{K}_{h}$, the right hand side in (3.2) is computable using $\Pi_{K}^{\Delta}v_{h}$, and thus $\Pi_{K}^{0}v_{h}$ depends only on the values of the degrees of freedom for $v_{h}$ and $\nabla v_{h}$. Actually, it is easy to check that on the space $W^{K}_{h}$ the projectors $\Pi_{K}^{0}v_{h}$ and $\Pi_{K}^{\Delta}v_{h}$ are the same operator. In fact:

[TABLE]

In what follows, we keep the notation $\Pi_{K}^{\Delta}$ for both operators.

We can now present the global virtual space: for every decomposition $\mathcal{T}_{h}$ of $\Omega$ into simple polygons $K$, we define

[TABLE]

A set of degrees of freedom for $W_{h}$ is given by all pointwise values of $v_{h}$ on all vertices of $\mathcal{T}_{h}$ together to all pointwise values of $\nabla v_{h}$ on all vertices of $\mathcal{T}_{h}$, excluding the vertices on $\Gamma$ (where the values vanishes). Thus, the dimension of $W_{h}$ is tree times the number of interior vertices.

On the other hand, let $s_{K}(\cdot,\cdot)$ and $s_{K}^{0}(\cdot,\cdot)$ be any symmetric positive definite bilinear forms to be chosen as to satisfy:

[TABLE]

Then, we set

[TABLE]

where $a_{h,K}(\cdot,\cdot)$ and $b_{h,K}(\cdot,\cdot)$ are the local bilinear forms defined on $W^{K}_{h}\times W^{K}_{h}$ by

[TABLE]

The construction of the bilinear forms $a_{h,K}(\cdot,\cdot)$ and $b_{h,K}(\cdot,\cdot)$ guarantees the usual consistency and stability properties of VEM, as noted in the Proposition below. Since the proof follows standard arguments in the Virtual Element literature (see [4, 6, 8]) it is omitted.

Proposition 3.1**.**

The local bilinear forms bilinear forms $a_{h,K}(\cdot,\cdot)$ and $b_{h,K}(\cdot,\cdot)$ on each element $K$ satisfy

Consistency: for all $h>0$ and for all $K\in\mathcal{T}_{h}$ we have that

[TABLE]

• 2.

Stability: There exist positive constants $\alpha_{i},i=1,2,3,4,$ independent of $K$, such that:

[TABLE]

Now, we are in a position to write the virtual element discretization of problem (2.2).

Find $(\lambda_{h},w_{h})\in\mathbb{R}\times W_{h}$, $w_{h}\neq 0$, such that

[TABLE]

We observe that by virtue of (3.10), the bilinear form $a_{h}(\cdot,\cdot)$ is bounded. Moreover, as shown in the following lemma, it is also uniformly elliptic.

Lemma 3.1**.**

There exists a constant $\beta>0$, independent of $h$, such that

[TABLE]

Proof.

Thanks to (3.10) and Lemma 2.1, it is easy to check that the above inequality holds with $\beta:=\alpha\min\left\{\alpha_{1},1\right\}$. ∎

The discrete version of the operator $T$ is then given by

[TABLE]

where $u_{h}\in W_{h}$ is the solution of the corresponding discrete source problem

[TABLE]

Because of Lemma 3.1, the linear operator $T_{h}$ is well defined and bounded uniformly with respect to $h$. Once more, as in the continuous case, $(\lambda_{h},w_{h})\in\mathbb{R}\times W_{h}$ solves problem (3.12) if and only if $T_{h}w_{h}=\mu_{h}w_{h}$ with $\mu_{h}\neq 0$ and $w_{h}\neq 0$, in which case $\mu_{h}:=\frac{1}{\lambda_{h}}$. Moreover, $T_{h}$ is self-adjoint with respect to $a_{h}(\cdot,\cdot)$. Because of this, it is easy to prove the following spectral characterization.

Theorem 3.1**.**

The spectrum of $T_{h}$ consists of $M_{h}:=\mathop{\mathrm{\,dim}}\nolimits(W_{h})$ eigenvalues, repeated according to their respective multiplicities. All of them are real and positive.

In order to prove that the solutions of the discrete problem (3.12) converge to those of the continuous problem (2.2), the standard procedure would be to show that $T_{h}$ converges in norm to $T$ as $h$ goes to zero. However, such a proof does not seem straightforward in our case. In fact, the operator $T_{h}$ is not well defined for any $f\in{H_{0}^{2}(\Omega)}$, since the definition of bilinear form $b_{h,K}(\cdot,\cdot)$ in (3.7) needs the degrees of freedom and in particular the pointwise values of $f$, but it is for any $f\in W_{h}$.

To circumvent this drawback, we will resort instead to the spectral theory from [25] and [26]. In spite of the fact that the main use of this theory is when $T$ is a non-compact operator, it can also be applied to compact $T$ and we will show that in our case it works.

With this aim, we first recall the following approximation result which is derived by interpolation between Sobolev spaces (see for instance [31, Theorem I.1.4] from the analogous result for integer values of $s$). In its turn, the result for integer values is stated in [6, Proposition 4.2] and follows from the classical Scott-Dupont theory (see [17] and [4, Proposition 3.1]):

Proposition 3.2**.**

Let $v\in H^{2+s}(K)$ with $s\in(1/2,1]$. Assume that assumption A2 is satisfied. Then, there exist $v_{\pi}\in\mathbb{P}_{2}(K)$ and $C>0$ such that

[TABLE]

For the analysis we will introduce the broken $H^{2}$-seminorm:

[TABLE]

which is well defined for every $v\in L^{2}(\Omega)$ such that $v|_{K}\in H^{2}(K)$ for all polygon $K\in\mathcal{T}_{h}$.

Now, for $v\in W_{h}$, let $\Pi_{h}$ be defined in $L^{2}(\Omega)$ by $(\Pi_{h}v)|_{K}:=\Pi_{K}^{\Delta}v$ for all $K\in\mathcal{T}_{h}$, where $\Pi_{K}^{\Delta}$ has been defined in (3.1a)-(3.1b).

Lemma 3.2**.**

Let $v\in W_{h}$. Then, there exists $C>0$ such that

[TABLE]

Proof.

Let $v\in W_{h}$. Now, let $\Pi_{K}^{\Delta}v\in\mathbb{P}_{2}(K)$ as defined in (3.1a)-(3.1b). We have for all $r\in\mathbb{P}_{2}(K)$ that

[TABLE]

Thus,

[TABLE]

where we have used (3.3) and [34, Proposition 4.1] and the result follows. ∎

Now, the remainder of this section is devoted to prove the following properties which will be used in the sequel:

Lemma 3.3**.**

There exists $C>0$ such that, for all $f_{h}\in W_{h}$, if $u=Tf_{h}$ and $u_{h}=T_{h}f_{h}$, then

[TABLE]

for all $u_{I}\in W_{h}$ and for all $u_{\pi}\in L^{2}(\Omega)$ such that $u_{\pi}|_{K}\in\mathbb{P}_{2}(K)\quad\forall K\in\mathcal{T}_{h}$.

Proof.

Let $f_{h}\in W_{h}$. For $u_{I}\in W_{h}$, we set $v_{h}:=u_{h}-u_{I}$. Thus

[TABLE]

Now, thanks to Lemma 3.1, the definition of $a_{h,{K}}(\cdot,\cdot)$ and those of $T$ and $T_{h}$, we have

[TABLE]

Then, we bound the first term on the right hand side of the previous inequality as follows

[TABLE]

where we have used the consistency, Cauchy-Schwartz inequality and stability of $b_{h,K}(\cdot,\cdot)$.

Thus, from (3.14), using the above bound together with the Cauchy-Schwartz and triangular inequalities, we obtain

[TABLE]

Therefore, the proof follows from (3.13) and the above inequality. ∎

The next step is to find appropriate term $u_{I}$ that can be used in the above lemma. Thus, we have the following result.

Proposition 3.3**.**

Assume A1–A2 are satisfied, let $v\in{H^{2+s}(\Omega)}$ with $s\in(1/2,1]$. Then, there exist $v_{I}\in W_{h}$ and $C>0$ such that

[TABLE]

Proof.

The proof follows repeating the arguments from [11, Proposition 4.4] (see also [4, Proposition 3.1]). ∎

As we mention before, to prove that $T_{h}$ provides a correct spectral approximation of $T$, we will resort to the theory developed in [25] for noncompact operators. To this end, we first introduce some notations. For any linear operator $S:{H_{0}^{2}(\Omega)}\longrightarrow{H_{0}^{2}(\Omega)}$, we define the norm

[TABLE]

Moreover, we recall the definition of the gap $\widehat{\delta}$ between two closed subspaces $\mathcal{X}$ and $\mathcal{Y}$ of ${H_{0}^{2}(\Omega)}$:

[TABLE]

where

[TABLE]

The theory from [25] guarantees approximation of the spectrum of $T$, provide the following two properties are satisfied:

(P1): $\|T-T_{h}\|_{h}\rightarrow 0,\quad\text{as }h\rightarrow 0$,

• 2.

(P2): $\displaystyle\forall\phi\in{H_{0}^{2}(\Omega)},\quad\lim_{h\rightarrow 0}\delta(\phi,W_{h})=0$.

Property (P2) follows immediately from the approximation property of the virtual element space (see Proposition 3.3) and the density of smooth functions in ${H_{0}^{2}(\Omega)}$. Property (P1) is a consequence of the following lemma.

Lemma 3.4**.**

There exist $C>0$ and $s\in(1/2,1]$, independent of h, such that

[TABLE]

Proof.

Given $f_{h}\in W_{h}$, we have that (see Lemma 3.3)

[TABLE]

now, using Lemma 3.2, Propositions 3.2 and 3.3 and Lemma 2.2, we have

[TABLE]

The proof is complete. ∎

4 Convergence and error estimates

In this section we will adapt the arguments from [25, 26] to prove convergence of our spectral approximation as well as to obtain error estimates for the approximate eigenvalues and eigenfunctions.

The following results are consequence of property (P1) (see [25]):

Lemma 4.1**.**

Suppose that (P1) holds true and let $F\subset\rho(T)$ be closed. Then, there exist positive constants $C$ and $h_{0}$ independent of $h$, such that for $h<h_{0}$

[TABLE]

Theorem 4.1**.**

Let $U\subset\mathbb{C}$ be an open set containing $\mathop{\mathrm{sp}}\nolimits(T)$. Then, there exists $h_{0}>0$ such that $\mathop{\mathrm{sp}}\nolimits(T_{h})\subset U$ for all $h<h_{0}$.

An immediate consequence of this theorem is that the proposed virtual element method does not introduce spurious modes with eigenvalues interspersed among those with a physical meaning.

By applying the results from [25] to our problem, we conclude the spectral convergence of $T_{h}$ to $T$ as $h\to 0$. More precisely, let $\mu\neq 0$ be an isolated eigenvalue of $T$ with multiplicity $m$ and let $\mathcal{C}$ be an open circle in the complex plane centered at $\mu$, such that $\mu$ is the only eigenvalue of $T$ lying in $\mathcal{C}$ and $\partial\mathcal{C}\cap\mathop{\mathrm{sp}}\nolimits(T)=\emptyset$. Then, according to Section 2 in [25] for $h$ small enough there exist $m$ eigenvalues $\mu_{h}^{(1)},\ldots,\mu_{h}^{(m)}$ of $T_{h}$ (repeated according to their respective multiplicities) which lie in $\mathcal{C}$. Therefore, these eigenvalues $\mu_{h}^{(1)},\ldots,\mu_{h}^{(m)}$ converge to $\mu$ as $h$ goes to zero.

The next step is to obtain error estimates for the spectral approximation. With this aim, we will use the theory from [26]. However, we cannot apply the results from this reference directly to our problem, because of the variational crimes in the bilinear forms used to define the operator $T_{h}$. Therefore, we need to extend the results from this reference to our case. With this purpose, we follow an approach recently presented in [12].

Consider the eigenspace $\mathcal{E}$ of $T$ corresponding to $\mu$ and the $T_{h}-$invariant subspace $\mathcal{E}_{h}$ spanned by the eigenspaces of $T_{h}$ corresponding to $\mu_{h}^{(1)},\ldots,\mu_{h}^{(m)}$. As a consequence of Lemma 4.1, we have that

[TABLE]

for $h$ small enough.

Let $\mathcal{P}_{h}:H_{0}^{2}(\Omega)\to W_{h}\hookrightarrow H_{0}^{2}(\Omega)$ be the projector with range $W_{h}$ defined by the relation

[TABLE]

Notice that $\mathcal{P}_{h}$ is bounded uniformly on $h$ (namely $||\mathcal{P}_{h}u||_{2,\Omega}\leq||u||_{2,\Omega}$) and

[TABLE]

Let us define

[TABLE]

Notice that $\mathop{\mathrm{sp}}\nolimits(\widehat{T}_{h})=\mathop{\mathrm{sp}}\nolimits(T_{h})\cup\{0\}$.

Next, we introduce the following spectral projectors (the second one, is well defined at least for $h$ small enough):

The spectral projector of $T$ relative to $\mu$: $F:=\frac{1}{2\pi i}\int_{\partial\mathcal{C}}R_{z}(T)dz$;

• 2.

The spectral projector of $\widehat{T}_{h}$ relative to $\mu_{h}^{(1)},\ldots,\mu_{h}^{(m)}$: $\widehat{F}_{h}:=\frac{1}{2\pi i}\int_{\partial\mathcal{C}}R_{z}(\widehat{T}_{h})dz$.

We have the following result (see Lemma 1 in [26]).

Lemma 4.2**.**

There exist $h_{0}>0$ and $C>0$ such that

[TABLE]

Proof.

It is identical to that of Lemma 11 from [12]. ∎

Consequently, for $h$ small enough, the spectral projector $\widehat{F}_{h}$ is bounded uniformly in $h$.

Now, we define

[TABLE]

From Lemmas 2.2 and 3.2 we have that

[TABLE]

where ${\tilde{s}}\in(1/2,1]$ is such that $\mathcal{E}\subset H^{2+\tilde{s}}(\Omega)$.

The following result establish an error estimate for the eigenfunctions.

Theorem 4.2**.**

If $\mathcal{E}\subset H^{2+\tilde{s}}(\Omega)$ with $\tilde{s}>(1/2,1]$, there exist positive constants $h_{0}$ and $C$ such that, for all $h<h_{0}$,

[TABLE]

Proof.

It follows by arguing exactly as in the proof of Theorem 1 from [26] and using (4.1). ∎

Finally, we have the following result that provides an error estimate for the eigenvalues.

Theorem 4.3**.**

There exist positive constants $C$ and $h_{0}$ independent of $h$, such that, for all $h<h_{0}$,

[TABLE]

where ${\tilde{s}}\in(1/2,1]$ is such that $\mathcal{E}\subset H^{2+\tilde{s}}(\Omega)$.

Proof.

Let $w_{h}$ be such that $(\lambda_{h}^{(i)},w_{h})$ is a solution of (3.12) with $\left\|w_{h}\right\|_{2,\Omega}=1$. According to Theorem 4.2, $\delta(w_{h},\boldsymbol{\mathcal{E}})\leq Ch^{\tilde{s}}$. It follows that there exists $(\lambda,w)$ eigenpair solution of (2.2) such that

[TABLE]

From the symmetry of the bilinear forms and the facts that $a(w,v)=\lambda b(w,v)$ for all $v\in H_{0}^{2}(\Omega)$ (cf. (2.2)) and $a_{h}(w_{h},v_{h})=\lambda_{h}^{(i)}b_{h}(w_{h},v_{h})$ for all $v_{h}\in W_{h}$ (cf. (3.12)), we have

[TABLE]

from which we obtain the following identity:

[TABLE]

The next step is to bound each term on the right hand side above. The first and the second ones are easily bounded using the Cauchy-Schwarz inequality and (4.2):

[TABLE]

For the third term, for $w\in\boldsymbol{\mathcal{E}}$, we consider $w_{\pi}\in L^{2}(\Omega)$ defined on each $K\in\mathcal{T}_{h}$ so that $w_{\pi}|_{K}\in\mathbb{P}_{2}(K)$ and the estimate of Proposition 3.2 holds true. Then, we use (3.8) and (3.10) to write

[TABLE]

Then, adding and subtracting $w$, using triangular inequality, Proposition 3.2, and (4.2), we obtain

[TABLE]

For the last term in (4.3), using that $\Pi_{K}^{\Delta}$ is also the $L^{2}$-projector (see (3.3)), we obtain

[TABLE]

where we have used again Proposition 3.2 and (4.2).

On the other hand,

[TABLE]

thus, the theorem follows from from (4.3)–(4.5) and the inequalities above. ∎

5 Numerical results

We report in this section a couple of tests which have allowed us to assess the theoretical results proved above. With this aim, we have implemented in a MATLAB code the proposed VEM on arbitrary polygonal meshes, by following the ideas presented in [7].

Now, to complete the choice of the VEM, we had to fix the bilinear forms $s_{K}(\cdot,\cdot)$ and $s_{K}^{0}(\cdot,\cdot)$ satisfying (3.4) and (3.5), respectively. Proceeding as in [6], a natural choice for $s_{K}(\cdot,\cdot)$ is given by

[TABLE]

where $P_{1},\ldots,P_{N_{K}}$ are the vertices of $K$, $h_{P_{i}}$ corresponds to the maximun diameter of the elements with $P_{i}$ as a vertex and $\sigma_{K}>0$ is a multiplicative factor to take into account the magnitude of the material parameter, for instance, in the numerical tests a possible choice could be to set $\sigma_{K}>0$ as the mean value of the eigenvalues of the local matrix a_{K}\big{(}\Pi_{K}^{\Delta}u_{h},\Pi_{K}^{\Delta}v_{h}\big{)}. This ensure that the stabilizing term scales as $a_{K}(v_{h},v_{h})$. Now, a choice for $s_{K}^{0}(\cdot,\cdot)$ is given by

[TABLE]

In this case, we have multiplied the stabilizing term by the parameter $\sigma_{K}^{0}>0$, where in this case, a possible choice could be to set $\sigma_{K}^{0}>0$ as the mean value of the eigenvalues of the local matrix b_{K}\big{(}\Pi_{K}^{\Delta}u_{h},\Pi_{K}^{\Delta}v_{h}\big{)} to ensure (3.5). A proof of (3.4) and (3.5) for the above (standard) choices could be derived following the arguments in [8] (see also [4]). Finally, we mention that the above definitions of the bilinear forms $s_{K}(\cdot,\cdot)$ and $s_{K}^{0}(\cdot,\cdot)$ are according with the analysis presented in [34] in order to avoid spectral pollution.

We have tested the method by using different families of meshes (see Figure 1):

$\mathcal{T}_{h}^{1}$: rectangular meshes;

• 2.

$\mathcal{T}_{h}^{2}$: hexagonal meshes;

• 3.

$\mathcal{T}_{h}^{3}$: non-structured hexagonal meshes made of convex hexagons;

• 4.

$\mathcal{T}_{h}^{4}$: trapezoidal meshes which consist of partitions of the domain into $N\times N$ congruent trapezoids, all similar to the trapezoid with vertices $(0,0)$, $(\dfrac{1}{2},0)$, $(\dfrac{1}{2},\dfrac{2}{3})$ and $(0,\dfrac{1}{3})$.

The refinement parameter $N$ used to label each mesh is the number of elements on each edge of the plate.

5.1 Simply supported plate

First, we have considered a simply supported plate, because analytical solutions are available in this case (see [2, 5]). Even though our theoretical analysis has been developed only for clamped plates, we think that the results of the previous sections should hold true for more general boundary conditions as well. The results that follow give some numerical evidence of this. For the computations we took $\Omega:=(0,1)^{2}$.

In Table 1 we report the four lowest eigenvalues ($\lambda_{i},$ $i=1,2,3,4$) computed by our method with two different family of meshes and $N=32,64,128$ for a simply supported plate. The table includes computed orders of convergence, as well as more accurate values extrapolated by means of a least-squares fitting. The last column shows the exact eigenvalues. It can be seen from Table 1 that the method converges to the exact values with an optimal quadratic order. Notice that, for the $\mathcal{T}_{h}^{1}$ meshes, the second computed eigenvalue is double, because the meshes preserve the symmetry of the domain leading to an eigenvalue of multiplicity 2 in the continuous problem.

5.2 Clamped plate

In this numerical test we took $\Omega:=(0,1)^{2}$ and considered clamped boundary condition on the whole of $\partial\Omega$. We present numerical experiments which confirm the theoretical results proved above.

Table 2 shows the four lowest vibration frequencies computed with successively refined meshes of each type for a clamped plate. The table includes orders of convergence, as well as accurate values extrapolated by means of a least-squares fitting. Moreover, we compare the performance of the proposed method with the one presented in [35] with a mixed formulation for solving the plate vibration problem and a Galerkin method based on piecewise linear and continuous finite elements. With this aim, we include in the last column of Table 2 the values obtained by extrapolating those computed with method in [35] on uniform triangular meshes as those shown in Figure 2, for the same problem.

It is clear that the eigenvalue approximation order of our method is quadratic and that the results obtained by the two methods agree perfectly well.

5.3 L-shaped plate

Finally, we present two numerical experiments which confirm the theoretical results proved above. We have computed the vibration frequencies of an L-shaped plate: $\Omega:=(0,1)\times(0,1)\setminus[0.5,1)\times[0.5,1)$.

In the first test we considered clamped boundary condition on the whole of $\partial\Omega$ and we have used uniform triangular meshes as those shown in Figure 3. Once again, we compare the performance of the proposed method with the one presented in [35].

We report in Table 3 the four lowest vibration frequencies computed with the method analyzed in this paper. The table includes orders of convergence, as well as accurate values extrapolated by means of a least-squares fitting. The last column shows the values obtained by extrapolating those computed with method in [35] on the same uniform triangular meshes.

In this case, for the first vibration frequency, the method converges with order close to $1.28$, which is the expected one because of the singularity of the solution (see [32]). Instead, the method converges with larger orders for the second, third and fourth vibration frequencies.

In this case, we mention the following advantages of the proposed VEM method: the computational cost of our method is smaller than the method studied in [35]. In fact, the number of unknowns for our VEM method is, 3$N_{v}$, where $N_{v}$ denotes the number of vertices, whereas in [35] is 4$N_{v}$. Moreover, in this case, the eigenvalue problem to be solved is much simpler than the one arising from the formulation studied in [35]. In fact, the latter leads to a degenerate generalized matrix eigenvalue problem, which is shown to be well posed in [35, Appendix] but that cannot be solved with standard eigensolvers.

We show in Figure 3 the eigenfunctions corresponding to the four lowest eigenvalues for an L-shaped clamped plate.

Finally, Table 4 shows the four lowest vibration frequencies computed with successively refined triangular meshes for an L-shaped clamped-free plate. The table includes orders of convergence, as well as accurate values extrapolated by means of a least-squares fitting. We observe from the results reported in Table 4 that the order of convergence is again quadratic in this case.

Finally, we show in Figure 4 the eigenfunctions corresponding to the four lowest eigenvalues.

6 Conclusions

The mathematical and numerical analysis for the vibration problem of Kirchhoff-Love plates approximation by virtual elements was addressed in this paper. The variational formulation is written in terms of the transverse displacement of the plate and a conforming $H^{2}(\Omega)$ discrete formulation was proposed to numerically approximate the eigenvalue problem. It is established that the resulting scheme provides a correct spectral approximation and that the error estimates are of the optimal order for the eigenfunctions and eigenvalues. The proposed method is new on triangular meshes, and in this case the computational cost is almost $3N_{v}$, where $N_{v}$ denotes the number of vertices, thus providing a very competitive alternative in comparison to other classical techniques based on finite elements. The theoretical results obtained were validated numerically. Even though our theoretical analysis has been developed only for clamped plates, additional examples have been considered and we evidenced the results of the previous sections hold true for more general boundary conditions as well.

Acknowledgments

The first author was partially supported by CONICYT (Chile) through FONDECYT project 1140791 and by DIUBB through project 151408 GI/VC Universidad del Bío-Bío, Chile. The second author was partially supported by BASAL Project PFB 03, CMM, Universidad de Chile, Chile. The third author was partially supported by a CONICYT (Chile) fellowship.