Polygonal Finite Element for Two-Dimensional Lid-Driven Cavity Flow

This paper investigates a polygonal finite element (PFE) to solve a two-dimensional (2D) incompressible steady fluid problem in a cavity square. It is a well-known standard benchmark (i.e., lid-driven cavity flow)-to evaluate the numerical methods in solving fluid problems controlled by the Navier– Stokes (N–S) equation system. The approximation solutions provided in this research are based on our developed equal-order mixed PFE, called Pe1Pe1. It is an exciting development based on constructing the mixed scheme method of two equal-order discretisation spaces for both fluid pressure and velocity fields of flows and our proposed stabilisation technique. In this research, to handle the nonlinear problem of N-S, the Picard iteration scheme is applied. Our proposed method’s performance and convergence are validated by several simulations coded by commercial software, i.e., MATLAB. For this research, the benchmark is executed with variousReynolds numbers up to themaximum Re= 1000. All results then numerously compared to available sources in the literature.


Introduction
This research, instead of using widely used numerical approaches such as the finite difference method (FDMfinite volume method (FVM)), finite element method (FEM), e.g., [1][2][3], etc., proposes an advanced PFE method (PFEM) to solve the 2D lid-driven cavity problem controlled by the incompressible steady N-S equations. As known, PFEM is a robust numerical method offering a wide range of distinguished advantages, especially flexibility and the benefit of Voronoi algorithms to generate meshes with arbitrary element shapes [4,5]. Furthermore, PFEM offers better accuracy without the need for a sizeable overall mesh scale compared to its triangular and quadrilateral counterparts [6,7]. It means that PFEs can provide better solutions than triangular and quadrilateral elements [7,8]. However, the fact is that developments of PFEM in the fluid field is still too modest compared to the enormous potential of the method. The most current research of PFEM for fluid analysis, hitherto, was only given by Talischi et al. [6] in 2014, see and then is our research in 2019, see [9][10][11][12]. However, such research only focuses on the performance, stability and convergence of the PFEs for incompressible steady Stokes problems. Therefore, this research aims to use our recently developed PFE to solve 2D incompressible steady N-S cavity problems.
As that goal, this study adopts the equal-order mixed PFE, named Pe 1 Pe 1 . The primary advantage of our developed element is the computational ability for fluid flows on all kinds of mesh families, e.g., triangular, quadrilateral, hexagonal, random and centroidal Voronoi meshes. It is constructed by the mixed scheme between two equal-order discretisation spaces for both fluid pressure and velocity fields of flows. In this research, Wachspress basis shape functions are utilised to represent the approximation spaces of velocity and pressure fields. Furthermore, this research executes our novel stabilisation technique to eliminate the instability of the equal-order mixed scheme [9][10][11][12]. It is an innovation of the local polynomial pressure projection method introduced by Bochev et al. [13] in 2004. It automatically adapts the local stabilisation term for each arbitrary shape of element on a polygonal mesh. Our advanced stabilisation method adds a term to penalise pressure deviations from the 'consistent' polynomial order. And it helps to avoid the residual terms of the penalty method to maintain the symmetry of the numerical system.
As mentioned, this research employs a well-known benchmark (i.e., lid-driven cavity flow)-a classical standard to evaluate the numerical methods in solving the N-S equations. This benchmark's main advantage is the simplicity of the geometry, leading to applying numerical methods on this flow in terms of coding is relatively easy and straightforward. Despite its simple geometry, the driven cavity flow occupies a rich fluid physics flow [14,15]. The cavity problem is early and widely utilised by many researchers, e.g., Ghia et al. [16]; Botella et al. [14]; Bruneau et al. [17], etc. This study, hence, applies the 2D lid-driven cavity benchmark with various Reynolds numbers (i.e., Re = 100, 400 and 1000) to assess the performance of our developed PFE (i.e., Pe 1 Pe 1 , [10]) in solving the incompressible steady N-S equations. The solutions of this study are compared with the highly accurate benchmark solutions found in the literature.
N-S equations, as known, are a system of the nonlinear term of convection problems [18]. Nonlinear equations cannot, in general, be solved analytically [19,20]. Thus, in this case, nonlinear equations must be handled by iteration progress. The concept behind such iteration techniques starts at an arbitrary point-the closest possible point to the solution sought-and gradually reach the solution through a series of sequent tests [19,21,22]. Because of the advantages of the Picard iteration method (i.e., a huge ball of convergence [23][24][25], efficiency and simplicity [26], etc.), Picard is chosen to handle the nonlinear term in N-S equations in this study. In addition, to take advantage of the PFEM, the mesh generation algorithms and Voronoi diagrams' properties [4][5][6] are utilised. Besides, the Wachspress coordinates [27,28] are also adopted. In addition, the advanced techniques developed in [5,7,29], which handle the bottleneck of generating quality polygonal meshes, are applied.
The paper is organised as follows: Section 2 presents the incompressible steady N-S flow problems. In Section 3, we illustrate the iteration progress to solve the N-S equation system's nonlinear problem. Then, Section 4 presents the mixed discretisation scheme. Section 5 shows the numerical tests' results. Finally, the conclusions and future works are given in Section 6.

Incompressible Navier-Stokes Equations
As known, the steady-state N-S equation system for incompressible fluid flow is generally written as following [25,30,31]: where ν is a given positive constant called the fluid kinematic viscosity; f is the forcing term; u and p are the fluid velocity in H 1 ( ) and the modified pressure (after dividing by water density ρ) of the fluid in L 2 ( ), respectively. Then, the given domain has the following boundary conditions ( = D ∩ N in which D is the Dirichlet boundary and N is the Neumann boundary): where n and ∂u/∂n are the outward-pointing normal vector and directional derivative in the normal vector direction, respectively. For more simplification, body force vector f is set to zero with a little loss of generality. In this case, the conservation term of body force (gravity force) is the gradient of a scalar field, that is, f = −∇ , and thus it can be incorporated into the system by redefining the pressure (p → p + ). Then, by choosing the test spaces: where H 1 is the standard Sobolev space. The standard weak formulations of Eqs. (1) and (2) become [25,30,31]: where v and q are test functions; ν ∇u : ∇vd is diffusive term; (u · ∇u) ·vd is the convective term what is the nonlinear term of convection, which is identified by a trilinear form ς : where z is a given function in V 0 that is the subspace of divergence-free velocities as: Then, determine the bilinear form a z (., .) : V 0 × V 0 → R through: The fundamental feature is the skew-symmetric of convection term: ς (z; u, v) = −ς(z; v, u) over V 0 . It means that ς(z; u, u) = 0 and so we have [25,32,33]: So, the continuity becomes: The experience shows that the well-posedness of the weak formulations of N-S systems is a complicated problem because of the nonlinearity. To explain it, we define u : = u 1 − u 2 and p : = p 1 − p 2 where (u 1 , p 1 ) and (u 2 , p 2 ) represent distinct solutions, that is: It means that it is impossible to apply the homogeneous problem to establish uniqueness for Eqs. (7) and (8). To ensure the uniqueness of the weak solution, it can base on well-posedness conditions between the forcing function f and boundary data w as [25,34]: Alternatively, a well-known condition for uniqueness is [33]: where * is the best potential constant from Eq. (13). It firms that a unique weak solution exists when the viscosity parameter is large enough.

Nonlinear Iterations
To solve N-S equations, we need a linearised process to deal with the nonlinear problem at every computation step. So, we first need an "initial guess" that commonly is the solution of the corresponding Stokes problem (u 0 , p 0 ) ∈ H 1 × L 2 ( ). A sequence of iterates (u 0 , p 0 ), (u 1 , p 1 ), (u 2 , p 2 ),…∈ H 1 × L 2 ( ), then, is computed to find the convergent solution. Newton method is the first natural approach for the iteration process. This method applies a given iterate (u k , p k ) into Eqs. (7) and (8). It thereby creates a pair of nonlinear residuals R k (v) and r k (q) [25,33,34]: The solution of Eqs. (7) and (8) at the k th iteration that is u = u k + δu k and p = p k + δp k leads to the correction terms δu k ∈ H 1 0 and δp k ∈ L 2 0 ( ) satisfying: where D(u, δu, v), is the difference in the nonlinear terms, and dropping the quadratic term ς(δu k ; δu k , v) in Eqs. (19) and (20) produces a linear problem. Then, the corrected terms δu k ∈ H 1 0 and δp k ∈ L 2 0 ( ) satisfy [25,33,34]: The solution of Eqs. (21) and (22) is the so-called Newton correction. The previous iterate is updated by u k+1 = u k + δu k , p k+1 = p k + δp k to find the sequent iterate. The consistency of the iteration process is ensured by substituting u k = u and p k = p into the update formula provided that the right-hand side of Eqs. (21) and (22) is zero only when δu k = 0 and δp k = 0.

Polygonal Mixed Finite Discretisation
A discretisation based on two finite spaces M h ⊂ L 2 0 ( ) and X h ⊂ H 1 0 to approximate fluid pressure and velocity, respectively, independently leads to the mixed finite element method terminology. Specifically, for the given spaces of velocity, X h , and pressure, M h , the discretisation now becomes finding p h ∈ M h and u h ∈ X h so that: Based on the Newton iteration, find the correction terms δu h ∈ X h and δp h ∈ M h to give the finite-dimensional analogue of Eqs. (21) and (22) as: Here, R k (v h ) and r k (q h ) are the nonlinear residuals associated with the discrete formulation of Eqs. (25) and (26), respectively. If the term ς(δu h ; u h , v h ) is dropped, it is the discrete analogue of the Picard update.
To state the approximated solution, we introduce a set of vector-valued basis functions {φ j }−for velocity and a set of pressure basis functions {ψ k } for pressure as: In which, as mentioned, the polygonal basis shape functions for both velocity and pressure field in this research are constructed by Wachspress coordinates as: and their gradients are: where, n ne is the number vertices of a polygon e in the mesh I h ( = n e e=1 e ≈ I h ); p i (x) is defined as: in which h i (x) denote the perpendicular distance of the interior point v to the edge e i of the polygon; n i is the outward unit normal vector to the edge e i = [x i , x i+1 ] with counterclockwise ordering vertices indexed cyclically x n+1 = x 1 of the polygon [9,10,28,32].
Substituting Eqs. (29) and (30) into Eqs. (27) and (28) gives a system of linear equations as: Here, A is the vector-Laplacian matrix, B is the divergence matrix. Besides, N is the vectorconvection matrix, and W is the Newton derivative matrix, defined as follow: The right-hand side vectors in Eq. (34) are the nonlinear residuals associated with the current discrete solution estimates u h and p h , expanded via Eqs. (29) and (30): For Picard iteration, by omitting the Newton derivative matrix, the mixed finite discretisation system of Eqs. (27) and (28) becomes: Because of the efficiency and simplicity, the Picard iterative scheme is applied in this research to handle the nonlinear term in N-S equations. In case of unstable problems, we need to eliminate the zero block in the system Eq. (34) as well as the system Eq. (41) by a stabilisation matrix to get stable results. Therefore, the stabilised analogue of the system Eq. (41) is [11,25]: Here, C is the stabilisation matrix provided by c (·, ·) : M h × M h → R. The local stabilisation matrix, c e , of elements becomes [9,10]: for all e ∈ I h . In which the basis functions ψ i and ψ j does not vanish on the element e . Thus, the global stabilisation matrix C can be assembled from contribution matrices c e by: where A is the assembly operator, e is the e th element, n e is the total number of the elements in the domain I h .

Numerical Tests
As mentioned, this research applies the equal-order PFE, e.g., Pe 1 Pe 1 , to analyse the extensive well-known lid-driven cavity flow (see Fig. 1) of incompressible steady N-S flows. The detail of this benchmark is that its domain is a unit square = (0, 1) 2 (m). The boundary condition of velocity along the estuary of the cavity (0 ≤ x ≤ 1 (m) , y = 1(m)) are u x = −1(m/s) and u y = 0. Besides, the no-flow (zero velocity) conditions are applied on the left, right and bottom boundaries, see Fig. 1. The fluid density is ρ = 1 kg/m 3 . For validation, various literature, i.e., [14][15][16][17][35][36][37][38][39], of the lid-driven cavity flow is collected to make detailed comparisons for our results. Then, to test our advanced method with different Reynolds numbers up to Re = 1000, six progressively finer hexagonal meshes with a total of 4046, 16491, 36977, 66102, 14841 and 262075 elements are created, respectively. The initial results of this paper are presented in Fig. 2, indicating the velocity contours at the different Re of the cavity problem. As the initial observation, at the higher Re the positive contour centre of u x tend to reach to the bottom of the cavity. For the vertical velocity, the tendency is that the centres try to get lower and closer to the cavity walls. The next focus is the generation of the secondary contours at the cavity corners of simulations at the high Re, i.e., 1000, see Fig. 2. As fluid viscosity reduces, there is less resistance to fluid movement, which may explain those tendencies. Alternatively, as Re increases, the convection term becomes more prominent, resulting in more contour vortexes (see Eq. (42)).
For more profound validation, the velocities components (i.e., u x and u y ) on the cavity centrelines (i.e., x c -vertical centreline and y c -horizontal one, respectively) are compared to the relevant literature, see Fig. 3. As seen in the figure, the entire agreement of our results to those of other studies, such as Ghia et al. [16], Botella et al. [14], Bruneau et al. [17], Erturk et al. [15], is apparent. For further details, see Tab. 1 of velocity values at reference points on the cavity centrelines, affirming our proposed method's excellent performance. Like the table, the result of our proposed technique is very comparable to those seen in the literature, especially Ghia et al. [16]. For instance, the velocities at the domain's centre point (0.5, 0.5) vary only by 1%-4% from the literature. Thus, an excellent or even marginally improved solution is provided. Following that, the extrema velocities on the centrelines are compared (see Fig. 4). Our recent results clearly mirrored those of Ghia et al. [16], Bottela et al. [14], Vanka [36] as well as Bruneau et al. [38].
Tab. 2 contains more information about the apparent proof for the convergence and precision of the new approach in this study. Moreover, our models computed on increasingly finer meshes steadily exceed the reference and convergence solution, as seen in Fig. 4 and Tab. 2. Furthermore, readers will see that our new results are greater than Ghia et al. [16] and Vanka [36], i.e., the results of u min y and u max x , respectively. The subsequent validation provided in this research is the results of fluid pressure. These results are firstly presented in Figs. 5 and 6 In which, Fig. 5 depicts the fluid pressures on the entire cavity using contour plots. These results indicate that as Re rises, the pressure primary contour's centre moves closer to the centre of cavity. As seen in Fig. 5, the rise also causes more secondary pressure contours to surround primary ones. Besides, the singularity of pressure at the cavity's top corners (0, 1) and (1, 1) is noteworthy. It occurs as a function of the abrupt change in velocity at those corners. References [40,41] explains this rule in terms of saddle-point theory. The indications of pressure on the cavity centrelines are then seen in Fig. 6 to allow a detailed comparison with the literature. As the figure, it is evident that the current pressure results, i.e., the pressure at Re = 1000 is completely close to the reference solution, i.e., Bottela et al. [14], Bruneau et al. [17]. Tab. 3 reveals seventeen points of pressure on the cavity centrelines for a more detailed comparison. The pressure results on several different progressively finer meshes are given in this table to show that our proposed solution is convergent. The results in the table convince the comparatively strong convergence of our recent method. It can be seen that our pressure results give an excellent fitting to the reference results in [14,17] (the disparity is about 3.7%). As known, the error of the other methods, for example, FDM to [16], typically is 4% [14]. Hence, it can be confirmed that our recent polygonal method provided an excellent solution for both velocity and pressure. It is crucial evidence of the new technique's performance in solving 2D incompressible steady N-S problems. The only concern is that there is not enough comparison data for pressures at other Re, such as 100 and 400, resulting in poor validity. However, it makes a significant contribution in terms of providing further reference data for future research. In other words, this study adds to the body of knowledge about fluid pressure by using previously unpublished evidence, such as the pressure of a lid-driven cavity flow at Re = 100 and 400, as in Tab. 3.
All the above results already exhibited a perfect convergence as well as the precision of velocity and pressure. Thus, the remaining part is to illustrate another result of streamlines. The first streamlines regarding the lid-driven cavity flow at Re = 100 are seen in Fig. 7a. As can be seen from this figure, the current streamlines are identical to Ghia et al. [16]. Furthermore, Figs. 7b and 7c provide a more in-depth analysis of the streamline effects by displaying all of the vortex centre properties (i.e., location, intensity). Besides, Tab. 4 displays the values of such data. It is concrete proof for this validation. In terms of the figure and table, it is easy to firm that the streamlines result completely match [16,[35][36][37][38]. In addition, the convergence of the streamlines is also verified (see Tab. 4).  Fig. 8 and Tab. 5. They validate our recent results' great matches to [16,35,36,42], as well as the enhanced polygonal method's good convergence and accuracy. The formation of another secondary vortex in the cavity's right bottom corner, as seen in Fig. 8, is then of particular interest. As known, the flow at the bottom corners of the cavity is almost fragile (close to zero). The generation of secondary vortexes is supported by reduced friction force as the Reynolds number (viscosity) increases. It also explains the appearance of the tertiary vortices of the following simulation at Re = 1000 (see Fig. 9). They are incredibly close to the bottom corners of the cavity where the intensity of vortices is entirely trivial (extremely close to zero, see Tab. 6). The intensities of left and right tertiary vortex, in particular, approximate 8.9 × 10 −8 and 5.02 × 10 −8 , respectively. The accuracy of our current results is also identified in Fig. 9, which shows the closeness of the current vortex centres to the reference ones. Tab. 6 then indicates how the intensities and positions of vortex centres are close to those in [14,16,17,[35][36][37][38].   In conclusion, all the above results reveal that the present polygonal method is entirely good enough to deal with the 2D incompressible steady flows of lid-driven cavity benchmark. It is reinforced by numerous extensive comparisons to previously published research that used various highly precise techniques. For example, Botella et al. [14] in 1998 used the highly accurate Chebyshev collocation method associated with the subtraction method of the leading terms of the asymptotic expansion to obtain a highly accurate spectral solution with a maximum of grid mesh of N = 160 (polynomial degree). Erturk et al. [15] 2005 applied FDM and spectral method to provide a high spatial accuracy solution of driven cavity flow. Bruneau et al. [17] in 2006 utilised the finite difference discretisation and the multigrid solver with a cell-by-cell relaxation procedure.

Conclusion
To summarise, this research's commitment to designing new PFEM to solve 2D incompressible steady fluid flows controlled by N-S equations is practical. This contribution is focused on the development of a mixed-order equal-order system. This idea then applies Wachspress basis shape functions associated with our advanced stabilisation technique to generate the equal-order mixed PFE, called Pe 1 Pe 1 [10]. Furthermore, in this article, we successfully address the Picard iteration technique for our proposed PFE to deal with the nonlinear convection term of N-S equations. This paper incorporates a widely mathematical benchmark to evaluate the accuracy and performance of our developed PFE. It is the lid-driven cavity benchmark executed to measure the efficiency of numerical methods in solving N-S problems. Our established method shows an excellent agreement with the highly accurate solutions found in the literature regarding current research solutions. It means that the efficacy of our proposed technique for solving incompressible steady N-S fluid flows has been proven without any reasonable doubt.
Furthermore, the production of higher-order PFEs is a high-potential direction for future research in this field. For example, Rand et al. [43] suggest a novel quadratic serendipity shape function in 2014, a solid foundation for this direction. Alternatively, the Taylor-Hood elements, as described in [44], could be a good start. Additionally, extending our proposed PFEs to computations for transient fluid flow problems is a significant strategic work of this study. The proposed method's application to 3D problems is, therefore, a priority task. Other promising ideas for future research are free surface and fluid-structure interaction problems. In addition, in recent years, Deep Neural Networks (DNNs) developments are becoming a mathematical option to solve the partial differential equations (PDEs) of different phenomena in science and engineering [45,46]. Thus, the efficiency of DNNs is an exciting direction for this research.
Funding Statement: This work was supported by the VLIR-UOS TEAM Project, VN2017TEA454A 103, 'An innovative solution to protect Vietnamese coastal riverbanks from floods and erosion' funded by the Flemish Government.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.