The Optimal L 2 Error Estimate of Stabilized Finite Volume Method for the Stationary Navier-Stokes Problem

A finite volume method based on stabilized finite element for the two-dimensional stationary Navier-Stokes equations is analyzed. For the element, we obtain the optimal error estimates of the finite volume solution and . We also provide some numerical examples to confirm the efficiency of the FVM. Furthermore, the effect of initial value for iterative method is analyzed carefully.


Introduction
In paper 1 , G.He and Y.He introduce a finite volume method FVM based on the stabilized finite element method for solving the stationary Navier-Stokes problem and obtain the optimal H 1 error estimates for discretization velocity, however, to our dismay, without the optimal L 2 error estimate.It is inspiring that the following further numerical examples tell us that it has nearly second-order convergence rate.So, in this paper, we introduce a new technique to prove the optimal L 2 error of a generalized bilinear form and then gain the optimal L 2 error estimates of the stabilized finite volume method for the stationary Navier-Stokes problem.
For the convenience of analysis, we introduce the following useful notations.Let Ω be a bounded domain in R 2 assumed to have a Lipschitz continuous boundary ∂Ω and to satisfy a further smooth condition to ensure the weak solution's existence and regularity of

Advances in Numerical Analysis
Stokes problem.For more information, see the A1 assumption stated in 1, 2 .We consider the stationary Navier-Stokes equations −νΔu u • ∇ u ∇p f, div u 0, x ∈ Ω, u| ∂Ω 0, 1.1 where u u 1 x , u 2 x represents the velocity vector, p p x the pressure, f f x the prescribed body force, and ν > 0 the viscosity.
For the mathematical setting of problem 1.1 , we introduce the following Hilbert spaces: The spaces L 2 Ω m m 1, 2, 4 are endowed with the usual L 2 -scalar product •, • and norm • 0 , as appropriate.The space X is equipped with the scalar product ∇u, ∇v and norm ∇u 0 .
Define Au −Δu, which is the operator associated with the Navier-Stokes equations.It is positive self-adjoint operator from D A H 2 Ω 2 ∩ X onto Y , so, for α ∈ R , the power A α of A is well defined.In particular, D A 1/2 X, D A 0 Y , and We also introduce the following continuous bilinear forms a •, • and d •, • on X × X and X × M, respectively, by Under the above notations, the variational formulation of the problem 1.1 reads as follows: find u, p ∈ X, M such that for all v, q ∈ X, M : The following existence and uniqueness results are classical see 3 .
Theorem 1.1.Assume that ν and f ∈ Y satisfy the following uniqueness condition: where Then the problem 1.7 admits a unique solution u, p

FVM Based on Stabilized Finite Element Approximation
In this section, we consider the FVM for two-dimensional stationary incompressible Navier-Stokes equations 1.1 .Let h > 0 be a real positive parameter.The finite element subspace X h , M h of X, M is characterized by T h T h Ω , a partition of Ω into triangles, assumed to be regular in the usual sense see 4-7 .The mesh parameter h is given by h max{h K }, and the set of all interelement boundaries will be denoted by Γ h .Besides, we also use the configuration based on barycenter of element K i ∈ T h to construct a dual partition T * h of T h , which is shown in Figure 1.
Finite element subspaces of interest in this paper are defined as follows: the continuous piecewise linear velocity subspace the piecewise constant pressure subspace and the dual space of velocity subspace Actually, this choice of X * h is the span of the charaicteristic functions of the volume K * .Note that this mixed finite element method is unstable in the standard Babuška-Brezzi sense 8 .
Define the interpolation operator where N h {P i : Vertices of triangles in T h }.
where n is the out normal vector.We also define the trilinear forms b for all u h , v h , w h ∈ X h , the right side of term and a generalized bilinear form on Based on the dual partition and bilinear forms defined above, this paper still introduces the norms and seminorms 1, 9 : where with S v the area of Δv i v j v k see Figure 1 .
Formally, there are some differences between u h 0,h , A 1/2 h u h 0 and u h 0 , A 1/2 h u h 0 , respectively, but we, actually, have the following results 9-12 .

Lemma 2.1.
There exist constants c 1 , c 2 , independent of h, such that

2.11
So, for simplicity, we also denote u h 0,h and u h 1,h by u h 0 and u h 1 , respectively, without confusion.Below c with or without a subscript is a generic positive constant.
For the above finite element spaces X h and M h , it is well known that the following approximation properties and inverse inequality 12 hold see 4, 13 , where I h : D A → X h is the interpolation operator and J h :

Advances in Numerical Analysis
In order to define a locally stabilized formulation of the stationary Navier-Stokes problem, we also need a macroelement partition Λ h as follows: Given any subdivision T h , a macroelement partition Λ h may be defined such that each macroelement K is a connected set of adjoining elements from T h .Every element K must lie in exactly one macroelement.For each K, the set of interelement edges which are strictly in the interior of K will be denoted by Γ K , and the length of an edge e ∈ Γ K is denoted by h e .For a macroelement K the restricted pressure space is given by With the above choices of the velocity-pressure finite element spaces X h , X * h , M h and these additional definitions, a locally stabilized formulation of the Navier-Stokes problem 1.1 can be stated as follows. where p e q e ds,

2.15
for all p, q in the algebraic sum H 1 Ω M h , and • e is the jump operator across e ∈ Γ K and β > 0 is the local stabilization parameter.
A general framework for analyzing the locally stabilized formulation 2.14 can be developed using the notion of equivalence class of macroelements.As in Stenberg 14 , each equivalence class, denoted by E K , contains macroelements which are topologically equivalent to a reference macroelement K. See 1, 2 to get some examples.
The following stability results of these mixed methods for the macroelement partition defined above were formally established by Kay and Silvester 6 and Kechkar and Silvester 7 .Throughout the paper we will assume that β ≥ β 0 .Theorem 2.3.Given a stabilization parameter β ≥ β 0 > 0, suppose that every macroelement K ∈ Λ h belongs to one of the equivalence classes E K , and that the following macroelement connectivity condition is valid: for any two neighboring macroelements K 1 and K 2 with K 1 ∩K 2 ds / 0, there exists v ∈ X h such that Advances in Numerical Analysis 7 Then, for all p, q ∈ H 1 Ω M h , and where c > 0 is a constant independent of h and β, and β 0 is some fixed positive constant.

Error Estimates
In order to derive error estimates of u h , p h in the FVM, we need the existence and some regularities of the variational problem 2.14 see 1 .
Lemma 3.1.Under the assumptions of Theorem 2.3, there exist constants γ and α > 0 such that hold for all u h , p h and v h , q h ∈ X h , M h .
For the trilinear terms b u, v h , I * h w h and b v h , u, I * h w h , the following properties are useful 1, 2 .Set for any u h , v h , w h ∈ X h .Lemma 3.3.Suppose the assumptions of Theorem 2.3 and 3.4 hold, and the body force f satisfies the following uniqueness condition: Then there exists a unique solution (u h , p h ) of problem 2.14 satisfying the following estimate: In order to derive error estimates of the stabilized finite volume solution u h , p h , we need the following Galerkin projection R h , Q h : X, M → X h , M h defined by 3.9 for each v, q ∈ X, M .Note that, due to Lemma 3.1, R h , Q h is well defined.Now, we derive the following optimal L 2 error estimate of u h and p h defined in 3.9 .Using an argument similar to ones used by Layton and Tobiska in 15 , the following approximate properties can be obtained.

Lemma 3.4. Under the assumptions of Lemma 3.3, the projection R h , Q h satisfies
Proof.Equations 3.10 and 3.11 is the directly from 1 .Next, let v, q ∈ D A , H 1 Ω ∩M and introduce the dual Stokes problem: find Φ, Ψ ∈ X, M such that

3.13
Using the regularity assumption of Stokes problem See the A1 assumption in 1, 2 , there holds

3.14
Now, setting w v − R h v, q , r q − Q h v, q , using 2.18 and 3.9 , we obtain that for

3.15
For any v h ∈ X h , we have 9-11

3.16
Since the dual partition is formed by the barycenter, similar calculation in 10, 11, 16 allows us to have

Advances in Numerical Analysis
Since C h p, q h 0, for all p ∈ H 1 Ω ∩ M , for all q h ∈ M h , similarly in 7 , we have
Next, we will derive the following error estimates of the finite element solution u h , p h defined in Section 2.

3.21
Proof.It is well known that the weak solutions u, p ∈ D A ∩ V, H 1 Ω ∩ M .Hence, we derive from 2.14 and 3.9 that for all v h , q h ∈ X h , M h where e h R h u, p − u h and η h Q h u, p − p h .Taking v, q e h , η h in 3.22 , we obtain

3.25
Combining the above estimates with 3.24 and using the uniqueness condition 3.4 yield

3.26
Finally, one finds from 3.12 and 3.26 that

3.27
Hence, combining the above estimates with 3.8 gives 3.21 .

Numerical Example
For stationary Navier-Stokes problem, the iteration scheme, in general, is

4.1
The submatrices occurring in 4.1 correspond to differential operators as • .The right-hand side f contains information from the source information.
In general, this problem can be solved by the following Newton method: where R, r are the so-called nonlinear residual.Actually the difference between 4.2 and 4.1 is that, in computing the corrections v mid and p mid from R, r, the quadratic term N v mid v mid deduced from 4.1 is dropped and gives the linear problem 4.2 .

Numeric Example I
Consider a unit square domain with an exact solution given by

4.3
f is determined by 1.1 .After some computation using stretched gird, we have the following results.
Figure 2 is the relative error and convergence rate of velocity and pressure when ν 0.005, β 10.Table 1 lists the different errors and convergence rates of numerical velocity and pressure for the same ν and β.From the figure and table, we can see that the almost second-order L2 convergence is obtained, which confirms our theoretical prediction.
Figure 3 is the L 2 relative error of numerical velocity versus the number of iterate steps for different ν.The figure tells us the numerical velocities, in general, converge very fast.Moreover, the figure also tells the bigger the ν, the faster the convergence speed, which is consistent with the really case.If ν is not too small, for example, ν ≥ 0.01, only several Newton iterations are needed.

Numeric Example II
The second example is the classical lid-driven flow governed by stationary Navier-Stokes equations in a square cavity.We impose watertight boundary conditions, that is, u x 0, 1 u x 1, 1 0 and u x 1, for 0 < x < 1.From the streamlines in Figures 4-7, we can see there are some different performances for different ν for the problem based on the stretched grid 16, 17 with 128 × 128 grid .
The left subplot in Figure 4 is the velocity solution of Stokes problem which serves as the initial guess of Newton method; the right subplot in Figure 4 is the nonconvergence numeric velocity with that initial guess and 9 times Newton iterations.From Figure 4, we can see that there is a different performance from Numerical Example I for the lid-driven flow; if the ν is smaller, the initial value needed in Newton iteration has to be nearer the exact solution.The nonconvergence indicates that we must have a good initiate value for Newton iteration.For the ν 0.001, the usual Stokes initial value is not sufficient and a better initial value is needed, which can be computed by the following Picard method:   The difference between Picard method and Newton method is that the linear term N v mid v old is also dropped from 4.2 , and thus the Picard method commonly referred to as the Oseen system.
The left subplot in Figure 5 is the initial velocity for the Newton iteration based on two Picard iterations without Newton iteration, and the right subplot in Figure 5 is the streamlines of the convergence numeric velocity evaluated at the 2 Picard iterations, using 4 times Newton iterations.
Figures 6 and 7 give the behavior of different iteration results for ν 0.001 and ν 0.00033.From these figures, we can see that if ν is smaller, the initial value needed in Newton method should be better.The initial value computed by one or two steps Picard method is already insufficient and thus more Picard iterations are needed.In addition, we Iteration steps  can also see that the convergence speed of Picard is not as fast as Newton: if the initial value is sufficient for Newton iteration, the convergence speed of Newton's method is faster than Picard's method.Actually, the Picard method corresponds to a simple fixed point iteration strategy for solving 2.14 whose convection coefficient is evaluated at the current velocity.Thus, the rate of convergence of Picard method is only linear in general; whereas, for the added more linear  term, if the initial value is sufficient close to a nonsingular solution, the Newton method is locally convergence quadratic For more information see 18 .
It is necessary to pay attention to the "finest" number of Picard iteration in the computation of initial value for Newton iteration.Since the convergence radius of the Newton method is proportional to Reynolds number namely, 1/ν in general, in these computations, we roughly choose the times of Picard iteration to increase proportionately with Reynolds number.Many numerical tests show that this strategy is good enough for the success of the ensuing Newton iteration.

Conclusions
The main work in this paper is the demonstration of the optimal order in the L 2 error of the velocity and emphasis on some aspects of its associated numerical computation.Both the theoretical analysis and numerical results indicate the efficiency of the FVM for stationary Navier-Stokes equations.Further, numerical computations show the convergence of Newton method is closely related to the viscosity ν.Thus, as it is decreased, better and better initial values are needed, whereas the advantage of Picard method is that, relative to Newton method, it has a much large region of trust of convergence.As a result, a good choice is to combine the Newton method with Picard method in computing, and thus more complicated problems can be solved efficiently.

Figure 1 :
Figure 1: The partition and dual partition of a triangular.

Figure 2 :
Figure 2:The relative errors and convergence rates of velocity and pressure.

Figure 3 :
Figure 3: L 2 relative error of numerical velocity versus iteration steps for different ν.