Next Article in Journal
Link Prediction in Complex Networks Using Recursive Feature Elimination and Stacking Ensemble Learning
Next Article in Special Issue
Linear Full Decoupling, Velocity Correction Method for Unsteady Thermally Coupled Incompressible Magneto-Hydrodynamic Equations
Previous Article in Journal
Asynchronous Stabilization for Two Classes of Stochastic Switching Systems with Applications on Servo Motors
Previous Article in Special Issue
Analysis of the Element-Free Galerkin Method with Penalty for Stokes Problems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

College of Mathematics and System Sciences, Xinjiang University, Urumqi 830046, China
*
Author to whom correspondence should be addressed.
Submission received: 10 July 2022 / Revised: 11 August 2022 / Accepted: 11 August 2022 / Published: 15 August 2022

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 i n f - s u p 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 α -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 i n f - s u p 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 i n f - s u p 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 i n f - s u p 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 Ω be a polygonal domain in R 3 . This paper considers the 3D incompressible Stokes equations to illustrate the penalty virtual element method:
Δ u + p = f , in Ω , × u = 0 , in Ω , u = 0 , on Ω ,
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]
Δ u ϵ + p ϵ = f , in Ω , × u ϵ + ϵ p ϵ = 0 , in Ω , u ϵ = 0 , on Ω ,
where the penalty parameter is 0 < ϵ < 1 .
Consider the Hilbert spaces
V : = [ H 0 1 ( Ω ) ] 3 , Q : = L 0 2 ( Ω ) = { q L 2 ( Ω ) : Ω q d Ω = 0 } ,
with norms | | v | | 1 : = | | v | | [ H 0 1 ( Ω ) ] 3 , | | q | | Q : = | | q | | L 2 ( Ω ) .
The Galerkin variational formula of Equation (1): find ( u , p ) V × Q , namely
a ( u , v ) b ( v , p ) = ( f , v ) , v V , b ( u , q ) = 0 , q ) Q ,
where the bilinear forms a ( × ; × ) : V × V R and b ( × ; × ) : V × Q R be defined by:
a ( u , v ) : = Ω u v d Ω , u , v V ,
b ( v , q ) : = Ω d i v v q d Ω , v V , q Q .
It is well known that (see for [2,35]):
  • a ( × , × ) and b ( × , × ) are continuous, i.e.,
    | a ( u , v ) | | | u | | 1 | | v | | 1 u , v V ,
    | b ( v , q ) | | | v | | 1 | | q | | Q v V , a n d q Q ,
  • a ( × , × ) is coercive, i.e., there exists a positive constant α such that
    a ( v , v ) α | | v | | 1 2 , v V ;
  • Moreover, the bilinear form b ( × , × ) satisfies the i n f - s u p condition: where there exists a constant β 0 > 0 such that
    sup v V b ( v , q ) | | v | | 1 β 0 | | q | | 0 , q Q .
The variational formula of Equation (2): find ( u ϵ , p ϵ ) V × Q such that
a ( u ϵ , v ) b ( v , p ϵ ) = ( f , v ) , v V b ( u ϵ , q ) + ϵ ( p ϵ , q ) = 0 , q Q .
Theorem 1.
There exists at least a solution pair ( u ϵ , p ϵ ) V × Q which satisfies (4) and
| | u ϵ | | 0 | | f | | 1 .
Moveover, if f L 2 ( Ω ) , and ϵ satisfies ϵ C 0 1 / 2 , then the solution ( u ϵ , p ϵ ) of (4) satisfies the following regularity:
| | u ϵ | | 2 + 1 / 2 | | p ϵ | | 0 C 0 f 0
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].
Lemma 1.
Assume that f L 2 ( Ω ) , ( u , p ) and ( u ϵ , p ϵ ) are the solution of variational formulas (3) and (4), respectively. Then
| | ( u u ϵ ) | | 0 + | | u u ϵ | | 0 + | | p p ϵ | | 0 c ϵ f 1 .
Please refer to [1,34] for details of certification.

3. 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 C h K ;
  • for every face f K we have h f C h K and f is star-shaped with respect to every point of a disk of radius C h f ;
  • for every edge e f , we have | e | 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 1 , the same discussion as in [10] reveals that the dimension of V h K is
N K V = 2 n k + k ( k 1 ) , f o r d = 2 , 3 2 n k ( k + 1 ) + 1 2 ( k 1 ) k ( k + 1 ) , f o r d = 3 .
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 α = ( x x K ) α h K , | α | k }
where x = ( x , y , z ) is a vector, α = ( α 1 , α 2 , α 3 ) is a triple index, h K is the element K diameter, and 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
W k ( f K ) : = { v H 1 ( f K ) : v | f K B k ( f K ) , Δ v P k 2 ( f K ) } ,
where B k ( f K ) : = { v C 0 ( f K ) : v | e P k ( f K ) , e f K } . Here W k ( f ) is the function space of a face f K defined in an element K, B k ( f K ) space is defined on the face f K , e is an edge of the face f K , and f K is all the edges of the face f K . P k ( f K ) is the set of polynomials on f K of degree k .
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
W k K ( K ) : = { v C 0 ( K ) : v | f K W k ( f K ) , Δ v P k 2 ( K ) } ,
here, K is the set of all faces of element K.
The definition of local finite element spaces is given by
W k ( K ) : = { v H 1 ( K ) : v | K V k K ( K ) , Δ v | K P k 2 ( K ) } .
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 ) :
  • the value of v h at the vertices of K;
  • on each edge e, the value of 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 v h in f:
    f v h m α , α = 1 , , n k 2 ,
    where the scaled monomials m α M k 2 ( f ) are defined by
    m α ( x ) : = ( x x f h f ) α ;
  • the moments up to order k 2 of v h in K:
    K v h μ α , α = 1 , , ν k 2 ,
    where the scaled monomials μ α M k 2 ( K ) are defined in (1.5), and ν k 2 by
    ν α ( x ) : = ( x x K h K ) α .
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
( w h , ( Π K v h v h ) ) 0 , K = 0 , w h W h , P 0 ( Π K v h v h ) = 0 ,
where, P 0 is a projection operator of constant functions, that we choose as
P 0 v h : = 1 N K V i = 1 N K V v h ( V i ) , f o r k = 1 ,
P 0 v h : = 1 | K | K v h , f o r k 2 .
Remark 1.
Note that, integrating by parts, we have again, for every p k P k ( K ) ,
( p k , v h ) 0 , K = K p k v h + K p k n v h .
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.
W ˜ k ( K ) : = { v h W k ( K ) : ( q , v h Π K v h ) 0 , K = 0 , q M k 1 M k } .
The L 2 ( K ) -orthogonal projection Π K 0 : W ˜ k ( K ) P 0 ( K ) , and for any w h V k ( K ) ,
( w h Π K 0 w h , m k ) 0 , K = 0 , m k M 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
V k ( K ) : = [ W ˜ k ( K ) ] 3 , Q k ( K ) : = W ˜ k ( K ) , K T h .
Then, we can define the global finite element space V h as
V h : = { v h [ H 0 1 ( Ω ) ] 3 : v h | K V k ( K ) for all K T h } , Q h : = { q h L 0 2 ( Ω ) : q | K Q k K T h } .

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 ( × , × ) is set as
a h K ( u h , v h ) = ( ( Π K u h ) , ( Π K v h ) ) K + S K ( u h , v h ) K
S K ( u h , v h ) K = h K ( ( u h Π K u h ) , ( v h Π K v h ) ) K
where S K is any symmetric, positive definite and computable bilinear form that guarantees [5]:
c a K ( v h , v h ) S K ( v h , v h ) c a K ( v h , v h ) , v h k e r ( Π K ) .
Then, it is easy to prove that this discrete local bilinear form satisfies the following consistency and stability assumptions [36]:
  • k-compatibility: if v V k ( K ) , s P k ( K ) , have
    a h K ( v , s ) = a K ( v , s ) ;
  • Stability: there are two positive constants α and α dependent on h K and K, have
    α a K ( v , v ) a h K ( v , v ) α a K ( v , v ) , v V k ( K ) ;
  • 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
a h ( u h , v h ) : = K T h a h K ( u h , v h ) , u h , v h V h .
We define projection Π K div ( × v h ) : ( Π K div ( × v h ) , m α ) K = ( × v h , m α ) K , m α M ( K ) . The bilinear form b h K ( p h , v h ) and global contributions b h ( p h , v h ) are given as
b h K ( p h , v h ) = ( Π K 0 p h , Π K div ( × v h ) ) , p h Q h , v h V h , b h ( p h , v h ) : = K T h b h K ( p h , v h ) .
The linear form f ( v ) on the right-hand side of the variational problem is discrete by f h such that
( f h , v h ) = K T h ( f , Π K 0 v h ) ) 0 , K .
The discrete problem (2) which we can solve is written as: find u h V h such that
a h ( u h , v h ) b h ( p h , v h ) + b h ( u h , q h ) + d h ( p h , q h ) = ( f h , v h ) , v h V h , q h Q h ,
where
d h ( p h , q h ) = K T h ϵ ( Π K 0 p h , Π K 0 q h ) + K T h ϵ ( ( I Π K 0 ) p h , ( I Π 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 τ h , then for all u [ H s + 1 ( K ) ] 3 with 0 s k , there exists a polynomial function u π [ P k ( K ) ] 3 , such that
| | u u π | | 0 , K + h K | u u π | 1 , K C h K s + 1 | u | s + 1 , K .
We have the following proposition
Proposition 1.
Let u V [ H s + 1 ( K ) ] 3 with 0 s k . Under the mesh assumptions, there exists a polynomial function u I V h , such that
| | u u I | | 0 , K + h K | u u I | 1 , K C h K s + 1 | u | s + 1 , D ( K ) ,
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.
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 ( × , × ) satisfies the discrete i n f - s u p condition, i.e.,
sup v h V h v h 0 b ( u h , q h ) | | v h | | 1 β | | q h | | Q q h Q h .
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 : V V h , satisfying
b ( π h v , q h ) = b ( v h , q h ) , v V , q h Q h , | | π h v | | C π | | v | | 1 , v V ,
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 construct v ¯ h V h such that
b ( v v ¯ h , q ¯ h ) = 0 , q h Q h ,
and
| | v v ¯ h | | 1 C | | v | | 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 T h and all smooth enough functions φ defined on K, it holds:
| | φ Π K φ | | C h K s m | φ | s , K , s , m N , m s k + 1 , s 1 ,
| | φ Π K 0 φ | | C h K s m | φ | s , K , s , m N , m s k + 1 ,
| | φ φ I | | m , K C h K s m | φ | s , E , m , s N , m s k + 1 , s 2 .
Define a norm as | | | ( u h , p h ) | | | h 2 : = | | u h | | 1 2 + | | p h | | 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 { h K } , K T h . Let ( u ϵ , p ϵ ) and ( u h , p h ) is solution of (4) and (10), respectively. Then, it holds
| | u ϵ u h | | 1 + | | p ϵ p h | | 0 C ( | | u ϵ u I | | 1 + | | p ϵ p I | | 0 + | | × u ϵ Π K div × u ϵ | | 0 + | u ϵ Π K 0 u ϵ | 1 , h + | | p ϵ Π K div p ϵ | | 0 + | | f f h | | 0 ) .
Proof. 
Combined with the definition of projection operator Π K , Π K 0 , Π K d i v . When k = 1 , Π K = Π K 0 , Π K d i v is the component of Π K (references [4,5]). Firstly, we simply deal with u ϵ u h and p ϵ p h ,
u ϵ u h = u ϵ u I + u I u h = η u + ξ u ,
p ϵ p h = p ϵ p I + p I p h = η p + ξ p .
Next, using the above notation and k-consistency, we obtain
a h ( ξ u , v h ) = K T h a h K ( u I u h , v h ) = K T h [ a h K ( u I Π K 0 u ϵ , v h ) + a h K ( Π K 0 u ϵ , v h ) ] a h ( u h , v h ) = K T h [ a h K ( u I Π K 0 u ϵ , v h ) + a K ( u ϵ Π K 0 u ϵ , v h ) ] + a ( u ϵ , v h ) a h ( u h , v h ) = K T h a h K ( η u , v h ) + K T h a h K ( u ϵ Π K 0 u ϵ , v h ) K T h a K ( u ϵ Π K 0 u ϵ , v h ) + a ( u ϵ , v h ) a h ( u h , v h ) .
According to Lemmas 3 and 4, we have
b h ( v h , ξ p ) = K T h ( Π K div × v h , Π K 0 ( p I p h ) ) 0 , K = K T h ( Π K div × v h , p I p h ) 0 , K = K T h ( Π K div × v h , η p ) 0 , K + K T h ( × v h Π K div × v h , p ϵ ) 0 , K ( × v h , p ϵ ) + K T h ( Π K div × v h , p h ) 0 , K = K T h ( Π K div × v h , η p ) 0 , K + K T h ( × v h Π K div × v h , p ϵ Π K div p ϵ ) 0 , K ( × v h , p ϵ ) + K T h ( Π K div × v h , Π K 0 p h ) 0 , K = K T h ( Π K div × v h , p ϵ ) 0 , K + K T h ( × v h , p ϵ Π K div p ϵ ) 0 , K b ( v h , p ϵ ) + b h ( v h , p h ) .
Combined with Equation (16), we can obtain the estimation formula of b h ( ξ u , q h )
b h ( ξ u , q h ) = K T h ( × η u , Π K div q h ) 0 , K + b ( u ϵ , q h ) K T h ( × u ϵ Π K div × u ϵ , q h ) 0 , K b h ( u h , q h ) .
Similarly, we can do the same with d h ( ξ p , q h )
d h ( ξ p , q h ) = d h ( p I , q h ) d h ( p h , q h ) .
Add Equations (15)–(18), combine Lemma 2 and Proposition 1 in [4], we have
a h ( ξ u , v h ) b h ( v h , ξ p ) + b h ( ξ u , q h ) + d h ( ξ p , q h ) = K T h a h K ( η u , v h ) + K T h a h K ( u ϵ Π K 0 u ϵ , v h ) K T h a K ( u ϵ Π K 0 u ϵ , v h ) + K T h ( Π K div × v h , p ϵ ) 0 , K + ( f f h , v h ) + K T h ( × v h , p ϵ Π K div p ϵ ) 0 , K K T h ( × η u , Π K div q h ) 0 , K K T h ( × u ϵ Π K div × u ϵ , q h ) 0 , K + d h ( p I , q h ) C ( | | η u | | 1 + | | η p | | 0 + | | × u ϵ Π K div × u ϵ | | 0 + | u ϵ Π K 0 u ϵ | 1 , h + | | p ϵ Π K div p ϵ | | 0 + | | f f h | | 0 ) | | | ( v h , q h ) | | | h .
Thirdly, apply the properties of the defined projection operator to obtain
| | | ( ξ u , ξ p ) | | | h C ( | | η u | | 1 + | | η p | | 0 + | | × u ϵ Π K div × u ϵ | | 0 + | u ϵ Π K 0 u ϵ | 1 , h + | | p ϵ Π K div p ϵ | | 0 + | | f f h | | 0 ) ,
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 Ω = [ 0 , 1 ] × [ 0 , 1 ] × [ 0 , 1 ] .
We use projectors Π K 0 and Π K to evaluate the error:
  • L 2 -norm: u L 2 e r = K T h | | u Π K 0 u h | | 0 , K 2 ;
  • H 1 -norm: u H 1 e r = K T h | u Π K u h | 1 , K 2 .

4.1. 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 , u 3 ) with ϵ = O ( h K )
u 1 ( x , y , z ) = sin ( π y ( y 1 ) ) sin ( π z ( z 1 ) ) , u 2 ( x , y , z ) = 2 ( 1 y 2 ) y 2 ( 1 2 z ) ( 1 z ) z sin ( π x ( x 1 ) ) , u 3 ( x , y , z ) = 2 ( 1 z 2 ) z 2 ( 1 2 y ) ( 1 y ) y sin ( π x ( x 1 ) ) , p ( x , y , z ) = 3 ( x 3 + y 3 + z 3 ) .
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 e r r and u L 2 e r r and p L 2 e r r when ϵ = O ( h K ) . It can be seen from the table that the u H 1 e r r and p L 2 e r r can reach the order 1, and the u L 2 e r r 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 Ω is set as [ 0 , 1 ] 3 . The exact solution u for the velocity is generated by curl Ψ , where
Ψ = 10 x 2 y 2 z 2 ( x 1 ) 2 ( y 1 ) 2 ( z 1 ) 2 e x + y + z 100 s i n ( x y z ) 10 ( x 2 + y 2 + z 2 ) ,
and the pressure is given by
p = s i n ( 2 π x ) s i n ( 2 π y ) s i n ( 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.

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 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, Figure 7 and Figure 8, respectively, show the velocity vectors of Stokes flow in the XY 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.

Funding

This research received no external funding.

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:
Ppressure
udimensionless velocity components
ϵ penalty parameter
Π K 0 L2-orthogonal projection
Π K H1 projection operator
P 0 projection operator of constant functions
Kpolyhedral element
fface of polyhedral element
N K dimension of shape function space
h K diameter of polyhedral element
eridge of polyhedral element
T h mesh generation
h f diameter of face of polyhedral element

References

  1. Carey, G.; Krishnan, R. Penalty finite element method for the Navier-Stokes equations. Comput. Methods Appl. Mech. Eng. 1984, 2, 183–224. [Google Scholar] [CrossRef]
  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, 714–727. [Google Scholar] [CrossRef]
  3. Brenner, S.C.; Scott, L.R. The Mathematical Theory of Finite Element Methods. Texts Appl. Math. 2002, 3, 263–291. [Google Scholar]
  4. Beirão da Veiga, L.; Lovadina, C.; Vacca, G. Divergence free Virtual Elements for the Stokes problem on polygonal meshes. Mathematics 2015, 2, 327–346. [Google Scholar]
  5. Brezzi, F.; Falk, R.S.; Marini, L.D. Basic principles of mixed Virtual Element Methods. Esaim Math. Model. Numer. Anal. 2014, 4, 1227–1240. [Google Scholar] [CrossRef]
  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, 1210–1242. [Google Scholar] [CrossRef]
  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, 296–331. [Google Scholar] [CrossRef]
  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, 707–743. [Google Scholar]
  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, 694–711. [Google Scholar] [CrossRef]
  10. De Dios, B.; Lipnikov, K.; Manzini, G. The nonconforming virtual element method. ESAIM Math. Model. Numer. Anal. 2016, 3, 879–904. [Google Scholar] [CrossRef]
  11. Cangiani, A.; Manzini, G.; Sutton, O. Conforming and nonconforming virtual element methods for elliptic problems. IMA J. Numer. Anal. 2017, 3, 1317–1354. [Google Scholar] [CrossRef]
  12. Da Veiga, L.; Brezzi, F.; Marini, L.; Russo, A. H(div) and H(curl)-conforming virtual element methods. Numer. Math. 2016, 133, 303–332. [Google Scholar] [CrossRef]
  13. Cascio, M.; Milazzo, A.; Benedetti, I. Virtual element method for computational homogenization of composite and heterogeneous materials. Compos. Struct. 2020, 232, 111523. [Google Scholar] [CrossRef]
  14. Brenner, S.C.; Sung, L.Y. Virtual element methods on meshes with small edges or faces. Math. Model. Methods Appl. Sci. 2018, 7, 1291–1336. [Google Scholar] [CrossRef]
  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. [Google Scholar] [CrossRef]
  16. Da Veiga, L.; Brezzi, F.; Marini, L. Virtual elements for linear elasticity problems. SIAM J. Numer. Anal. 2013, 2, 794–812. [Google Scholar] [CrossRef]
  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, 327–346. [Google Scholar] [CrossRef]
  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, 132–160. [Google Scholar] [CrossRef]
  19. Antonietti, P.; Manzini, G.; Verani, M. The fully nonconforming virtual element method for biharmonic problems. Math. Model. Methods Appl. Sci. 2018, 2, 387–407. [Google Scholar] [CrossRef]
  20. Brezzi, F.; Marini, L. Virtual element methods for plate bending problems. Comput. Methods Appl. Mech. Eng. 2013, 253, 455–462. [Google Scholar] [CrossRef]
  21. Zhao, J.; Chen, S.; Zhang, B. The nonconforming virtual element method for plate bending problems. Math. Model. Methods Appl. Sci. 2016, 9, 1671–1687. [Google Scholar] [CrossRef]
  22. Zhao, J.; Zhang, B.; Chen, S. The Morley-type virtual element for plate bending problems. J. Sci. Comput. 2018, 1, 610–629. [Google Scholar] [CrossRef]
  23. Antonietti, P.; Da Veiga, L.; Scacchi, S.; Verani, M. A C1 virtual element method for the Cahn–Hilliard equation with polygonal meshes. SIAM J. Numer. Anal. 2016, 1, 34–56. [Google Scholar] [CrossRef]
  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, 386–404. [Google Scholar] [CrossRef]
  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, 509–535. [Google Scholar] [CrossRef]
  26. Cangiani, A.; Gyrya, V.; Manzini, G. The nonconforming virtual element method for the Stokes equations. SIAM J. Numer. Anal. 2016, 6, 3411–3435. [Google Scholar] [CrossRef]
  27. da Veiga, L.; Manzini, G. A virtual element method with arbitrary regularity. IMA J. Numer. Anal. 2014, 2, 759–781. [Google Scholar] [CrossRef]
  28. Vacca, G. Virtual element methods for hyperbolic problems on polygonal meshes. Comput. Math. Appl. 2017, 5, 882–898. [Google Scholar] [CrossRef]
  29. Vacca, G.; da Veiga, L. Virtual element methods for parabolic problems on polygonal meshes. Numer. Methods Partial. Differ. Equations 2015, 6, 2110–2134. [Google Scholar] [CrossRef]
  30. Perugia, I.; Pietra, P.; Russo, A. A plane wave virtual element method for the Helmholtz problem. ESAIM Math. Model. Numer. Anal. 2016, 3, 783–808. [Google Scholar] [CrossRef]
  31. Mascotto, L. Ill-conditioning in the virtual element method: Stabilizations and bases. Numer. Methods Partial. Differ. Equ. 2018, 4, 1258–1281. [Google Scholar] [CrossRef]
  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, 18–40. [Google Scholar] [CrossRef]
  33. Da Veiga, L.; Chernov, A.; Mascotto, L.; Russo, A. Basic principles of hp virtual elements on quasiuniform meshes. Math. Model. Methods Appl. Sci. 2016, 8, 1567–1598. [Google Scholar] [CrossRef]
  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, 708–725. [Google Scholar] [CrossRef]
  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, 1503–1514. [Google Scholar] [CrossRef]
  36. Chen, L.; Huang, J. Some error analysis on virtual element methods. Calcolo 2018, 1, 1–23. [Google Scholar] [CrossRef]
  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, 75–89. [Google Scholar]
  38. Kashiwabara, T.; Oikawa, I.; Zhou, G. Penalty method with P1/P1 finite element approximation for the Stokes equations under the slip boundary condition. Numer. Math. 2016, 4, 705–740. [Google Scholar] [CrossRef]
Figure 1. Polyhedral mesh.
Figure 1. Polyhedral mesh.
Entropy 24 01129 g001
Figure 2. Polyhedral mesh.
Figure 2. Polyhedral mesh.
Entropy 24 01129 g002
Figure 3. Cubic cavity flow problem with boundary conditions.
Figure 3. Cubic cavity flow problem with boundary conditions.
Entropy 24 01129 g003
Figure 4. Polyhedral mesh.
Figure 4. Polyhedral mesh.
Entropy 24 01129 g004
Figure 5. Velocity vectors of x–z plane at y = 0.5 for Stokes flow in a cubic cavity.
Figure 5. Velocity vectors of x–z plane at y = 0.5 for Stokes flow in a cubic cavity.
Entropy 24 01129 g005
Figure 6. 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.
Entropy 24 01129 g006
Figure 7. Velocity vectors of x–y plane at z = 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.
Entropy 24 01129 g007
Figure 8. 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.
Entropy 24 01129 g008
Table 1. Error results with arbitrary hexahedron mesh and Kershaw mesh.
Table 1. Error results with arbitrary hexahedron mesh and Kershaw mesh.
DOF u H 1 e r u H 1 e r rate u L 2 e r u L 2 e r rate p L 2 e r p L 2 e r rate
7291.2090 × 10 1  - 8.0410 × 10 3  - 1.1038 × 10 1  - 
49136.3291 × 10 2 1.01772.2149 × 10 3 2.02735.6913 × 10 2 1.0471
35,9373.1722 × 10 2 1.04145.7288 × 10 4 2.03872.7656 × 10 2 1.0827
7296.3351 × 10 1  - 2.3690 × 10 2  - 8.2720 × 10 2  - 
49133.0385 × 10 1 1.15536.9287 × 10 3 1.93304.8572 × 10 2 0.8371
35,9371.4716 × 10 1 1.09301.8020 × 10 3 2.03042.5461 × 10 2 0.9737
Table 2. Error results with checkerboard mesh and twist polyhedral mesh.
Table 2. Error results with checkerboard mesh and twist polyhedral mesh.
DOF u H 1 e r u H 1 e r rate u L 2 e r u L 2 e r rate p L 2 e r p L 2 e r rate
6252.1395× 10 1  - 1.0443× 10 1  - 3.3572× 10 1  - 
44171.1618× 10 1 0.93683.1281× 10 2 1.84941.8848× 10 2 0.8857
33,0256.1371× 10 2 0.95173.0917× 10 3 1.97418.9848× 10 3 0.9878
30806.6406× 10 2  - 2.1750× 10 2  - 1.6897× 10 2  - 
20,1603.3853× 10 2 1.07596.0897× 10 3 2.03278.6744× 10 3 1.0647
63,2402.2962× 10 2 1.01872.9105× 10 3 1.93735.8729× 10 4 1.0235
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, L.; Su, H.; He, Y. Penalty Virtual Element Method for the 3D Incompressible Flow on Polyhedron Mesh. Entropy 2022, 24, 1129. https://doi.org/10.3390/e24081129

AMA Style

Li L, Su H, He Y. Penalty Virtual Element Method for the 3D Incompressible Flow on Polyhedron Mesh. Entropy. 2022; 24(8):1129. https://doi.org/10.3390/e24081129

Chicago/Turabian Style

Li, Lulu, Haiyan Su, and Yinnian He. 2022. "Penalty Virtual Element Method for the 3D Incompressible Flow on Polyhedron Mesh" Entropy 24, no. 8: 1129. https://doi.org/10.3390/e24081129

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop