A Subgrid Stabilized Method for Navier-Stokes Equations with Nonlinear Slip Boundary Conditions

. In this paper, we consider a subgrid stabilized Oseen iterative method for the Navier-Stokes equations with nonlinear slip boundary conditions and high Reynolds number. We provide one-level and two-level schemes based on this stability algorithm. The two-level schemes involve solving a subgrid stabilized nonlinear coarse mesh inequality system by applying m Oseen iterations, and a standard one-step Newton linearization problems without stabilization on the ﬁne mesh. We analyze the stability of the proposed algorithm and provide error estimates and parameter scalings. Numerical examples are given to conﬁrm our theoretical ﬁndings.


Introduction
We consider the approximations for steady, incompressible Navier-Stokes equations: −ν∆u + (u · ∇)u + ∇p = f in Ω, divu = 0 in Ω, (1.1) as Re := U L/ν, where U is a characteristic velocity and L is a characteristic length. The equations (1.1) are completely specified with an appropriate boundary condition for the velocity u. In this paper, we consider the following nonlinear slip boundary conditions: u = 0 on Γ 0 , u n = 0, |σ τ (u)| ≤ g on Γ 1 , (1.2) where ∂Ω = Γ 0 ∪ Γ 1 . We set the usual Dirichlet boundary condition on Γ 0 and slip, non-leak boundary condition on Γ 1 , respectively. For simplicity, we suppose meas(Γ 0 ) = 0, meas(Γ 1 ) = 0 and Γ 0 ∩ Γ 1 = ∅. Here and what follows, the unit outer normal of the boundary is expressed by n. If a is a vector defined on the boundary, then a n = a · n is the normal component of a, and a τ = a · τ is the tangential component of a; σ τ (u) = ν ∂uτ ∂n is the tangential component of stress vector σ; g is a given positive function on Γ 1 . These boundary conditions can be found in [8,9,10] where some problems in hydrodynamics are investigated.
We notice that the velocity-pressure variational formulation (1.1)-(1.2) is the variational inequality of the second kind. Some numerical results have been obtained for these variational inequalities [12,13,19,20,21,28]. There are two main challenges in solving this problem. First, the velocity u and pressure p are coupled by the incompressible constraint divu = 0, which brings a saddle problem [2,5,6]. A popular approach to solve this situation is to relax the incompressible constraint to obtain a pseudo-compressible system. There are many ways of implementing this technique, such as penalty method [4,17], pressure stabilization method [14] and the projection method [11]. Second, the problem involves nonlinear terms, including (u · ∇)u and slip boundary conditions. For (u · ∇)u, we usually use linearization method such as Newton method [3], twolevel method [16,20] and other methods. In [18], the Uzawa iteration is used to handle the slip boundary conditions, which is also transformed into the saddle point problem by using the augmented Lagrangian approach [5].
The main idea of the two-level method is first to solve the fully nonlinear problem on a coarse grid, and then solve a linear problem on fine grid to correct the solution. Both theoretical analysis and numerical examples show that this method is very efficient [16,20]. However, the nonlinear system may fail to converge on the coarse mesh when we use standard two-level method to solve high Reynolds number problems [7,22]. Therefore, the stabilization strategies need to be combined.
There are many stabilization methods for solving high Reynolds number problems, including the subgrid stabilization methods, the defect-correct methods, the variational multiscale methods and so on. The subgrid method first divides the approximation space into resolved scales and subgrid scales and then modifies the Galerkin approximation, using an asymptotically consistent artificial diffusion term on the subgrid scales [15,22,26]. The defect-correct method achieves stabilization by several correction steps after the initial defect step [20,21,25]. The key idea of the variational multiscale method is to decompose the flow scale and design the large scale by projection into approximate subspaces [27]. Shang [22] proposed a two-level subgrid stabilized Oseen iter-ative method for steady Naiver-Stokes. In this method, a subgrid scale model based on an elliptic projection of the velocity is utilized to stabilize the numerical form of the incompressible Navier-Stokes equations on a coarse grid, and then a linear problem that the convection term is fixed by the coarse grid solution which is solved on the fine grid. Both theoretical analysis and numerical examples show that the method can efficiently simulate the high Reynolds number flows.
In this paper, based on the one-level scheme we develop a two-level subgrid stabilized Oseen iterative method for the Navier-Stokes equations with nonlinear slip boundary conditions. The first step solves a small variational inequality based on an elliptic projection of the velocity on the coarse mesh. The second step solves a large linearization correction problem with nonlinear slip boundary conditions on the fine mesh in term of Newton iteration. We obtain the following error estimate for two-level solution (u h , p h ): where H is the coarse grid, h is the fine grid, and θ is the stabilization parameter. Compared to the well-known two-level methods, our method can simulate high Reynolds number flows with nonlinear slip boundary conditions. Compared to the two-level defect-correction stabilization method [20], there are two main differences between them. Firstly, our method only applies the subgrid stabilization method on coarse grid, while the method [20] applies the artificial viscosity stabilized method on both coarse and fine grids. Secondly, our method uses an elliptic projection to establish the subgrid stabilization model, while the method [20] adopts L 2 -projection and an artificial viscosity term to achieve stabilization. Compared to the two-level subgrid stabilized Oseen iterative method for steady Naiver-Stokes [22], our method can work well for flows with nonlinear slip boundary conditions. The rest of the paper is organized as follows. In Section 2, we will introduce some notations and recall some results for Navier-Stokes equations with nonlinear slip boundary conditions. In Section 3, we give the concrete schemes of one-level subgrid stabilized Oseen iterative method based on an elliptic projection of the velocity for Navier-Stokes equations and obtain the stability and error estimates. Two-level subgrid stabilized Oseen iterative method is presented and analyzed in Section 4. The numerical results are presented in Section 5. Finally, we give some conclusions in the last section.

Preliminaries
We first introduce some notations in this section. Let H 1 0 (Ω) be the standard Sobolev space (see [1]) equipped with the usual norm · 1 . We use | · | 1 for its semi-norm. Let V , Q and V 0 be defined by The scalar product and norm in Q are denoted by the usual L 2 (Ω) inner product and · respectively. Next we introduce the bilinear forms: and for the treatment of the convective term in Equation (1.1), the following trilinear form is also considered: This form is well defined and continuous on these spaces(see [24]). The following estimates will be used repeatedly in the sequel, see [23,24] and the references therein: Based on the above expression, the velocity-pressure variational formulation for (1.1)-(1.2) can be converted to the variational inequality of the second kind: where j(v) = Γ1 g|v|ds. In order to guarantee the well-posedness of the problems (2.3), we should assume the LBB condition holds: Then the existence and uniqueness theorem of the solution to the problem (2.3) is as follows(see [17]): then the problem (2.3) admits a unique solution u such that

4)
Remark 1. From the condition of Theorem 1, we can see that when the Reynolds number is high, i.e. ν is very small, it will bring some difficulties to the solution. However, as long as the body force f , the function g and boundary Γ 1 are properly adjusted, the existence of the solution can be guaranteed. In our last numerical example, when the Reynolds number is high, we indeed encounter non-convergence. This non-convergence may be due to the need for stabilization, but the other may be caused by the nonexistence of the solution. However, we adjust the function g and Γ 1 , and finally solve the problem successfully.
To present the subgrid stabilization method, we first consider the usual finite element approximation to (1.1)-(1.2), which can be described as follows. Assume that V µ and Q µ (here µ = H, h with H > h) are the finite element subspaces of V and Q. We make some assumptions. For each (u, p) ∈ V × Q, there exists an approximation (π µ u, ρ µ p) ∈ V µ × Q µ such that We shall always assume that the following discrete LBB condition holds: Then the usual finite element approximation (u µ , p µ ) is calculated by solving the nonlinear system:

One-level subgrid stabilization method
In this section, we give the concrete schemes of one-level subgrid stabilized Oseen iterative method based on an elliptic projection of the velocity for Navier-Stokes equations with nonlinear slip boundary. We assume that T µ (Ω) = {K} is a shape-regular triangulation of Ω with the mesh size µ. Then we introduce the elliptic projection operator Π µ : [26]: and has the following property here P 1 is the space of polynomials of degree less than or equal to one. Based on this projection operator, we can define the subgrid stabilization term as where θ > 0 is a user-defined stabilization parameter whose value will be mentioned in the following. By combining the stabilization term (3.2) and applying the Oseen iterative method to deal with nonlinear term, we give the following one-level subgrid stabilization method for Navier-Stokes equations with nonlinear slip boundary: for n = 1, 2, · · · , and u 0 µ = 0. Next we give the stability of one-level subgrid method (3.3).
then, we get Setting v µ = 0, q µ = p n µ in (3.3), using the In what follows, c denotes a generic constant, possibly different at different occurrences, which may depend on the data f, g, ν, the domain Ω and the continuous solution u, but is independent of the mesh size µ (µ = H, h) and stabilization parameters θ. The next theorem provides the error estimate for the solution produced by one-level subgrid stabilized Oseen iterative method. Theorem 3. Let (u, p) ∈ V × Q and (u n µ , p n µ ) ∈ V µ × Q µ be the solutions of (2.3) and (3.3), and let the LBB condition (2.6) hold, then where c > 0 is independent of θ and µ.
Proof. Set v = u n µ and v = 2u − v µ in the first equation of (2.3), respectively. We have Summing of the two aforementioned inequalities, we get Setting (e n , η n ) = (u − u n µ , p − p n µ ) and from (2.3), (3.3) and (3.6), we obtain the following relation: Denote the u I , p I as the projections onto space V µ , Q µ , respectively. Then, we get Taking into account the property of trilinear form (2.2) and the boundedness of u (2.4)and u n µ (3.4), we obtain For the stabilization term, considering the (3.1) we have In addition, we have Then, considering we can easily get LetṼ µ be the finite element subset of V 0 . Set v = u + w µ in the first equation of (2.3), where w µ ∈Ṽ µ , which yields Similarly, from the first equation of (3.3), we can obtain (ν+θ)(∇u n µ , ∇w µ )+b(u n−1 µ , u n µ , w µ )+d(w µ , p n µ )=(f, w µ ) + θ(∇Π µ u n−1 µ , ∇w µ ).
Therefore, we get Taking into account the establishment of LBB condition (2.6) and inequality (2.2), we get the following inequalities: where Taking (3.8) into (3.7) gives Then we estimate the right-hand side using the triangle and Young's inequalities: At last, using the inequality |u − u n µ | 1 ≤ |u − u I | 1 + |u I − u n µ | 1 and omitting the high-order term, we get the error estimate for velocity For pressure term, it is easy to get the following formula when we take into account the inequality p − p n µ ≤ p − p I + p I − p n µ and (3.8) which concludes the proof of Theorem 3.
From the Theorem 3, we can see that one-level subgrid method can get optimal solutions if the stabilization parameter θ is small enough and the iteration number n tends to infinity. Of course, the stability parameters θ can not be taken to be too small. Otherwise, it will not have a stabilization effect. In fact, we can choose θ = O(µ) to get optimal solutions. As for the Oseen iteration error, we only need to impose a certain stop condition to control the iteration number n. From numerical examples in Section 5, we can see that the one-level subgrid stabilized method solves the high Reynolds number problems very well. But at the same time, we also see that the computing time increases rapidly when mesh grid µ tends to zero. In Section 4, we try to improve computational efficiency while maintaining the same level of accuracy.

Two-level subgrid stabilization method
In order to improve computational efficiency, we combine two-grid discretization strategy with the subgrid stabilized Oseen iterative method. We develop our two-level subgrid stabilized Oseen iterative method as follows: for n = 1, 2, · · · , m, and u 0 H = 0.

Find a find grid solution
Stabilization measures are only used for the coarse grid in the first step, while the fine grid linear problem is a standard one-step Newton linearization. The second step can also be based on Oseen iteration. But Newton's method is more suitable as we have had a relatively accurate solution. The next Theorem provides the error estimate of the solution produced by the two-level subgrid stabilized Oseen iterative method. Theorem 4. Let (u, p) ∈ V × Q and (u h , p h ) ∈ V h × Q h be the solutions of (2.3) and (4.1), and let the LBB condition (2.6) hold, then where c > 0 is independent of H, h and θ.
Proof. Set v µ = v h and u n µ = u h in (3.6), added by (4.1), we obtain where we set (e h , η h ) = (u − u h , p − p h ).
Denote the u I h , p I h as the projections onto space V h , Q h , respectively. Then we have . Similar to the proof of Theorem 3, we can get where we use the LBB condition (2.6). Considering the inequality For pressure terms, we can easily get the conclusion by using the triangle inequality |η h | 1 ≤ |p − p I h | 1 + |p I h − p h | 1 and the formula (4.4).
A direct consequence of Theorem 4 is that This follows by inserting the estimate (3.5) for the coarse grid solution and the approximation assumption (2.5) of the finite element spaces into the right-hand of estimate (4.2). The last term related to the Oseen iterations on the coarse grid indicates the influence of the nonlinear term. It shows that the nonlinear solver should be taken to be small enough to obtain a good approximate solution. Moreover, (4.5) leads to the typical choices of parameters: which ensures the optimal order of the solution.
The meshes which consist of triangular elements and the P 2 − P 1 finite element pair are used for discretization. In all experiments, the projection Π µ u n−1 µ (µ = H, h) is computed by solving Laplace equations with piecewise linear finite element functions. The iterative tolerance for Oseen iterations is 10 −6 . The maximum number of Oseen iterations is 5000. When the number of iterations exceeds 5000, we declare that the method is not convergent. We first consider a simple problem with known analytical solutions. Let Ω = (0, 1) × (0, 1), Γ 1 = S 1 S 2 (see Figure 1) and exact solution (u, p) of Navier-Stokes equations (5.1) be Then the body force f is given by (1.1). By simple calculation we can know that Then the position function g can be chosen by g = σ τ . Let the initial value λ 0 = 1 and the parameter ρ = ν.
First, we investigate the effect of stabilization parameter θ on the errors of one-level subgrid method. For the case 1, we fix H = 1/64, ν = 0.05 and take parameter θ for five different values. For the case 2, we fix H = 1/64, ν = 0.005 and take parameter θ for six different values. m is the Oseen iterations count satisfying the stopping condition. Table 1 shows errors and numbers of Oseen iterations, from which we can see that the error and number become smaller when the parameter θ gets smaller. This means in order to improve the accuracy and save time, the parameter values should be taken small. On the other hand, we know that the parameters can't be very small because of the instability of the high Reynolds number problems. Therefore, combining these two cases, we intend to take the value of the parameter as 0.01H in the following examples. Next we show the errors of one-level subgrid stabilized Oseen iterative method with mesh size h = 1/49, 1/64, . . . , 1/144 and θ = 0.01h in Table 2, from which we can see the convergence rates are consistent with the estimate (3.5). Table 3 shows the errors of the two-level subgrid stabilized Oseen iterative method, where the fine mesh size is set as h = 1/49, 1/64, . . . , 1/144, corresponding H satisfying H = h 1/2 and the stabilization parameter is set as θ = 0.01H. Comparing Table 2 and Table 3, we find that two-level subgrid stabilized method takes much less computational time while accuracies of the two methods are comparable to each other.
Then we study the necessity of the stabilization in our two-level method for high Reynolds number problems. We set ν = 0.01, θ = 0.0005, H = 1/16, h = 1/256 and θ = 0.02h, and then we use our two-level subgrid stabilization method and the standard two-level method based on Newton iteration to solve the problem. The experimental results show that our two-level stabilization method works well, but the standard two-level method can not obtain the solution of equations since no convergence solution is obtained on the coarse grid. Figure 2 shows the computed streamlines and pressure contours, which confirms the necessity of subgrid stabilization. Second, we consider the 2D lid-driven cavity flow problem which is defined in the unit square with nonlinear slip boundary conditions only in the boundary {0 < x < 1, y = 1}. This is also investigated in [20]. The Reynolds number for this problem can be defined as Re = U L/ν, where U is a characteristic velocity and L is a characteristic length. We need to know the position function g in this problem. Therefore, we consider the lid-driven cavity flow on a unit square with no-slip conditions only in the boundary {0 < x < 1, y = 1} imposing the velocity u = (1, 0). We solve this problem on the fairly finer mesh(h = 1/257) to get a solution as the exact solution and we get g = σ τ . Then, the nonlinear slip conditions g = σ τ is imposed on the top one while zero Dirichlet condition is imposed on the rest of boundary.
We compute approximate solutions for the cases at Re=1000, 5000 and 10000 for the lid-driven cavity flow. Table 4 shows the numbers of Oseen iterations by one-level subgrid stabilized method with H = 1/64, θ = H, 0.5H, . . . , 0.01H. We can see from Table 4 at Re=1000, the iterative method converges for all values of stabilized parameter, while at Re=5000 and Re=10000, the method fails to converge for small values. Figures 3-5 depict the numerical streamlines and pressure contours computed by the two-level subgrid stabilized method with H = 1/64, h = 1/128 and θ = 0.1H. For this problem, we also tested two other larger Reynolds numbers, Re=30000 and Re=50000. Here, we calculate g = σ τ at grid h = 1/513. As mentioned above, the small size of the coarse mesh will lead to non convergence, so we still choose H = 1/64. Figure 6 shows the numerical result with h = 1/256 and θ = 0.08H. We find that with the increase of Reynolds number, in order to ensure the convergence of the algorithm, it will be more difficult to find a suitable value of θ. As mentioned in [9], the Navier-Stokes equations with nonlinear slip boundary conditions can simulate the blood flow of patients with arteriosclerosis. So in the last example, we simulate a simplified blood flow problem of this type. It should be noted that the parameters we set for the convenience of calcula- tion are still somewhat different from the real blood flow parameters, but they can basically explain some facts. We first introduce Dirichlet boundaries with parabolic shape normal to the fixed boundary, which is given by where l is the length of the boundary part where the flow velocity is presented andm is the magnitude of the flow velocity at the centre of this boundary part. The detailed boundary conditions and regional distribution are shown in Figure 7 which in our example represents a bifurcated blood vessel. The inlet is set to Dirichlet boundary conditions with the maximum flow velocitym = 1 and the length l = 1.2, and the outlet is set to open boundary, i.e.g 1 = 0, where g 1 = (ν∇u − pI) · n. A part of the top of the blood vessel is set to nonlinear slip boundary conditions, and the remaining boundary is set to zero Dirichlet condition.
Like the previous example, we also need to know the position function g. Similarly, we also calculate a problem on a relatively finer grid mesh(h = 1/200) with no-slip conditions to obtain a solution as the exact solution and obtain g = σ τ . We compute approximate solutions for the cases at ν = 0.1 and ν = 0.001 for the blood flow with H = 1/64, h = 1/128 and θ = 0.1H. The influence of the stable term is hardly reflected when ν = 0.1. In the twolevel algorithm, regardless of whether there is a stable term, the number of Oseen iterations on the coarse mesh is m = 6. But for the case of ν = 0.001, stabilization is necessary. If there is no stabilizing term, then the method will not converge on the coarse mesh, which will lead to the failure of calculation. When we add stabilizing term and set θ = 0.1H, the number of Oseen iterations on the coarse mesh is m = 99. Figures 8-9 depict the numerical streamlines and pressure contours computed by the two-level subgrid stabilized method.

Conclusions
In this paper, we propose a two-level subgrid stabilized method for Navier-Stokes equations with nonlinear slip boundary conditions and high Reynolds number. In this way, we use an elliptic projection to stabilize problems only on the coarse grid. Stability and error estimates of the method were analyzed. Numerical results verified the theoretical results and demonstrated the efficiency of the method. We will go on working on the highly efficient methods of numerical solutions to Navier-Stokes equations.