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

T.Vu-Huu,C.Le-Thanh,H.Nguyen-Xuan and M.Abdel-Wahab

1Faculty of Civil Engineering,Vietnam Maritime University,484 Lach Tray Str.,Hai Phong,Vietnam

2Faculty of Civil Engineering and Electricity,Ho Chi Minh City Open University,Vietnam

3CIRTech Institute,Ho Chi Minh City University of Technology(HUTECH),Ho Chi Minh City,Vietnam

4Department of Architectural Engineering,Sejong University,209 Neungdongro,Gwangjingu,Seoul,05006,Republic of Korea

5Soete Laboratory,Faculty of Engineering and Architecture,Ghent University,Technologiepark Zwijnaarde,903,B-9052,Zwijnaarde,Belgium

Abstract: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 various Reynolds numbers up to the maximum Re=1000.All results then numerously compared to available sources in the literature.

Keywords: Lid-driven cavity; incompressible; steady; Navier-Stokes equations; polygonal finite element method

1 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-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-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 Pe1Pe1.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-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., Pe1Pe1, [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-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-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.

2 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 andpare the fluid velocity in H1(Ω)and the modified pressure (after dividing by water densityρ) of the fluid inL2(Ω), respectively.Then, the given domainΩhas the following boundary conditions(Γ=ΓD∩ΓNin whichΓDis the Dirichlet boundary andΓNis the Neumann boundary):

where n and∂u/∂nare 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:

whereH1is the standard Sobolev space.The standard weak formulations of Eqs.(1) and (2)become [25,30,31]:

where v andqare test functions;ν∇vdΩis diffusive term;vdΩis the convective term what is the nonlinear term of convection, which is identified by a trilinear form→R as:

where z is a given function in V0that is the subspace of divergence-free velocities as:

Then, determine the bilinear formaz(.,.): V0×V0→R through:

The fundamental feature is the skew-symmetric of convection term:ς(z;u,v)=-ς(z;v,u)over V0.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: =u1-u2andp: =p1-p2where(u1,p1)and(u2,p2)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.

3 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(u0,p0)∈H1×L2(Ω).A sequence of iterates(u0,p0),(u1,p1),(u2,p2),…∈H1×L2(Ω), 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(uk,pk)into Eqs.(7) and (8).It thereby creates a pair of nonlinear residualsRk(v)andrk(q)[25,33,34]:

The solution of Eqs.(7) and (8) at thekthiteration that is u=uk+δukandp=pk+δpkleads to the correction termsδuk∈andδpk∈satisfying:

whereD(u,δu,v), is the difference in the nonlinear terms,D(u,δu,v)=ς(δu;δu,v)+ς(δu;u,v)+ς(u;δu,v).ExpandingD(uk,δuk,v)and dropping the quadratic termς(δuk;δuk,v)in Eqs.(19) and(20) produces a linear problem.Then, the corrected termsδuk∈andδpk∈satisfy [25,33,34]:

The solution of Eqs.(21) and (22) is the so-called Newton correction.The previous iterate is updated by uk+1=uk+δuk,pk+1=pk+δpkto find the sequent iterate.The consistency of the iteration process is ensured by substituting uk=u andpk=pinto the update formula provided that the right-hand side of Eqs.(21) and (22) is zero only whenδuk=0 andδpk=0.

The second approach of linearisation is the Picard method.It bases on the dropping of the quadratic termς(δuk;δuk,v)and the linear termς(δuk;uk,v)in Eqs.(19) and (20).Thus, instead of Eqs.(21) and (22), we have an alternative linear system to findδuk∈andδpk∈as [25,33,34]:

4 Polygonal Mixed Finite Discretisation

A discretisation based on two finite spacesMh⊂and Xh⊂to approximate fluid pressure and velocity, respectively, independently leads to the mixed finite element method terminology.Specifically, for the given spaces of velocity, Xh, and pressure,Mh, the discretisation now becomes findingph∈Mhand uh∈Xhso that:

Based on the Newton iteration, find the correction termsδuh∈Xhandδph∈Mhto give the finite-dimensional analogue of Eqs.(21) and (22) as:

Here,Rk(vh)andrk(qh)are the nonlinear residuals associated with the discrete formulation of Eqs.(25) and (26), respectively.If the termς(δuh;uh,vh)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:

in whichhi(x)denote the perpendicular distance of the interior point v to the edge eiof the polygon; niis the outward unit normal vector to the edge ei=[xi,xi+1] with counterclockwise ordering vertices indexed cyclically xn+1=x1of 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 uhandph, 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 byc(·,·):Mh×Mh→R.The local stabilisation matrix,ce, of elements becomes [9,10]:

for allΩe∈ℑh.In which the basis functionsψiandψjdoes not vanish on the elementΩe.Thus,the global stabilisation matrix C can be assembled from contribution matricesceby:

where A is the assembly operator,eis theethelement,neis the total number of the elements in the domain ℑh.

5 Numerical Tests

As mentioned, this research applies the equal-order PFE, e.g., Pe1Pe1, 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))areux=-1(m/s)anduy=0.Besides, the no-flow (zero velocity) conditions are applied on the left, right and bottom boundaries, see Fig.1.The fluid density isρ= 1For validation, various literature,i.e., [14-17,35-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 toRe=1000, six progressively finer hexagonal meshes with a total of 4046, 16491, 36977, 66102, 14841 and 262075 elements are created, respectively.

Figure 1: Left-lid-driven Cavity Domain; Right-A sample of hexagonal mesh

The initial results of this paper are presented in Fig.2, indicating the velocity contours at the differentReof the cavity problem.As the initial observation, at the higherRethe positive contour centre ofuxtend 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 highRe, i.e., 1000, see Fig.2.As fluid viscosity reduces, there is less resistance to fluid movement, which may explain those tendencies.Alternatively, asReincreases, the convection term becomes more prominent, resulting in more contour vortexes (see Eq.(42)).

Figure 2: Contours of velocity fields (Left-Re=100; Middle - Re=400; Right-Re=1000).(a)Horizontal velocity field, ux.(b) Vertical velocity field, uy

For more profound validation, the velocities components (i.e.,uxanduy) on the cavity centrelines (i.e.,xc-vertical centreline andyc-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].

Figure 3: Velocity on the cavity centerlines.(a) ux-vertical velocity.(b) uy-horizontal velocity

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 ofandrespectively.

Figure 4: Locations of extrema velocities on the cavity centerlines (left-(m/s); middle-and right (a) Re=100 (b) Re=400 (c) Re=1000

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 asRerises, 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.

Figure 5: Fluid pressure contours (a) Re=100.(b) Re=400.(c) Re=1000.

Figure 6: Fluid pressure on cavity centrelines.(a) Fluid pressure on the vertical centerline.(b)Fluid pressure on the horizontal centerline

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 atRe=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 otherRe, 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 atRe= 100 and 400, as in Tab.3.

Figure 7: Streamline results of the cavity flow at Re=100.(a) Streamlines of the cavity flow.(b)Centre locations of the primary vortex (PV).(c) Centre locations of the left secondary vortex(LSV)

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 atRe= 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-38].In addition, the convergence of the streamlines is also verified (see Tab.4).

Besides, the streamlined results of the lid-driven cavity simulation at Re = 400 are the subject of the following evaluation.They are illustrated in 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 atRe=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-8and 5.02×10-,8respectively.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-38].

Figure 8: Streamline results of the cavity flow at Re=400.(a) Streamlines of the cavity flow.(b)Centre locations of the primary vortex (PV).(c) Centre locations of the left secondary vortex(LSV).(d) Centre locations of the right secondary vortex (RSV)

Figure 9: Streamline results of the cavity flow at Re=1000.(a) Streamlines of the cavity flow.(b)Centre locations of PV.(c) Centre locations of LSV.(d) Centre locations of the right secondary vortex (RSV).(e) Centre locations of the left tertiary vortex (LTV).(f) Centre locations of the right tertiary vortex (RTV)

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 ofN=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.

6 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 Pe1Pe1[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.

FundingStatement: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.

Appendix A.

Table 1: The velocity on centrelines of the cavity

Table 2: Extrema velocities on cavity centrelines

Table 3: The pressure on centrelines of the cavity

Table 4: Vortex intensity at Re=100;(x,y) refers to the centroid vortex location

Table 5: Vortex intensity at Re=400; (x,y) refers to the centroid vortex location

Table 6: Vortex intensity at Re=1000; (x,y) refers to the centroid vortex location

(Continued)

Table 6: Continued