AN EMA-CONSERVING, PRESSURE-ROBUST AND RE-SEMI-ROBUST METHOD WITH A ROBUST RECONSTRUCTION METHOD FOR NAVIER–STOKES ⋆

. Proper EMA-balance (balance of kinetic energy, linear momentum and angular momentum), pressure-robustness and Re-semi-robustness ( Re : Reynolds number) are three important properties of Navier–Stokes simulations with exactly divergence-free elements. This EMA-balance makes a method conserve kinetic energy, linear momentum and angular momentum in an appropriate sense; pressure-robustness means that the velocity errors are independent of the pressure; Re-semi-robustness means that the constants appearing in the error bounds of kinetic and dissipation energies do not explicitly depend on inverse powers of the viscosity. In this paper, based on the pressure-robust reconstruction framework and certain suggested reconstruction operators in Linke and Merdon [ Comput. Methods Appl. Mech. Eng. 311 (2016) 304–326], we propose a reconstruction method for a class of non-divergence-free simplicial elements which admits almost all the above properties. The only exception is the energy balance, where kinetic energy should be replaced by a suitably redefined discrete energy. The lowest order case is the Bernardi–Raugel element on general shape-regular meshes. Some numerical comparisons with exactly divergence-free methods, the original pressure-robust reconstruction methods and the EMAC method are provided to confirm our theoretical results.


Introduction
In this paper, we study the finite element methods for the unsteady Navier-Stokes equations (NSEs): ) ) where  = (0,  ] is a bounded time interval, Ω ⊂ R  ( = 2, 3) is a bounded domain with Lipschitz-continuous polyhedral boundary Γ,  and  represent the velocity and pressure unknowns,  > 0 is the constant kinematic viscosity,  represents the external force, and  0 is the initial velocity.We assume  () ∈  2 (Ω) for all  ∈ .
For designing accurate numerical schemes, it is widely believed that preserving the fundamental (physical or mathematical) properties of the continuous problem is of great importance.For NSEs, these fundamental properties include the divergence constraint (1.1b), the balance laws for some physical quantities (e.g., kinetic energy, linear momentum, angular momentum, vorticity, helicity, and enstrophy) [12,19], and an invariance property for the velocity with respect to the gradient field in  [34,41].Among these properties, the divergence constraint is of central importance.The papers [19,34] showed that exactly divergence-free mixed methods preserve the balance laws and the invariance property mentioned above.The latter means that these methods are pressure-robust; namely the velocity errors are independent of the pressure.Moreover, it was demonstrated in [55] that the constants in error estimates, including the Gronwall constant, did not depend on  −1 explicitly for divergence-free finite element methods.This property was called (Re-)semi-robustness, uniform, or quasiuniform estimates [35] sometimes.
Due to these fascinating properties, constructing exactly divergence-free elements has been an increasingly hot topic in recent years [14, 28-30, 46, 47, 58].However, compared to classical non-divergence-free elements, the construction of these elements is not trivial in most cases.For example, the well-known -th order Scott-Vogelius elements usually require a barycentric-refined mesh.
Another popular idea is relaxing the continuity condition but enforcing the divergence constraint strongly, which results in the so-called nonconforming (div)-conforming methods [15,27,36,53,56,57].In this paper, we focus on the conforming mixed methods for the Navier-Stokes equations.
The method introduced in this paper is also included in this topic.Pressure-robustness plays an important role on the accuracy of a method when "gradient forces dominate the momentum balance" [44]; the velocity errors of the methods which are not Re-semi-robust might grow quickly with respect to time for higher Reynolds number flows [55]; the EMAC scheme is one of the "enhanced-physics" based schemes which have a long history such as [1,4,20,51,52], and the paper [49] showed that an improper treatment of energy, linear momentum and angular momentum produced lower bounds for  2 velocity error.It is worth mentioning that, the properties mentioned above are usually not mutually independent.For example, in the paper [49], Olshanskii and Rebholz proved that the Gronwall constants in EMAC estimates did not depend on the viscosity explicitly, which was exactly Re-semi-robustness except that the constant related to a pressure-induced error term polynomially depended on the inverse of viscosity.Another example is the popular grad-div stabilization [11,48,50].It can not only weaken the impact of the pressure on velocity errors [34,48] but also make the usual skew-symmetric scheme Re-semi-robust [16].Finally, we also refer the readers to the review articles [22,34,35] for more details.
In this paper, we propose a reconstruction formulation which is EMA-conserving for the reconstructed discrete velocity (here the energy should be redefined as (2.15) below), pressure-robust, and Re-semi-robust.For simplicity, we shall refer to this reconstruction method as the "EMAPR" method throughout this paper.Our method is based on the pressure-robust reconstruction formulation in [43].The main difference lies on the discretization of the convective term.For the convective term, two pressure-robust discrete forms were proposed in [43]: the convective form and the rotational form.Following a similar strategy as in [12], one can check that these two forms do not conserve linear momentum and angular momentum (the latter conserves kinetic energy).
Here we propose an EMA-conserving form, i.e., it does not produce any extra energy, momentum and angular momentum in an appropriate sense.Then we give a pressure-robust and Re-semi-robust error estimate for the continuous-in-time case, provided that the continuous solution  is in  2 (︀ ;  1,∞ (Ω) )︀ . To obtain such an estimate, we also need to slightly modify the discretization of the evolutionary term by introducing a stabilization.
Finally, we shall prove that our formulation can be easily applied to a class of simplicial locally mass-conserving elements whose pressures are discontinuous.Moreover, the reconstruction also lowers the requirement of the quadrature rules.Note that in a practical view a drawback of this class of elements is that it uses bubble functions to enrich velocity spaces which require higher order quadrature rules.To the best of our knowledge, the EMAPR method is the first method on conforming non-divergence-free elements which is EMA-conserving, pressure-robust, and Re-semi-robust simultaneously.For nonconforming and non-divergence-free methods, a reconstructed hybrid discontinuous Galerkin method in [39] (see formulas "(6.3d)" and "(6.5)" therein) probably admits most of these properties also.The lowest number of DoFs (degrees of freedom) of our method is the same as of the lowest order divergence-free methods proposed in [14,30], while our method uses simpler bubble functions which do not require any special split of meshes.
We also conjecture that the EMAC formulation with grad-div stabilization is a simple and efficient strategy in practice, as the grad-div stabilization can dramatically lower the effect of the pressure and also make a method Re-semi-robust.Moreover, in the original paper [48], it was shown that the grad-div stabilization could also improve the performance of linear algebraic solvers.However, for a much more complicated or very large pressure (which can happen in practice; cf.[21,42]), the grad-div stabilization might not be enough since it can not remove the the effect of the pressure totally.Some numerical comparisons on this aspect are also done in this paper.
The remainder of this paper is organized as follows.In Section 2 we discuss the EMAPR methods and some balance laws.Section 3 is devoted to a pressure-robust and Re-semi-robust error estimate for the EMAPR method.We show that a class of non-divergence-free elements can be easily incorporated into our framework in Section 4. Finally we carry out some numerical experiments in Section 5 and do some conclusions in Section 6.
In what follows we will use , with or without a subscript, to denote a generic positive constant.The standard inner product for With the convention the subscripts ,  and  will be omitted for  = 0,  = 2 and  = Ω, respectively.  () coincides with  ,2 () and  , () coincides with [ , ()]  .

The divergence-free reconstruction operator
Let  ℎ denote a partition of Ω.We define the mesh size ℎ := max ∈ ℎ ℎ  with ℎ  the diameter of elements .Denote by   the diameter of the biggest ball inscribed in .Here we assume that  ℎ is shape-regular, i.e., there exists a positive constant  such that We define and where  is the unit external normal vector on Γ. Furthermore we define the bilinear form  :  ×  → R by (, ) := (∇ • , ) for any (, ) ∈  × .
Let ( ℎ ,  ℎ ,  ℎ ) ⊂ ( , ,  ) denote a triple of finite element spaces satisfying inf Denote by the spaces of divergence-free velocity functions and discretely divergence-free velocity functions, respectively.Note that if ∇ •  ℎ ⊆  ℎ , we have  0 ℎ ⊂  0 , which means that the functions in  0 ℎ are exactly divergence-free.For most classical elements, this relationship does not hold.
Let   () denote the space of polynomials on  of degree no more than .We also define for all  ∈  ℎ }︁ .
We suppose that the velocity space  ℎ is of order  ( ≥ 1), i.e., there exists a non-negative integer  such that   ⊂  ℎ and  +1 ̸ ⊂  ℎ .We assume the existence of a divergence-free reconstruction operator Π ℎ :  ℎ →  ℎ (we do not give the detailed definition here) which satisfies that ) and (, Note that (2.4) can be derived from (2.3) and (2.5).We also assume that Π ℎ satisfies the following properties.
Assumption 1.There exist two operators Π 1 ℎ :  ℎ →  ℎ and Π R ℎ :  ℎ →  ℎ and two positive constants Remark 1.The fundamental requirement of the projection Π ℎ is: For any  ℎ ∈  ℎ , Π ℎ  ℎ can be decomposed into a sufficiently approximate  1 -conforming component and a "small" (div)-conforming component (consider Π 1 ℎ  ℎ and Π R ℎ  ℎ ).This is the prerequisite to constructing our methods.Regarding the discretization of the convection term (see (2.14) below), a rough description of our basic idea is: Apply the  1 -conforming part to guarantee the accuracy and "abuse" the (div)-conforming part to guarantee the conservation of energy, momentum and angular momentum.In Section 4, we shall show that the reconstruction operators in Remark 4.2 of [43] and their higher order versions on a class of non-divergence-free simplicial elements satisfy all the above properties.

The EMAPR method for classical elements
and (0) =  0 .A straightforward semi-discrete analog of (2.9) is to find ( ℎ ,  ℎ ) :  →  ℎ ×  ℎ satisfying  ℎ (0) =  0 ℎ ∈  ℎ with  0 ℎ being some approximation of  0 and ) for all ( ℎ ,  ℎ ) ∈  ℎ ×  ℎ .However, it is well-known that the above scheme is not energy-stable and pressurerobust unless  ℎ is exactly divergence-free (or equivalently, ∇ •  ℎ ⊆  ℎ ).To obtain a pressure-robust velocity for classical elements, in [43] Linke and Merdon proposed a novel finite element formulation which reads Here ( ℎ ,  ℎ ) should be interpreted as (Π ℎ  ℎ ,  ℎ ) via property (2.5).By using divergence-free reconstructions, the above formulation restores the  2 -orthogonality between discretely divergence-free test functions and gradient fields, and thus remove the effect of the continuous pressure in velocity errors.There is some consistency error arising from the diffusion term.
To make the method energy-stable, a pressure-robust and energy-conserving rotational form (ROT) of the nonlinear term was also proposed in [43]: where ∇× is the curl operator [25].However, following the strategy in [12] or Section 2.3 below, one can prove that the rotational form does not preserve momentum and angular momentum (see Rem. 3 below).To resolve this issue, we propose an EMA-conserving form, which results in the EMAPR reconstruction: ) for all  ℎ ∈  ℎ ,  ℎ ∈  ℎ and  ℎ (0) =  0 ℎ .Here  ℎ is given by where  is a positive parameter.The trilinear form  ℎ is defined by Remark 2. The bilinear form  ℎ with  = 0 is commonly used in pressure-robust reconstructions for the (, )-like term [2,43].However, for the case  = 0 we can not obtain a Re-semi-robust estimate theoretically.
In practice, we find that the stabilization term )︀ is of importance for high order elements ( ≥ 2) in the case that  is very small or equal to zero (the Euler equation).In this case, without this term, the  1 error of the discrete velocity might be large.

EMA-balance in semi-discrete schemes
Now we are in the position to analyze the balance laws with the EMAPR reconstruction.We define kinetic energy  :  → R, linear momentum  :  → R  , and angular momentum   :  → R 3 by for any  * ∈ .For any two-dimensional vector  * = ( * 1 ,  * 2 ), to compute angular momentum or cross product, one can always embed it into three dimensions by setting  * = ( * 1 ,  * 2 , 0).Let  be the solution of (2.9) which satisfies the following balance laws [12,49]: where the balance laws of momentum and angular momentum are based on some appropriate assumptions.
Similarly to the assumptions in [12,13,49], here (only for the analysis of momentum and angular momentum) we assume that (, ) is compactly supported in Ω (e.g., consider an isolated vortex).We also define a discrete energy for any  * ∈  0 ℎ , which is asymptotically convergent to the real energy.Note that we have (Π ℎ  * ) ≤   ( * ).The following lemma is essential for EMA analysis.
Lemma 2. For any Using a chain rule like 12) and applying (2.12b) and Lemma 2, one immediately obtains the following theorem.
Theorem 1.Let  ℎ be the solution of (2.12).Then it satisfies the following balance of energy: Boundary conditions may influence the balance of momentum and angular momentum.For simplicity, we ignore the effect of boundary here by doing some extra assumptions which are similar to those of the continuous case.Similar assumptions were also used for the analysis of the EMAC formulation in [12].
Assumption 2. The finite element solution ( ℎ ,  ℎ ), Π ℎ  ℎ and the external force  are only supported on a sub-mesh Tℎ ⊂  ℎ such that there exists an operator  : ) Here   ∈ R  is the unit vector whose -th component equals 1.
In fact, for discontinuous pressure elements (like the elements we used herein), the support of Π ℎ  ℎ is the same as that of  ℎ , since the reconstruction can be locally performed on each element.Furthermore, note that ∇ •   = ∇ • ( ×   ) = 0 and   and  ×   are linear polynomials.So the equalities for Π ℎ and Π 1 ℎ in (2.17) are not hard to satisfy.The reconstruction operators in Section 4 fulfill these equalities.
For the case  = 0 (the Euler equations), if we apply the no-penetration boundary condition ( •  = 0 on  × Γ), the skew-symmetry of  ℎ still holds by Lemma 1. Thus Theorem 1 implies that the method (2.12) conserves a discrete energy for  = 0 and  = 0, with only no-penetration boundary condition strongly imposed.Theorem 2 implies that (2.12) conserves linear momentum and angular momentum (of Π ℎ  ℎ ) for  with zero momentum and zero angular momentum under Assumption 2, respectively.Remark 3 (Momentum analysis for the rotational form (2.11)).We use linear momentum as an example.By equation ( 8) of [12], we have Then taking  ℎ = (  ) gives Thus linear momentum is not preserved by  rot .
Split the error  −  ℎ as We also define Next, we introduce the dual norm ‖•‖ ( 0 ℎ ) ′ for any linear functional  on  0 ℎ : or for any  ∈  2 (Ω): We also define a mesh-dependent norm ||| • ||| * on  ℎ by The following two inequalities will be used to estimate the nonlinear terms.
Lemma 3. Let  inv denote a generic constant in inverse inequalities.Then one has ) Proof.It follows from the Schwarz's inequality and the inverse inequality that For the second inequality, from (2.8) similarly we have This completes the proof.
Theorem 3. Let  be the solution of (2.9) and  ℎ be the solution of (2.12).Under Assumptions 1, 3, and the assumptions that , and  0 ℎ = Π S ℎ  0 , with  independent of ℎ and  the following estimate holds: where Proof.Substituting (3.5) into (3.2) and taking  ℎ =  ℎ give that Let us estimate each term in (3.11).For the evolutionary term we have For  ℎ term the estimate can be found in [43]: Now, let us estimate the convective terms.We use a similar decomposition with [49,55]: It follows from the Schwarz's inequality, Young's inequality, Lemma 3, Assumptions 1, and 3 that )︂ .
Based on the above estimates for the nonlinear terms, one can obtain that where Then substituting (3.12)-(3.14)into (3.11)provides ) by Lemma 3 of [9] with ""= Π 1 ℎ  ℎ and ""= Π R ℎ  ℎ therein.Integrating over , and applying the fact  ℎ (0) = Π S ℎ (0) and the Gronwall inequality, we can get Then (3.9) follows immediately from a combination of the above inequality and the Young's inequality.
Remark 4 (The choice of ).For the error estimate (3.9), on the one hand, a larger  results in a smaller Gronwall constant (,  ), but on the other hand, the existence of the terms does not allow a too large .The optimal choice of  is still an open problem.By observing the Gronwall constant, it seems that  ≪ 1 is not a good choice since the -related part can dominate the other part with such an .Some numerical study on  will be given in Section 5.

The reconstruction operators for a class of simplicial non-divergence-free elements
In this section, we focus on a class of simplicial locally mass-conserving elements which satisfy the inf-sup condition (2.2), and give the corresponding divergence-free reconstruction operators.First, let us recall their Table 1.Local velocity spaces on an element .For  = 1, one usually calls it as Bernardi-Raugel space.
Consider an arbitrary element  with vertices   , 1 ≤  ≤  + 1. Denote by   the face opposite to   and   the unit outward normal vector corresponding to   , 1 ≤  ≤  + 1.Further,   , 1 ≤  ≤  + 1, denote the corresponding barycentric coordinates.Then the (normal-weighted) face bubbles are defined by We also define , and P () := span Then the local finite element spaces for velocity on an element  are defined as Table 1.
For -th order velocity spaces, the matching pressure space is the space of discontinuous piecewise polynomials of degree no more than  − 1, whatever the dimension is.In what follows,  ℎ will denote a velocity space of order  mentioned above and  ℎ is the corresponding pressure space.From Table 1 one can see there exists a space of bubble functions   ℎ ∼ =  ℎ /  such that  ℎ =   ⊕   ℎ .For any  ℎ ∈  ℎ , it is natural to split it into two parts: We consider a class of divergence-free reconstruction operators which were also discussed in Remark 4.2 of [43].These operators are defined as follows.Let  ℎ :  0 (︀ Ω)︀ →   denote the standard interpolation operator to   (such as the interpolations defined in [25], pp.99, 100).Then Π 1 ℎ is defined by Let Π RT ℎ be the common Raviart-Thomas interpolation of order  − 1 [7].The operator Π R ℎ is defined by Hence, for any  ℎ ∈  ℎ one has At this time the space  ℎ could be chosen as for all  ∈  ℎ }︁ .
In other words, these reconstruction operators only change the bubble part of the elements.Thus the reconstruction is low-cost and has not changed much compared to the previous classical formulation.This is especially the case for the first order element and the second order element in two/three dimensions, since it is not hard to find that, for any  ℎ belonging to these spaces, Remark 6.Another advantage of EMAPR is that it can lower the requirement of the quadrature rules dramatically.For classical discretizations, usually the convective term requires higher order quadrature rules.The use of bubble functions exacerbates the situations further since these functions are polynomials of higher degree.This limits the practical application of the simplicial locally mass-conserving elements.One can easily check that, to compute the convective term (•, •, •) exactly, the quadrature rules with algebraic precision of degree 5 ( 8) are needed for two-dimensional (three-dimensional, respectively) Bernardi-Raugel element.However, in any dimensions for Bernardi-Raugel elements one only needs the quadrature rules with algebraic precision of degree 2 to compute  ℎ (•, •, •) exactly.

On the implementation of Π ℎ for Bernardi-Raugel elements
This subsection discusses a detailed construction of Π ℎ for Bernardi-Raugel elements.The reconstructed shape functions are explicitly provided, which can be used to assemble coefficient matrices for our methods in a standard way.Recall that   ∈ R  denotes the -th vertex of .On an element , the shape functions for the lowest order Raviart-Thomas space have the form The relationship between the shape functions and reconstructed shape functions for the Bernardi-Raugel elements is listed in Table 2. Indeed, simple calculations give Then Table 2 follows immediately from the definition of Π ℎ .One can use Table 2 and standard Raviart-Thomas basis functions to build the basis functions for Π ℎ  ℎ easily and locally.Then the matrix assembly of our method is very similar to that of the classical conforming methods, since the former also consists of volume integrals.
Lemma 4. The operators  ℎ and Π RT ℎ satisfy and for all  ∈  ℎ .
Proof.We refer the readers to Theorem 4.4.
where we repeatedly use the local estimates in Lemma 4.5.3 of [8] and the interpolation error of Π RT ℎ (4.5).Then (4.6) follows immediately from the triangle inequality.Note that (4.6) does not hold if  ℎ is an arbitrary function in  .
Up to now, we have proven that a class of divergence-free reconstruction operators satisfies the assumptions in Section 2.1.Thus it admits an a priori error estimate in Theorem 3. The only question is whether the right-hand side of (3.9) is an  (︀ ℎ  )︀ quantity if  is sufficiently smooth.The non-trivial terms are the ones corresponding to   ,  1  ,  R  , and ‖∆ ∘ (1 − Π ℎ )‖ ( 0 ℎ ) ′ .The following lemmas are to answer this question.
Finally, based on the results in Theorem 3, Lemmas 6 and 7, as well as the approximation properties of  ℎ and Π S ℎ [8,25], we get the convergence rates of the kinetic and dissipation energy errors of  ℎ (or Π ℎ  ℎ ) for (2.12), with the elements and reconstruction operators in this section.

Numerical experiments
The meaning of some abbreviations used in the legends below is listed in Table 3.
Table 3.Some abbreviations used in legends for numerical experiments and their meaning.Divergence-free second order Scott-Vogelius elements [5] CONV L-M 2016 Pressure-robust reconstruction methods in [43] with convective form ROT L-M 2016 Pressure-robust reconstruction methods in [43] with rotation form GD Grad-div stabilization [50] 5.1.Example 1: the lattice vortex problem In the first example, we consider the lattice vortex problem [49,55] on Ω = (0, 1) T with an appropriate ,  fulfills an exact unsteady NSE with  = 0. We choose a small  and a large  :  = 1 × 10 −5 and  = 10.We do some numerical experiments by Bernardi-Raugel elements on an unstructured mesh with 14 112 triangles and 49 778 DoFs.To make a comparison, the EMAC formulation and the pressure-robust reconstruction method in [43] with the rotation formulation are also tested.Further, considering that the low order Taylor-Hood element is much more popular in practice, we also show some results from this element with EMAC and/or the grad-div stabilization.The well-known second order Scott-Vogelius (SV2) element,  2 / disc 1 [5,34], with the convective formulation, is also used to make a comparison, which is a kind of divergence-free method.We try to make all the methods have almost the same number of DoFs.The methods with Taylor-Hood element and SV2 element are run on a mesh with 11 018 triangles (50 154 DoFs) and a barycentric-refined mesh with 7068 triangles (49 686 DoFs), respectively.For the time discretizations, we use the BDF2 scheme with ∆ = 0.001.At each time step we solve a nonlinear system by the Newton iteration method with a tolerance of 10 −8 in the  1 norm.
We also do some comparison between the results from different .And for  = 1 we try two cases: the full related stabilization and the only left-hand one.The latter means that the -related stabilization is only added to the coefficient matrix, that is, is not for the time step   .In the legends we use " = 1 (L)" to represent this case.For other  we only use the full one.
Some results are shown in Figures 1-3.From Figure 1 one can see the errors of EMAC and pressure-robust ROT formulations with the Bernardi-Raugel element grow quickly at an early time.We conjecture that the reason for EMAC is that it is not pressure-robust and the reason for pressure-robust ROT is that it is not Re-semi-robust.The former can be considerably improved by using the low order Taylor-Hood element, since its pressure approximation order is one order higher than Bernardi-Raugel pressure (but the number of DoFs of the pressure is instead less due to the continuity).For all the methods in Figure 1, EMAPR with  = 1 (only left-hand) shows best performance.And the cost of this method is also the second least (see Tab. 4).
In Figure 2 we show some comparison with the grad-div stabilized EMAC (the stabilization coefficient is always taken to be 0.1) and Taylor-Hood elements.We also consider a more complicated pressure there, by setting  = 10 2 ∇(), 10 3 ∇().Note that this changes the pressure but not the velocity.Without a modified continuous pressure, the grad-div stabilized EMAC gives a better performance than our method.However, as the pressure becomes more complicated, its performance becomes worse.For EMAPR method, more complicated pressures do not affect the discrete velocity.From Figure 3 one can see the results from  = 0 and  = 10 −3 are virtually the same, which further verify that a small  is not a meaningful choice.The full stabilization with  = 10 3 gives the best accuracy before  = 8, while the only left-hand stabilization with  = 1 performs best at the final time.

Example 2: potential flow
For the second example we consider the potential flow in Example 6 of [43].On Ω = (0, 1) 2 the velocity is prescribed as  = min{, 1}∇ with  =  3  −  3 .We set  = 0 such that the pressure gradient exactly balances the gradient filed produced by the velocity terms.Due to the quadratic convective term, the pressure is much more complicated than the velocity [24].The pressure-robustness will play a key role on the accuracy of the simulations of this problem.  .We set ∆ = 0.01 and  = 2.We also give some results from the classical scheme and the pressure-robust reconstruction scheme in [43] with the convective form for the nonlinear term.In each time step, we solve a linear problem by replacing the advective velocity in trilinear forms with a linear extrapolation of the previous step velocities ( ).Some results are shown in Tables 5 and 6.All the pressure-robust schemes have a much better accuracy than the classical scheme.

Example 3: flow passing a cylinder in two dimensions
This example describes a time-dependent flow passing a cylinder in two dimensions [31,54].The computational domain is set as Ω = (0, 2.2) × (0, 0. such that the mean inflow/outflow speed is  () = sin(/8) and max 0≤≤8  () = 1.On all other boundaries the no-slip boundary condition ( = 0) is enforced.The external force is  = 0 and the viscosity is  = 10 −3 such that the Reynolds number Re satisfies 0 ≤ Re ≤ 100.The flow starts with static ( 0 = 0).In this example a vortex street develops behind the cylinder around  = 4 [31].
The geometry setting and the mesh we used can be found in Figure 4. We employee the Bernardi-Raugel element to check the performance of our method, whose number of DoFs is 41368.The second order Scott-Vogelius pair (with the standard convective form for the nonlinear discretization) is used as a reference and comparison, which runs on the barycentric refinement of the mesh in Figure 4, producing totally 245 898 DoFs.For time discretizaton, the BDF2 method is used with ∆ = 0.0025.On each time step a linearized system is solved by replacing the advective velocity  ,  ℎ ) and  ⋆ =  for the Scott-Vogelius element.Some contours of speed can be found in Figure 5.The results from our method is very close to those from the Scott-Vogelius element method, and they are all very close to the plots in Figure 6 from [18].Some values corresponding to maximum drag and lift are listed in Table 7.We try different  from 0 to 10 3 .In the paper [54], the reference intervals of maximum drag and lift are given as [2.93, 2.97] and [0.47, 0.49], respectively.It seems that the -related stabilization is not necessary for this example.One possible reason is that this is a low Reynolds number flow.Our methods with  = 0, 0.1, 1 give very similar results.On the other hand, we also find that for a larger  ( = 10, 10 3 ) the maximal lift coefficient can fall outside of the reference interval [0.47, 0.49], which implies that the choice of  does have an influence on the performance of our method, although this influence seems to be not large here.We conjecture that for low Reynolds number flows a small  ( ≤ 1) is sufficient and a larger  is not recommended.

Conclusions
We have developed an EMA-conserving, pressure-robust and Re-semi-robust reconstruction method for a class of finite element methods of the Navier-Stokes equations.The lowest order case is the Bernardi-Raugel element.In practice, we conjecture grad-div stabilized EMAC is a very simple and efficient method in many situations.However, when the pressure is much larger or more complicated than the velocity, besides divergence-free elements, we believe that the method proposed herein can also be a good choice due to its easier construction and mild requirement of the meshes.Further research on stabilizations is an important part of future work.

Figure 1 .
Figure 1.Example 1.  2 errors (left) and  1 errors (right) over time for different variants/modifications of the Bernardi-Raugel finite element method, Taylor-Hood finite element method with EMAC form, and Scott-Vogelius finite element method with convective form.

Figure 3 .
Figure 3. Example 1.  2 errors (left) and  1 errors (right) over time for the Bernardi-Raugel finite element method with EMAPR reconstruction and different parameters/choices of the stabilization.

Table 4 .
The CPU time cost of some methods in Example 1.The "Total iterations" term refers to the total number of Newton iteration steps.

Table 5 .
Example 2. Errors for different variants/modifications of the Bernardi-Raugel finite element method (classical CONV/EMAPR  = 0/CONV L-M 2016).We consider the case of  = 5 × 10 −4 and apply the Bernardi-Raugel element and the second order element in

Table 6 .
Example 2. Errors for different variants/modifications of  bubble

Table 7 .
Example 3. Maximal values of drag and lift for the Bernardi-Raugel finite element method with EMAPR reconstruction and different choices of  and Scott-Vogelius finite element method with convective form. ,max )  ,max ( ,max )  ,max