# Longer time simulation of the unsteady Navier-Stokes equations based on a modified convective formulation

For the discretization of the convective term in the Navier-Stokes equations (NSEs), the commonly used convective formulation (CONV) does not preserve the energy if the divergence constraint is only weakly enforced. In this paper, we apply the skew-symmetrization technique in [B. Cockburn, G. Kanschat and D. Schötzau, Math. Comp., 74 (2005), pp. 1067-1095] to conforming finite element methods, which restores energy conservation for CONV. The crucial idea is to replace the discrete advective velocity with its a H(div)-conforming divergence-free approximation in CONV. We prove that the modified convective formulation also conserves linear momentum, helicity, 2D enstrophy and total vorticity under some appropriate senses. Its a Picard-type linearization form also conserves them. Under the assumption u∈ L^2(0,T;W^1,∞(Ω)), it can be shown that the Gronwall constant does not explicitly depend on the Reynolds number in the error estimates. The long time numerical simulations show that the linearized and modified convective formulation has a similar performance with the EMAC formulation and outperforms the usual skew-symmetric formulation (SKEW).

## Authors

• 42 publications
• 3 publications
05/10/2022

### Improved long time accuracy for projection methods for Navier-Stokes equations using EMAC formulation

We consider a pressure correction temporal discretization for the incomp...
02/04/2020

### Longer time accuracy for incompressible Navier-Stokes simulations with the EMAC formulation

In this paper, we consider the recently introduced EMAC formulation for ...
10/09/2020

### An energy, momentum and angular momentum conserving scheme for a regularization model of incompressible flow

We introduce a new regularization model for incompressible fluid flow, w...
06/12/2020

### Continuous data assimilation applied to a velocity-vorticity formulation of the 2D Navier-Stokes equations

We study a continuous data assimilation (CDA) algorithm for a velocity-v...
09/01/2020

### Pressure-robust error estimate of optimal order for the Stokes equations on domains with edges

The velocity solution of the incompressible Stokes equations is not affe...
11/06/2019

### Variational multiscale modeling with discretely divergence-free subscales

We introduce a residual-based stabilized formulation for incompressible ...
07/15/2020

### Versatile Mixed Methods for the Incompressible Navier-Stokes Equations

In the spirit of the "Principle of Equipresence" introduced by Truesdell...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

This paper is concerned with the finite element discretizations of the unsteady Navier-Stokes equations equationparentequation

 (1.1a) ut−νΔu+(u⋅∇)u+∇p =f in (0,T]×Ω, (1.1b) ∇⋅u =0 in (0,T]×Ω, (1.1c) u(0) =u0 in Ω, (1.1d) u =0 on (0,T]×Γ,

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 [7]. Since at least Arakawa devised an energy and enstrophy conservation method for two-dimensional NSEs in [2], 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 [7] that this formulation was not EMA-conserving unless exactly divergence-free elements were used. We note that most classical elements are not divergence-free [24]. Also in [7], 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 [36], 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 [8], 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. [10]. 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 [10]. 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 [7]). 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 [7]. 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 [9] 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

 Find u:(0,T]→V and p:(0,T]→W such that (2.1a) (∂u∂t,v)+νa(u,v)+c(u,u,v)−b(v,p) =(f,v)∀ v∈V, (2.1b) b(u,q) =0 ∀ q∈W,

where

 a(u,v)=(∇u,∇v),
 c(u,v,w)=((u⋅∇)v,w),

and

 b(v,q)=(∇⋅v,q),

for any and . Here the trilinear form is the so-called convective formulation (CONV).

A semi-discrete analog of eq. 2.1 reads: equationparentequation

 Find (uh,ph):(0,T]→Vh×Wh such that (2.2a) (∂uh∂t,vh)+νa(uh,vh)+c(uh,uh,vh)−b(vh,ph) =(f,vh)∀ vh∈Vh, (2.2b) b(uh,qh) =0  ∀ qh∈Wh.

Let be the velocity solution of eq. 2.1 or eq. 2.2. Introduce the kinetic energy:

 Kinetic energy E:=12∫Ω|u∗|2 dx.

The following identity is widely used in energy-stability analysis:

 (2.3) c(u,v,w)=−c(u,w,v)−((∇⋅u)v,w)∀ u,w,v∈V.

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 [4] and the Taylor-Hood elements [18, 5]. Then the skew-symmetric structure loses. The discretely divergence-free subspace of is defined by

 (2.4) V0h={vh∈Vh:b(vh,qh)=0∀ qh∈Wh}.

To fix the lack of skew-symmetry, one commonly used method is to modify the CONV formulation into the following SKEW formulation:

 (2.5) cskew(u,v,w)=c(u,v,w)+12((∇⋅u)v,w).

However, Olshanskii and Rebholz [36] 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.

We define

 H0(div;Ω)={v∈H(div;Ω):v⋅n|Γ=0},

with

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:

###### Assumption 1.

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)

 (2.6) ch(uh,vh,wh):=c(Πhuh,vh,wh).

We consider the reconstructed scheme: equationparentequation

 (2.7a) (∂uh∂t,vh)+νa(uh,vh)+ch(uh,uh,vh)−b(vh,ph) =(f,vh)∀ vh∈Vh, (2.7b) b(uh,qh) =0 ∀ qh∈Wh.

The following lemma is very essential for our analysis.

###### Lemma 1.

For any we have

 c(uh,vh,wh)=−c(uh,wh,vh)

if

• ;

• or or on .

###### Proof.

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

 (2.8) c(uh,vh,wh)=((∇vh)uh,wh)=∑1≤i,j≤d∑T∈Th∫Tui∂vj∂xiwj dx.

The integration by parts on each element gives

 (2.9) ∫Tui∂vj∂xiwj dx =−∫T∂(uiwj)∂xivj dx+∫∂Tuinivjwj ds =−∫T(ui∂wj∂xivj+∂ui∂xivjwj) dx+∫∂Tuinivjwj ds,

where denotes the -th component of the unit external normal vector of , . A combination of the continuity and boundary conditions of and implies that

 (2.10) ∑1≤i,j≤d∑T∈Th∫∂Tuinivjwj ds=∑T∈Th∫∂T(uh⋅n)(vh⋅wh) ds=0,

and the first condition in Lemma 1 gives

 (2.11) ∑1≤i,j≤d∑T∈Th∫Tui∂wj∂xivj+∂ui∂xivjwj dx =c(uh,wh,vh)+((∇⋅uh)vh,wh) =c(uh,wh,vh).

Then Lemma 1 follows immediately from eqs. 2.11, 2.10, 2.9 and 2.8. ∎

By Lemma 1 and Assumption 1 one can obtain

 (2.12) ch(uh,vh,wh)=−ch(uh,wh,vh) for all uh∈V0h,vh,wh∈Vh.

### 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 [18]. Then we define

 Linear momentum M:=∫Ωu∗ dx;
 Angular momentum Mx:=∫Ωu∗×x dx;
 Helicity H:=∫Ωu∗⋅w∗ dx;
 Enstrophy HE:=12∫Ωw∗⋅w∗ dx;
 Total vorticity W:=∫Ωw∗ dx.
###### Remark 1.

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.

###### Theorem 2.1.

Under Assumption 1 ➀, the modified method eq. 2.7 properly balances the kinetic energy:

 (2.13) 12ddt||uh||2+ν∥∇uh∥2=(f,uh).
###### Proof.

Equation 2.12 implies that . Then eq. 2.13 follows immediately from taking in eq. 2.7a. ∎

###### Remark 2.

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 [7].

###### Assumption 2.

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 .

###### Theorem 2.2.

Under Assumption 1 ➀ and Assumption 2, the modified method eq. 2.7 conserves momentum (for with zero linear momentum), helicity (for and ), 2D enstrophy (for and ), and total vorticity.

We divide the proof of the above theorem into several subsections.

#### 2.1.1. Momentum

Testing with in eq. 2.7 and applying eq. 2.12 give that

 (f,ei)=ddt(uh,ei)+ch(uh,uh,ei)=ddt(uh,ei)−ch(uh,ei,uh)=ddt(uh,ei).

Then the conservation of momentum follows.

###### Remark 3.

For the analysis of angular momentum, note that . Then substituting into eq. 2.7, one could similarly obtain that

 ddt(uh,x×ei)−ch(uh,x×ei,uh)=(f,x×ei).

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.

#### 2.1.2. Helicity

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:

 (2.14) wt+(u⋅∇)w−(w⋅∇)u−νΔw=∇×f.

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 [7], ‘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.15) (∂wh∂t,vh)+ch(uh,wh,vh) −ch(wh,uh,vh)+ ν(∇wh,∇vh)− b(vh,ηh)+b(wh,qh)=(∇×f,vh),

for any and , where is the solution of eq. 2.7. Note that eq. 2.15 implies that for all and further .

Testing with in eq. 2.7 gives

 (2.16) (∂uh∂t,wh)+ν(∇uh,∇wh)+ch(uh,uh,wh)=(f,wh).

Meanwhile, testing with in eq. 2.15 gives

 (2.17) (∂wh∂t,uh)+ch(uh,wh,uh)−ch(wh,uh,uh)+ν(∇wh,∇uh)=(∇×f,uh).

Adding eq. 2.16 and eq. 2.17 and setting and , one obtain

 ddt(uh,wh)+ch(uh,wh,uh)−ch(wh,uh,uh)+ch(uh,uh,wh)=ddt(uh,wh)=0,

where we use the property eq. 2.12. Thus helicity is also conserved.

#### 2.1.3. 2D Enstrophy

In two dimensions the vorticity satisfies

 wt+(u⋅∇)w−νΔw=∇×f

The corresponding finite element scheme reads:

 (2.18) (∂wh∂t,v)+((Πhuh⋅∇)wh,vh)+ν(∇wh,∇vh)=(∇×f,vh)

Similarly, we have . Thus, taking and for the case we have

 (2.19) 12ddt∥wh∥2=0.

Thus we prove the conservative property of the 2D enstrophy.

For 3D flows, taking in eq. 2.15 we arrive at

 (2.20) 12ddt∥wh∥2−ch(wh,uh,wh)=0,

in the case where and . Note that the exact vorticity satisfies

 (2.21) 12ddt∥w∥2−c(w,u,w)=0.

There is a little difference between the continuous case eq. 2.21 and the discrete case eq. 2.20 formally.

#### 2.1.4. Vorticity

Note that under Assumption 2. Set in eq. 2.15 and we obtain

 (2.22) (∂wh∂t,ei)+ch(uh,wh,ei)−ch(wh,uh,ei)=0.

Equation 2.12 implies that

 (2.23) ch(uh,wh,ei)=−ch(uh,ei,wh)=0,

and

 (2.24) ch(wh,uh,ei)=−ch(wh,ei,uh)=0.

Substituting eq. 2.23 and eq. 2.24 into eq. 2.22 gives that

 (∂wh∂t,ei)=0,

which implies that the total vorticity is conserved.

### 2.2. Semi-discrete error analysis

Denote by the Stokes projection which satisfies that

 (2.25) (∇ΠShz,∇vh)=(∇z,∇vh)∀ z∈H10(Ω),vh∈V0h.
###### Assumption 3.

The Stokes projection satisfies the following estimate

 (2.26) ∥∥∇ΠShw∥∥p⩽C∥∇w∥p,1≤p≤∞.

To satisfy the estimate eq. 2.26, some extra conditions might be required. For example, the domain is convex. For the concrete analysis of eq. 2.26, we refer the readers to the paper [17]. Let

 e=u−uh=u−ΠShu+ΠShu−uh=η+ϕh,

and

 eπ=u−Πhuh=u−ΠhΠShu+Πh(ΠShu−uh)=ηπ+Πhϕh.

The following error equation will be used in error analysis.

 (2.27) (et,vh)+c(u,u,vh) −c(Πhuh,uh,vh) +ν(∇e,∇vh)−b(vh,p−qh)=0∀ vh∈V0h,qh∈Wh.

The above equation can be further rewritten as

 (2.28) ((ϕh)t,vh) +ν(∇ϕh,∇vh)=−(ηt,vh)−(c(u,u,vh)−c(Πhuh,uh,vh)) +b(vh,p−qh)=0∀ vh∈V0h,qh∈Wh,

where we use the fact that .

###### Theorem 2.3 (Semi-discrete estimate).

Let solve eq. 2.7 and solve eq. 2.1 with , and Under Assumption 1 and Assumption 3 it holds

 (2.29) ∥e(T)∥2 +ν∫T0∥∇e∥2 dt≤∥η(T)∥2+ν∥∇η∥2L2(0,T;L2)+K(u)(∥ϕh(0)∥2 +ν−1∥ηt∥2L2(0,T;H−1)+ν−1infqh∈L2(0,T;Wh)∥p−qh∥2L2(0,T;L2) +B(u)(∥∇η∥2L4(0,T;L2)+∥ηπ∥2L4(0,T;L2)))

with

 K(u)=exp(C(∥u∥L1(0,T;L∞)+∥∇u∥L1(0,T;L∞))),

and

 B(u)=C(∥u∥L2(0,T;L∞)+∥∇u∥L2(0,T;L∞)),

independent of

###### Proof.

Set in eq. 2.28 and we arrive at

 (2.30) 12ddt||ϕh||2+ ν||∇ϕh||2=−(ηt,ϕh)− (c(u,u,ϕh)−c(Πhuh,uh,ϕh))+b(ϕh,p−qh).

Applying the Cauchy-Schwarz inequality and the weighted Young’s inequality one obtains

 (2.31) |(ηt,ϕh)|≤ν−1∥ηt∥2H−1+ν4∥∇ϕh∥2,
 (2.32) |b(ϕh,p−qh)|≤Cν−1∥p−qh∥2+ν4∥∇ϕh∥2.

To analyze the nonlinear terms, here we use a similar splitting way as [40, 36]. We split the difference of the two nonlinear terms as

 (2.33) c(u,u,ϕh)−c(Πhuh,uh,ϕh)=c(u,η,ϕh)+c(u,ΠShu,ϕh)−c(Πhuh,uh,ϕh).

Further,

 (2.34) c(u,ΠShu,ϕh)−c(Πhuh,uh,ϕh) =c(eπ,ΠShu,ϕh)+c(Πhuh,ΠShu,ϕh) −c(Πhuh,uh,ϕh) =c(eπ,ΠShu,ϕh)+c(Πhuh,ϕh,ϕh) =c(eπ,ΠShu,ϕh) =c(ηπ,ΠShu,ϕh)+c(Πhϕh,ΠShu,ϕh).

Then the Cauchy-Schwarz inequality, Equation 2.26 and the Young’s inequality give that

 (2.35) |c(u,η,ϕh)|≤||u||∞||||∇η||||ϕh||≤C||u||∞(||∇η||2+||ϕh||2),
 (2.36) |c(ηπ,ΠShu,ϕh)|≤||ηπ||||∇ΠShu||∞||ϕh||≤C||∇u||∞(||ηπ||2+||ϕh||2),

and

 (2.37) |c(Πhϕh,ΠShu,ϕh)|≤||Πhϕh||||∇ΠShu||∞||ϕh||≤C||∇u||∞||ϕh||2,

where in the last inequality we also use Assumption 1. Substituting eqs. 2.37, 2.36 and 2.35 into eqs. 2.34 and 2.33 provides

 |c(u,u,ϕh)−c(Πhuh,uh,ϕh)|≤C(∥u∥∞∥∇η∥2+∥∇u∥∞∥ηπ∥2+[∥u∥∞+∥∇u∥∞]∥ϕh∥2),

which, together with eqs. 2.32 and 2.31, gives that

 (2.38)

Then integrating over and using the Gronwall inequality and Hölder inequality we finally obtain

 (2.39) ∥ϕh(T)∥2+ ν∫T0∥∇ϕh∥2 dt≤K(u)(∥ϕh(0)∥2+ν−1∥ηt∥2L2(0,T;H−1) +ν−1infqh∈L2(0,T;Wh)∥p−qh∥2L2(0,T;L2) +B(u)(∥∇η∥2L4(0,T;L2)+∥ηπ∥2L4(0,T;L2))).

Finally, the estimate eq. 2.29 follows immediately from eq. 2.39 and the triangle inequality. ∎

## 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

 (3.1) Given upreh,unh∈Vh and pnh∈Wh, find (un+1h,pn+1h)∈Vh×Wh such that (un+1h−unhΔt,vh)+νa(un+12h,vh) +ch(upreh,un+12h,vh) −b(vh,pn+12h) +b(un+1h,qh)=(fn+12,vh)

for all , where