Convergence results for the vector penalty-projection and two-step artificial compressibility methods

In this paper, we propose and analyze a new artificial compressibility splitting method which is issued from the recent vector penalty-projection method for the numerical solution of unsteady incompressible viscous flows introduced in [1], [2] and [3]. This method may be viewed as an hybrid two-step prediction-correction method combining an artificial compressibility method and an augmented Lagrangian method without inner iteration. The perturbed system can be viewed as a new approximation to the incompressible Navier-Stokes equations. In the main result, we establish the convergence of solutions to the weak solutions of the Navier-Stokes equations when the penalty parameter tends to zero.

1. Introduction and setting of the problem. The artificial compressibility method was introduced by Chorin [6] and Temam [17] for the solution of the unsteady incompressible Stokes or Navier-Stokes equations; see also [20] for the theoretical analysis. Then, some other numerical schemes to efficiently compute the solutions of Navier-Stokes problems can be viewed as discretizations of perturbed systems of the type of penalization [14] or pseudo-compressibility. This is the case with the famous projection methods from Chorin [7] and Temam [18,19] and their variants [10], see e.g. [15].
Here, we present a new approximation method for the Navier-Stokes equations modeling incompressible viscous flows in a bounded regular open set Ω endowed with Dirichlet boundary conditions on Γ = ∂Ω (Lipschitz-continuous). With a given source term f , the Navier-Stokes system reads: where R e denotes the Reynolds number.
According to the identity, −∆ϕ = curl curl ϕ−∇ div ϕ, we consider the following approximate method to obtain a solution of the above Navier-Stokes system, with the parameters r ≥ 0, γ > 0 and, ε > 0 We associate to the previous system the following boundary conditions and initial data, v ε (0) = v 0 , v ε (0) = 0, p ε (0) = p 0 , v ε · ν where ν denotes the outward unit normal vector on Γ. To vanish, at the limit process, the two tangential component of the velocity fields, v ε ∧ ν and v ε ∧ ν, we use a penalization method which will be detailed below.
This method is close to the artificial compressibility method of Chorin [6] and Temam [17], but presents one important difference. It is a two-step splitting method. The first equation of the previous system gives a predicted velocity v ε and the second one is the approximate projection of v ε on the free-divergence vector fields. This equation may be seen as an approximate method to solve the well-posed problem (see appendix A) : div v ε = −div v ε , Remark 1. This approximate method is issued from the Vector Penalty-Projection (VPP r,ε ) methods for the numerical solution of unsteady incompressible viscous flows introduced in [1] and [3]. A fast version of these methods, the so-called (VPP ε ) method, is recently proposed also for the numerical solution of the nonhomogeneous Navier-Stokes equations in [2,4]. It is shown to be very efficient to compute multiphase flows, i.e. fast, cheap, and robust whatever the density, viscosity or permeability jumps.
Even for r = 0, the resulting method which corresponds to a two-step pseudocompressibility method, is different from the original artificial compressibility method of Chorin [6] and Temam [17,20].
The new important point is the penalty term 1 ε (∇div v ε + ∇div v ε ) that appears in the velocity correction step which allows us a direct estimate on the divergence of the velocity. Moreover, this system is quite easy to solve and presents good stability properties, see [1,2,3]. The velocity v ε and the pressure p ε satisfy the equations: The vanishing of the tangential component at the limit process, is fullfilled by a penalization method, which implies that this boundary condition is satisfied at the order ε for the approximate solution.
Proposition 1. Under the previous hypothesis, one has the following properties: Moreover, there exists one constant C > 0 depending only on Ω such that: (1) Besides, if we suppose that the open set Ω is simply-connected, there exist two constants λ 0 and λ 1 depending only on Ω such that: (Ω) and we have: Proof. All these results may be found in [9] and [8].
For a Banach space E we introduce the Nikolskii space defined for 1 ≤ q < +∞, 0 < σ < 1: endowed with the following norm: Let us recall the following property see for example [5] page 105,

Proposition 2. Let H be an Hilbert space and f a function given in
where M σ is a constant depending only on σ.
2. Main result. We associate to the previous approximate system, the variational problem where the tangential components of the velocities v ε and v ε are penalized. This problem is studied in the next section.
Then the velocity v ε and the pressure p ε satisfy in Remark 2. In order to establish the strong convergence of the sequence (v ε ) ε>0 when ε → 0, we use in Section 4 the Leray's orthogonal decomposition in the bounded domain. The curl-free component vanishes with the penalty term introduced by our method, whereas the divergence-free component strongly converges thanks to an estimate of a fractional derivative in time, see [20]. However, this requires to consider velocity fields having only their normal component which is zero on the boundary. Since at the limit process, we aim at solving the Navier-Stokes problem with homogeneous Dirichlet boundary condition, we also penalize the tangential part of the velocity fields.
We prove in section 3 the following results.
Lemma 2.1. Let us suppose that f belongs to L 2 (]0, T [; L 2 (Ω)). Then, there exists at least a solution to the system (2). This solution is unique in two space dimension. For the dimension d ≤ 3, this solution satisfies the following energy inequality: For the dimension d = 2, one has the following energy equality: This result is quite classical and we only give the sketch of proof in the section 3.
In fact, we can precise the previous energy inequality if we suppose that the data f belongs to L ∞ (R + ; L 2 (Ω)). This shows the absolute stability of the approximate method.
then, there exists a constant α independent of the data, such that The main goal of this paper is to prove the following convergence theorem: (3) that converges to a solution (v, p) to the system of Navier-Stokes equations with homogeneous Dirichlet boundary conditions.
. We now give an interpretation of the pressure and precise its convergence. Let us define The scalar function q ε appears to be the effective approximate pressure, and we have These convergence results for both velocity and pressure are proved in Section 4.
3. Energy estimates. We first establish the following existence result.
If d = 2, the unique solution of (2) satisfies the following regularity results: Remark 3. In the three-dimensional case, the equalities valid in the trace sense.

PHILIPPE ANGOT AND PIERRE FABRIE
Proof. For fixed parameters ε > 0, r ≥ 0 and γ > 0, we build approximate solutions by a classical Galerkin process.
Let us introduce the self-adjoint operator A = curl curl − ∇ div defined on the domain H 1 ν (Ω) ∩ (H 2 (Ω)) d . Then, for the approximation of the two fields of velocity, we use as special basis the eigenfunctions of this operator associated with the following boundary conditions: For the pressure, one can use as special basis the eigenfunctions of the self-adjoint operator A = −∆ with domain H 2 (Ω) associated to the Neumann boundary conditions. This approximate finite dimensional system is then a classical ordinary differential equation which has a unique solution. Next, to perform the limit we use the same strategy as for the classical Navier-Stokes equations i.e. a priori estimates and compactness results using an estimate on the temporal derivative, see for example [13], [20], [5].
Now, we will focus our attention on the estimates on the time derivative according to the dimension d. Let us begin with the three-dimensional case. We have to estimate the two nonlinear terms v ε · ∇w, w div v ε , with either w = v ε or w = v ε and the pressure p ε .
Suppose first that d = 3. By Sobolev embedding, the two nonlinear terms of the form v ε · ∇w and w div v ε belong to L : The bounds of the linear terms are straightforward and we have These estimates show that the velocities ( v ε , v ε ) are equal almost everywhere to continuous functions with values in (H 1 ν (Ω)) ′ . Besides, the pressure p ε is equal almost everywhere to a continuous function with value in L 2 0 (Ω). For the two-dimensional case, the situation is quite different. We observe first that the velocity fields v ε and v ε belong to L 4 (]0, T [; L 4 (Ω)), so that: So it follows from the equation (2) that We now observe that the two velocity fields v ε and v ε belong to Thus the functions ( v ε , v ε ) are equal almost everywhere to continuous functions with values in L 2 (Ω).
This ends the proof of proposition 3.

3.1.
Stability. In the case of three-dimensional vector spaces, we do not have an equality for the conservation of the energy, we have only an inequality. Nevertheless for two-dimensional vector spaces, the weak solutions satisfy the energy equality.
Proof. Through classical computations one obtains, with equations (2) and (3): Multiplying (4) by r ε and (5) by ε and summing with (6) and (7), one obtains: Let us now give some bounds of the right-hand side terms. The According to the estimate of the L 2 norm in H 1 ν (Ω) given by equation (1), we bound the source terms in the following way: Using these bounds, we get from the previous equation the following fundamental estimate: After integration in time, we deduce from the previous estimate that there exists a continuous function g defined on [0, T ] such that: for all t > 0, with This inequality is the key point to establish the convergence result.
To improve the convergence result in the two-dimensional case, we use the following energy equality derived as above without using (5).
This concludes the proof of lemma 2.1.

3.2.
Uniform stability for the approximate solution. In this section, we deal with the stability of the proposed approximation method. We notice that this property is valid for all solutions satisfying the energy inequality (8), as it is the case when they are built by a finite dimensional approximation method such as the Galerkin method for example.
4. Convergence analysis and compactness results.

4.1.
Compactness results for the velocity. Let us introduce the Leray projection w ε of a velocity field v ε (t) ∈ H 1 ν (Ω) defined as follows By the estimate (9), we see that the irrotational part of v ε goes to zero with ε. Thus it remains to bring to the fore the behavior of the free divergence part w ε and to obtain an estimate on a fractional time derivative of this term. We detail the different steps of this strategy.
From the regularity of the Leray projector (see R. Temam [20] page 18), one has: Moreover, we can easily prove the following lemma.
Lemma 4.1. There exists two constants depending only on T and Ω such that: Proof: The function q ε belongs to H 2 (Ω) and satisfies This implies using the estimate (9) Besides, we have so that, with Poincaré-Neumann inequality, we get The regularity properties of the Neumann problem give Moreover, by orthogonality of the Leray projector in L 2 , one has This concludes the proof of the lemma 4.1.
So by interpolation and using estimates (13)- (14), we have proved the result below.
Corollary 1. The function q ε satisfies: ∇q ε strongly converges to 0 in L p (]0, T [; L 2 (Ω)) , ∀p, 1 ≤ p < +∞. Now we have to write the equation satisfied by w ε . As the Leray projection is orthogonal in L 2 (Ω), this equation reads Now we introduce the extension by 0 of w ε (resp. v ε ) outside [0, T ] denoted, only in this part, by w ε (resp. v ε ) and we take the Fourier transform in time of the equation (15) Following Boyer-Fabrie [5, page 253], we take ϕ = F( w ε )(τ ) as test function in the previous equation to obtain for all τ ∈ R: As we look for an estimate independent of ε, we have to pay a special attention to the imaginary part of the penalty term: By writing w ε = v ε − ∇q ε , we have: So, the imaginary part of A ε is bounded as follows: From estimates (12) and (9), we have So there exists a function f 4 ε (τ ) ∈ L 1 (R) bounded independently on ε such that: Now we can derive the estimate of |τ | Ω |F( w ε )(τ )| 2 dω, and we have We now estimate each term of the right-hand side of the previous inequality for d ≤ 3.
We also have the inequality: which implies, according to (9), that v ε · ∇v ε is bounded in L (Ω)). So, by Hausdorff-Young theorem, the family of functions F( v ε · ∇) v ε is bounded in L 6 (R; L 6 5 (Ω)). Then with Hölder inequality, The same arguments show that According to the regularity of the Leray projection recalled above and estimate (9), one has: As we have seen by the estimate (18), By hypothesis,f is a given function in L 2 (R; L 2 (Ω)) and from estimate (9), we get that v ε is bounded in L 2 (R; L 2 (Ω)), so as w ε is the Leray projection of v ε , the function w ε is also bounded in L 2 (R; L 2 (Ω)). Finally, we obtain These two terms come from the Dirac measure when we derive discontinuous functions. Let us consider f 7 ε .
We treat in the same way the term f 6 ε , and thus: We are now able to show that the set of functions ( w ε ) ε is bounded in an appropriate Nikolskii space.
For all γ < 1, there exists a constant d such that: Let us denote f 0 ε (τ ) = F( w ε )(τ ) 2 L 2 , which belongs to L 1 (R), the previous inequality reads with (19): If we suppose that the function τ → 1 1 + |τ | γ belongs to L ∞ (R) ∩ L 2 (R), then the function τ → h(τ ) belongs to L 1 (R). This condition is satisfied for γ ∈] 2 3 , 1[. So we have proved : Then, from lemma 4.1 and 4.2 we deduce the following key result: There exists a sequence (ε k ) k which converges to zero and a function v ∈ L 2 (]0, T [; L 2 (Ω)) satisfying div v = 0 such that: Proof. The function v ε is the sum of two terms ∇q ε and w ε . From corollary 1 the first term converges strongly to 0 in L 2 (]0, T [; L 2 (Ω)). Now, from Aubin-Lions-Simon Theorem, it follows from lemma 4.2, that there exists a sequence (ε k ) k such that: Moreover, since div w ε k = 0, we have div v = 0.

4.2.
Convergence of the method. We first give a general convergence theorem for a subsequence solution of the approximate scheme (3), to a weak solution of the initial Navier-Stokes problem, in the case d ≤ 3. For the two-dimensional case, since the weak solution of the Navier-Stokes equation is unique, the whole sequence of approximate solution v ε converges to v. Moreover, in this case, we prove that the convergence is strong.
According to estimate (9), there exists a sequence ε k such that 2 is a norm equivalent to the H 1 -norm on H 1 ν (Ω), we have ∇v ε k ⇀ ∇v in L 2 (0, T ; L 2 (Ω)) weakly. and so, we can take the limit on the term From estimate (9), we have for the tangential traces and since for any function v in H 1 (Ω) d , This implies (v ∧ ν) |Γ = 0, and so, since by construction (v · ν) |Γ = 0, v belongs to Finally, at the limit process we obtain From the identity which is valid for all functions (v, ϕ) ∈ (H 1 0 (Ω)) d × (H 1 0 (Ω)) d , the previous equality shows that the limit function v satisfies the classical Navier-Stokes equations in a weak sense.

4.2.2.
The special case d = 2. The key point to establish the strong convergence in the two-dimensional case, lies on an idea of R. Temam [20]. It is based on the fact that in this case, the solution of the approximate problem and the solution of the Navier-Stokes equation verify the equality of energy. The idea is to bring to the fore an energy equation satisfied by the difference between the approximate solution and the exact solution.
We first observe that, according to the equality −∆ = curl curl − ∇ div the classical weak solution v of the Navier-Stokes equation satisfies for any test function in (H 1 0 (Ω)) d , with free-divergence: The equation satisfied by the error v ε −v = u ε with a free-divergence test function After integration in time, this equation gives curl u ε · curl ϕ dω dτ = Ω u ε (0) · ϕ dω.
Taking the limit when ε goes to 0, one obtains with the convergences properties stated in the previous section: Following R. Temam [20], we introduce Here v is the unique solution of the two-dimensional Navier-Stokes equations. Due to the energy equality (10), χ ε (t) satisfies By weak convergence in L 2 (]0, t[; H 1 0 (Ω)) of the sequence (v ε ) ε and from the lemma 4.4, we observe that : In the two-dimensional case, the unique solution of the Navier-Stokes equation satisfies the following energy equality Moreover, from estimate (10) We are now able to precise the convergence for the effective pressure and establish the theorem 2.4 Let us write the equation satisfied by the velocity v ε and the pressure p ε in the distribution sense. We have ∂v ε ∂t + (v ε · ∇)v ε + 1 2 (div v ε )v ε + 1 R e curl curl v ε + ∇ p ε − 1 + ε εR e + r div v ε = f.
Introducing the effective pressure q ε = p ε − 1 + ε εR e + r div v ε , this equation reads and the proof follows, from the previous steps.

Remark 4.
1. The function u belongs to H 1 ν (Ω), which is not the case for v ε or u 1 , without some additional regularity hypotheses on the function g. Nevertheless, the function u ε = v ε − u − u 1 belongs to H 1 ν (Ω). 2. In the case where g = 0, the function v ε belongs to H 1 ν (Ω).