Two-step algorithms for the stationary incompressible Navier–Stokes equations with friction boundary conditions

Accepted Manuscript Two-step algorithms for the stationary incompressible Navier–Stokes equations with friction boundary conditions

Hailong Qiu, Rong An, Liquan Mei, Liang Xue

PII: DOI: Reference:

S0168-9274(17)30118-6 http://dx.doi.org/10.1016/j.apnum.2017.05.003 APNUM 3209

To appear in:

Applied Numerical Mathematics

Received date: Revised date: Accepted date:

17 June 2016 12 February 2017 8 May 2017

Please cite this article in press as: H. Qiu et al., Two-step algorithms for the stationary incompressible Navier–Stokes equations with friction boundary conditions, Appl. Numer. Math. (2017), http://dx.doi.org/10.1016/j.apnum.2017.05.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Two-step algorithms for the stationary incompressible Navier-Stokes equations with friction boundary conditions Hailong Qiua,∗, Rong Anb , Liquan Meic , Liang Xuea a

School of Mathematics and Physics, Yancheng Institute of Technology, Yancheng, 224051, China b College of Mathematics and Information Science, Wenzhou University, 325035 Wenzhou, China c School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, 710049, China

Abstract Two-step algorithms for the stationary incompressible Navier-Stokes equations with friction boundary conditions are considered in this paper. Our algorithms consist of solving one Navier-Stokes variational inequality problem used the linear equal-order finite element pair (i.e., P1 − P1 ) and then solving a linearization variational inequality problem used the quadratic equal-order finite element pair (i.e., P2 − P2 ). Moreover, the stability and convergence of our two-step algorithms are derived. Finally, numerical tests are presented to check theoretical results. Keywords: Navier-Stokes equations, Friction boundary conditions, Linear equal-order pair, Quadratic equal-order pair, Two-step strategy, Error estimate. 1. Introduction Incompressible viscous flow with friction boundary conditions has been successfully applied to some flow phenomena in environmental and medical problems such as oil flow over or beneath sand layer and blood flow in the thoracic aorta [7]. The governing equations are the following incompressible ∗

Corresponding author. E-mail address: [email protected] (H. Qiu).

Preprint submitted to Applied Numerical Mathematics

May 10, 2017

Navier-Stokes ones:  −μΔu + (u · ∇)u + ∇p = f, divu = 0,

in Ω, in Ω,

with the friction boundary conditions:  u = 0, on Γ, un = 0, −στ (u) ∈ g∂|uτ |, on S,

(1)

(2)

where Ω is a bounded convex domain in R2 with Lipschitz continuous boundary ∂Ω partitioned into nonoverlapping parts, i.e., Γ ∩ S = ∅, Γ ∪ S = ∂Ω. Here the vector u = (u1 (x), u2 (x)) and the scalar p = p(x) denote the velocity and the pressure of the incompressible flow at the space point x, respectively, and f = (f1 (x), f2 (x)) represents some body force vector, g a scalar function and μ > 0 the viscous coefficient. We assume that un = u·n and uτ = u−un n are the normal and tangential components of the velocity, where n stands for the unit vector of the external normal to S. στ (u) = σ − σn n, independent of p, is the tangential component of the stress vector σ which is defined by ∂uj ∂ui σi = σi (u, p) = (μeij (u) − pδij )nj , where eij (u) = ∂x j + ∂xi , i, j = 1, 2. We assume that the set ∂ψ(a) denotes a subdifferential of the function ψ at the point a. From the definition of the subdifferential in the following Section, we can know that the variational formulation of (1)-(2) is the variational inequality problem of the second kind. For the friction boundary conditions (2), Fujita in [8] showed well-posedness of weak solution to the Stokes problem. Subsequently, Saito in [32] studied regularity of the weak solution by Yoshida’s regularized method and difference quotients for the Stokes problem. Meanwhile, the well-posedness of the weak solution to the Navier-Stokes problems under this typical boundary conditions were obtained in [23, 27, 15]. For numerical methods for the stationary Navier-Stokes equations, Li and An in [24] studied pressure stabilized finite element method for this kind friction boundary conditions. Kashiwabara in [16, 17] studied discrete variational inequality problem for the Stokes equations with slip/leak boundary conditions of friction type. Other theoretical and numerical results for the stationary and unstationary Stokes/Navier-Stokes problems with such boundary conditions can be found in [25, 26, 28, 3, 29, 30, 4, 31]. Despite the available computing power is developing very fast nowdays, computing the Navier-Stokes equations numerically is still remain some important but challenging problems. As we know, using finite element methods 2

to solve the Navier-Stokes equations may bring out some shortcomings: the oscillations due to the dominating convection and the violation of the discrete inf-sup condition and so on. Therefore, to increase the efficiency of numerical algorithms, to the authors’ knowledge, a good idea is two-grid method. The main idea of two-grid discretization method is to capture the low modes, or global solution envelope by solving an initial approximation problem on a coarse mesh. The fine structures are captured by computing one linearized problem on a fine mesh. Some details of two-grid strategy can be found in the works of Layton et al. [18, 19], Ervin et al. [5], He et al. [10, 11, 13]. Although some stable mixed finite element pairs have been studied over the years [33, 9], while some finite element pairs not satisfying the discrete inf-sup condition are also popular. Recently, a new family of lower-older stabilized method based on polynomial pressure projection was studied in [1, 2], which has some prominent features: without higher-order derivatives and edge-based data structures, and stabilization being completely local at the element grid. Hence, this stabilized method gains more and more popularity in computational fluid dynamics. Moreover, based on the work of Bochev et al. [1, 2], Li et al. in [21, 22] studied on stabilization of the linear equal-order stabilized finite element method for the Stokes/NavierStokes equations. Subsequently, Zheng et al. in [34] studied the quadratic equal-order stabilized finite element method for the Stokes equations. these methods obtain that the discrete inf-sup condition holds by adding to a bilinear term. The focus of this paper is to consider two-step algorithms for the stationary incompressible Navier-Stokes equations with friction boundary conditions and give the corresponding numerical analysis. The main idea of our twostep algorithms for computing Navier-Stokes variational inequality problem is to solve an initial solution based on the linear equal-order element P1 − P1 , then to solve a linearized problem (Stokes linearization, Oseen linearization, Newton linearization) based on the quadratic equal-order element P2 − P2 . More specifically, our method includes three algorithms: Stokes two-step algorithm, Oseen two-step algorithm and Newton two-step algorithm. Denote (uh , ph ) the Stokes/Oseen/Newton two-step algorithm solution to the problem (1)-(2), then we derive that our two-step algorithm solution (uh , ph ) satisfies the following error estimate u − uh V + p − ph  ≤ ch2 . This paper is organized as follows. In next Section, we will introduce 3

some notations, and give the variational formulation of the problem (1)(2). In Section 3, we will recall two stabilized finite element methods. The stability and convergence of two-step algorithms are established in Section 4. The numerical results is presented in last Section. Throughout this paper, the letter c always denotes some positive constant independent of the mesh parameter h and may be different even in the same formulation. 2. Preliminaries As in [32], we introduce the definition of the subdifferential property [6]. Let ψ : R2 → R = (−∞, +∞] be a given function possessing the properties of convexity and weak semi-continuity from below (ψ is not identical with +∞). We say that the set ∂ψ(a) is a subdifferential of the function ψ at the point a if and only if ∂ψ(a) = {b ∈ R2 : ψ(h) − ψ(a) ≥ b · (h − a),

∀h ∈ R2 }.

Introduce the following Hilbert spaces: V = {u ∈ H 1 (Ω)2 , u|Γ = 0, u · n|S = 0}, Vσ = {u ∈ V,

∇ · u = 0},

V0 = H01 (Ω)2 ,  2 2 M = L0 (Ω) = {q ∈ L (Ω), qdx = 0}. Ω

k

2

Let  · k be the norm in Hilbert space H (Ω) . Let (·, ·) and  ·  be the inner product and the norm in L2 (Ω)2 . We equip the inner product and the norm in V by (∇·, ∇·) and  · V = ∇ · , respectively, because ∇ ·  is equivalent to  · 1 . We define the following bilinear forms and trilinear form: ⎧ ∀ u, v ∈ V, ⎨ a(u, v) = μ(eij (u), eij (v)), d(v, p) = (p,∇ · v), ∀ v ∈ V, p ∈ M, ⎩ b(u, v, w) = Ω (u · ∇)v · wdx, ∀ u, v, w ∈ V. Using Korn’s inequality [14], i.e., there exists a constant μ0 > 0 such that a(v, v) ≥ μ0 v2V ,

v ∈ V.

Moreover, if ∇ · u = 0, then trilinear form b(·, ·, ·) satisfies b(u, v, w) = ((u · ∇)v, w) + 12 ((∇ · u)v, w) = 12 ((u · ∇)v, w) − 12 ((u · ∇)w, v), 4

∀ u, v, w ∈ V.

Therefore b(u, v, w) = −b(u, w, v), Denote N = sup u,v,w∈V

∀ u, v, w ∈ V.

b(u, v, w) , uV vV wV

then it is not difficult to prove that the following estimates ⎧ ⎨ b(u, v, v) = 0, b(u, v, w) ≤ N uV vV wV , ⎩ |b(u, v, w)| + |b(v, u, w)| + |b(w, u, v)| ≤ N uV v2 w,

hold [9, 11, 20]:

∀ u, v ∈ V, ∀ u v w ∈ V, ∀u ∈ V, v ∈ H 2 (Ω)2 , w ∈ L2 (Ω)2 . (3) With the above notations, the variational formulation of the problem (1)-(2) reads as follows: ⎧ ⎨ Find (u, p) ∈ V × M such that a(u, v − u) + b(u, u, v − u) + j(vτ ) − j(uτ ) − d(v − u, p) ≥ (f, v − u), ∀ v ∈ V, ⎩ d(u, q) = 0, ∀ q ∈ M, (4)  1/2 2 where j(η) = S g|η|ds for η ∈ H (S) . In [32], Saito obtained the following inf-sup condition holds, i.e., there exists some positive β > 0 such that βq ≤ sup v∈V



d(v, q) . vV

Hence, the variational inequality (4) is equivalent to Find u ∈ Vσ such that a(u, v − u) + b(u, u, v − u) + j(vτ ) − j(uτ ) ≥ (f, v − u),

∀ v ∈ Vσ . (5) The existence and uniqueness of the solution of the variational problem (5) has been proved in [23, 24], we borrow Theorem 2.1. Given f ∈ L2 (Ω)2 and g ∈ L2 (S). If 4κ1 N (f  + gL2 (S) ) < 1, μ20

(6)

then the variational inequality problem (5) admits a unique solution u ∈ μ0 1 Kσ = {v ∈ Vσ : vV ≤ 2κ (f  + gL2 (S) ) < 2N }, where κ1 > 0 satisfies μ0 |(f, v) − j(vτ )| ≤ κ1 (f  + gL2 (S) )vV , 5

∀ v ∈ Vσ .

3. Stabilized finite element methods We assume that h be a real positive parameter tending to 0. The finite element subspaces Vh × Mh , Vh × Mh of V × M is characterized by Kh , a partitioning of Ω into triangles K with the mesh size h, assumed to be uniformly regular in the usual sense. Let e = ∂K ∩ S. In addition, we denote Σ1,h = set of all vertices of triangles in Kh , Σ2,h = set of all midpoints of sides of triangles in Kh , and Σh S1,h S2,h S˙ 1,h S˙ 2,h

= Σ1,h ∪ Σ2,h , = S ∩ Σ1,h = {Q1 , Q2 , · · · , Qm , Qm+1 }, = S ∩ Σh = {Q1 , Q3/2 , Q2 , · · · , Qm , Qm+1/2 , Qm+1 }, = S ∩ Σ1,h = {Q2 , · · · , Qm }, = S ∩ Σh = {Q3/2 , Q2 , · · · , Qm , Qm+1/2 }.

Then we introduce finite element spaces: = {u ∈ C 0 (Ω)2 ∩ V : u|K ∈ P1 (K)2 , ∀K ∈ Kh }, = {q ∈ C 0 (Ω) ∩ M : q|K ∈ P1 (K), ∀K ∈ Kh }, = {u ∈ C 0 (Ω)2 ∩ V : u|K ∈ P2 (K)2 , ∀K ∈ Kh }, = {q ∈ C 0 (Ω) ∩ M : q|K ∈ P2 (K), ∀K ∈ Kh },

Vh Mh Vh Mh

where Pi (K) is the space of all polynomials on K of degree less than i (i = 1 or i = 2). As note, the pair Vh × Mh or Vh × Mh does not satisfy the following discrete inf-sup conditions: d(v h , q h ) β0 q  ≤ sup , h v h ∈Vh v V

∀q h ∈ Mh ,

d(v h , q h ) , v h V

∀q h ∈ Mh ,

h

and 

β0 q h  ≤ sup

v h ∈Vh



where β0 > 0, β0 > 0 are independent of h. In order to get discrete inf-sup conditions, we need two stabilized generalized bilinear forms [22, 34]. Let B(uh , ph ; v h , q h ) = a(uh , v h ) − b(v h , ph ) + b(uh , q h ) + G(ph , q h ), where G(ph , q h ) is defined by G(ph , q h ) = ((I − Π)ph , (I − Π)q h ), 6

and let ¯ h , q h ), B(uh , ph ; v h , q h ) = a(uh , v h ) − b(v h , ph ) + b(uh , q h ) + G(p ¯ h , q h ) is defined by where G(p ¯ h , q h ) = h2 ((I − Π)∇ph , (I − Π)∇q h ), G(p Then the projection operator Π satisfies the following properties [22, 34]: Πp ≤ cp, p − Πp ≤ chp,

p ∈ L2 (Ω), p ∈ H 1 (Ω).

As in [22, 34], we define the stabilized term G(ph , q h ) based on local Gauss integration as follows:    h h h h G(p , q ) = { p q dx − ph q h dx}, ph , q h ∈ Mh , K∈Kh

K,2

K,1

 where K,i g(x)dx denotes a local Guass integral over K that is exact for polynomials of degree i (i = 1, 2). And we also define the stabilized term ¯ h , q h ) as follows: G(p    h h 2 h h ¯ G(p , q ) = h { ∇p ∇q dx − ∇ph ∇q h dx}, ph , q h ∈ Mh . K∈Kh

K,2

K,1

Moveover, the inverse inequalities are recalled v h 1 ≤ ch−1 v h , ph 1 ≤ ch−1 ph . We are now in a position to consider the following stabilized problems: Find (uh , ph ) ∈ Vh × Mh such that for all (v, q) ∈ Vh × Mh B(uh , ph ; v − uh , q − ph ) + b(uh , uh , v − uh ) + j(vτ ) − j(uhτ ) ≥ (f, v − uh ). (7) And find (uh , ph ) ∈ Vh × Mh such that for all (v, q) ∈ Vh × Mh B(uh , ph ; v − uh , q − ph ) + b(uh , uh , v − uh ) + j(vτ ) − j(uhτ ) ≥ (f, v − uh ). (8) The following theorem gives the continuity property and weak coercivity of the bilinear forms B(·, ·; ·, ·) for the finite element spaces (Vh , Mh )×(Vh , Mh ) ¯ ·; ·, ·) for the finite element spaces (Vh , Mh ) × (Vh , Mh ) [22, 34]: and B(·, 7

Theorem 3.1. For all (uh , ph ), (v h , ph ) ∈ Vh × Mh , B(·, ·; ·, ·) satisfies the continuity property |B(uh , ph ; v h , q h )| ≤ β1 (uh V +ph )(v h V +q h ),

∀ (uh , ph ), (v h , q h ) ∈ (Vh , Mh ),

and the weak coercivity property β2 (uh V + ph ) ≤

B(uh , ph ; v h , q h ) , h  + q h  v V  (v h ,q h )∈(V ,M ) h h sup

h , Mh ), ∀ (uh , ph ) ∈ (V

h = {v ∈ V0 : v|K ∈ P1 (K), ∀ K ∈ Kh } be the finite element where V subspace of V0 , β1 > 0, β2 > 0 are two constants independent of h. ¯ ·; ·, ·) satisfies the conMoreover, for all (uh , ph ), (v h , ph ) ∈ Vh × Mh , B(·, tinuity property ¯ h , ph ; v h , q h )| ≤ β3 (uh V +ph )(v h V +q h ), |B(u

∀ (uh , ph ), (v h , q h ) ∈ (Vh , Mh ),

and the weak coercivity property β4 (uh V + ph ) ≤

sup  ,M (v h ,q h )∈(V h h

¯ h , ph ; v h , q h ) B(u , h  + q h  v V )

, M ), ∀ (uh , ph ) ∈ (V h h

= {v ∈ V : v| ∈ P (K), ∀ K ∈ K } be the finite element where V h 0 K 2 h subspace of V0 , β3 > 0, β4 > 0 are two constants independent of h. As in [24, 28, 31], we give the following results. Theorem 3.2. Under the assumption of the Theorem 2.1, If the exact solution (u, p) be in V ∩ H 2 (Ω) ∩ H 2 (∂Ω) × M ∩ H 1 (Ω). Then (uh , ph ) ∈ Vh × Mh of problem (7) satisfies the following stability and error estimate: μ0 uh V ≤ κ1 (f  + gL2 (S) ), u − uh  + h(u − uh V + p − ph ) ≤ ch2 . Moreover, if the exact solution (u, p) be in V ∩ H 3 (Ω) ∩ H 3 (∂Ω) × M ∩ H 2 (Ω). Then (uh , ph ) ∈ Vh ×Mh of problem (8) satisfies the following stability and error estimate: μ0 uh V ≤ κ1 (f  + gL2 (S) ), u − uh  + h(u − uh V + p − ph ) ≤ ch3 . 8

4. Two-step algorithms Two-step stabilized algorithms based on local Gauss integration we study are as follows. As in [34], we introduce ˜ ¯ ¯ q), B((u, p); (v, q)) = B((u, p); (v, q)) − G(p,

∀(v, q) ∈ Vh × Mh ,

and define the projection operator (Rh , Qh ) : V × M → Vh × Mh through ¯ ˜ B((R h u, Qh p); (v, q)) = B(((u, p); (v, q))),

∀(v, q) ∈ Vh × Mh ,

which is well defined by Theorem 3.1. A similar proof to that applied in [22, 34] obtains the following approximation result: u − Rh u + h(u − Rh uV + p − Qh p) ≤ ch3 . Algorithm 4.1 (Stokes two-step algorithm) Step I: Solve a Navier-Stokes variational inequality problem by linear equal-order finite element pair (P1 − P1 ), i.e., find (uh , ph ) ∈ Vh × Mh such that for all (v, q) ∈ Vh × Mh , B(uh , ph ; v − uh , q − ph ) + b(uh , uh , v − uh ) + j(vτ ) − j(uhτ ) ≥ (f, v − uh ). (9) Step II: Solve a general Stokes variational inequality problem by quadratic equal-order finite element pair (P2 − P2 ), i.e., find (uh , ph ) ∈ Vh × Mh such that for all (v, q) ∈ Vh × Mh , ¯ h , ph ; v − uh , q − ph ) + b(uh , uh , v − uh ) + j(vτ ) − j(uhτ ) ≥ (f, v − uh ). (10) B(u Under the condition of Theorem 2.1, the discrete problems given in Algorithms 4.1-4.3 (below) admit unique solution solutions (uh , ph ) and (uh , ph ). The proofs are similar to the theorem 3.2 in the [24] (or [28]), thus we skep it. Now, we state and derive the stability and error estimate of Algorithm 4.1. Theorem 4.1. Suppose that (uh , ph ) ∈ (Vh , Mh ) is the solution to the problem (10), then satisfies 2 4 2 ¯ h , ph ) ≤ 2κ1 (f +gL2 (S) )2 + 2κ1 N (f +gL2 (S) )4 . (11) μ0 uh 2V +2G(p μ0 μ50

9

Proof: Letting v = 0, q = ph in (9), using (3) and Young’s inequality, we easily find that μ0 uh 2V + G(ph , ph ) ≤ (f, uh ) − j(uhτ ) ≤ κ1 (f  + gL2 (S) )uh V μ0 κ2 ≤ uh 2V + 1 (f  + gL2 (S) )2 , 2 2μ0 Thus we derive μ0 uh 2V + 2G(ph , ph ) ≤

κ21 (f  + gL2 (S) )2 . μ0

(12)

Letting v = 0, q = ph in (10), it yields ¯ h , ph ) μ0 uh 2V + G(p ≤ (f, uh ) − j(uhτ ) + b(uh , uh , uh ) ≤ κ1 (f  + gL2 (S) )uh V + N uh 2V uh V μ0 κ2 κ4 N 2 ≤ uh 2V + 1 (f  + gL2 (S) )2 + 1 5 (f  + gL2 (S) )4 . 2 μ0 μ0 Then we obtain (11).  Theorem 4.2. Under the assumption of the Theorem 2.1, If the exact solution (u, p) is in V ∩ H 3 (Ω) ∩ H 3 (∂Ω) × M ∩ H 2 (Ω), then (uh , ph ) defined by scheme (10) satisfies the following bound: u − uh V + p − ph  ≤ ch2 ,

(13)

where c > 0 is independent of h. Proof: Taking v = Rh u, q = ph − Qh p in (10), Thanks to the definition of ¯ we have the bilinear form B, ¯ h − Qh p, ph − Qh p) μuh − Rh u2V + G(p ¯ = B(uh − Rh u, ph − Qh p; uh − Rh u, ph − Qh p) ¯ h u, Qh p; uh − Rh u, ph − Qh p) ¯ h , ph ; uh − Rh u, ph − Qh p) − B(R = B(u = a(uh , uh − Rh u) − d(uh − Rh u, ph ) + d(uh , ph − Qh p) ¯ h u, Qh p; uh − Rh u, ph − Qh p) ¯ h , ph − Qh p) − B(R +G(p h ≤ (f, uh − Rh u) + b(u , uh , uh − Rh u) ¯ h u, Qh p; uh − Rh u, ph − Qh p). +j(Rh uτ ) − j(uhτ ) − B(R (14) 10

Setting v = uh and v = 2u − Rh u in (4), one has a(u, uh − u) + b(u, u, uh − u) − d(uh − u, p) + j(uhτ ) − j(uτ ) ≥ (f, uh − u) and a(u, u−Rh u)+b(u, u, u−Rh u)−d(u−Rh u, p)+j((2u − Rh u)τ )−j(uτ ) ≥ (f, u−Rh u), which gives a(u, uh − Rh u) + b(u, u, uh − Rh u) − d(uh − Rh u, p) +j((2u − Rh u)τ ) + j(uhτ ) − 2j(uτ ) ≥ (f, uh − Rh u). Substituting the above inequality into (14) yields ¯ h − Qh p, ph − Qh p) μuh − Rh u2V + G(p ≤ a(u, uh − Rh u) + b(u, u, uh − Rh u) + j(Rh uτ ) − 2j(uτ ) + j(2uτ − Rh uτ ) −d(uh − Rh u, p) + b(uh , uh , Rh u − uh ) − a(Rh u, uh − Rh u) ¯ h p, ph − Qh p) +d(uh − Rh u, Qh p) − d(Rh u, ph − Qh p) − G(Q = a(u − Rh u, uh − Rh u) − d(uh − Rh u, p − Qh p) + d(u − Rh u, ph − Qh p) +b(u, u, uh − Rh u) − b(uh , uh , uh − Rh u) + j(Rh uτ ) ¯ h p, ph − Qh p) −2j(uτ ) + j(2uτ − Rh uτ ) − G(Q ≤ |a(u − Rh u, uh − Rh u)| + |d(u − Rh u, ph − Qh p) − d(uh − Rh u, p − Qh p)| +|b(u, u, uh − Rh u) − b(uh , uh , uh − Rh u)| +|j(Rh uτ ) − 2j(uτ ) + j(2uτ − Rh uτ )| ¯ h p, ph − Qh p) −G(Q = J1 + · · · + J 5 , (15) where we use d(u, q) = 0 for all q ∈ M . Now, we estimate the five terms of the right side of (15). Due to Young inequality, we estimate J1 , J2 as follows: J1 ≤ μu − Rh uV uh − Rh uV μ ≤ uh − Rh u2V + cu − Rh u2V 9 μ ≤ ch4 + uh − Rh u2V . 9 and J2 ≤ u − Rh uV ph − Qh p + uh − Rh uV p − Qh p μ ≤ δph − Qh p2 + cu − Rh u2V + uh − Rh u2V + cp − Qh p2 9 μ ≤ ch4 + δph − Qh p2 + uh − Rh u2V . 9 11

where δ > 0 is a sufficiently small constant. By using (3), we estimate J3 as follows: J3 = |b(u, u, Rh u − uh ) − b(uh , uh , Rh u − uh )| = |b(uh − u, u, Rh u − uh ) − b(uh , u − uh , Rh u − uh )| = |b(uh − u, u, Rh u − uh ) + b(u, uh − u, Rh u − uh ) +b(uh − u, uh − u, Rh u − uh )| ≤ 2N u2 Rh u − uh V uh − u + N u − uh 2V Rh u − uh V μ ≤ c(u − uh 2 + u − uh 4V ) + Rh u − uh 2V 9 μ ≤ ch4 + Rh u − uh 2V . 9 For J4 , we have

J4 ≤ cRh u − uL2 (S) .

Estimate J5 by J5 = |G(Qh p − p, ph − Qh p) + G(p, ph − Qh p)| ≤ h2 Qh ∇p − ∇p · ∇(ph − Qh p) + h2 ∇p − Π∇p · ∇(ph − Qh p) 1 1 ≤ 2δh2 ∇(ph − Qh p)2 + h2 Qh ∇p − ∇p2 + h2 ∇p − Π∇p2 4δ 4δ ≤ ch4 + 2δc0 ph − Qh p2 , where we apply inverse inequality. Then substituting J1 , · · · , J5 into (15), dropping G(ph − Qh p, ph − Qh p), we can obtain uh − Rh u2V ≤ ch4 +

3δc0 ph − Qh p2 . 2μ

According to triangular inequality, we have u − uh 2V

≤ 2u − Rh u2V + 2uh − Rh u2V 3δc0 ph − Qh p2 . ≤ ch4 + μ

(16)

Next we estimate ph − Qh p. By using Theorem 3.1, we have ¯ h − Rh u, ph − Qh p; wh , qh ) B(u wh V + qh   ,M ) (wh ,qh )∈(V h h ¯ h − u, ph − p; wh , qh ) + B(u ¯ − Rh u, p − Qh p; wh , qh ) B(u = sup . wh V + qh   ,M ) (w ,q )∈(V

β4 ph − Qh p ≤

sup

h

h

h

h

(17) 12

⊂ V , taking v = u ± w in (4) yields Let wh ∈ V h 0 h a(u, wh ) + b(u, u, wh ) − d(wh , p) ≥ (f, wh ),

. ∀ wh ∈ V h

and . ∀ wh ∈ V h

a(u, −wh ) + b(u, u, −wh ) − d(−wh , p) ≥ (f, −wh ), Thus a(u, wh ) + b(u, u, wh ) − d(wh , p) = (f, wh ),

. ∀ wh ∈ V h

Similarly, taking v = uh ± wh in (10), one has a(uh , wh ) + b(uh , uh , wh ) − d(wh , ph ) = (f, wh ),

. ∀ wh ∈ V h

Hence a(u − uh , wh ) + b(u, u, wh ) − b(uh , uh , wh ) = d(wh , p − ph ),

. ∀ wh ∈ V h

¯ we have By the definition of B, ¯ h − u, ph − p; wh , qh ) B(u ¯ h − p, qh ) = a(uh − u, wh ) − d(wh , ph − p) + d(uh − u, qh ) + G(p h h ¯ = b(u , u , wh ) − b(u, u, wh ) − G(p, qh ) ¯ qh ) = b(uh − u, u, wh ) + b(u, uh − u, wh ) − G(p, h h h ¯ qh ) = b(u − u, u, wh ) + b(u − u, u − u, wh ) + b(u, uh − u, wh ) − G(p, h h 2 2 ≤ 2N u2 u − u wh V + N u − uV wh V + h ∇p − Π∇p∇qh  ≤ ch2 wh V + ch2 qh , (18) ¯ where we use d(uh −u, qh )+ G(ph , qh ) = 0 and inverse inequality. Substituting (18) into (17), we get that

Thus we obtain

ph − Qh p ≤ ch2 .

(19)

u − uh V ≤ ch2 .

(20)

Then according to the triangular inequality and (19), we have p − ph  ≤ ch2 . 13



Algorithm 4.2 (Oseen two-step algorithm) Step I: Solve a Navier-Stokes variational inequality problem by linear equal-order finite element pair (P1 − P1 ), i.e., find (uh , ph ) ∈ (Vh , Mh ) by (9). Step II: Solve a Oseen variational inequality problem by quadratic equalorder finite element pair (P2 − P2 ), i.e., find (uh , ph ) ∈ Vh × Mh such that for all (v, q) ∈ Vh × Mh , ¯ h , ph ; v − uh , q − ph ) + b(uh , uh , v − uh ) + j(vτ ) − j(uhτ ) ≥ (f, v − uh ). (21) B(u We now in a position to show and prove the stability and error estimate of Algorithm 4.2. Theorem 4.3. Suppose that (uh , ph ) ∈ (Vh , Mh ) is the solution to the problem (21), then satisfies 2 ¯ h , ph ) ≤ κ1 (f  + gL2 (S) )2 . μ0 uh 2V + 2G(p μ0

Proof:

(22)

Letting v = 0, q = ph in (21) and using (4), it yields ¯ h , ph ) μ0 uh 2V + G(p ≤ (f, uh ) − j(uhτ ) ≤ κ1 (f  + gL2 (S) )uh V μ0 κ2 ≤ uh 2V + 1 (f  + gL2 (S) )2 . 2 2μ0

Then we obtain (22).  Theorem 4.4. Under the assumption of the Theorem 2.1, If the exact solution (u, p) is in V ∩ H 3 (Ω) ∩ H 3 (∂Ω) × M ∩ H 2 (Ω), then (uh , ph ) defined by scheme (21) satisfies the following bound: u − uh V + p − ph  ≤ ch2 , where c > 0 is independent of h.

14

(23)

Proof:

Proceeding as the in proof of Theorem 4.2, we obtain

μuh − Rh u2V + G(ph − Qh p, ph − Qh p) ≤ a(u, uh − Rh u) + b(u, u, uh − Rh u) − b(uh , uh , uh − Rh u) +j(Rh uτ ) − 2j(uτ ) + j(2uτ − Rh uτ ) − d(uh − Rh u, p) − G(Qh p, ph − Qh p) −a(Rh u, uh − Rh u) + d(uh − Rh u, Qh p) − d(Rh u, ph − Qh p) = a(u − Rh u, uh − Rh u) − d(uh − Rh u, ph − Qh p) +d(u − Rh u, ph − Qh p) + b(u, u, uh − Rh u) − b(uh , uh , uh − Rh u) +j(Rh uτ ) − 2j(uτ ) + j(2uτ − Rh uτ ) − G(Qh p, ph − Qh p) ≤ |a(u − Rh u, uh − Rh u)| + |d(u − Rh u, ph − Qh p) −d(uh − Rh u, p − Qh p)| + |b(u, u, uh − Rh u) − b(uh , uh , uh − Rh u)| +|j(Rh uτ ) − 2j(uτ ) + j(2uτ − Rh uτ )| − G(Qh p, ph − Qh p) = J6 + · · · + J10 . (24) In the above equation, J6 , J7 , J9 , J10 have been estimated in the proof of Theorem 4.2. Here, we only estimate J8 as follows: J8 = |b(u, u, uh − Rh u) − b(uh , uh , uh − Rh u)| = |b(uh , Rh u − uh , Rh u − uh ) + b(uh , u − Rh u, Rh u − uh ) + b(u − uh , u, Rh u − uh )| = |b(uh , Rh u − u, Rh u − uh ) + b(uh − u, u, Rh u − uh )| ≤ N uh V Rh u − uV Rh u − uh V + N u2 Rh u − uh V uh − u μ ≤ c(Rh u − u2V + uh − u2 ) + Rh u − uh 2V 9 μ 4 2 ≤ ch + Rh u − uh V . 9 ¯ h − Qh p, ph − Qh p) again, We substitute J6 , · · · , J10 into (24), dropping G(p according to triangular inequality, we have u − uh 2V

≤ 2u − Rh u2V + 2uh − Rh u2V 3δc0 Qh p − ph 2 . ≤ ch4 + μ

15

(25)

¯ we derive By the definition of B, ¯ h − u, ph − p; wh , qh ) B(u = a(uh − u, wh ) − d(wh , ph − p) + d(uh − u, qh ) + G(ph − p, qh ) = b(uh , uh , wh ) − b(u, u, wh ) − G(p, qh ) = b(uh − u, uh , wh ) + b(u, uh − u, wh ) − G(p, qh ) = b(uh − u, uh − u, wh ) + b(uh − u, u, wh ) + b(u, uh − u, wh ) − G(p, qh ) ≤ N uh − uV uh − uV wh V + uh − uu2 wh V +N uV u − uh V wh V − G(p, qh ) ≤ (c1 hu − uh V + ch2 + c2 u − uh V )wh V + ch2 qh , (26) where we use d(uh −u, qh )+G(ph , qh ) = 0 and inverse inequality. Substituting (26) into (17) and according to Theorem 3.1, we obtain that β4 Qh p − ph  ≤ (c1 h + c2 )u − uh V + ch2 .

(27)

Then for sufficiently small h, δ such that

1 3δc0 c1 h + c2 · < . μ β4 2 Thus we derive

u − uh V ≤ ch2 .

(28)

Finally, using the triangular inequality and (27)-(28), we get p − ph  ≤ ch2 .  Algorithm 4.3 (Newton two-step algorithm) Step I: Solve a Navier-Stokes variational inequality problem by linear equal-order finite element pair (P1 − P1 ), i.e., find (uh , ph ) ∈ (Vh , Mh ) by (9). Step II: Solve a Newton linearized variational inequality problem by quadratic equal-order finite element pair (P2 −P2 ), i.e., find (uh , ph ) ∈ Vh ×Mh such that for all (v, q) ∈ Vh × Mh , ¯ h , ph ; v − uh , q − ph ) + b(uh , uh , v − uh ) + b(uh , uh , v − uh ) B(u +j(vτ ) − j(uhτ ) ≥ (f, v − uh ) + b(uh , uh , v − uh ).

(29)

We are now in a position to give and prove the stability and error estimate of Algorithm 4.3. 16

Theorem 4.5. Suppose that (uh , ph ) ∈ (Vh , Mh ) is the solution to the problem (29), then satisfies 4 2 ¯ h , ph ) ≤ 4κ1 (f +gL2 (S) )2 + 4κ1 N (f +gL2 (S) )4 . (30) μ0 uh 2V +4G(p μ0 μ50

Proof: Letting v = 0, q = ph in (29), using (3) and Young’s inequality, we easily find that ¯ h , ph ) μ0 uh 2V − |b(uh , uh , uh )| + G(p ≤ (f, uh ) − j(uhτ ) + b(uh , uh , uh ) ≤ κ1 (f  + gL2 (S) )uh V + N uh 2V uh V μ0 κ2 N2 h 4 ≤ uh 2V + 1 (f  + gL2 (S) )2 + u V . 2 μ0 μ0 Since μ0 uh 2V − |b(uh , uh , uh )| ≥ μ0 uh 2V − N uh 2V uh V ≥

3μ0 uh 2V . 4

Then we have 2 2 3μ0 ¯ h , ph ) ≤ μ0 uh 2 + κ1 (f  + gL2 (S) )2 + N uh 4 . uh 2V + G(p V V 4 2 μ0 μ0

Using (12), we obtain (30).  Theorem 4.6. Under the assumption of the Theorem 2.1, If the exact solution (u, p) is in V ∩ H 3 (Ω) ∩ H 3 (∂Ω) × M ∩ H 2 (Ω), then (uh , ph ) defined by scheme (29) satisfies the following bound: u − uh V + p − ph  ≤ ch2 , where c > 0 is independent of h.

17

(31)

Proof:

Proceeding as the in proof of Theorem 4.2, we obtain

μuh − Rh u2V + G(ph − Qh p, ph − Qh p) ≤ a(u, uh − Rh u) + b(u, u, uh − Rh u) + b(uh , uh , uh − Rh u) −b(uh , uh , uh − Rh u) − b(uh , uh , uh − Rh u) + j(Rh uτ ) −2j(uτ ) + j(2uτ − Rh uτ ) − d(uh − Rh u, p) − G(Qh p, ph − Qh p) −a(Rh u, uh − Rh u) + d(uh − Rh u, q h ) − d(Rh u, p − Qh p) = a(u − Rh u, uh − Rh u) − d(uh − Rh u, ph − Qh p) +d(u − Rh u, ph − Qh p) + b(u, u, uh − Rh u) + b(uh , uh , uh − Rh u) −b(uh , uh , uh − Rh u) − b(uh , uh , uh − Rh u) + j(Rh uτ ) − 2j(uτ ) +j(2uτ − Rh uτ ) − G(Qh p, ph − Qh p) ≤ |a(u − Rh u, uh − Rh u)| + |d(u − Rh u, ph − Qh p) −d(uh − Rh u, p − Qh p)| + |b(u, u, uh − Rh u) + b(uh , uh , uh − Rh u) −b(uh , uh , uh − Rh u) − b(uh , uh , uh − Rh u)| +|j(Rh uτ ) − 2j(uτ ) + j(2uτ − Rh uτ )| − G(Qh p, ph − Qh p) = J11 + · · · + J15 . (32) In the above equation, J11 , J12 , J14 , J15 have been estimated in the proof of Theorem 4.2. Here, we only estimate J13 as follows: J13 = |b(uh , uh , Rh u − uh ) + b(uh , uh , Rh u − uh ) − b(u, u, Rh u − uh ) − b(uh , uh , Rh u − uh )| = |b(uh − u, u, Rh u − uh ) + b(uh , uh − u, Rh u − uh ) − b(uh − uh , uh − uh , Rh u − uh )| = |b(uh − Rh u, u, Rh u − uh ) + b(Rh u − u, u, Rh u − uh ) + b(uh , Rh u − u, Rh u − uh ) +b(Rh u − uh , Rh u − uh , Rh u − u) + b(uh − Rh u, Rh u − uh , Rh u − uh )| ≤ N uV Rh u − uh 2V + N uV Rh u − uV Rh u − uh V +|b(uh − Rh u, Rh u − u, Rh u − uh ) + b(Rh u, Rh u − u, Rh u − uh ) +b(Rh u − uh , Rh u − uh , Rh u − u) + b(uh − Rh u, Rh u − uh , Rh u − uh )| ≤ N uV Rh u − uh 2V + N uV Rh u − uV Rh u − uh V +N uh − Rh u2V Rh u − uV + Rh uV Rh u − uV Rh u − uh V +N Rh u − uh V Rh u − uh V Rh u − uh V + N uh − Rh u2V Rh u − uh V μ ≤ ch2 Rh u − uh V + c3 (h2 + h)Rh u − uh 2V + Rh u − uh 2V 2 5μ 4 2 2 Rh u − uh 2V . ≤ ch + c3 (h + h)Rh u − uh V + 9 μ For sufficiently small h such that c3 (h2 + h) = , We substitute J11 , · · · , J15 6 ¯ h − Qh p, ph − Qh p) again, according to triangular into (32), dropping G(p

18

inequality, we have u − uh 2V

≤ 2u − Rh u2V + 2uh − Rh u2V 36δc0 ≤ ch4 + Qh p − ph 2 . μ

(33)

¯ we have By the definition of B, ¯ h − u, ph − p; wh , qh ) B(u

= a(uh − u, wh ) − d(wh , ph − p) + d(uh − u, qh ) + G(ph − p, qh ) = b(uh , uh , wh ) + b(uh , uh , wh ) − b(u, u, wh ) − b(uh , uh , wh ) − G(p, qh ) = b(uh − u, u, wh ) + b(u, Rh u − u, wh ) + b(u − uh , u − Rh u, wh ) +b(uh − uh , Rh u − uh , wh ) + b(uh , uh − Rh u, wh ) − G(p, qh ) ≤ (N uV uh − uV + N uV u − Rh uV + N u − uh V u − Rh uV )wh V +N uh − uh V uh − Rh uV wh V + N uh V uh − Rh uV wh V +h2 ∇p − Π∇p∇qh  μ ≤ ( u − uh V + ch2 + ch2 u − uh V )wh V 2 +N (uh − uV + u − uh V )(u − uh V + u − Rh uV )wh V +N uh V (uh − uV + u − Rh uV )wh V + ch2 qh  μ ≤ ( u − uh V + ch2 + ch2 u − uh V )wh V 2 +cN (uh − uV + h)(h + h2 )wh V +N uh V (uh − uV + h2 )wh V + ch2 qh  ≤ (μ + c4 h2 + c5 N (h + h2 ))u − uh V wh V + ch2 wh V + ch2 qh , (34) where we use d(uh −u, qh )+G(ph , qh ) = 0 and inverse inequality. Substituting (34) into (17) and according to Theorem 3.1, we obtain that β4 Qh p − ph  ≤ (μ + c4 h2 + c5 N (h + h2 ))u − uh V + ch2 .

(35)

Then, for sufficiently small h, δ such that

36δc0 1 (μ + c4 h2 + c5 N (h + h2 )) < . 2 μβ4 2 Thus we show

u − uh V ≤ ch2 .

(36)

Then according to the triangular inequality and (35)-(36), we get p − ph  ≤ ch2 .  19

5. Numerical results In this section, we present the performance of our two-step algorithms described in the previous sections for Navier-Stokes equations with friction boundary conditions. We apply the Uzawa iteration algorithm introduced by [26] to solve variational inequality problem (4), which is equivalent to the following variational formulation: ⎧  ⎨ a(u, v) + b(u, u, v) − d(v, p) + S λgvτ ds = (f, v), ∀ v ∈ V, d(u, q) = 0, ∀ q ∈ M, (37) ⎩ λuτ = |uτ |, a.e. on S, where λ ∈ Λ = {γ ∈ L2 (S) : |γ| ≤ 1 a.e. on S}. The discretized formulations corresponding to (37) is ⎧  ⎨ a(uh , vh ) + b(uh , uh , vh ) − d(vh , ph ) + S λh gvhτ ds = (f, vh ), ∀ vh ∈ Vh or Vh , ∀ qh ∈ Mh or Mh , d(uh , qh ) = 0, ⎩ λh uhτ = |uhτ |, a.e. on S, (38) where ˜ h | |γh (Q)| ≤ 1, ∀Q ∈ S˙ 1,h (orS˙ 2,h )}, λh ∈ Λh = {γh ∈ Λ and ˜ h = {γh ∈ C 0 (S) | γh |e ∈ Pk (e), γh (Q1 ) = γh (Qm+1 ) = 0, k = 1, 2}. Λ Now we can compute the problem (37) by the following Uzawa iteration algorithm introduced in [26]: λ0 ∈ Λ is given ,

(39)

then λn is known, we solve (un , pn ) and λn+1 by   a(un , v) + b(un , un , v) − b(v, pn ) = (f, v) − S λn gvτ ds, b(un , q) = 0, ∀ q ∈ M,

∀ v ∈ V, (40)

and λn+1 = PΛ (λn + ρgunτ ), where PΛ (γ) = sup(−1, inf(1, γ)), 20

ρ > 0, ∀ γ ∈ L2 (S).

(41)

(a) Figure 1: Domain Ω

The discretized formulations corresponding to (39)-(41) are λh0 ∈ Λh

is given,

(42)

then λhn is known, we solve (uhn , phn ) and λh(n+1) by   a(uhn , vh ) + b(uhn , uhn , vh ) − b(vh , phn ) = (f, vh ) − S λhn gvhτ ds, ∀ vh ∈ Vh or Vh , ∀ qh ∈ Mh or Mh , b(uhn , qh ) = 0, (43) and ρ > 0, (44) λh(n+1) = PΛh (λ(hn) + ρguhnτ ), ˜ h → Λh is explicitly expressed as [16]: where the projection operator PΛh : Λ ⎧ if γh (Q) > 1, ⎨ +1 γh (Q) if |γh (Q)| ≤ 1, ∀Q ∈ S1,h (orS2,h ), PΛh (γh )(Q) = ⎩ −1 if γh (Q) < −1. ˜ h. for each γh ∈ Λ 5.1. 2D exact solution problem As is known, the problem (39) is a standard variational formulation of Navier-Stokes equations. We consider the problem (1)-(2) in the fixed square domain (0, 1) × (0, 1) (see Fig.1). The exact solution is given by u(x, y) = (u1 (x, y), u2 (x, y)), p(x, y) = (2x − 1)(2y − 1), u1 (x, y) = −x2 y(x − 1)(3y − 2), u2 (x, y) = xy 2 (y − 1)(3x − 2) 21

and f is defined by (1) with μ = 1/Re. It is easy to check the exact solution u satisfies u = 0 on Γ, u · n = u1 = 0, u2 = 0 on S1 and u1 = 0, u · n = u2 = 0 on S2 . Moreover, the tangential vector τ on S1 and S2 are (0, 1) and (−1, 0), thus,  στ = 4μy 2 (y − 1) on S1 , στ = 4μx2 (x − 1) on S2 . Moreover, from boundary conditions (2), we get |στ | ≤ g, thus we can chose the function g such that g = −στ ≥ 0 on S1 and S2 . Let the iteration initial value λh0 = 1 and the parameter ρ = μ. To access the theoretical results, numerical results obtained by P1 − P1 stabilized algorithm, P2 − P2 stabilized algorithm and two-step algorithms are compared. From Tables 1-5, we can see that two-step algorithms keep the convergence rates just like the theoretical analysis and the CPU-time of two-step algorithms is between the P1 − P1 and P2 − P2 stabilized algorithms. Especially, it is easy to find that the errors of two-step algorithms and P2 − P2 stabilized algorithm are the same older, but two-step algorithms are more efficient and cost less CPU-time. In addition, from table 6, we can know that Newton twostep algorithm is better than other algorithms in the precision for different μ. Thus we recommend algorithm 4.3 for actual use. h 1/8 1/16 1/24 1/32 1/40

u − uh V uV 0.31321 0.151205 0.0991352 0.0736784 0.0586082

Order / 1.0506 1.0412 1.0316 1.0255

p − ph  p 0.0552333 0.0205971 0.0113004 0.00734181 0.00524753

Order

CPU

/ 1.4231 1.4806 1.4991 1.5050

0.109 0.258 0.458 0.752 1.123

Table 1 Using the P1 − P1 stabilized finite element algorithm with μ = 1/Re = 0.1. h 1/8 1/16 1/24 1/32 1/40

u − uh V uV 0.0353693 0.00809836 0.00347295 0.001916 0.00121139

Order / 2.1268 2.0881 2.0674 2.0546

p − ph  p 0.0121835 0.00301414 0.00134455 7.5637e-04 4.78096e-04 22

Order

CPU

/ 0.189 2.0151 0.841 2.0462 1.642 2.0469 2.494 2.0557 3.523

Table 2 Using the P2 − P2 stabilized finite element algorithm with μ = 1/Re = 0.1. h 1/8 1/16 1/24 1/32 1/40

u − uh V uV 0.0353663 0.00809596 0.0034763 0.00191689 0.00121283

Order / 2.1271 2.0850 2.0692 2.0514

p − ph  p 0.0121281 0.00302714 0.00134149 7.52362e-04 4.79845e-04

Order

CPU

/ 0.165 2.0023 0.487 2.0071 0.978 2.0103 1.569 2.0155 2.643

Table 3 Using the algorithm 4.1 with μ = 1/Re = 0.1. h 1/8 1/16 1/24 1/32 1/40

u − uh V uV 0.0353663 0.00809596 0.0034763 0.00191689 0.00121283

Order / 2.1271 2.0850 2.0692 2.0514

p − ph  p 0.0121289 0.00302725 0.00134132 7.52336e-04 4.79825e-04

Order

CPU

/ 0.166 2.0024 0.491 2.0076 0.982 2.0099 1.574 2.0156 2.693

Table 4 Using the algorithm 4.2 with μ = 1/Re = 0.1. h 1/8 1/16 1/24 1/32 1/40

u − uh V uV 0.035322 0.00809442 0.00347209 0.0019157 0.00121135

Order / 2.1256 2.0875 2.0671 2.0540

p − ph  p 0.0121279 0.00302721 0.00134117 7.52325e-04 4.79811e-04

Order

CPU

/ 2.0023 2.0078 2.0096 2.0156

0.165 0.492 0.983 1.573 2.691

Table 5 Using the algorithm 4.3 with μ = 1/Re = 0.1. μ

Algorithm

0.1 0.1 0.1 0.005 0.005 0.005

4.1 4.2 4.3 4.1 4.2 4.3

u − uh V uV 0.00121283 0.00121283 0.00121135 0.00438431 0.00438415 0.00438373

p − ph  p 4.79845e-04 4.79825e-04 4.79811e-04 6.53021e-03 6.53007e-03 6.52991e-03

Table 6 Using algorithms 4.1-4.3 with h = 1/40. 23

1 0.8 0.6 0.4 0.2 0 1 1 0.8

0.5

0.6 0.4 0

0.2 0

(a) Figure 2: Domain Ω

5.2. 3D exact solution problem We consider the problem (1)-(2) in the fixed cube domain (0, 1)3 (see Fig.2), the boundary of which consists of two portions S and Γ given by S = (0, 1) × (0, 1) × {1}, and Γ = (0, 1) × {1} × (0, 1) ∪ {1} × (0, 1) × (0, 1) ∪ (0, 1) × (0, 1) × {0} ∪(0, 1) × {0} × (0, 1) ∪ {0} × (0, 1) × (0, 1). The exact solution is given by u(x, y, z) = (u1 (x, y, z), u2 (x, y, z), u3 (x, y, z)), u1 (x, y, z) = −(x3 − x2 )((3y 2 − 2y)(z 3 − z 2 ) − (y 3 − y 2 )(3z 2 − 2z)), u2 (x, y, z) = (y 3 − y 2 )((3x2 − 2x)(z 3 − z 2 ) − (x3 − x2 )(3z 2 − 2z)), u3 (x, y, z) = (z 3 − z 2 )((3y 2 − 2y)(x3 − x2 ) − (3x2 − 2y 2 )(y 3 − y 2 )), p(x, y, z) = (2x − 1)(2y − 1)(2z − 1), and f is defined by (1) with μ = 1/Re. By direct computation, we have max |στ | = S

max

0≤x≤1,0≤y≤1,

|μ(x3 − x2 )(4y 3 − 7y 2 + 2y)| = 24

8μ . 27

then we can chose the function g = 8μ on S. Let the iteration initial value 27 λh0 = 1 and the parameter ρ = μ. We display the CPU time and the convergence orders of the new method in Tables 7-9 for the 3D case. From these tables, we observe that the H 1 -error for the velocity and the L2 -error for the pressure have convergent property O(h2 ) as h → 0,in agreement with the theoretical analysis. Meanwhile, we have an interesting observation that the error of the velocity has better convergence rate than analytical results. h 1/8 1/9 1/10

u − uh V uV 0.16541 0.126453 0.098234

Order / 2.2801 2.3967

p − ph  p 0.0325045 0.0256832 0.0208029

Order

CPU

/ 84.625 1.9998 156.708 2.0002 274.639

Table 7 Using the algorithm 4.1 with μ = 1/Re = 0.01. h 1/8 1/9 1/10

u − uh V uV 0.16541 0.126450 0.098228

Order / 2.2803 2.3971

p − ph  p 0.0325041 0.0256825 0.0208021

Order

CPU

/ 85.774 1.9999 158.221 2.0003 276.335

Table 8 Using the algorithm 4.2 with μ = 1/Re = 0.01. h 1/8 1/9 1/10

u − uh V uV 0.16539 0.126441 0.098221

Order / 2.2799 2.3971

p − ph  p 0.0325038 0.0256821 0.0208018

Order

CPU

/ 86.643 2.0000 159.312 2.0003 278.001

Table 9 Using the algorithm 4.3 with μ = 1/Re = 0.01. 5.3. Lid-driven cavity flow In this example, we consider the lid-driven cavity flow on a unit square with friction boundary conditions only in the boundary {(x, 1) : 0 < x < 1} as follows. First, we consider the lid-driven cavity flow on a unit square with on-slip boundary conditions only in the boundary {(x, 1) : 0 < x < 1} with the velocity u = (1, 0) in the fairly finer mesh width as exact solution and we obtain g = |στ |. Then, the friction boundary conditions g = |στ | is imposed 25

IsoValue -0.698299 -0.585212 -0.472125 -0.359038 -0.245952 -0.132865 -0.0197782 0.0933086 0.206395 0.319482 0.432569 0.545655 0.658742 0.771829 0.884916 0.998002 1.11109 1.22418 1.33726 1.45035

Vec Value 0 0.0526343 0.105269 0.157903 0.210537 0.263172 0.315806 0.36844 0.421075 0.473709 0.526343 0.578978 0.631612 0.684247 0.736881 0.789515 0.84215 0.894784 0.947418 1.00005

(a)

(b)

Figure 3: The pressure level lines and the velocity field with P2 − P2 . IsoValue -0.698299 -0.585212 -0.472125 -0.359038 -0.245952 -0.132865 -0.0197782 0.0933086 0.206395 0.319482 0.432569 0.545655 0.658742 0.771829 0.884916 0.998002 1.11109 1.22418 1.33726 1.45035

Vec Value 0 0.0526343 0.105269 0.157903 0.210537 0.263172 0.315806 0.36844 0.421075 0.473709 0.526343 0.578978 0.631612 0.684247 0.736881 0.789515 0.84215 0.894784 0.947418 1.00005

(a)

(b)

Figure 4: The pressure level lines and the velocity field with algorithm 4.3.

on the top side {(x, 1) : 0 < x < 1} and zero Dirichlet condition is imposed on the rest of the boundary. Comparison is made between the P2 − P2 stabilized method and the algorithm 4.3, with Re = 200 and h = 1/25 for the lid-driven cavity flow. In Figures 3-4, we show the pressure level lines and the velocity field of cavity flow. However, the CPU-time of the algorithm 4.3 is 3.674(s), while the one of the P2 − P2 stabilized method is 14.398(s), which is nearly three times. Thus our algorithm can take less time to get the same result as the P2 − P2 stabilized method.

26

Conclusions In this paper, we have provided theoretical analysis with two-step algorithms for the incompressible flows with friction boundary conditions. Error estimations of these numerical algorithms are proved. Numerical examples show that our two-step algorithms are valid for the incompressible NavierStokes equation with friction boundary conditions, and the numerical results are consistent with the theoretical analysis. [1] P. Bochev, M. Gunzburger, An absolutely stable pressure-poisson stabilized finite element method for the Stokes equations, SIAM J. Numer. Anal. 42 (2004) 1189-1207. [2] P. Bochev, C. Dohrmann, M. Gunzburger, Stabilization of low-order mixed finite elements for the Stokes equations, SIAM J. Numer. Anal. 44 (2006) 82-101. [3] J. Djoko, On the time approximation of the Stokes equations with nonlinear slip boundary conditions, Int. J. Numer. Anal. Model. 11 (1) (2014) 34-53. [4] J. Djoko, J. Kokob, Numerical methods for the Stokes and Navier-Stokes equations driven by threshold slip boundary conditions, Comput. Methods Appl. Mech. Engrg. 305 (2016) 936-958. [5] V. Ervin, W. Layton, J. Maubach, A posteriori error estimators for a two-level finite element method for the Navier-Stokes equations, Numer. Methods Part. Diff. Eq. 12 (1996) 333-346. [6] L. Evans, Partial Differential Equations, American Mathematical Society, 1998. [7] H. Fujita, Flow Problems with Unilateral Boundary Conditions, Lecons, Coll`ege de France 1993. [8] H. Fujita, A mathematical analysis of motions of viscous incompressible fluid under leak or slip boundary conditions, RIMS Kokyuroku 888 (1994) 199-216. [9] V. Girault, P. A. Raviart, Finite Element Method for Navier-Stokes Equations: Theory and Algorithms, Springer-Verlag, Berlin, Herdelberg, 1986. 27

[10] Y. He, K. Li, Two-level stabilized finite element methods for the steady Navier-Stokes problem, Computing 74 (2005) 337-351. [11] Y. He, A. Wang, A simplified two-level method for the steady NavierStokes equations, Comput. Methods Appl. Mech. Eng. 197 (2008) 15681576. [12] Y. He, J. Li, A stabilized finite element method based on local polynomial pressure projection for the stationary Navier-Stokes equations, Appl. Numer. Math. 58 (2008) 1503-1514. [13] P. Huang, X. Feng, Y. He, An efficient two-step algorithm for the incompressible flow problem, Adv. Comput. Math. 41 (2015) 1059-1077. [14] N. Kikuchi, J. T. Oden, Contact Problems in Elasticity, SIAM, Philadelphia, 1988. [15] T. Kashiwabara, On a strong solution of the non-stationary NavierStokes equations under slip or leak boundary conditions of friction type, J. Differential Equations 254 (2013) 756-778. [16] T. Kashiwabara, On a finite element approximation of the Stokes equations under a slip boundary condition of the friction type, Japan J. Indust. Appl. Math. 30 (2013) 227-261. [17] T. Kashiwabara, Finite element method for Stokes equations under leak boundary condition of friction type, SIAM J. Numer. Anal. 51 (2013) 2448-2469. [18] W. Layton, W. Lenferink, Two-level picard and modified picard methods for the Navier-Stokes equations, Appl. Math. Comput. 69 (1995) 263274. [19] W. Layton, L. Tobiska, A two-level method with backtracking for the Navier-Stokes equations, SIAM J. Numer. Anal. 35 (1998) 2035-2054. [20] W. Layton, Introduction to Finite Element Methods for Incompressible Viscous Flows, SIAM, Philadelphia, 2008. [21] J. Li, Y. He, Z. Chen, A new stabilized finite element method for the transient Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 197 (2007) 22-35. 28

[22] J. Li, Y. He, A stabilized finite element method based on two local Gauss integrations for the Stokes equations, J. Comp. Appl. Math. 214 (2008) 58-65. [23] Y. Li, K. Li, Existence of the solution to stationary Navier-Stokes equations with nonlinear slip boundary conditions, J. Math. Anal. Appl. 381 (2011) 1-9. [24] Y. Li, K. Li, Pressure projection stabilized finite element method for Navier-Stokes equations with nonlinear slip boundary conditions, Computing 87 (2010) 113-133. [25] Y. Li, R. An, Semi-discrete stabilized finite element methods for NavierStokes equations with nonlinear slip boundary conditions based on regularization procedure, Numer. Math. 117 (1) (2011) 1-36. [26] Y. Li, K. Li, Uzawa iteration method for Stokes type variational inequality of the second kind, Acta Math. Appl. Sin-E 27 (2011) 302-315. [27] Y. Li, K. Li, Global strong solutions of two-dimensional Navier-Stokes equations with nonlinear slip boundary conditions, J. Math. Anal. Appl. 393 (2012) 1-13. [28] Y. Li, R. An, Penalty finite element method for Navier-Stokes equations with nonlinear slip boundary conditions, Int. J. Numer. Meth. Fluids 69 (2012) 550-566. [29] H. Qiu, L. Mei, H. Liu, S. Cartwright, Defect-correction stabilized finite element method for Navier-Stokes equations with friction boundary conditions, Appl. Numer. Math. 90 (2015) 9-21. [30] H. Qiu, L. Mei, Two-level defect-correction stabilized finite element method for Navier-Stokes equations with friction boundary conditions, J. Comp. Appl. Math. 280 (2015) 80-93. [31] H. Qiu, L. Mei, Y. Zhang, A two-grid stabilized finite element method for Navier-Stokes equations with friction boundary conditions, Submitted (2016). [32] N. Saito, On the Stokes equations with the leak and slip boundary conditions of friction type: regularity of solutions, Publ. RIMS, Kyoto Univ. 40 (2004) 345-383. 29

[33] R. Temma, Navier-Stokes Equation: Theory and Numerical Analysis, North-Holland, Ansterdam, New York, Oxford, 1984. [34] H. Zheng, L. Shan, Y. Hou, A quadratic equal-order stabilized method for Stokes problem based on two local Gauss integrations, Numer. Meth. Part. Differ. Eq. 26 (2010) 1180-1190.

30