Penalty Virtual Element Method for the 3D Incompressible Flow on Polyhedron Mesh

Abstract

In this paper, a penalty virtual element method (VEM) on polyhedral mesh for solving the 3D incompressible flow is proposed and analyzed. The remarkable feature of VEM is that it does not require an explicit computation of the trial and test space, thereby bypassing the obstacle of standard finite element discretizations on arbitrary mesh. The velocity and pressure are approximated by the practical significative lowest equal-order virtual element space pair $(X_h,Q_h)$ which does not satisfy the discrete $inf$-$sup$ condition. Combined with the penalty method, the error estimation is proved rigorously. Numerical results on the 3D polygonal mesh illustrate the theoretical results and effectiveness of the proposed method.

1. Introduction

In computational fluid dynamics, in order to simulate fluid flow, the incompressible fluid equation is numerically solved. There are classical finite element methods [1] and finite volume methods [2]. It is well known that the traditional finite element method or finite volume method are based on the known interpolation. They are mostly used in structural meshes or quadrilateral meshes, the requirements of meshes are coordinated in the implementation process, so the smoothness of interpolation functions is difficult to improve, which greatly limits the numerical accuracy. The newly developed virtual element method (VEM) [3,4,5] benefits from not requiring a construction of the basis function explicitly to avoid these difficulties. The virtual element space can be constructed by the degrees of freedom of the interior and boundary of the element, and the function space does not need to be expressed explicitly. Therefore, it can be applied to any polygon (polyhedron) with convex or non-convex vertices, and the change of the number of nodes is not needed to recalculate the basis function. Thus, the mesh selection bears great flexibility in the calculation [6,7,8,9].

Due to its flexibility in mesh processing and avoidance of explicit construction of shape functions, VEMs have been successfully applied to a large number of problems. For example, conforming and non-conforming VEMs are presented for elliptic problems in [10,11]. Furthermore, H(div)-conforming and H(curl)-conforming virtual elements are introduced in [12].

In the application of practical problems, the virtual meta method has also made outstanding contributions [13,14,15]. These were developed in [16,17,18] for conforming virtual elements for elasticity problems. Conforming and non-conforming VEMs for fourth-order problems were presented in [19,20,21,22]. $C^1$-continuous VEMs for the Cahn–Hilliard equation were presented in [23]. For the Stokes problem, several VEMs emerged, such as flow VEM formulations [24], scatter-free virtual elements [25] and non-conforming-required VEMs [26]. When constructing the virtual element space, these methods use more complex degrees of freedom to solve the Stokes equation. For example, in [25], in addition to the degrees of freedom we use, the degrees of freedom also need the moment of the normal vector outside the unit on the edge, and our method is relatively simple. By combining with the penalty method, we can solve the problem of velocity pressure compatibility encountered in the solution of Stokes equation. For the discretization of higher order problems, a class of $H^{\alpha }$-conforming ($\alpha >0$) elements was proposed in [27] to satisfy arbitrary regularity requirements. VEMs for parabolic and hyperbolic problems were also developed. Moreover, parabolic and hyperbolic problems were developed in [28,29], respectively. A plane-wave virtual element method for the Helmholtz problem was proposed in the literature [30]. In [31], Mascotto studied the behavior of the stiffness matrix resulting from the approximation of the two-dimensional Poisson problem via the virtual element method. In addition, the SUPG stabilization of the virtual element formulation for advection–diffusion problems was presented in [32], and the $h_p$ virtual element was studied in [33]. In these references, most of the work focused on conforming VEMs with discrete spaces as subspaces of the original space.

In the finite element approximation, a variety of hybrid finite element algorithms which satisfy the velocity pressure $inf$-$sup$ condition have been deeply studied. This condition makes more ideal finite element space unavailable, such as the lowest equal-order $P_1$-$P_1$ finite element pairs. In practical scientific calculation, the equal-order velocity–pressure pairs are very practical owing to their simple data structure, small amount of calculation and high accuracy. In addition, they violate the discrete $inf$-$sup$ condition, that is, the compatibility between velocity and pressure space, leading to pressure oscillation. In order to compensate for the lack of stability of the $inf$-$sup$ condition, several stabilized finite element methods for incompressible flows are proposed, such as the penalty method [1], regularization method [34], rich multi-scale method [18], local Gauss integral method [2], and so on. Therefore, it is natural to combine the virtual element method with the penalty method to solve the incompressible flow problem.

The outline of this paper is as follows. In the next section, we provide the theoretical results and related continuous forms of the penalty Stokes equations. Section 3 introduces the virtual element space and the virtual element method. Then, the stability and error results of the penalty virtual element method are given. Finally, the numerical results are provided to verify the theoretical analysis and the effectiveness of the proposed method.

2. The Stationary Penalty Stokes Equations

Let $\Omega$ be a polygonal domain in ${\mathbb{R}}^3$. This paper considers the 3D incompressible Stokes equations to illustrate the penalty virtual element method:

(1)$\begin{array}{c}\hfill \left\{\begin{array}{cc}-\Delta \mathbf{u}+\nabla p=\mathbf{f},\hfill & \mathrm{in}\Omega ,\hfill \\ \nabla \times \mathbf{u}=0,\hfill & \mathrm{in}\Omega ,\hfill \\ \mathbf{u}=\mathbf{0},\hfill & \mathrm{on}\partial \Omega ,\hfill \end{array}\right.\end{array}$

where $\mathbf{u}$ and p are the velocity field and the pressure field, respectively. $\mathbf{f}$ represents the external force [35].

The penalty method applied to (1) is to approximate the solution $(\mathbf{u},p)$ by $({\mathbf{u}}_{ϵ},p_{ϵ})$ satisfying the following steady-state penalty Stokes Equation [34]

(2)$\begin{array}{c}\hfill \left\{\begin{array}{cc}-\Delta {\mathbf{u}}_{ϵ}+\nabla p_{ϵ}=\mathbf{f},\hfill & \mathrm{in}\Omega ,\hfill \\ \nabla \times {\mathbf{u}}_{ϵ}+ϵp_{ϵ}=0,\hfill & \mathrm{in}\Omega ,\hfill \\ {\mathbf{u}}_{ϵ}=\mathbf{0},\hfill & \mathrm{on}\partial \Omega ,\hfill \end{array}\right.\end{array}$

where the penalty parameter is $0<ϵ<1$.

Consider the Hilbert spaces

$\begin{array}{c}\mathbf{V}:={\left[H_0^1(\Omega )\right]}^3,Q:=L_0^2(\Omega )=\{q\in L^2(\Omega ):{\int }_{\Omega }qd\Omega =0\},\hfill \end{array}$

with norms ${\left|\right|\mathbf{v}\left|\right|}_1:={\left|\right|\mathbf{v}\left|\right|}_{{\left[H_0^1(\Omega )\right]}^3}$, ${\left|\right|q\left|\right|}_Q:={\left|\right|q\left|\right|}_{L^2(\Omega )}$.

The Galerkin variational formula of Equation (1): find $(\mathbf{u},p)\in \mathbf{V}\times Q$, namely

(3)$\begin{array}{c}\hfill \left\{\begin{array}{c}a(\mathbf{u},\mathbf{v})-b(\mathbf{v},p)=(\mathbf{f},\mathbf{v}),\forall \mathbf{v}\in \mathbf{V},\hfill \\ b(\mathbf{u},q)=\mathbf{0},\forall q)\in Q,\hfill \end{array}\right.\end{array}$

where the bilinear forms $a(\times ;\times ):\mathbf{V}\times \mathbf{V}\to \mathbb{R}$ and $b(\times ;\times ):\mathbf{V}\times Q\to \mathbb{R}$ be defined by:

$a(\mathbf{u},\mathbf{v}):={\int }_{\Omega }\nabla \mathbf{u}\nabla \mathbf{v}d\Omega ,\forall \mathbf{u},\mathbf{v}\in \mathbf{V},$

$b(\mathbf{v},q):={\int }_{\Omega }div\mathbf{v}qd\Omega ,\forall \mathbf{v}\in \mathbf{V},q\in Q.$

It is well known that (see for [2,35]):

• $a(\times ,\times )$ and $b(\times ,\times )$ are continuous, i.e.,

  $\left|a(\mathbf{u},\mathbf{v})\right|\le {\left|\right|\mathbf{u}\left|\right|}_1{\left|\right|\mathbf{v}\left|\right|}_1\forall \mathbf{u},\mathbf{v}\in \mathbf{V},$

  $\left|b(\mathbf{v},q)\right|\le {\left|\right|\mathbf{v}\left|\right|}_1{\left|\right|q\left|\right|}_Q\forall \mathbf{v}\in \mathbf{V},and\forall q\in Q,$

• $a(\times ,\times )$ is coercive, i.e., there exists a positive constant $\alpha$ such that

  $a(\mathbf{v},\mathbf{v})\ge {\alpha \left|\right|\mathbf{v}\left|\right|}_1^2,\forall \mathbf{v}\in \mathbf{V};$

• Moreover, the bilinear form $b(\times ,\times )$ satisfies the $inf$-$sup$ condition: where there exists a constant ${\beta }_0>0$ such that

  $\begin{array}{c}\hfill \underset{\mathbf{v}\in \mathbf{V}}{sup}\frac{b(\mathbf{v},q)}{{\left|\right|\mathbf{v}\left|\right|}_1}\ge {\beta }_0{\left|\right|q\left|\right|}_0,\forall q\in Q.\end{array}$

The variational formula of Equation (2): find $({\mathbf{u}}_{ϵ},p_{ϵ})\in \mathbf{V}\times Q$ such that

(4)$\begin{array}{c}\hfill \left\{\begin{array}{c}a({\mathbf{u}}_{ϵ},\mathbf{v})-b(\mathbf{v},p_{ϵ})=(\mathbf{f},\mathbf{v}),\forall \mathbf{v}\in \mathbf{V}\hfill \\ b({\mathbf{u}}_{ϵ},q)+ϵ(p_{ϵ},q)=0,\forall q\in Q.\hfill \end{array}\right.\end{array}$

Theorem 1.There exists at least a solution pair $({\mathit{u}}_{ϵ},p_{ϵ})\in \mathit{V}\times Q$ which satisfies (4) and

(5)$\begin{array}{c}\hfill \left|\right|\nabla {\mathit{u}}_{ϵ}{\left|\right|}_0\le {\left|\right|f\left|\right|}_{-1}.\end{array}$

Moveover, if $f\in L^2(\Omega )$, and ϵ satisfies $ϵC_0\le 1/2$, then the solution $({\mathit{u}}_{ϵ},p_{ϵ})$ of (4) satisfies the following regularity:

(6)$\begin{array}{c}\hfill \left|\right|{\mathit{u}}_{ϵ}{\left|\right|}_2+1/2\left|\right|p_{ϵ}{\left|\right|}_0\le C_0{\parallel f\parallel }_0\end{array}$

where $C_0$ denotes some generic constant depending on the data $(\Omega ,f)$, which may stand for different values at its different occurrences.

There are the following error estimates [1,34].

Lemma 1.

Assume that $\mathbf{f}\in L^2(\Omega )$, $(\mathit{u},p)$ and $({\mathit{u}}_{ϵ},p_{ϵ})$ are the solution of variational formulas (3) and (4), respectively. Then

$\begin{array}{c}\hfill \left|\right|\nabla (\mathit{u}-{\mathit{u}}_{ϵ}){\left|\right|}_0+\left|\right|\mathit{u}-{\mathit{u}}_{ϵ}{\left|\right|}_0+\left|\right|p-p_{ϵ}{\left|\right|}_0\le cϵ{\parallel f\parallel }_{-1}.\end{array}$

Please refer to [1,34] for details of certification.

3. The Penalty Virtual Element Method for Stokes Equations

The model problem (1) where $\Omega$ is (for simplicity) a convex polyhedron. We suppose that we are given a decomposition ${\mathcal{T}}_h$ of $\Omega$ in rather general polyhedra K. We can assume that ${\mathcal{T}}_h$ satisfies the regularity assumption [5]: ∃ a constant $C>0$ assume that

• each polyhedron $K\in {\mathcal{T}}_h$ is star-shaped with respect to every point of a ball of radius $\ge C{h}_K$;

• for every face $f\in \partial K$ we have $h_f\ge C{h}_K$ and f is star-shaped with respect to every point of a disk of radius $\ge C{h}_f$;

• for every edge $e\in \partial f$, we have $\left|e\right|\ge C{h}_K$

where $h_K$ is the diameter of element K, that is, the maximum distance between two nodes of element K; f is a face of element K; $h_f$ is the diameter of the face f; e is the edge length of the element K.

For $k\ge 1$, the same discussion as in [10] reveals that the dimension of ${\mathbf{V}}_h^K$ is

$\begin{array}{c}\hfill N_K^{\mathbf{V}}=\left\{\begin{array}{cc}2nk+k(k-1),\hfill & ford=2,\hfill \\ \frac{3}{2}nk(k+1)+\frac{1}{2}(k-1)k(k+1),\hfill & ford=3.\hfill \end{array}\right.\end{array}$

One of the advantages of the virtual element method is that it does not need to display the construction unit basis function, so one of our most important calculation tools is the definition of scale monomials. The form of scale monomials is given below

$M_k:=\{m_{\alpha }=\frac{{(\mathbf{x}-{\mathbf{x}}_K)}^{\mathbf{\alpha }}}{h_K},\left|\mathbf{\alpha }\right|\le k\}$

where $\mathbf{x}=(x,y,z)$ is a vector, $\mathbf{\alpha }=({\alpha }_1,{\alpha }_2,{\alpha }_3)$ is a triple index, $h_K$ is the element K diameter, and ${\mathbf{x}}_K=(x_K,y_K,z_K)$ is the element K centroid.

Different from the construction of the two-dimensional (2D) virtual element space, in the case of 2D, the lower order polynomial of polygon cannot be obtained. Therefore, we calculate on the edge of the element and extend it to the whole space by using the boundary value problem. Therefore, the 2D virtual element space mainly considers the nodes, edges and interior of the element. In the construction of the 3D virtual element space, we should not only consider the nodes and interior of the element, but also the nodes, ridge and interior of the face on each face

${\mathbf{W}}_k\left(f_K\right):=\{\mathbf{v}\in H^1\left(f_K\right):\mathbf{v}{|}_{\partial f_K}\in {\mathbb{B}}_k(\partial f_K),\Delta \mathbf{v}\in {\mathbb{P}}_{k-2}\left(f_K\right)\},$

where ${\mathbb{B}}_k(\partial f_K):=\{\mathbf{v}\in C^0(\partial f_K):\mathbf{v}{|}_e\in {\mathbb{P}}_k\left(f_K\right),\forall e\subset \partial f_K\}$. Here ${\mathbf{W}}_k\left(f\right)$ is the function space of a face $f_K$ defined in an element K, ${\mathbb{B}}_k(\partial f_K)$ space is defined on the face $f_K$, e is an edge of the face $f_K$, and $\partial f_K$ is all the edges of the face $f_K$. ${\mathbb{P}}_k\left(f_K\right)$ is the set of polynomials on $f_K$ of degree $\le k$.

With the definition of the space of a single face ${\mathbf{W}}_k\left(f\right)$, we give the definition of the virtual meta space of all faces

${\mathbf{W}}_k^{\partial K}\left(K\right):=\{\mathbf{v}\in C^0(\partial K):\mathbf{v}{|}_{f_K}\in {\mathbf{W}}_k\left(f_K\right),\Delta \mathbf{v}\in {\mathbb{P}}_{k-2}\left(K\right)\},$

here, $\partial K$ is the set of all faces of element K.

The definition of local finite element spaces is given by

${\mathbf{W}}_k\left(K\right):=\{\mathbf{v}\in H^1{\left(K\right):\mathbf{v}|}_{\partial K}\in {\mathbf{V}}_k^{\partial K}\left(K\right),\Delta \mathbf{v}{|}_K\in {\mathbb{P}}_{k-2}\left(K\right)\}.$

For each polyhedron K, we define a local finite element space ${\mathbf{V}}_k\left(K\right)$. Roughly speaking, ${\mathbf{V}}_k\left(K\right)$ contains all polynomials of degree k (which is essential for convergence) plus other functions whose restriction on a face is still a polynomial of degree k.

Note that even here a polynomial of degree k satisfies the conditions above, so that ${\mathbb{P}}_k\left(K\right)$ is a subspace of ${\mathbf{V}}_k\left(K\right)$. We can take the following degrees of freedom in ${\mathbf{V}}_k\left(K\right)$:

• the value of ${\mathbf{v}}_h$ at the vertices of K;

• on each edge e, the value of ${\mathbf{v}}_h$ at the $k-1$ internal points of the $(k+1)$-points Gauss–Lobatto quadrature rule on e;

• for each face f the moments up to order $k-2$ of ${\mathbf{v}}_h$ in f:

  ${\int }_f{\mathbf{v}}_h{m}_{\alpha },\alpha =1,\dots ,n_{k-2},$

  where the scaled monomials $m_{\alpha }\in M_{k-2}\left(f\right)$ are defined by

  $m_{\alpha }\left(\mathbf{x}\right):={\left(\frac{\mathbf{x}-{\mathbf{x}}_f}{h_f}\right)}^{\alpha };$

• the moments up to order $k-2$ of ${\mathbf{v}}_h$ in K:

  ${\int }_K{\mathbf{v}}_h{\mu }_{\alpha },\alpha =1,\dots ,{\nu }_{k-2},$

  where the scaled monomials ${\mu }_{\alpha }\in M_{k-2}\left(K\right)$ are defined in (1.5), and ${\nu }_{k-2}$ by

  ${\nu }_{\alpha }\left(\mathbf{x}\right):={\left(\frac{\mathbf{x}-{\mathbf{x}}_K}{h_K}\right)}^{\alpha }.$

For each polyhedron $K\in {\mathcal{T}}_h$ and ${\mathbf{v}}_h\in {\mathbf{W}}_h\left(K\right)$, we can construct ${\Pi }_K^{\nabla }{\mathbf{v}}_h:{\mathbf{W}}_k\left(K\right)\to {\mathbb{P}}_k\left(K\right)$ according to the method in [6], which is defined as follows

$\begin{array}{c}\hfill \left\{\begin{array}{c}{(\nabla {\mathbf{w}}_h,\nabla ({\Pi }_K^{\nabla }{\mathbf{v}}_h-{\mathbf{v}}_h))}_{0,K}=0,\forall {\mathbf{w}}_h\in {\mathbf{W}}_h,\hfill \\ P_0({\Pi }_K^{\nabla }{\mathbf{v}}_h-{\mathbf{v}}_h)=0,\hfill \end{array}\right.\end{array}$

where, $P_0$ is a projection operator of constant functions, that we choose as

(7)$\begin{array}{c}\hfill \begin{array}{c}{P}_0{\mathbf{v}}_h:={\displaystyle \frac{1}{N_K^{\mathbf{V}}}}{\displaystyle \sum _{i=1}^{N_K^{\mathbf{V}}}}{\mathbf{v}}_h\left({\mathbf{V}}_i\right),fork=1,\hfill \end{array}\end{array}$

(8)$\begin{array}{c}\hfill \begin{array}{c}{P}_0{\mathbf{v}}_h:={\displaystyle \frac{1}{\left|K\right|}}{\int }_K{\mathbf{v}}_h,fork\ge 2.\hfill \end{array}\end{array}$

Remark 1.

Note that, integrating by parts, we have again, for every $p_k\in {\mathbb{P}}_k\left(K\right)$,

(9) $\begin{array}{c}\hfill \begin{array}{c}{(\nabla p_k,\nabla {\mathit{v}}_h)}_{0,K}=-{\int }_K▵p_k{\mathit{v}}_h+{\int }_{\partial K}{\displaystyle \frac{\partial p_k}{\partial n}}{\mathit{v}}_h.\hfill \end{array}\end{array}$

Since $▵p_k\in {\mathbb{P}}_{k-2}\left(K\right)$, the first term can again (for $k>1$) be computed using the degrees of freedom (iv). The second term, instead, cannot be computed directly from the degrees of freedom (iii), since on each face f of $\partial K$, ${\displaystyle \frac{\partial p_k}{\partial n}}$ is in ${\mathbb{P}}_{k-1}\left(f\right)$, but the choice of using $v_h{|}_f\in W_k\left(f\right)$ allows us to compute the moments of order k-1 and k as well. It is known in [5]: When constructing $L^2$- projection operator, the shape function on every polygon K is in the space $W_k$. For any $k\ge 1$ and the degree of freedom of given shape function $w_h$, the moments up to order k can be calculated accurately.

However, there only exists the $L^2$-projection on all polyhedrons of degree no more than $(k-2)$ in $W_k\left(K\right)$, which is not enough to deal with more complex situations, e.g., the $L^2$-error estimates and lower-order terms. To this end, the modified virtual element spaces on every face element K.

${\tilde{\mathbf{W}}}_k\left(K\right):=\{v_h\in {\mathbf{W}}_k\left(K\right):{(q^{\ast },v_h-{\Pi }_K^{\nabla }{v}_h)}_{0,K}=0,\forall q^{\ast }\in M_{k-1}^{\ast }\cup M_k^{\ast }\}.$

The $L^2\left(K\right)$-orthogonal projection ${\Pi }_K^0:{\tilde{\mathbf{W}}}_k\left(K\right)\to {\mathbb{P}}_0\left(K\right)$, and for any ${\mathbf{w}}_h\in {\mathbf{V}}_k\left(K\right)$,

${({\mathbf{w}}_h-{\Pi }_K^0{\mathbf{w}}_h,m_k)}_{0,K}=0,\forall m_k\in M_k\left(K\right).$

For each polyhedron K, a local finite element space ${\mathbf{V}}_k\left(K\right)$ is defined. Now, the local virtual element spaces for the velocity and pressure are, respectively, introduced by

$V_k\left(K\right):={\left[{\tilde{\mathbf{W}}}_k\left(K\right)\right]}^3,Q_k\left(K\right):={\tilde{\mathbf{W}}}_k\left(K\right),\forall K\in {\mathcal{T}}_h.$

Then, we can define the global finite element space ${\mathbf{V}}_h$ as

$\begin{array}{cc}\hfill {\mathbf{V}}_h& :=\{{\mathbf{v}}_h\in {\left[H_0^1(\Omega )\right]}^3:{\mathbf{v}}_{{h|}_K}\in {\mathbf{V}}_k\left(K\right)\mathrm{for}\mathrm{all}K\in {\mathcal{T}}_h\},\hfill \\ \hfill Q_h& :=\{q_h\in L_0^2(\Omega ):q{|}_K\in Q_k\forall K\in {\mathcal{T}}_h\}.\hfill \end{array}$

3.1. Constructing the Discrete Matrix

With the definition of $H^1$-projection operator, we first define the local discrete version for (4). The discrete version of $a_h^K(\times ,\times )$ is set as

$a_h^K({\mathbf{u}}_h,{\mathbf{v}}_h)={(\nabla \left({\Pi }_K^{\nabla }{\mathbf{u}}_h\right),\nabla \left({\Pi }_K^{\nabla }{\mathbf{v}}_h\right))}_K+S_K{({\mathbf{u}}_h,{\mathbf{v}}_h)}_K$

$S_K{({\mathbf{u}}_h,{\mathbf{v}}_h)}_K=h_K{(\nabla ({\mathbf{u}}_h-{\Pi }_K^{\nabla }{\mathbf{u}}_h),\nabla ({\mathbf{v}}_h-{\Pi }_K^{\nabla }{\mathbf{v}}_h))}_K$

where $S_K$ is any symmetric, positive definite and computable bilinear form that guarantees [5]:

$c_{\ast }{a}_K({\mathbf{v}}_h,{\mathbf{v}}_h)\le S_K({\mathbf{v}}_h,{\mathbf{v}}_h)\le c^{\ast }{a}_K({\mathbf{v}}_h,{\mathbf{v}}_h),\forall {\mathbf{v}}_h\in ker\left({\Pi }_K^{\nabla }\right).$

Then, it is easy to prove that this discrete local bilinear form satisfies the following consistency and stability assumptions [36]:

• k-compatibility: if $\mathbf{v}\in {\mathbf{V}}_k\left(K\right),\mathbf{s}\in {\mathbb{P}}_k\left(K\right)$, have

  $a_h^K(\mathbf{v},\mathbf{s})=a_K(\mathbf{v},\mathbf{s});$

• Stability: there are two positive constants ${\alpha }_{\ast }$ and ${\alpha }^{\ast }$ dependent on $h_K$ and K, have

  ${\alpha }_{\ast }{a}_K(\mathbf{v},\mathbf{v})\le a_h^K(\mathbf{v},\mathbf{v})\le {\alpha }^{\ast }{a}_K(\mathbf{v},\mathbf{v}),\forall \mathbf{v}\in {\mathbf{V}}_k\left(K\right);$

• Computability: we can know the computability of the discrete bilinear form from Remark 1.

Finally, we define the global approximated bilinear form $a_h(\times ,\times )$ by simply summing the local contributions

$a_h({\mathbf{u}}_h,{\mathbf{v}}_h):=\sum _{K\in {\mathcal{T}}_h}{a}_h^K({\mathbf{u}}_h,{\mathbf{v}}_h),\forall {\mathbf{u}}_h,{\mathbf{v}}_h\in {\mathbf{V}}_h.$

We define projection ${\Pi }_K^{\mathrm{div}}(\nabla \times {\mathbf{v}}_h):{({\Pi }_K^{\mathrm{div}}(\nabla \times {\mathbf{v}}_h),m_{\alpha })}_K={(\nabla \times {\mathbf{v}}_h,m_{\alpha })}_K,\forall m_{\alpha }\in M\left(K\right)$. The bilinear form $b_h^K(p_h,{\mathbf{v}}_h)$ and global contributions $b_h(p_h,{\mathbf{v}}_h)$ are given as

$\begin{array}{cc}\hfill b_h^K(p_h,{\mathbf{v}}_h)& =({\Pi }_K^0{p}_h,{\Pi }_K^{\mathrm{div}}(\nabla \times {\mathbf{v}}_h)),\forall p_h\in Q_h,{\mathbf{v}}_h\in {\mathbf{V}}_h,\hfill \\ \hfill b_h(p_h,{\mathbf{v}}_h)& :=\sum _{K\in {\mathcal{T}}_h}{b}_h^K(p_h,{\mathbf{v}}_h).\hfill \end{array}$

The linear form $\mathbf{f}\left(\mathbf{v}\right)$ on the right-hand side of the variational problem is discrete by ${\mathbf{f}}_h$ such that

$({\mathbf{f}}_h,{\mathbf{v}}_h)=\sum _{K\in {\mathcal{T}}_h}(\mathbf{f},{\Pi }_K^0{\mathbf{v}}_h){)}_{0,K}.$

The discrete problem (2) which we can solve is written as: find ${\mathbf{u}}_h\in {\mathbf{V}}_h$ such that

(10)$\begin{array}{c}\hfill \begin{array}{c}{a}_h({\mathbf{u}}_h,{\mathbf{v}}_h)-b_h(p_h,{\mathbf{v}}_h)+b_h({\mathbf{u}}_h,q_h)+d_h(p_h,q_h)=({\mathbf{f}}_h,{\mathbf{v}}_h),\forall {\mathbf{v}}_h\in {\mathbf{V}}_h,q_h\in Q_h,\hfill \end{array}\end{array}$

where

$d_h(p_h,q_h)=\sum _{K\in {\mathcal{T}}_h}ϵ({\Pi }_K^0{p}_h,{\Pi }_K^0{q}_h)+\sum _{K\in {\mathcal{T}}_h}ϵ((I-{\Pi }_K^0)p_h,(I-{\Pi }_K^0)q_h).$

The first term in $d_h(p_h,q_h)$ is the penalty term, and the second term is the stable term.

3.2. Theoretical Analysis

We begin by proving an approximation result for the virtual local space $V_h$. First of all, let us recall a classical result by Scott–Dupont [37].

Lemma 2.

Let $K\in {\tau }_h$, then for all $\mathit{u}\in {\left[H^{s+1}\left(K\right)\right]}^3$ with $0\le s\le k$, there exists a polynomial function ${\mathit{u}}_{\pi }\in {\left[{\mathbb{P}}_k\left(K\right)\right]}^3$, such that

(11) $\begin{array}{c}\hfill \left|\right|\mathit{u}-{\mathit{u}}_{\pi }{\left|\right|}_{0,K}+h_K|\mathit{u}-{\mathit{u}}_{\pi }{|}_{1,K}\le C{h}_K^{s+1}{\left|\mathit{u}\right|}_{s+1,K}.\end{array}$

We have the following proposition

Proposition 1.

Let $\mathit{u}\in \mathit{V}\cap {\left[H^{s+1}\left(K\right)\right]}^3$ with $0\le s\le k$. Under the mesh assumptions, there exists a polynomial function ${\mathit{u}}_I\in {\mathit{V}}_h$, such that

(12) $\begin{array}{c}\hfill \left|\right|\mathit{u}-{\mathit{u}}_I{\left|\right|}_{0,K}+h_K|\mathit{u}-{\mathit{u}}_I{|}_{1,K}\le C{h}_K^{s+1}{\left|\mathit{u}\right|}_{s+1,D\left(K\right)},\end{array}$

where C is a constant independent of h, and $D\left(K\right)$ denotes the “diamond” of K, i.e., the union of the polygon in ${\mathcal{T}}_h$ intersecting K.

Lemma 3.

Given the discrete spaces $V_h$ and $Q_h$ defined in (2.3) and (2.4), there exists a positive β, independent of h, such that: the bilinear form $b_h(\times ,\times )$ satisfies the discrete $inf$-$sup$ condition, i.e.,

(13) $\begin{array}{c}\hfill \underset{{\mathit{v}}_h\in {\mathit{V}}_h{\mathit{v}}_h\ne 0}{sup}\frac{b({\mathit{u}}_h,q_h)}{\left|\right|{\mathit{v}}_h{\left|\right|}_1}\ge \beta \left|\right|q_h{\left|\right|}_Q\forall q_h\in Q_h.\end{array}$

Proof. 

We only sketch the proof, because it essentially follows the guidelines of Theorem 3.1 in [38]. Since the continuous inf-sup condition (2.6) is fulfilled, it is sufficient to construct a linear operator ${\pi }_h:V\to V_h$, satisfying

(14)$\begin{array}{c}\hfill \left\{\begin{array}{c}b({\pi }_h\mathbf{v},q_h)=b({\mathbf{v}}_h,q_h),\forall \mathbf{v}\in \mathbf{V},\forall q_h\in Q_h,\hfill \\ \left|\right|{\pi }_h\mathbf{v}\left|\right|\le C_{\pi }{\left|\right|\mathbf{v}\left|\right|}_1,\forall \mathbf{v}\in \mathbf{V},\hfill \end{array}\right.\end{array}$

where $C_{\pi }$ is a positive h-independent constant. Given $v\in V$, using arguments borrowed from [38] and considering the VEM interpolant $v_I$ presented in Proposition (4.1), we first construct ${\overline{\mathbf{v}}}_h\in V_h$ such that

$b(\mathbf{v}-{\overline{\mathbf{v}}}_h,{\overline{q}}_h)=0,\forall q_h\in Q_h,$

and

$\left|\right|\mathbf{v}-{\overline{\mathbf{v}}}_h{\left|\right|}_1\le {C\left|\right|\mathbf{v}\left|\right|}_1.$

□

Thus, we can directly have the following two lemmas with analogous proof as in [37], which are omitted here.

Lemma 4.

There exists a positive constant C such that, for all $E\in {\mathcal{T}}_h$ and all smooth enough functions φ defined on K, it holds:

$\left|\right|\phi -{\Pi }_K^{\nabla }\phi \left|\right|\le C{h}_K^{s-m}{\left|\phi \right|}_{s,K},s,m\in \mathbb{N},m\le s\le k+1,s\ge 1,$

$\left|\right|\phi -{\Pi }_K^0\phi \left|\right|\le C{h}_K^{s-m}{\left|\phi \right|}_{s,K},s,m\in \mathbb{N},m\le s\le k+1,$

$\left|\right|\phi -{\phi }_I{\left|\right|}_{m,K}\le C{h}_K^{s-m}{\left|\phi \right|}_{s,E},m,s\in N,m\le s\le k+1,s\ge 2.$

Define a norm as $\left|\right||({\mathbf{u}}_h,p_h){\left|\right||}_h^2:=\left|\right|{\mathbf{u}}_h{\left|\right|}_1^2+\left|\right|p_h{\left|\right|}_0^2+d_h(p_h,p_h)$. We have the following theorem.

Theorem 2.

Assume that the mesh assumption is satisfied and $h=max\left\{h_K\right\}$, $\forall K\in {\mathcal{T}}_h$. Let $({\mathit{u}}_{ϵ},p_{ϵ})$ and $({\mathit{u}}_h,p_h)$ is solution of (4) and (10), respectively. Then, it holds

$\begin{array}{cc}\hfill \left|\right|{\mathit{u}}_{ϵ}-{\mathit{u}}_h{\left|\right|}_1+\left|\right|p_{ϵ}-p_h{\left|\right|}_0& \le C\left(\right||{\mathit{u}}_{ϵ}-{\mathit{u}}_I|{|}_1+\left|\right|p_{ϵ}-p_I|{|}_0+||\nabla \times {\mathit{u}}_{ϵ}-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathit{u}}_{ϵ}|{|}_0\hfill \\ \hfill & +|{\mathit{u}}_{ϵ}-{\Pi }_K^0{\mathit{u}}_{ϵ}{|}_{1,h}+\left|\right|p_{ϵ}-{\Pi }_K^{\mathrm{div}}{p}_{ϵ}{\left|\right|}_0+\left|\right|\mathit{f}-{\mathit{f}}_h{\left|\right|}_0).\hfill \end{array}$

Proof. 

Combined with the definition of projection operator ${\Pi }_K^{\nabla }$, ${\Pi }_K^0$, ${\Pi }_K^{div}$. When $k=1$, ${\Pi }_K^{\nabla }={\Pi }_K^0$, ${\Pi }_K^{div}$ is the component of ${\Pi }_K^{\nabla }$ (references [4,5]). Firstly, we simply deal with ${\mathbf{u}}_{ϵ}-{\mathbf{u}}_h$ and $p_{ϵ}-p_h$,

${\mathbf{u}}_{ϵ}-{\mathbf{u}}_h={\mathbf{u}}_{ϵ}-{\mathbf{u}}_I+{\mathbf{u}}_I-{\mathbf{u}}_h={\eta }_{\mathbf{u}}+{\xi }_{\mathbf{u}},$

$p_{ϵ}-p_h=p_{ϵ}-p_I+p_I-p_h={\eta }_p+{\xi }_p.$

Next, using the above notation and k-consistency, we obtain

(15)$\begin{array}{cc}\hfill a_h({\xi }_{\mathbf{u}},{\mathbf{v}}_h)& =\sum _{K\in {\mathcal{T}}_h}{a}_h^K({\mathbf{u}}_I-{\mathbf{u}}_h,{\mathbf{v}}_h)\hfill \\ \hfill & =\sum _{K\in {\mathcal{T}}_h}[a_h^K({\mathbf{u}}_I-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)+a_h^K({\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)]-a_h({\mathbf{u}}_h,{\mathbf{v}}_h)\hfill \\ \hfill & =\sum _{K\in {\mathcal{T}}_h}[a_h^K({\mathbf{u}}_I-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)+a^K({\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)]+a({\mathbf{u}}_{ϵ},{\mathbf{v}}_h)-a_h({\mathbf{u}}_h,{\mathbf{v}}_h)\hfill \\ \hfill & =-\sum _{K\in {\mathcal{T}}_h}{a}_h^K({\eta }_{\mathbf{u}},{\mathbf{v}}_h)+\sum _{K\in {\mathcal{T}}_h}{a}_h^K({\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)\hfill \\ \hfill & -\sum _{K\in {\mathcal{T}}_h}{a}^K({\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)+a({\mathbf{u}}_{ϵ},{\mathbf{v}}_h)-a_h({\mathbf{u}}_h,{\mathbf{v}}_h).\hfill \end{array}$

According to Lemmas 3 and 4, we have

(16)$\begin{array}{cc}\hfill -b_h({\mathbf{v}}_h,{\xi }_p)& =-\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,{\Pi }_K^0(p_I-p_h))}_{0,K}\hfill \\ \hfill & =-\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,p_I-p_h)}_{0,K}\hfill \\ \hfill & =-\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,{\eta }_p)}_{0,K}+\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\mathbf{v}}_h-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,p_{ϵ})}_{0,K}\hfill \\ \hfill & -(\nabla \times {\mathbf{v}}_h,p_{ϵ})+\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,p_h)}_{0,K}\hfill \\ \hfill & =-\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,{\eta }_p)}_{0,K}+\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\mathbf{v}}_h-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,p_{ϵ}-{\Pi }_K^{\mathrm{div}}{p}_{ϵ})}_{0,K}\hfill \\ \hfill & -(\nabla \times {\mathbf{v}}_h,p_{ϵ})+\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,{\Pi }_K^0{p}_h)}_{0,K}\hfill \\ \hfill & =\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,p_{ϵ})}_{0,K}+\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\mathbf{v}}_h,p_{ϵ}-{\Pi }_K^{\mathrm{div}}{p}_{ϵ})}_{0,K}\hfill \\ \hfill & -b({\mathbf{v}}_h,p_{ϵ})+b_h({\mathbf{v}}_h,p_h).\hfill \\ \hfill \end{array}$

Combined with Equation (16), we can obtain the estimation formula of $b_h({\xi }_{\mathbf{u}},q_h)$

(17)$\begin{array}{cc}\hfill b_h({\xi }_{\mathbf{u}},q_h)& =-\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\eta }_{\mathbf{u}},{\Pi }_K^{\mathrm{div}}{q}_h)}_{0,K}+b({\mathbf{u}}_{ϵ},q_h)\hfill \\ \hfill & -\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\mathbf{u}}_{ϵ}-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{u}}_{ϵ},q_h)}_{0,K}-b_h({\mathbf{u}}_h,q_h).\hfill \\ \hfill \end{array}$

Similarly, we can do the same with $d_h({\xi }_p,q_h)$

(18)$\begin{array}{c}\hfill d_h({\xi }_p,q_h)=d_h(p_I,q_h)-d_h(p_h,q_h).\end{array}$

Add Equations (15)–(18), combine Lemma 2 and Proposition 1 in [4], we have

$\begin{array}{cc}\hfill & a_h({\xi }_{\mathbf{u}},{\mathbf{v}}_h)-b_h({\mathbf{v}}_h,{\xi }_p)+b_h({\xi }_{\mathbf{u}},q_h)+d_h({\xi }_p,q_h)\hfill \\ \hfill & =-\sum _{K\in {\mathcal{T}}_h}{a}_h^K({\eta }_{\mathbf{u}},{\mathbf{v}}_h)+\sum _{K\in {\mathcal{T}}_h}{a}_h^K({\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)-\sum _{K\in {\mathcal{T}}_h}{a}^K({\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ},{\mathbf{v}}_h)\hfill \\ \hfill & +\sum _{K\in {\mathcal{T}}_h}{({\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{v}}_h,p_{ϵ})}_{0,K}+(\mathbf{f}-{\mathbf{f}}_h,{\mathbf{v}}_h)+\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\mathbf{v}}_h,p_{ϵ}-{\Pi }_K^{\mathrm{div}}{p}_{ϵ})}_{0,K}\hfill \\ \hfill & -\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\eta }_{\mathbf{u}},{\Pi }_K^{\mathrm{div}}{q}_h)}_{0,K}-\sum _{K\in {\mathcal{T}}_h}{(\nabla \times {\mathbf{u}}_{ϵ}-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{u}}_{ϵ},q_h)}_{0,K}+d_h(p_I,q_h)\hfill \\ \hfill & \le C\left(\right|\left|{\eta }_{\mathbf{u}}\right|{|}_1+\left|\right|{\eta }_p|{|}_0+||\nabla \times {\mathbf{u}}_{ϵ}-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{u}}_{ϵ}|{|}_0+|{\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ}{|}_{1,h}\hfill \\ \hfill & +\left|\right|p_{ϵ}-{\Pi }_K^{\mathrm{div}}{p}_{ϵ}{\left|\right|}_0+\left|\right|\mathbf{f}-{\mathbf{f}}_h{\left|\right|}_0\left)\right|\left|\right|({\mathbf{v}}_h,q_h){\left|\right||}_h.\hfill \end{array}$

Thirdly, apply the properties of the defined projection operator to obtain

(19)$\begin{array}{cc}\hfill \left|\right||({\xi }_{\mathbf{u}},{\xi }_p){\left|\right||}_h\le & C\left(\right|\left|{\eta }_{\mathbf{u}}\right|{|}_1+\left|\right|{\eta }_p|{|}_0+||\nabla \times {\mathbf{u}}_{ϵ}-{\Pi }_K^{\mathrm{div}}\nabla \times {\mathbf{u}}_{ϵ}|{|}_0\hfill \\ \hfill & +|{\mathbf{u}}_{ϵ}-{\Pi }_K^0{\mathbf{u}}_{ϵ}{|}_{1,h}+\left|\right|p_{ϵ}-{\Pi }_K^{\mathrm{div}}{p}_{ϵ}{\left|\right|}_0+\left|\right|\mathbf{f}-{\mathbf{f}}_h{\left|\right|}_0),\hfill \end{array}$

combining (19) and applications of the triangle inequality, this proof is completed. □

4. Numerical Experiments

In this section, we will present the results of three numerical experiments to illustrate some of the characteristics of the methods discussed in the previous sections. In the first experiment, we used arbitrary hexahedral mesh to verify the effectiveness and convergence of the method. In the second numerical experiment, we used mesh generation with suspension points and polyhedral mesh. In the third experiment, we considered square cavity flow without true solution. In addition, the meshes we use are divided in the cube area of $\Omega =[0,1]\times [0,1]\times [0,1]$.

We use projectors ${\Pi }_K^0$ and ${\Pi }_K^{\nabla }$ to evaluate the error:

• $L^2$-norm: $u_{L^2}{e}_r=\sqrt{{\sum }_{K\in {\mathcal{T}}_h}\left|\right|\mathbf{u}-{\Pi }_K^0{\mathbf{u}}_h{\left|\right|}_{0,K}^2}$;

• $H^1$-norm: $u_{H^1}{e}_r=\sqrt{{\sum }_{K\in {\mathcal{T}}_h}{|\mathbf{u}-{\Pi }_K^{\nabla }{\mathbf{u}}_h|}_{1,K}^2}$.

4.1. Smooth Solution

We consider the parabolic Equation (1), where the load term $\mathbf{f}$ and the Dirichlet boundary are chosen according to the exact solution $\mathbf{u}(x,y,z)=(u_1,u_2,u_3)$ with $ϵ=O\left(h_K\right)$

$\begin{array}{c}\hfill \left\{\begin{array}{c}{u}_1(x,y,z)=sin\left(\pi y(y-1)\right)sin\left(\pi z(z-1)\right),\hfill \\ u_2(x,y,z)=2(1-y^2)y^2(1-2z)(1-z)zsin\left(\pi x(x-1)\right),\hfill \\ u_3(x,y,z)=-2(1-z^2)z^2(1-2y)(1-y)ysin\left(\pi x(x-1)\right),\hfill \\ p(x,y,z)=3-(x^3+y^3+z^3).\hfill \end{array}\right.\end{array}$

The corresponding results are shown in Table 1, where we use arbitrary hexahedron mesh and Kershaw mesh in Figure 1. Table 1 shows the ${\mathbf{u}}_{H_{1}err}$ and ${\mathbf{u}}_{L_{2}err}$ and $p_{L_{2}err}$ when $ϵ=O\left(h_K\right)$. It can be seen from the table that the ${\mathbf{u}}_{H_{1}err}$ and $p_{L_{2}err}$ can reach the order 1, and the ${\mathbf{u}}_{L_{2}err}$ can reach the order 2, which is consistent with our theoretical results.

4.2. True Solution

Let us provide a numerical example. The solution domain $\Omega$ is set as ${[0,1]}^3$. The exact solution $\mathbf{u}$ for the velocity is generated by $\mathbf{curl}\Psi$, where

$\Psi =10x^2{y}^2{z}^2{(x-1)}^2{(y-1)}^2{(z-1)}^2\left(\begin{array}{c}{e}^{x+y+z}\\ 100sin\left(xyz\right)\\ -10(x^2+y^2+z^2)\end{array}\right),$

and the pressure is given by

$p=-sin\left(2\pi x\right)sin\left(2\pi y\right)sin\left(2\pi z\right).$

Through this numerical experiment, we can see from Table 2 that our method is still valid for mesh with hanging points Figure 2a. In addition, for twist polyhedral mesh Figure 2b, our method can reach the theoretical order.

4.3. Driven Cavity Flow

We employ a cubic cavity flow to illustrate the feasibility of the 3D Stokes flows. Consider a cube per unit volume here. In this numerical example, the unit tangential velocity in the x-direction is prescribed on the top surface ($z=1$), and $\mathbf{u}=\mathbf{0}$ is prescribed on the remaining bounding surfaces, see Figure 3. The experiments of driven cavity flow are mainly carried out on locally refinement mesh to verify that our method is effective for mesh with hanging points (see Figure 4).

Figure 5 shows the velocity cross-section of the first locally refined mesh when the degrees of freedom are 1881 and 13,073, and $y=0.5$. Figure 6 shows the velocity cross-section of the second locally refined mesh when $y=0.5$ with 1333 and 9097 degrees of freedom. On the other hand, Figure 7 and Figure 8, respectively, show the velocity vectors of Stokes flow in the X–Y plane with $z=0.5$ in the two grid cubic cavities. The velocity affects the distribution of the intensity of the vorticity, so we look more closely at the vorticity profile for more cross-sections.

It can be seen that when the degrees of freedom are similar, the cavity flow phenomena under different grids are similar. In the same grid, the phenomena of different degrees of freedom are stable, which shows that our method is effective.

5. Conclusions

In this work, an interesting combination of VEM and the lowest order element with practical significance is proposed for the three-dimensional incompressible flow on polyhedral meshes. In order to better illustrate the method, the optimal error estimate is strictly given by taking the Stokes equation as an example. Numerical experiments verify the theoretical analysis and the effectiveness of the method for arbitrary polyhedral meshes, meshes with hanging points and distorted polyhedral meshes. In our ongoing work, we will consider the adaptive scheme in the VEM framework.

Author Contributions

Conceptualization, Y.H. and H.S.; methodology, H.S. and L.L.; software, L.L.; validation, Y.H. and H.S.; formal analysis, H.S. and L.L.; investigation, H.S. and L.L.; resources, Y.H. and H.S.; data curation, L.L.; writing—original draft preparation, L.L.; writing—review and editing, Y.H., H.S. and L.L.; visualization, L.L.; supervision, Y.H.; project administration, Y.H. and H.S.; funding acquisition, Y.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the editor and referees for their valuable comments and suggestions, which helped us to improve the results of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

The following abbreviations are used in this manuscript:

P	pressure
u	dimensionless velocity components
$ϵ$	penalty parameter
$\Pi$${}_K^0$	L2-orthogonal projection
$\Pi$${}_K^{\nabla }$	H1 projection operator
P${}_0$	projection operator of constant functions
K	polyhedral element
f	face of polyhedral element
N${}^K$	dimension of shape function space
h${}_K$	diameter of polyhedral element
e	ridge of polyhedral element
$\mathcal{T}$${}_h$	mesh generation
h${}_f$	diameter of face of polyhedral element

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Polyhedral mesh.

Figure 2. Polyhedral mesh.

Figure 3. Cubic cavity flow problem with boundary conditions.

Figure 4. Polyhedral mesh.

Figure 5. Velocity vectors of x–z plane at y = 0.5 for Stokes flow in a cubic cavity.

Figure 6. Velocity vectors of x–z plane at y = 0.5 for Stokes flow in a cubic cavity.

Figure 7. Velocity vectors of x–y plane at z = 0.5 for Stokes flow in a cubic cavity.

Figure 8. Velocity vectors of x–y plane at z = 0.5 for Stokes flow in a cubic cavity.

Table 1

Error results with arbitrary hexahedron mesh and Kershaw mesh.

$\mathit{DOF}$	${\mathit{u}}_{{\mathit{H}}_1}{\mathit{e}}_{\mathit{r}}$	${\mathit{u}}_{{\mathit{H}}_1}{\mathit{e}}_{\mathit{r}}\mathit{rate}$	${\mathit{u}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}$	${\mathit{u}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}\mathit{rate}$	${\mathit{p}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}$	${\mathit{p}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}\mathit{rate}$
729	1.2090 × 10${}^{-1}$	-	8.0410 × 10${}^{-3}$	-	1.1038 × 10${}^{-1}$	-
4913	6.3291 × 10${}^{-2}$	1.0177	2.2149 × 10${}^{-3}$	2.0273	5.6913 × 10${}^{-2}$	1.0471
35,937	3.1722 × 10${}^{-2}$	1.0414	5.7288 × 10${}^{-4}$	2.0387	2.7656 × 10${}^{-2}$	1.0827
729	6.3351 × 10${}^{-1}$	-	2.3690 × 10${}^{-2}$	-	8.2720 × 10${}^{-2}$	-
4913	3.0385 × 10${}^{-1}$	1.1553	6.9287 × 10${}^{-3}$	1.9330	4.8572 × 10${}^{-2}$	0.8371
35,937	1.4716 × 10${}^{-1}$	1.0930	1.8020 × 10${}^{-3}$	2.0304	2.5461 × 10${}^{-2}$	0.9737

Table 2

Error results with checkerboard mesh and twist polyhedral mesh.

$\mathit{DOF}$	${\mathit{u}}_{{\mathit{H}}_1}{\mathit{e}}_{\mathit{r}}$	${\mathit{u}}_{{\mathit{H}}_1}{\mathit{e}}_{\mathit{r}}\mathit{rate}$	${\mathit{u}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}$	${\mathit{u}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}\mathit{rate}$	${\mathit{p}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}$	${\mathit{p}}_{{\mathit{L}}_2}{\mathit{e}}_{\mathit{r}}\mathit{rate}$
625	2.1395× 10${}^{-1}$	-	1.0443× 10${}^{-1}$	-	3.3572× 10${}^{-1}$	-
4417	1.1618× 10${}^{-1}$	0.9368	3.1281× 10${}^{-2}$	1.8494	1.8848× 10${}^{-2}$	0.8857
33,025	6.1371× 10${}^{-2}$	0.9517	3.0917× 10${}^{-3}$	1.9741	8.9848× 10${}^{-3}$	0.9878
3080	6.6406× 10${}^{-2}$	-	2.1750× 10${}^{-2}$	-	1.6897× 10${}^{-2}$	-
20,160	3.3853× 10${}^{-2}$	1.0759	6.0897× 10${}^{-3}$	2.0327	8.6744× 10${}^{-3}$	1.0647
63,240	2.2962× 10${}^{-2}$	1.0187	2.9105× 10${}^{-3}$	1.9373	5.8729× 10${}^{-4}$	1.0235

References

1. Carey, G.; Krishnan, R. Penalty finite element method for the Navier-Stokes equations. Comput. Methods Appl. Mech. Eng.; 1984; 2, pp. 183-224. [DOI: https://dx.doi.org/10.1016/0045-7825(84)90025-2]

2. Feng, X.; Kim, I.; Nam, H. Locally stabilized P1-nonconforming quadrilateral and hexahedral finite element methods for the Stokes equations. J. Comput. Appl. Math.; 2011; 5, pp. 714-727. [DOI: https://dx.doi.org/10.1016/j.cam.2011.06.009]

3. Brenner, S.C.; Scott, L.R. The Mathematical Theory of Finite Element Methods. Texts Appl. Math.; 2002; 3, pp. 263-291.

4. Beirão da Veiga, L.; Lovadina, C.; Vacca, G. Divergence free Virtual Elements for the Stokes problem on polygonal meshes. Mathematics; 2015; 2, pp. 327-346.

5. Brezzi, F.; Falk, R.S.; Marini, L.D. Basic principles of mixed Virtual Element Methods. Esaim Math. Model. Numer. Anal.; 2014; 4, pp. 1227-1240. [DOI: https://dx.doi.org/10.1051/m2an/2013138]

6. Da Veiga, L.B.; Lovadina, C.; Vacca, G. Virtual Elements for the Navier-Stokes Problem on Polygonal Meshes. SIAM J. Numer. Anal.; 2018; 3, pp. 1210-1242. [DOI: https://dx.doi.org/10.1137/17M1132811]

7. Cáceres, E.; Gatica, G.N. A mixed virtual element method for the pseudostress—Velocity formulation of the Stokes problem. IMA J. Numer. Anal.; 2017; 1, pp. 296-331. [DOI: https://dx.doi.org/10.1093/imanum/drw002]

8. Cáceres, E.; Gatica, G.N.; Sequeira, F.A. A mixed virtual element method for the Brinkman problem. Math. Model. Methods Appl. Sci.; 2017; 4, pp. 707-743.

9. Liu, X.; Li, J.; Chen, Z. A nonconforming virtual element method for the Stokes problem on general meshes. Comput. Methods Appl. Mech. Eng.; 2017; 15, pp. 694-711. [DOI: https://dx.doi.org/10.1016/j.cma.2017.03.027]

10. De Dios, B.; Lipnikov, K.; Manzini, G. The nonconforming virtual element method. ESAIM Math. Model. Numer. Anal.; 2016; 3, pp. 879-904. [DOI: https://dx.doi.org/10.1051/m2an/2015090]

11. Cangiani, A.; Manzini, G.; Sutton, O. Conforming and nonconforming virtual element methods for elliptic problems. IMA J. Numer. Anal.; 2017; 3, pp. 1317-1354. [DOI: https://dx.doi.org/10.1093/imanum/drw036]

12. Da Veiga, L.; Brezzi, F.; Marini, L.; Russo, A. H(div) and H(curl)-conforming virtual element methods. Numer. Math.; 2016; 133, pp. 303-332. [DOI: https://dx.doi.org/10.1007/s00211-015-0746-1]

13. Cascio, M.; Milazzo, A.; Benedetti, I. Virtual element method for computational homogenization of composite and heterogeneous materials. Compos. Struct.; 2020; 232, 111523. [DOI: https://dx.doi.org/10.1016/j.compstruct.2019.111523]

14. Brenner, S.C.; Sung, L.Y. Virtual element methods on meshes with small edges or faces. Math. Model. Methods Appl. Sci.; 2018; 7, pp. 1291-1336. [DOI: https://dx.doi.org/10.1142/S0218202518500355]

15. Zhang, X.; Chi, H.; Paulino, G. Adaptive multi-material topology optimization with hyperelastic materials under large deformations: A virtual element approach. Comput. Methods Appl. Mech. Eng.; 2020; 370, 112976. [DOI: https://dx.doi.org/10.1016/j.cma.2020.112976]

16. Da Veiga, L.; Brezzi, F.; Marini, L. Virtual elements for linear elasticity problems. SIAM J. Numer. Anal.; 2013; 2, pp. 794-812. [DOI: https://dx.doi.org/10.1137/120874746]

17. Da Veiga, L.; Lovadina, C.; Mora, D. A virtual element method for elastic and inelastic problems on polytope meshes. Comput. Methods Appl. Mech. Eng.; 2015; 295, pp. 327-346. [DOI: https://dx.doi.org/10.1016/j.cma.2015.07.013]

18. Gain, A.; Talischi, C.; Paulino, G. On the virtual element method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes. Comput. Methods Appl. Mech. Eng.; 2014; 282, pp. 132-160. [DOI: https://dx.doi.org/10.1016/j.cma.2014.05.005]

19. Antonietti, P.; Manzini, G.; Verani, M. The fully nonconforming virtual element method for biharmonic problems. Math. Model. Methods Appl. Sci.; 2018; 2, pp. 387-407. [DOI: https://dx.doi.org/10.1142/S0218202518500100]

20. Brezzi, F.; Marini, L. Virtual element methods for plate bending problems. Comput. Methods Appl. Mech. Eng.; 2013; 253, pp. 455-462. [DOI: https://dx.doi.org/10.1016/j.cma.2012.09.012]

21. Zhao, J.; Chen, S.; Zhang, B. The nonconforming virtual element method for plate bending problems. Math. Model. Methods Appl. Sci.; 2016; 9, pp. 1671-1687. [DOI: https://dx.doi.org/10.1142/S021820251650041X]

22. Zhao, J.; Zhang, B.; Chen, S. The Morley-type virtual element for plate bending problems. J. Sci. Comput.; 2018; 1, pp. 610-629. [DOI: https://dx.doi.org/10.1007/s10915-017-0632-3]

23. Antonietti, P.; Da Veiga, L.; Scacchi, S.; Verani, M. A C¹ virtual element method for the Cahn–Hilliard equation with polygonal meshes. SIAM J. Numer. Anal.; 2016; 1, pp. 34-56. [DOI: https://dx.doi.org/10.1137/15M1008117]

24. Antonietti, P.; Da Veiga, L.; Mora, D.; Verani, M. A stream virtual element formulation of the Stokes problem on polygonal meshes. SIAM J. Numer. Anal.; 2014; 1, pp. 386-404. [DOI: https://dx.doi.org/10.1137/13091141X]

25. Da Veiga, L.; Lovadina, C.; Vacca, G. Divergence free virtual elements for the Stokes problem on polygonal meshes. ESAIM Math. Model. Numer. Anal.; 2017; 2, pp. 509-535. [DOI: https://dx.doi.org/10.1051/m2an/2016032]

26. Cangiani, A.; Gyrya, V.; Manzini, G. The nonconforming virtual element method for the Stokes equations. SIAM J. Numer. Anal.; 2016; 6, pp. 3411-3435. [DOI: https://dx.doi.org/10.1137/15M1049531]

27. da Veiga, L.; Manzini, G. A virtual element method with arbitrary regularity. IMA J. Numer. Anal.; 2014; 2, pp. 759-781. [DOI: https://dx.doi.org/10.1093/imanum/drt018]

28. Vacca, G. Virtual element methods for hyperbolic problems on polygonal meshes. Comput. Math. Appl.; 2017; 5, pp. 882-898. [DOI: https://dx.doi.org/10.1016/j.camwa.2016.04.029]

29. Vacca, G.; da Veiga, L. Virtual element methods for parabolic problems on polygonal meshes. Numer. Methods Partial. Differ. Equations; 2015; 6, pp. 2110-2134. [DOI: https://dx.doi.org/10.1002/num.21982]

30. Perugia, I.; Pietra, P.; Russo, A. A plane wave virtual element method for the Helmholtz problem. ESAIM Math. Model. Numer. Anal.; 2016; 3, pp. 783-808. [DOI: https://dx.doi.org/10.1051/m2an/2015066]

31. Mascotto, L. Ill-conditioning in the virtual element method: Stabilizations and bases. Numer. Methods Partial. Differ. Equ.; 2018; 4, pp. 1258-1281. [DOI: https://dx.doi.org/10.1002/num.22257]

32. Benedetto, M.; Berrone, S.; Borio, A.; Pieraccini, S.; Scialo, S. Order preserving SUPG stabilization for the virtual element formulation of advection–diffusion problems. Comput. Methods Appl. Mech. Eng.; 2016; 311, pp. 18-40. [DOI: https://dx.doi.org/10.1016/j.cma.2016.07.043]

33. Da Veiga, L.; Chernov, A.; Mascotto, L.; Russo, A. Basic principles of h_p virtual elements on quasiuniform meshes. Math. Model. Methods Appl. Sci.; 2016; 8, pp. 1567-1598. [DOI: https://dx.doi.org/10.1142/S021820251650038X]

34. He, Y.; Li, J. A penalty finite element method based on the Euler implicit/explicit scheme for the time-dependent Navier-Stokes equations. J. Comput. Appl. Math.; 2010; 3, pp. 708-725. [DOI: https://dx.doi.org/10.1016/j.cam.2010.06.025]

35. He, Y.; Li, J. A stabilized finite element method based on local polynomial pressure projection for the stationary Navier-Stokes equations. Appl. Numer. Math.; 2008; 10, pp. 1503-1514. [DOI: https://dx.doi.org/10.1016/j.apnum.2007.08.005]

36. Chen, L.; Huang, J. Some error analysis on virtual element methods. Calcolo; 2018; 1, pp. 1-23. [DOI: https://dx.doi.org/10.1007/s10092-018-0249-4]

37. Veiga, L.; Brezzi, F.; Marini, L. Virtual Element Methods for general second order elliptic problems on polygonal meshes. Math. Model. Methods Appl. Sci.; 2014; 4, pp. 75-89.

38. Kashiwabara, T.; Oikawa, I.; Zhou, G. Penalty method with P₁/P₁ finite element approximation for the Stokes equations under the slip boundary condition. Numer. Math.; 2016; 4, pp. 705-740. [DOI: https://dx.doi.org/10.1007/s00211-016-0790-5]