This paper is concerned with the finite element discretizations of the unsteady Navier-Stokes equations equationparentequation
where () is a bounded domain with Lipschitz and polyhedral boundary ; the velocity and pressure are two unknowns; is the constant kinematic viscosity; represents the external body force at ; and denotes the initial velocity.
Physically, there are several balance laws implied in the Navier-Stokes equations such as the balances of kinetic energy, linear momentum and angular momentum (EMA) [14, 7]. In some cases (e.g., vanishing external force and viscosity), these quantities are conservative under some appropriate assumptions . Since at least Arakawa devised an energy and enstrophy conservation method for two-dimensional NSEs in , scientists have recognised that these conservation (balance) laws are important references for unlocking efficient and stable numerical methods for NSEs [15, 39, 38, 1, 14, 36]. A violation of these laws may bring unexpected instability.
In this paper, we consider inf-sup stable mixed methods for NSEs. The crucial points for these balance laws are the discretization of the nonlinear term and the finite elements one uses. The most common formulation for this term might be the so-called convective formulation (CONV). However, it was shown in  that this formulation was not EMA-conserving unless exactly divergence-free elements were used. We note that most classical elements are not divergence-free . Also in , Charnyi et al. proposed an EMA-conserving formulation (”EMAC”) for the common elements. The numerical experiments in [7, 29, 28, 37] demonstrated that this formulation was always among the best discretizations of the nonlinear term. In addition, the EMAC formulation was also shown to preserve the 2D enstrophy, helicity and total vorticity under certain appropriate assumptions. In the paper , Olshanskii and Rebholz proved that the Gronwall constant in the EMAC error estimates does not explicitly depend on the Reynolds number, which is contrast to other formulations such as SKEW. We note that, although the EMAC formulation has many fascinating properties, it is not easy to find a fully satisfactory linearization scheme for it. In , Charnyi et al. proposed two linearization schemes for this method: a skew-symmetric linearization and the Newton linearization. The former did not conserve momentum and angular momentum, which gave a poor performance for some problems, the latter conserved momentum and angular momentum but was not energy-stable. The objective of this paper is to design an alternative formulation that is very suited to Picard linearization.
The idea for this paper goes back to a class of locally discontinuous Galerkin (DG) methods for NSEs, cf. . By replacing one of the velocity with its a divergence-free approximation in the nonlinear term, a class of energy-stable DG methods were developed in . Similar techniques were also applied in [21, 26] and to coupled flow-transport problems [34, 19, 24]. Meanwhile, seeking an exactly divergence-free approximation to some discretely divergence-free velocity is also an objective of the pressure-robust reconstruction methods [32, 27, 31], where a large class of divergence-free reconstruction operators were proposed, almost covered all the classical conforming elements. In this paper, we will use these operators to reconstruct CONV for conforming elements.
One motivation for our methods is that, the CONV formulation not only has the easiest form but also usually shows good performance before the scheme blows up due to the energy instability (e.g., see the comparison experiments in ). In our opinion, it has the potential to simulate a problem accurately and the divergence-free reconstruction is one of the keys to unlock this potential. Here we focus on unsteady NSEs and give some theoretical analysis of the reconstructed CONV formulation. In what follows, we shall refer to it simply as “modified CONV” sometimes. We prove that modified CONV conserves kinetic energy, linear momentum, total vorticity, 2D enstrophy and helicity under the same assumptions as EMAC . Further, we also prove that the Picard linearization conserves them as well. Compared to EMAC, our formulation does not conserve angular momentum theoretically. However, the numerical experiments show that modified CONV preserves this quantity well. For the semi-discrete scheme, we prove that the Gronwall constant of the error bounds does not depend on the Reynolds number explicitly under the assumption . In this aspect a related concept is called (Re-)semi-robustness [40, 22], which is a type of robustness that the constants (including Gronwall constant) in velocity error estimates do not explicitly depend on the inverse of the viscosity. For Re-semi-robust methods/analysis, we refer to [42, 41, 40, 6, 11, 12, 16].
The rest of this paper is organized as follows. In Section 2 we state the reconstructed scheme and discuss the conservative properties and error estimates of the semi-discrete scheme. Section 3 is devoted to investigating the conservative properties of Picard linearization, matching the Crank-Nicolson time stepping method. Section 4 studies some numerical experiments.
Throughout the paper we use with or without subscript, to denote a generic positive constant. The norm (seminorm) of the Sobolev space or () is denoted by (, respectively). The standard inner product of or is denoted by . When (, ), with the convention that the index (, , respectively) is omitted.
2. A modified convective formulation
Let be shape-regular partitions of  and with the diameter of element . Let , and denotes a pair of inf-sup stable finite element spaces on . Here we do not give them a concrete definition first.
A continuous variational form for eq. 1.1 is: equationparentequation
for any and . Here the trilinear form is the so-called convective formulation (CONV).
A semi-discrete analog of eq. 2.1 reads: equationparentequation
The following identity is widely used in energy-stability analysis:
Equation 2.3 gives a skew-symmetric structure for and in when is divergence-free. Based on this fact, one could further prove that the continuous velocity is energy-stable in eq. 2.1. In the discrete level, however, the divergence constraint is usually relaxed for most of the classical elements such as the Bernardi-Raugel elements  and the Taylor-Hood elements [18, 5]. Then the skew-symmetric structure loses. The discretely divergence-free subspace of is defined by
To fix the lack of skew-symmetry, one commonly used method is to modify the CONV formulation into the following SKEW formulation:
However, Olshanskii and Rebholz  showed that the SKEW formulation might give a poor accuracy for long time simulations with higher Reynolds number. In this paper, we consider another skew-symmetrization way for CONV, where we shall use the divergence-free projection operators.
the unit external normal vector on. Denote by a generic divergence-free projection operator, which may be different for different space pairs. Introduce a corresponding finite element space such that . Here we assume that has the following properties:
➀ for all ; ➁ there exists a constant independent of such that for all .
For most cases, the inequality in Assumption 1 ➁ can be obtained by a combination of the approximation property of , the inverse inequality and triangle inequality.
Introduce the modified convective formulation (modified CONV)
We consider the reconstructed scheme: equationparentequation
The following lemma is very essential for our analysis.
For any we have
or or on .
Lemma 1 can be regarded as a special case of the skew-symmetric discontinuous Galerkin formulation (see, e.g., [13, Lemma 6.39]). Here for completeness we reprove this special case. To make the proof process more clear, we will derive it in the component form. Let denotes the component form of any . We have
The integration by parts on each element gives
where denotes the -th component of the unit external normal vector of , . A combination of the continuity and boundary conditions of and implies that
and the first condition in Lemma 1 gives
2.1. Conservation laws
In this subsection, we focus on the conservative properties of the modified method eq. 2.7. First, let us give the definitions of momentum, angular momentum, helicity, enstrophy and vorticity. Denote by the exact velocity in eq. 1.1 or the discrete velocity in eq. 2.7. Let when represents the continuous solution; let be the solution of some finite element discretization of the Navier-Stokes vorticity equation (see eq. 2.14 and eq. 2.15 below) when is the discrete velocity. Note that the operator curl in two dimensions denotes a scalar operator . Then we define
In what follows we shall discuss some conservative properties for the case (the Euler equations) sometimes. For an Euler equation, the boundary conditions should be altered in eq. 1.1. Here we replace the no-slip boundary condition () with the no-penetration boundary condition () for it. In this case the method eq. 2.7 should only strongly enforce the no-penetration boundary condition also. By abuse of notation, we shall use eq. 1.1 and eq. 2.7 to denote the Euler equation and the corresponding finite element methods as well, respectively.
When and , one can similarly obtain . In this case, kinetic energy is conserved by our method.
To analysis the conservative properties of the other quantities, we need some extra assumptions, which are similar to the EMAC analysis .
The exact solution , the discrete solution and the source term are only supported on a sub-domain such that there exists restriction for with in and on . Here represents the -th usual basis of .
We divide the proof of the above theorem into several subsections.
For the analysis of angular momentum, note that . Then substituting into eq. 2.7, one could similarly obtain that
Here we apply the fact that and . Assuming that has zero angular momentum, i.e., , the above equation gives that since in general. Thus angular momentum is not preserved by our methods.
Let be the curl of the exact solution : . Note that . In two dimensions, represents a scalar and we use the symbol to denote it. Taking the curl operator on eq. 1.1a gives the following Navier-Stokes vorticity equations:
For the discrete case, we apply the strategy in [35, 7]. Here we do not analyze the quantity for the discrete solution . We consider a sightly altered virticity which is a solution of some direct discretization of eq. 2.14. As it was said in , ‘this discrete vorticity still depends on the computed velocity , but more implicitly, through the equation coefficients’. We further assume vorticity also vanishes on and near the boundary due to Assumption 2. The corresponding finite element methods to eq. 2.14 reads: Find and the Lagrange multifier such that
2.1.3. 2D Enstrophy
In two dimensions the vorticity satisfies
The corresponding finite element scheme reads:
Similarly, we have . Thus, taking and for the case we have
Thus we prove the conservative property of the 2D enstrophy.
2.2. Semi-discrete error analysis
Denote by the Stokes projection which satisfies that
The Stokes projection satisfies the following estimate
The following error equation will be used in error analysis.
The above equation can be further rewritten as
where we use the fact that .
Theorem 2.3 (Semi-discrete estimate).
Set in eq. 2.28 and we arrive at
Applying the Cauchy-Schwarz inequality and the weighted Young’s inequality one obtains
Then the Cauchy-Schwarz inequality, Equation 2.26 and the Young’s inequality give that
Then integrating over and using the Gronwall inequality and Hölder inequality we finally obtain
3. The Picard linearization scheme
In practice, a commonly used linearization way is replacing one of the velocity solutions of the nonlinear term with last time step solutions. In this section we prove that this way preserves all the conservative properties from the semi-discrete version when matching the Crank-Nicolson time discretizations. The linearized Crank-Nicolson scheme is that
for all , where