# The OptimalL2 Error Estimate of Stabilized Finite Volume Method for the Stationary Navier-Stokes Problem

Authors: Guoliang He; Jian Su; Wenqiang Dai
Journal: Advances in Numerical Analysis
 (2012)

---

## Abstract

A finite volume method based on stabilized finite element for the two-dimensional
stationary Navier-Stokes equations is analyzed. For theP1–P0
element, we obtain the optimal L2 error estimates of the finite volume solution uh and ph. 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.


---

## Body

## 1. 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 H1 error estimates for discretization velocity, however, to our dismay, without the optimal L2 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 L2 error of a generalized bilinear form and then gain the optimal L2 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 ℝ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 Stokes problem. (For more information, see the A1 assumption stated in [1, 2].) We consider the stationary Navier-Stokes equations
(1.1)-νΔu+(u·∇)u+∇p=f,divu=0,x∈Ω,u∣∂Ω=0,
where u=(u1(x),u2(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:
(1.2)X=(H01(Ω))2,Y=(L2(Ω))2,M={q∈L2(Ω):∫Ω‍qdx=0},H={v∈L2(Ω)2;divv=0inΩ,v·n∣∂Ω=0}.
The spaces (L2(Ω))m (m=1,2,4) are endowed with the usual L2-scalar product (·,·) and norm ∥·∥0, as appropriate. The space X is equipped with the scalar product (∇u,∇v) and norm ∥∇u∥0.DefineAu=-Δu, which is the operator associated with the Navier-Stokes equations. It is positive self-adjoint operator from D(A)=(H2(Ω))2∩X onto Y, so, for α∈ℝ, the power Aα of A is well defined. In particular, D(A1/2)=X, D(A0)=Y, and
(1.3)(A1/2u,A1/2v)=(∇u,∇v),∥A1/2u∥0=(∇u,∇u)1/2,
for all u,v∈X.We also introduce the following continuous bilinear formsa(·,·) and d(·,·) on X×X and X×M, respectively, by
(1.4)a(u,v)=ν((u,v))∀u,v∈X,d(v,q)=-(v,∇q)=(q,divv)∀v∈X,q∈M,
a generalized bilinear form on (X,M)×(X,M) by
(1.5)ℬ((u,p);(v,q))=a(u,v)-d(v,p)+d(u,q),
and a trilinear form on X×X×X by
(1.6)b(u,v,w)=((u·∇)v,w)+12((divu)v,w).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):
(1.7)ℬ((u,p);(v,q))+b(u,u,v)=(f,v).
The following existence and uniqueness results are classical (see [3]).Theorem 1.1.
Assume thatν and f∈Y satisfy the following uniqueness condition:
(1.8)1-N1ν2∥f∥-1>0,
where
(1.9)N1=supu,v,w∈H01(Ω)b(u,v,w)∥A1/2u∥0∥A1/2v∥0∥A1/2w∥0.
Then the problem (1.7) admits a unique solution (u,p)∈(D(A)∩X,H1(Ω)∩M) such that
(1.10)∥A1/2u∥1≤ν-1∥f∥-1,∥u∥2+∥p∥1≤c∥f∥0.


## 2. 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 (Xh,Mh) of (X,M) is characterized by Th=Th(Ω), 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{hK}, and the set of all interelement boundaries will be denoted by Γh. Besides, we also use the configuration based on barycenter of element Ki∈Th to construct a dual partition Th* of Th, which is shown in Figure 1.Figure 1
The partition and dual partition of a triangular.Finite element subspaces of interest in this paper are defined as follows: the continuous piecewise linear velocity subspace(2.1)Xh={v∈X:v∣K∈(P1(K))2,∀K∈Th},
the piecewise constant pressure subspace
(2.2)Mh={q∈M:q∣K∈P0(K),∀K∈Th},
and the dual space of velocity subspace Xh*(2.3)Xh*={v∈(L2(Ω))2:v∣K*∈(P0(K*))2,∀K*∈Th*}.
Actually, this choice of Xh* 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 operatorIh*:Xh→Xh*,
(2.4)Ih*uh=∑xi∈Nh‍uh(xi)χi(x),
where Nh={Pi:VerticesoftrianglesinTh}.Let us introduce the continuous bilinear formsa~(·,·), d~(·,·), and d(·,·) on Xh×Xh, Xh×Mh and Xh×Mh as follows:
(2.5)a~(uh,Ih*vh)=ν((uh,Ih*vh))=-ν∑Ki*∈Th*‍∫∂Ki*‍vh(xi)∂uh∂nds,∀uh,vh∈Xh,d~(Ih*vh,ph)=(Ih*vh,∇ph)=∑Ki*∈Th*‍∫∂Ki*‍vh(xi)ph·nds,∀uh∈Xh,ph∈Mh,d(uh,qh)=-(uh,∇qh)=(qh,divuh),∀uh∈Xh,qh∈Mh,
where n is the out normal vector. We also define the trilinear forms b~(·,·,·) on Xh×Xh×Xh by
(2.6)b~(uh,vh,Ih*wh)=((uh·∇)vh,Ih*wh),
for all uh,vh,wh∈Xh, the right side of term
(2.7)(f,Ih*vh)=∑Ki*∈Th*‍∫Ki*‍vh(xi)fdx,∀vh∈Xh,
and a generalized bilinear form on
(2.8)ℬ~((uh,ph);(Ih*vh,qh))=a~(uh,Ih*vh)-d~(Ih*vh,ph)+d(uh,qh).Based on the dual partition and bilinear forms defined above, this paper still introduces the norms and seminorms [1, 9]:
(2.9)∥uh∥0,h=(∑K∈Th‍∥uh∥0,h,K2)1/2,∥A~h1/2uh∥0=(∑K∈Th‍∥A~h1/2uh∥0,h,K2)1/2,∥uh∥1,h=(∥uh∥0,h2+∥A~h1/2uh∥02)1/2,
where
(2.10)∥uh∥0,h,K=[Sv3(uPi2+uPj2+uPk2)]1/2,∥A~h1/2uh∥0,h,K={[(∂uh(p)∂x)2+(∂uh(p)∂y)2]Sv}1/2,
with Sv the area of Δvivjvk (see Figure 1).Formally, there are some differences between∥uh∥0,h,∥A~h1/2uh∥0 and ∥uh∥0,∥Ah1/2uh∥0, respectively, but we, actually, have the following results [9–12].Lemma 2.1.
There exist constantsc1,c2, independent of h, such that
(2.11)c1∥uh∥0,h≤∥uh∥0≤c2∥uh∥0,h,∀uh∈Xh,∥A~h1/2uh∥0=∥Ah1/2uh∥0,∀uh∈Xh.So, for simplicity, we also denote∥uh∥0,h and ∥uh∥1,h by ∥uh∥0 and ∥uh∥1, respectively, without confusion. Below c (with or without a subscript) is a generic positive constant.For the above finite element spacesXh and Mh, it is well known that the following approximation properties and inverse inequality
(2.12)∥Ah1/2v∥0≤ch-1∥v∥0,∀v∈Xh,∥v-Ihv∥0+h∥Ah1/2(v-Ihv)∥0≤ch2∥Ahv∥0,∀v∈D(A),∥v-Ih*v∥0≤ch∥Ah1/2v∥0,∀v∈Xh,∥q-Jhq∥0≤ch∥q∥1,∀q∈H1(Ω)∩M
hold (see [4, 13]), where Ih:D(A)→Xh is the interpolation operator and Jh:H1(Ω)∩M→Mh is the L2-orthogonal projection.In order to define a locally stabilized formulation of the stationary Navier-Stokes problem, we also needa macroelement partitionΛh as follows: Given any subdivision Th, a macroelement partition Λh may be defined such that each macroelement 𝒦 is a connected set of adjoining elements from Th. Every element K must lie in exactly one macroelement. For each 𝒦, the set of interelement edges which are strictly in the interior of 𝒦 will be denoted by Γ𝒦, and the length of an edge e∈Γ𝒦 is denoted by he. For a macroelement 𝒦 the restricted pressure space is given by
(2.13)M0,h={q∈L02(𝒦):q∣K∈P0(K),∀K∈𝒦}.With the above choices of the velocity-pressure finite element spacesXh, Xh*, Mh and these additional definitions, a locally stabilized formulation of the Navier-Stokes problem (1.1) can be stated as follows.Definition 2.2 (locally stabilized FVM formulation).
Find(uh,ph)∈(Xh,Mh), such that for all (v,q)∈(Xh,Mh)(2.14)ℬ~h((uh,ph);(Ih*v,q))+b~(uh,uh,Ih*vh)=(f,Ih*v),
where
(2.15)ℬ~h((uh,ph);(Ih*v,q))=ℬ~((uh,ph);(Ih*v,q))+β𝒞h(ph,q),∀(uh,ph),(v,q)∈(Xh,Mh),𝒞h(p,q)=∑𝒦∈Λh‍∑e∈Γ𝒦‍he∫e[p]e[q]eds,
for all p,q in the algebraic sum H1(Ω)+Mh, and [·]e is the jump operator across e∈Γ𝒦 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 ℰ𝒦^, contains macroelements which are topologically equivalent to a reference macroelement 𝒦^. 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 𝒦∈Λh belongs to one of the equivalence classes ℰ𝒦^, and that the following macroelement connectivity condition is valid: for any two neighboring macroelements 𝒦1 and 𝒦2 with ∫𝒦1∩𝒦2‍ds≠0, there exists v∈Xh such that
(2.16)suppv⊂𝒦1∪𝒦2,∫𝒦1∩𝒦2‍v·nds≠0.
Then,
(2.17)∥𝒞h(p,q)∥0≤c∑K∈Th‍(∫K‍(∥p∥02+h2∥∇p∥02)dx)1/2(∫K‍(∥q∥02+h2∥∇q∥02)dx)1/2,
for all p,q∈H1(Ω)+Mh, and
(2.18)𝒞h(p,qh)=0,𝒞h(ph,q)=0,𝒞h(p,q)=0∀p,q∈H1(Ω),ph,qh∈Mh,
where c>0 is a constant independent of h and β, and β0 is some fixed positive constant.


## 3. Error Estimates

In order to derive error estimates of(uh,ph) 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 Theorem2.3, there exist constants γ and α>0 such that
(3.1)ν∥A~h1/2uh∥02+β𝒞h(ph,ph)=ℬ~h((uh,ph);(Ih*uh,ph)),(3.2)|ℬ~h((uh,ph);(Ih*vh,qh))|≤γ(∥A~h1/2Auh∥0+∥ph∥0)(∥A~h1/2Avh∥0+∥qh∥0),(3.3)α(∥A~h1/2uh∥0+∥ph∥0)≤sup(vh,qh)∈(Xh,Mh)ℬ~h((uh,ph);(Ih*vh,qh))∥A~h1/2vh∥0+∥qh∥0
hold for all (uh,ph) and (vh,qh)∈(Xh,Mh).For the trilinear termsb~(u,vh,Ih*wh) and b~(vh,u,Ih*wh), the following properties are useful [1, 2]. Set
(3.4)N2=supu,vh,wh∈H01(Ω)b~(u,vh,Ih*w)∥A1/2u∥0∥A~h1/2v∥0∥A~h1/2w∥0,N3=supu,vh,wh∈H01(Ω)b~(vh,u,Ih*w)∥A~h1/2v∥0∥A1/2u∥0∥A~h1/2w∥0,N=max{N1,N2,N3}.Lemma 3.2.
The trilinear formb~ satisfies
(3.5)|b~(uh,vh,Ih*wh)|≤c∥A~huh∥01/2∥A~h1/2uh∥01/2∥A~h1/2vh∥0∥A~hwh∥01/2∥A~h1/2wh∥01/2,
for any uh,vh,wh∈Xh.Lemma 3.3.
Suppose the assumptions of Theorem2.3 and (3.4) hold, and the body force f satisfies the following uniqueness condition:
(3.6)1-4Nν2∥f∥-1>0.
Then there exists a unique solution (uh,ph) of problem (2.14) satisfying the following estimate:
(3.7)ν∥A~h1/2uh∥02+∥ph∥02≤κ,(3.8)∥A~h1/2(u-uh)∥0+∥p-ph∥0≤κh.In order to derive error estimates of the stabilized finite volume solution(uh,ph), we need the following Galerkin projection (R~h,Q~h):(X,M)→(Xh,Mh) defined by
(3.9)ℬ~h((R~h(v,q),Q~h(v,q));(Ih*vh,qh))=ℬ~((v,q);(Ih*vh,qh))∀(vh,qh)∈(Xh,Mh),
for each (v,q)∈(X,M).Note that, due to Lemma3.1, (R~h,Q~h) is well defined. Now, we derive the following optimal L2 error estimate of uh and ph 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 Lemma3.3, the projection (Rh,Qh) satisfies
(3.10)∥A~h1/2(v-R~h(v,q))∥0+∥q-Q~h(v,q)∥0≤c(∥A1/2v∥0+∥q∥0),
for all (v,q)∈(X,M),
(3.11)∥A~h1/2(v-R~h(v,q))∥0+∥q-Q~h(v,q)∥0≤ch(∥Av∥0+∥A1/2q∥0),
for all (v,q)∈(D(A),H1(Ω)∩M), and
(3.12)∥v-R~h(v,q)∥0+h∥A~h1/2(v-R~h(v,q))∥+h∥q-Q~h(v,q)∥0≤ch2(|Av|0+∥q∥1),
for all (v,q)∈(D(A),H1(Ω)∩M).Proof.
Equations (3.10) and (3.11) is the directly from [1]. Next, let (v,q)∈(D(A),H1(Ω)∩M) and introduce the dual Stokes problem: find (Φ,Ψ)∈(X,M) such that
(3.13)ℬ~((w,r);(Φ,Ψ))=(w,Ih*(v-R~h(v,q))),∀(w,r)∈(X,M).
Using the regularity assumption of Stokes problem (See the A1 assumption in [1, 2]), there holds
(3.14)∥Φ∥2+∥Ψ∥1≤c∥v-R~h(v,q)∥0.
Now, setting w=v-R~h(v,q),r=q-Q~h(v,q), using (2.18) and (3.9), we obtain that for (Φh,Ψh)=(Ih*Φ,JhΨ)∈(Xh,Mh),
(3.15)∥v-R~h(v,q)∥02=ℬ~((v-R~h(v,q),q-Q~h(v,q));(Φ,Ψ))=ℬ~h((v-R~h(v,q),q-Q~h(v,q));(Φ,Ψ))=ℬ~h((v-R~h(v,q),q-Q~h(v,q));(Φ-Φh,Ψ-Ψh))=a~(v-R~h(v,q),Ih*(Φ-Φh))-d~(Ih*(Φ-Φh),q-Q~h(v,q))+d(v-R~h(v,q),Ψ-Ψh).
For any vh∈Xh, we have [9–11]
(3.16)a~(v-R~h(v,q),Ih*(Φ-Φh))=(Δ((v-R~h(v,q))-Ih(v-R~h(v,q))),Ih*(Φ-Φh)),∥d~(Ih*(Φ-Φh),q-Q~h(v,q))∥0≤c∥d((Φ-Φh),q-Q~h(v,q))∥0.
Since the dual partition is formed by the barycenter, similar calculation in [10, 11, 16] allows us to have
(3.17)∥a~(v-R~h(v,q),Ih*(Φ-Φh))∥0≤κh2∥v-R~h(v,q)∥0×(∥A~h1/2uh∥0+∥u∥2+∥f∥1+∥A~h1/2(v-R~h(v,q))∥0),(3.18)∥d~(Ih*(Φ-Φh),q-Q~h(v,q))∥0≤c∥d((Φ-Φh),q-Q~h(v,q))∥0≤c∥∇(Φ-Φh)∥0∥q-Q~h(v,q)∥0.
By (3.10), (3.11), and (3.14), we have
(3.19)∥d~(Ih*(Φ-Φh),q-Q~h(v,q))∥0≤ch2∥AΦ∥0(∥Av∥0+∥A1/2q∥0)≤ch2∥v-R~h(v,q)∥0(∥Av∥0+∥A1/2q∥0).
Since Ch(p,qh)=0, forallp∈(H1(Ω)∩M), forallqh∈Mh, similarly in [7], we have
(3.20)∥(div(v-R~h(v,q)),Ψ-Ψh)∥0≤ch∥∇(v-R~h(v,q))∥0∥Ψ∥1≤κh2(∥u∥2+∥p∥1)∥(v-R~h(v,q))∥0.
Finally, combining (3.10) with (3.17), (3.19), and (3.20) yields (3.12).Next, we will derive the following error estimates of the finite element solution(uh,ph) defined in Section 2.Theorem 3.5.
Assume that the assumptions of Lemma3.3 hold. Then the stabilized finite element solution (uh,ph) satisfies the error estimates
(3.21)∥u-uh∥0+h(∥A~h1/2(u-uh)∥0+∥p-ph∥0)≤ch2.Proof.
It is well known that the weak solutions(u,p)∈(D(A)∩V,H1(Ω)∩M). Hence, we derive from (2.14) and (3.9) that for all (vh,qh)∈(Xh,Mh)(3.22)ℬh((eh,ηh);(Ih*vh,qh))+b~(u-uh,u,Ih*vh)+b~(uh,u-uh,Ih*vh)=0,
where eh=R~h(u,p)-uh and ηh=Q~h(u,p)-ph. Taking (v,q)=(eh,ηh) in (3.22), we obtain
(3.23)ν∥A~h1/2eh∥2+β0𝒞h(ηh,ηh)+b~(eh,u,Ih*eh)+b~(u,eh,Ih*eh)≤|b(u-R~h(u,p),u,eh)|+|b(uh,u-R~h(u,p),eh)|.
We find from (1.10), (3.1), (3.4), and (3.6) that
(3.24)ν∥A~h1/2eh∥02-|b~(eh,u,Ih*eh)|-|b~(u,eh,Ih*eh)|≥ν∥A~h1/2eh∥02-2N∥Ah1/2u∥0∥A~h1/2eh∥02≥ν(1-2∥f∥0ν-2)∥A~h1/2eh∥02.
Moreover, By (3.5), (3.11), and Poincaré's estimate, we have
(3.25)∥b~(uh,u-R~h(u,p),eh)∥0+∥b~(u-R~h(u,p),u,eh)∥0≤∥b~(u,u-R~h(u,p),eh)∥0+∥b~(u-R~h(u,p),u,eh)∥0+∥b~(u-R~h(u,p),u-R~h(u,p),eh)∥0+∥b~(eh,u-R~h(u,p),eh)∥0≤c1∥Au∥0∥A~h1/2(u-R~h(u,p))∥0∥A~h1/2eh∥0+c2(∥A~h1/2(u-R~h(u,p))∥0+∥A~h1/2eh∥0)∥A~h1/2(u-R~h(u,p))∥0∥A~h1/2eh∥0≤ch2∥A~h1/2eh∥0.
Combining the above estimates with (3.24) and using the uniqueness condition (3.4) yield
(3.26)∥A~h1/2eh∥0≤ch2.
Finally, one finds from (3.12) and (3.26) that
(3.27)∥u-uh∥0≤∥eh∥0+∥u-Rh(u,p)∥0≤c3∥A~h1/2eh∥0+ch2(∥Au∥0+∥p∥1)≤ch2.
Hence, combining the above estimates with (3.8) gives (3.21).


## 4. Numerical Example

For stationary Navier-Stokes problem, the iteration scheme, in general, is(4.1)νAv+N(v)v+Bp=f-BTv+βC=0.
The submatrices occurring in (4.1) correspond to differential operators as A~-diag(νΔ), N(v)~v·∇, B~∇, -BT~div, and C~𝒞h(·,·). The right-hand side f contains information from the source information.In general, this problem can be solved by the following Newton method:(4.2)(1)R=f-(νAvold+N(vold))vold-Bpold,r=-BTvold;(2)(νAvmid+N(vold))vmid+N(vmid)vold+Bpmid=R,BTvmid=r;(3)vnew=vold+vmid,pnew=pold+pmid,
where R,r are the so-called nonlinear residual. Actually the difference between (4.2) and (4.1) is that, in computing the corrections vmid and pmid from R,r, the quadratic term N(vmid)vmid deduced from (4.1) is dropped and gives the linear problem (4.2).


### 4.1.  Numeric Example I

Consider a unit square domain with an exact solution given by(4.3)u(x,y)=(u1(x,y),u2(x,y)),p(x,y)=10(2x-1)(2y-1),u1(x,y)=10x2(x-1)2y(y-1)(2y-1),u2(x,y)=-10x(x-1)(2x-1)y2(y-1)2.f is determined by (1.1). After some computation using stretched gird, we have the following results.Figure2 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.Table 1
Numerical results of the FVM (ν=0.005, β=10).
h
∥
A
~
h
1
/
2
(
u
-
u
h
)
∥
0
/
∥
A
~
h
1
/
2
u
∥
0
Con. rate
∥
u
-
u
h
∥
0
/
∥
u
∥
0
Con. rate
∥
p
-
p
h
∥
0
/
∥
p
∥
0
Con. rate
1/16
0.06316279
—
0.00391173
—
0.07127911
—
1/32
0.03204341
0.97904713
0.00099767
1.97116812
0.03572585
0.99651095
1/64
0.01627933
0.97698618
0.00025462
1.97021278
0.01790786
0.99637536
1/128
0.00827253
0.97664010
0.00006499
1.96996167
0.00897672
0.99633198Figure 2
The relative errors and convergence rates of velocity and pressure.Figure3 is the L2 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.Figure 3
L
2 relative error of numerical velocity versus iteration steps for different ν.


### 4.2.  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,ux(0,1)=ux(1,1)=0 and ux=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).Nonconvergence streamlines forν=0.001,β=10 ((a) Stokes, (b) N9).
(a)
(b)Streamlines forν=0.001,β=10 ((a) P2N0, (b) P2N4).
(a)
(b)Nonlinear residual versus iteration steps ((a)ν=0.001,β=10 and (b) ν=0.00033,β=20).
(a)
(b)Streamlines forν=0.000333,β=30 ((a) P5N10 (wrong result), (b) P10N5 (right result)).
(a)
(b)The left subplot in Figure4 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:
(4.4)(1)R=f-(νAvold+N(vold))vold-Bpold,r=-BTvold;(2)(νAvmid+N(vold))vmid+Bpmid=R,BTvmid=r;(3)vnew=vold+vmid,pnew=pold+pmid.The difference between Picard method and Newton method is that the linear termN(vmid)vold is also dropped from (4.2), and thus the Picard method commonly referred to as the Oseen system.The left subplot in Figure5 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.Figures6 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 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.


## 4.1.  Numeric Example I

Consider a unit square domain with an exact solution given by(4.3)u(x,y)=(u1(x,y),u2(x,y)),p(x,y)=10(2x-1)(2y-1),u1(x,y)=10x2(x-1)2y(y-1)(2y-1),u2(x,y)=-10x(x-1)(2x-1)y2(y-1)2.f is determined by (1.1). After some computation using stretched gird, we have the following results.Figure2 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.Table 1
Numerical results of the FVM (ν=0.005, β=10).
h
∥
A
~
h
1
/
2
(
u
-
u
h
)
∥
0
/
∥
A
~
h
1
/
2
u
∥
0
Con. rate
∥
u
-
u
h
∥
0
/
∥
u
∥
0
Con. rate
∥
p
-
p
h
∥
0
/
∥
p
∥
0
Con. rate
1/16
0.06316279
—
0.00391173
—
0.07127911
—
1/32
0.03204341
0.97904713
0.00099767
1.97116812
0.03572585
0.99651095
1/64
0.01627933
0.97698618
0.00025462
1.97021278
0.01790786
0.99637536
1/128
0.00827253
0.97664010
0.00006499
1.96996167
0.00897672
0.99633198Figure 2
The relative errors and convergence rates of velocity and pressure.Figure3 is the L2 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.Figure 3
L
2 relative error of numerical velocity versus iteration steps for different ν.


## 4.2.  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,ux(0,1)=ux(1,1)=0 and ux=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).Nonconvergence streamlines forν=0.001,β=10 ((a) Stokes, (b) N9).
(a)
(b)Streamlines forν=0.001,β=10 ((a) P2N0, (b) P2N4).
(a)
(b)Nonlinear residual versus iteration steps ((a)ν=0.001,β=10 and (b) ν=0.00033,β=20).
(a)
(b)Streamlines forν=0.000333,β=30 ((a) P5N10 (wrong result), (b) P10N5 (right result)).
(a)
(b)The left subplot in Figure4 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:
(4.4)(1)R=f-(νAvold+N(vold))vold-Bpold,r=-BTvold;(2)(νAvmid+N(vold))vmid+Bpmid=R,BTvmid=r;(3)vnew=vold+vmid,pnew=pold+pmid.The difference between Picard method and Newton method is that the linear termN(vmid)vold is also dropped from (4.2), and thus the Picard method commonly referred to as the Oseen system.The left subplot in Figure5 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.Figures6 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 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.


## 5. Conclusions

The main work in this paper is the demonstration of the optimal order in theL2 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.



---

*Source: 251908-2012-11-11.xml*