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

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 (Xh,Qh) 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.


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 α -conforming (α > 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 in f -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 in f -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 in f -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.

The Stationary Penalty Stokes Equations
Let Ω be a polygonal domain in R 3 . This paper considers the 3D incompressible Stokes equations to illustrate the penalty virtual element method: in Ω, where u and p are the velocity field and the pressure field, respectively. f represents the external force [35]. The penalty method applied to (1) is to approximate the solution (u, p) by (u , p ) satisfying the following steady-state penalty Stokes Equation [34] where the penalty parameter is 0 < < 1. The Galerkin variational formula of Equation (1): find (u, p) ∈ V × Q, namely where the bilinear forms a(×; ×) : It is well known that (see for [2,35]): • a(×, ×) and b(×, ×) are continuous, i.e., • a(×, ×) is coercive, i.e., there exists a positive constant α such that • Moreover, the bilinear form b(×, ×) satisfies the in f -sup condition: where there exists a constant β 0 > 0 such that The variational formula of Equation (2): Theorem 1. There exists at least a solution pair (u , p ) ∈ V × Q which satisfies (4) and Moveover, if f ∈ L 2 (Ω), and satisfies C 0 ≤ 1/2, then the solution (u , p ) of (4) satisfies the following regularity: where C 0 denotes some generic constant depending on the data (Ω, f ), which may stand for different values at its different occurrences. There are the following error estimates [1,34].

The Penalty Virtual Element Method for Stokes Equations
The model problem (1) where Ω is (for simplicity) a convex polyhedron. We suppose that we are given a decomposition T h of Ω in rather general polyhedra K. We can assume that T h satisfies the regularity assumption [5]: ∃ a constant C > 0 assume that • each polyhedron K ∈ T h is star-shaped with respect to every point of a ball of radius ≥ Ch K ; • for every face f ∈ ∂K we have h f ≥ Ch K and f is star-shaped with respect to every point of a disk of radius ≥ Ch f ; • for every edge e ∈ ∂ f , we have |e| ≥ Ch 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 ≥ 1, the same discussion as in [10] reveals that the dimension of V K h is 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 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 With the definition of the space of a single face W k ( f ), we give the definition of the virtual meta space of all faces here, ∂K is the set of all faces of element K. The definition of local finite element spaces is given by For each polyhedron K, we define a local finite element space V k (K). Roughly speaking, V k (K) 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 P k (K) is a subspace of V k (K). We can take the following degrees of freedom in V k (K): where the scaled monomials m α ∈ M k−2 ( f ) are defined by • the moments up to order k − 2 of v h in K: where the scaled monomials µ α ∈ M k−2 (K) are defined in (1.5), and ν k−2 by For each polyhedron K ∈ T h and v h ∈ W h (K), we can construct Π ∇ K v h : W k (K) → P k (K) according to the method in [6], which is defined as follows where, P 0 is a projection operator of constant functions, that we choose as Remark 1. Note that, integrating by parts , we have again, for every p k ∈ P k (K), Since p k ∈ P k−2 (K), 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 ∂K, ∂p k ∂n is in P k−1 ( f ), but the choice of using v h | f ∈ W k ( f ) 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 ≥ 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 (K), 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.
The L 2 (K)-orthogonal projection Π 0 K :W k (K) → P 0 (K), and for any w h ∈ V k (K), For each polyhedron K, a local finite element space V k (K) is defined. Now, the local virtual element spaces for the velocity and pressure are, respectively, introduced by Then, we can define the global finite element space V h as

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 K h (×, ×) is set as is any symmetric, positive definite and computable bilinear form that guarantees [5]: Then, it is easy to prove that this discrete local bilinear form satisfies the following consistency and stability assumptions [36]: • Stability: there are two positive constants α * and α * dependent on h K and K, have • Computability: we can know the computability of the discrete bilinear form from Remark 1.
Finally, we define the global approximated bilinear form a h (×, ×) by simply summing the local contributions The linear form f(v) on the right-hand side of the variational problem is discrete by f h such that The discrete problem (2) which we can solve is written as: find u h ∈ V h such that where The first term in d h (p h , q h ) is the penalty term, and the second term is the stable term.

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].
We have the following proposition Under the mesh assumptions, there exists a polynomial function u I ∈ V h , such that where C is a constant independent of h, and D(K) denotes the "diamond" of K, i.e., the union of the polygon in T h intersecting K.
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 π h : where C π is a positive h-independent constant. Given v ∈ V, using arguments borrowed from [38] and considering the VEM interpolant v I presented in Proposition (4.1), we first 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 ∈ T h and all smooth enough Combined with Equation (16), we can obtain the estimation formula of Similarly, we can do the same with d h (ξ p , q h ) Add Equations (15)-(18), combine Lemma 2 and Proposition 1 in [4], we have Thirdly, apply the properties of the defined projection operator to obtain combining (19) and applications of the triangle inequality, this proof is completed.

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 We use projectors Π 0 K and Π ∇ K to evaluate the error:

Smooth Solution
We consider the parabolic Equation (1), where the load term f and the Dirichlet boundary are chosen according to the exact solution u(x, y, z) = (u 1 , u 2 , The corresponding results are shown in Table 1, where we use arbitrary hexahedron mesh and Kershaw mesh in Figure 1. Table 1 shows the u H 1 err and u L 2 err and p L 2 err when = O(h K ). It can be seen from the table that the u H 1 err and p L 2 err can reach the order 1, and the u L 2 err can reach the order 2, which is consistent with our theoretical results.

True Solution
Let us provide a numerical example. The solution domain Ω is set as [0, 1] 3 . The exact solution u for the velocity is generated by curlΨ, where and the pressure is given by p = −sin(2πx)sin(2πy)sin(2πz).
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.

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 u = 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, Figures 7 and 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.

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.  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.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.