Djoko Kamdem
email: jules.djokokamdem@gmail.com

Jonas Koko
email: jonas.koko@uca.fr

Mohamed Mbehou
email: mbehoumoh@gmail.com

Toni Sayah
email: toni.sayah@usj.edu.lb

J K Djoko

Stokes and Navier-Stokes equations under power law slip boundary condition: Numerical Analysis

Keywords: power law slip boundary condition, Stokes equations, Navier-Stokes equations, finite element method, monotonicity, error estimates, rate of convergence

Introduction

Let Ω ⊂ R d with d = 2, 3 a bounded domain occupied by a fluid with a boundary ∂Ω. We would like to study the incompressible flows governed by the Stokes equations

-2ν div Du + ∇p = f in Ω, (1.1) div u = 0 in Ω , (1.2)
where; ν is the positive viscosity of the fluid velocity, u(x), the velocity, the pressure is p(x) and f is the external force acting on the fluid, Du is the symmetric part of the velocity gradient defined through 2Du = ∇u+(∇u) T . The equations (1.1) and (1.2) are supplemented by boundary conditions that will be described next. We are interested in (1.1), (1.2) when the position and the direction of the slip boundary condition are taken into account (see [START_REF] Djoko | Flows of incompressible viscous liquids with anisotropic wall slip[END_REF][START_REF] Roux | On the Navier Stokes equations with anisotrpic wall slip conditions[END_REF]). We then assume that the boundary ∂Ω of Ω is made of two components S and Γ, such that ∂Ω = S ∪ Γ, with S ∩ Γ = ∅. We assume the homogeneous Dirichlet condition on S, that is u = 0 on S .

(1.3)

Thus S is the porous or artificial boundary where the fluid is prescribed. On Γ, we assume the impermeability condition u • n = 0 on Γ , (1.4) where n : Γ -→ R d is the normal outward unit vector to Γ. Γ is an impermeable solid surface along which the fluid may slip. The force within the fluid is the Cauchy stress tensor T given through T = 2νDu -pI on Ω , with I the d-dimensional identity matrix. Taking the scalar product of u and (1.1), one has

Ω T : Dudx + Γ (-T n) τ • u τ dσ = Ω f • udx (1.5)
with dσ being the surface measure associated to Γ. We recall that for any vector w define on ∂Ω, we set w τ = w -(w • n)n. Thus (T n) τ denotes the projection of the normal stress into the corresponding tangent plane. Now, from a simple interpretation, (1.5) expresses the balance of energy in the system. (-T n) τ • u τ dσ stands for dissipative processes on the boundary. We need a functional relation between (T n) τ and u τ that can render the dissipative energy on Γ non negative. For that purpose, the most general relation between u τ and (T n) τ is the implicit constitutive relation [START_REF]On the implicit constitutive theories[END_REF] ψ(u τ , (T n) τ ) = 0 (1.6) where ψ is function. The simplest form of (1.6) that ensure the non-negativity of

Γ (-T n) τ •u τ dσ
is the choice (-T n) τ = αu τ dσ , α > 0 .

This is the Navier's slip boundary conditions. If (T n) τ = 0, then one gets a perfect slip boundary condition, while if u τ = 0, then there is no slip. But in this work we are interested in power law slip boundary condition given as follows

(T n) τ + |Ku τ | s-2 K 2 u τ = 0 on Γ , (1.7)
where |v| 2 = v • v is the Euclidean norm. K is an anisotropic tensor, assumed to be uniformly positive definite, symmetric, and bounded. s is a real, strictly positive number representing the flow behaviour index. The tangential shear is a power law function of the tangential velocity. Such a boundary condition arises when the contact surface is lubricated with a thin layer of a non-Newtonian fluid. It is manifest that for s = 2 and K = I, (1.7) is reduced to the classical Navier's slip condition. In this work we will assume that s = 2. We are also interested in the Navier-Stokes equations, where (1.1) is replaced by with g a tangential surface traction and H a bounded, uniformly Lipschitz continuous function are examined by C. Le Roux in [START_REF] Djoko | Flows of incompressible viscous liquids with anisotropic wall slip[END_REF][START_REF] Roux | On the Navier Stokes equations with anisotrpic wall slip conditions[END_REF]. The properties on the function H had helped Le Roux to control the H 1 -norm of the velocity with only the slip boundary condition u • n = 0 on the entire boundary because the Korn's inequality in that context reads (see [START_REF] Ciarlet | Three dimensional elasticity[END_REF])

-2ν div D(u) + (u • ∇)u + ∇p = f in Ω, ( 1
Ω |Dv| 2 dx + ∂Ω |v τ | 2 dσ ≥ c Ω |∇v| 2 dx .
The situation in our hands is different. Indeed, for s = 2, the function u → Γ |Ku τ | s-2 dσ does not have the properties enumerated for the function H. Hence the results obtained in [START_REF] Djoko | Flows of incompressible viscous liquids with anisotropic wall slip[END_REF][START_REF] Roux | On the Navier Stokes equations with anisotrpic wall slip conditions[END_REF] are not applicable here. It is therefore clear that the problems described need to be analysed. First, one wishes to find a way to describe/construct the solution(s) of the resulting boundary value problems, and secondly, how can we compute the solutions of the resulting boundary value problems. Our interest to study this type of boundary conditions is mainly motivated by the desire of extending some results known when classical boundary conditions are used. Thus this work also participates in a better mathematical comprehension of fluids with non classical boundary conditions. From the modelling perspective, flows of incompressible Newtonian/Non Newtonian fluids are most of the times analysed with Dirichlet type (or periodic ) boundary conditions on the velocity field. But, in some important situations such as oil ducts, water supply, microfluids channels or blood in arteries, some formulations with boundary conditions involving components of the velocity field, stress or pressure are of interest. A recent review of some of the formulations presented in the literature, and their associated boundary conditions, with a focus on applications to air and blood flows can be found in [START_REF] Fouchet-Incaux | Arificial boundaries and formulations for the incompressible Navier Stokes equations: applications to air and blood flows[END_REF]. It is important to observe that most papers dealing with the Stokes or Navier-Stokes equations are written by means of conservation of the momentum in terms of the Laplacian of the velocity, whereas in this text we prefer the equivalent formulation given in terms of the divergence of the symmetric gradient which is useful in the fluid-structure interaction problems. This formulation has many important implications. At the continuous level, the wellposedness of the problems are proved in the same functional spaces as for the standard boundary conditions. Moreover, the discrete analysis allows us to obtain convergence (and even optimal one) with standard inf-sup stable finite element spaces without requiring stabilized formulations. In this work, one also observes that the mathematical and computational models share this important feature: the fluid model is described in terms of mechanical stress tensor, which is more appropriate from the modelling viewpoint, and it is able to properly take into account the fact that the flow continues beyond the boundaries, thanks to the specific form of the boundary conditions. If one considers that ∂Ω is not partitioned, with (1.7) and u • n = 0 on ∂Ω, then one deduces that

Ω |Du| 2 dx + ∂Ω |Ku τ | s dσ = Ω f • udx .
Hence we will have some difficulties to control the velocity with H 1 -norm. To overcome this deficiency, we have assumed that the boundary ∂Ω is divided into a porous or artificial boundary where the fluid is prescribed, and the other part of the boundary is impermeable which may allow the fluid to slip. Ultimately, we would like to address the problem when the distance between Γ and S is zero, but this would introduce complications that go beyond the scope of this study. The goal of this work is to analyse from the mathematical and numerical point of view both (1.1)-(1.4), (1.7) and (1.2)-(1.4), (1.7) and (1.8). Firstly, for both problems, the existence of solution is constructed without restriction on the data by using monotone theory for the problem related to the Stokes equations, while for the problem associated with the Navier-Stokes equations, we use Galerkin's approximation, Brouwer's fixed point and some compactness results. Secondly, both variational problems are discretized by finite element schemes in a polygonal or polyhedral domain. We derive existence, uniqueness (conditional for the Navier-Stokes problem), and a priori error estimates for the respective problems. Thirdly, for the implementation of the nonlinear finite element problems, we formulate a solution technique based on the linearization of nonlinear terms. Next, we studied its feasibility, and the convergence of the iterative scheme associated with the problem with the Navier-Stokes equations. Finally, we present some numerical experiments that confirm the theoretical findings. We shall mentioned at at this juncture that, we decide to focus on the implementation of the nonlinear problem related to the Navier-Stokes equations because it has more difficulties/nonlinearities than the problem related to the Stokes equations. We note that partial differential equations with the nonlinear boundary condition (1.7) (without (1.4)) have been used in the past by [START_REF] Bodart | A fictitious domain decomposition method for a nonlinear bonded structure[END_REF][START_REF] Geymonat | A domain decomposition method for bonded structures[END_REF][START_REF] Koko | An optimization based domain decomposition method for a bonded structure[END_REF][START_REF] Koko | Convergence analysis of optimization-based domain decomposition methods for a bonded structure[END_REF] when studying bonded elastic materials in solid mechanics. A nonlinear boundary condition similar to power law slip boundary condition (1.7) were studied by [START_REF] Bresch | An optimization-based domain decomposition method for nonlinear wall laws in coupled systems[END_REF][START_REF] Koko | Uzawa Conjugate Gradient Domain Decomposition Methods for Coupled Stokes Flows[END_REF] for two coupled Stokes fluids with the nonlinear boundary condition representing the transmission condition. The interface condition in the later representing many situations; adhesion, friction, coupling of free-surface and exchange with ground water flows.

The study of flows whose motion is described by the Stokes or Navier-Stokes equations driven by anisotropic slip boundary condition have been the subject of some publications recently (see [START_REF] Djoko | Flows of incompressible viscous liquids with anisotropic wall slip[END_REF][START_REF] Roux | On the Navier Stokes equations with anisotrpic wall slip conditions[END_REF][START_REF] Charrault | Interfacial slip on rough, patterned and soft surfaces: a review of experiments and simulations[END_REF][START_REF] Lu | Characterization on the anisotropic slip for flows over unidirectional fibrous porous media for advanced composites manufacturing[END_REF][START_REF] Busse | Influence of an anisotropic slip-length boundary condition on turbulent channel flow[END_REF][START_REF] Najeeb Alam Khan | Effect of anisotropic slip and magnetic field on the flow and heat transfer of eyring powell fluid over an inifinite rotating disk[END_REF][START_REF] Hsiang-Chin | A Lubrication Theory for Anisotropic Slips and Flow Rheology[END_REF][START_REF] Cooper | The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer[END_REF]). But it is important to mention that the numerical analysis (solvability and a priori error estimates) of these models has never been analysed in the literature. Although this work takes into account the vast knowledge on the Stokes/Navier -Stokes analysis, it should be mentioned that the nonlinear slip boundary condition (1.7) adds some difficulties in the treatment of the resulting variational problems. The new challenge in the finite element analysis being the derivation of the uniform a priori error estimate with respect to the parameter in the power law constitutive relation. It should be noted that the power law boundary condition (1.7) we analyse in the work is not the same as the power law fluid analysed in [START_REF] Baranger | Analyse numerique des ecoulements quasi-Newtoniens dont la viscosite obeit a la loi puissance ou la loi de carreau[END_REF][START_REF] Glowinski | Sur l approximation par elements finis d ordre un et la resolution par penalisation-dualite d une classe de problemes de Dirirchlet nonlineaires[END_REF][START_REF] Sandri | Sur l approximation des ecoulements numeriques quasi-Newtoniens dont la viscosite obeit a la loi de puissance ou de Carreau[END_REF][START_REF] Djoko | Power law Stokes equations with threshold slip boundary conditions: Numerical methods and implementation[END_REF]. The later being defined on the flow domain Ω, and the constitutive relation bringing together the Cauchy stress tensor and the symmetric part of the velocity gradient. Finally, we observe that (1.7) is a generalisation of the Navier's slip boundary condition, which is an active research direction (see [START_REF] Verfurth | Finite element approximation of incompressible Navier Stokes equations with slip boundary condition II[END_REF][START_REF] Verfurth | Finite element approximation of incompressible Navier Stokes equations with slip boundary condition II[END_REF][START_REF] Dione | Stokes equations with penalized slip boundary conditions[END_REF][START_REF] Dione | Penalty finite element approximation of Stokes equations with slip boundary conditions[END_REF][START_REF] Zhou | The CrouzeixRaviart element for the Stokes equations with the slip boundary condition on a curved boundary[END_REF][START_REF] Zhou | Penalty method with Crouzeix-Raviart approximation for the Stokes equations under slip boundary condition[END_REF]).

The rest of the paper is organised in the following way. In Section 2, we present the variational problems and study their solvability. Section 3 is devoted to the description and construction of the finite element solution. Section 4 is devoted to the convergence of the finite element solution. Here we establish the convergence of the finite element solution by showing that difference in appropriate norm of the continuous and discrete solution tends to zero with the discretization parameter. In Section 5, we formulate the iterative schemes associated to the finite element approximations. The focus is mainly on the problem associated to the Navier-Stokes equations. We establish some qualitative properties and perform numerical experiments with FreeFem code [START_REF] Hecht | New development in FreeFem++[END_REF]. We complete this study with some concluding remarks.

Analysis of the continuous problems

In this section, we formulate the weak problems associated respectively with (1.1)-(1.4), (1.7), and (1.2)-(1.4), (1.7)-(1.8), and study their well posedness.

Variational analysis of the Stokes system

We first introduce some notations that will be used throughout this work. Next, we propose and discuss the weak solution of (1.1)-(1.4), and (1.7).

Notations

We adopt the standard definitions of Lebesgue and Sobolev spaces (see [START_REF] Brezis | Functional Analysis, Sobolev Spaces and partial differential equations[END_REF]). For 1 ≤ p < ∞ L p (Ω) = {v : Ω → R; v is measurable and

Ω |v(x)| p dx < ∞}, with associated norm v p L p (Ω) = Ω |v(x)| p dx
for which it is a Banach space, and for the special case p = 2, it is a Hilbert space with the norm that will be denoted by • . Finally for p = ∞, one has

L ∞ (Ω) = v : Ω → R, v is measurable
v L ∞ (Ω) = v ∞ = inf{c, |v(x)| ≤ c a.e on Ω} .
Now, for any non-negative integer m and a real number p ≥ 1, we defined the Sobolev space

W m,p (Ω) = {v ∈ L p (Ω), ∂ α v(x) ∈ L p (Ω) for all |α| ≤ m}, with α = (α 1 , ..., α d ) the multi-index with |α| = α 1 + • • • + α d , and ∂ α v(x) the distributional derivative of v given by ∂ α v = ∂ |α| v ∂x α1 1 • • • ∂x α d d .
The space W m,p (Ω) is equipped with the semi-norm and norm

|v| p W m,p (Ω) = |α|=m Ω |∂ α v(x)| p dx , v p m,p = 0≤|α|≤m Ω |∂ α v(x)| p dx .
The duality between, say, E and its dual E is denoted as • . The bold characters denote vector quantities and L p (Ω) = L p (Ω) d , W m,p (Ω) = W m,p (Ω) d , etc.... We note that W m,2 (Ω) = H m (Ω) and H 0 (Ω) coincides with L 2 (Ω), for which the norm and inner product are denoted as • and (•, •), respectively. In order to introduce functions spaces related to the equations (1.1)-(1.4), (1.7), we recall that (1.5) is re-written as follows 2ν

Ω |Du| 2 dx + Γ |Ku τ | s dσ -p, div u = f , u . (2.1)
Hence, it is manifest from (2.1) that for the velocity one needs u ∈ H 1 (Ω), and u τ ∈ L s (Γ). Then the pressure must belong to

L 2 (Ω) because div u ∈ L 2 (Ω). It should be noted that for v ∈ H 1 (Ω), v τ is an element of H 1/2 (S) ⊂ L 2 (Γ) which is not necessarily a subspace of L s (Γ). Hence we have two possibilities. If 1 ≤ s < 2, we have u τ L 1 (Γ) ≤ c u τ L s (Γ) ≤ c u τ L 2 (Γ)
. Hence the velocity is described in

V = {v ∈ H 1 (Ω) : with v| S = 0 and v • n| Γ = 0} ,
equipped with the regular H 1 -norm, while the pressure is in the space

M = {p ∈ L 2 (Ω) : (p, 1) = 0} = L 2 0 (Ω) .
We mentioned that the zero mean value is added on the pressure to avoid the pressure given in (1.1), (1.2) to be determined up to a constant. For 2 < s, we have u τ L 2 (Γ) ≤ c u τ L s (Γ) . Then the velocity is described in

V = {v ∈ H 1 (Ω) , v τ ∈ L s (Γ) : with v| S = 0 and v • n| Γ = 0}
equipped with the norm

|||v||| 2 = v 2 1 + v τ 2 L s (Γ)
. The pressure is described in the same function space L 2 0 (Ω), because we only need p, div u to make sense with div u ∈ L 2 (Ω). We shall often use the following Sobolev imbedding (see [START_REF] Brezis | Functional Analysis, Sobolev Spaces and partial differential equations[END_REF]):

for any real number p ≥ 1 when d = 2, or 1 ≤ p ≤ 2d d -2
when d ≥ 3, there exist constants c p and

c 0 p such that for all v ∈ H 1 (Ω), v L p (Ω) ≤ c p v 1 , for all v ∈ V , v L p (Ω) ≤ c 0 p ∇v . (2.2)
We note that when p = 2, the second inequality in (2.2) is the Poincaré's inequality. Also of importance is the Korn's inequality which reads (see [START_REF] Ciarlet | Three dimensional elasticity[END_REF]); there exists a constant c such that

for u ∈ V , Ω |Du| 2 dx ≥ c Ω |∇u| 2 dx . (2.3)

Variational formulation

In this part, we formulate the weak problem associated with (1.1)-(1.4), (1.7) and construct the weak solution.

We assume once and for all that for all that f ∈ L 2 (Ω).

Let q ∈ M , we multiply (1.2) by q, integrate the resulting equation over Ω. Next, we take the dot product between (1.1) and v ∈ V and integrate the resulting equation over Ω. After utilization of the Green's formula, (1.3), (1.4) and (1.7), one obtains:

         Find (u, p) ∈ V × M such that for all (v, q) ∈ V × M, a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ -b(v, p) = (v) , b(u, q) = 0 , (2.4) with a(u, v) = 2ν Ω Du : Dvdx , b(v, q) = Ω q div vdx , (v) = Ω f • vdx, where A : B = 1≤i,j≤d A ij B ij .
In the derivation of (2.4), we have also used the identity 1≤i,j≤d

D ij u ∂v i ∂x j = 1≤i,j≤d D ij uD ij v which
is the consequence of the symmetry of Du. Hence we have the definition.

Definition 2.1 A couple (u, p) is called weak solution of (1.1)-(1.4), (1.7) if (u, p) ∈ V × M and satisfy (2.4).
One of the key step in the study of mixed variational problem like (2.4) is the definition of the kernel of the bilinear form b(•, •) which permits the elimination of the pressure in the problem. Thus we consider the space

V div = {v ∈ V, for all q ∈ M, b(v, q) = 0} = {v ∈ V, div v = 0} .
With the space V div , the variational problem (2.4) is equivalent to

     Find u ∈ V div such that for all v ∈ V div , a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ = (v) . (2.5)
In the study of mixed problem like (2.4), the compatibility or inf-sup condition on b(•, •) is very important. It reads; there exists β > 0 such that

inf q∈M sup 0 =v∈V b(v, q) v V q ≥ β . (2.6)
Actually for 1 ≤ s < 2, one has V = {v ∈ H 1 (Ω) : with v| S = 0 and v • n| Γ = 0} and M = L 2 0 (Ω) . Eq. (2.6) is obtained by observing that H 1 0 (Ω) ⊂ V and the pair (H 1 0 (Ω), L 2 0 (Ω)) is inf-sup stable (see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF][START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]), hence (2.6) holds. For 2 < s, V = {v ∈ H 1 (Ω) , v τ ∈ L s (S) : with v| S = 0 and v • n| Γ = 0} and M = L 2 0 (Ω). We were not able to obtain (2.6) which is a key point in obtaining the pressure when studying (2.4).

In the rest of this analysis, we consider 1 < s < 2. This corresponds to the tangential shear thinning with same physical interpretation as the power law fluid (see [START_REF] Baranger | Analyse numerique des ecoulements quasi-Newtoniens dont la viscosite obeit a la loi puissance ou la loi de carreau[END_REF][START_REF] Glowinski | Sur l approximation par elements finis d ordre un et la resolution par penalisation-dualite d une classe de problemes de Dirirchlet nonlineaires[END_REF][START_REF] Sandri | Sur l approximation des ecoulements numeriques quasi-Newtoniens dont la viscosite obeit a la loi de puissance ou de Carreau[END_REF]).

To study (2.5) or (2.4) it is convenient to recall the following monotonicity and continuity properties (see [START_REF] Glowinski | Sur l approximation par elements finis d ordre un et la resolution par penalisation-dualite d une classe de problemes de Dirirchlet nonlineaires[END_REF][START_REF] Sandri | Sur l approximation des ecoulements numeriques quasi-Newtoniens dont la viscosite obeit a la loi de puissance ou de Carreau[END_REF]): there exists a constant c independent of x, y elements of R n such that for 1 ≤ s < 2;

(|y| + |x|) 2-s |y| s-2 y -|x| s-2 x, y -x ≥ c|x -y| 2 , (2.7)
and

|x| s-2 x -|y| s-2 y ≤ c|x -y| s-1 . (2.8)
Now we introduce the mapping v -→ A(v) defined as follows

for all w ∈ V div , A(v)w = a(v, w) + Γ |Kv τ | s-2 Kv τ • Kw τ dσ .
(2.9)

We have the following properties for A.

Lemma 2.2 Let 1 < s < 2, then (a)
A maps V div into its dual V div , and is bounded on all bounded subsets of V div .

(b) For all v, u elements of

V div A(v) -A(u) V ≤2ν v -u 1 + c K s L ∞ (Γ) u -v s-1 1
.

Proof. Let u, v, w elements of V div . From Hölder's and Cauchy-Schwarz's inequalities, we have

A(v)w ≤2ν v 1 w 1 + K s L ∞ (Γ) Γ |v τ | s-1 |w τ | dσ ≤2ν v 1 w 1 + c K s L ∞ (Γ) v s-1 1 w 1 , so A(v) V ≤ 2ν v 1 + c K s L ∞ (Γ) v s-1 1 .
Next, from (2.8), Hölder's and Cauchy-Schwarz's inequalities and since K is bounded, one has

(A(v) -A(u), w) ≤2ν v -u 1 w 1 + Γ |Ku τ | s-2 Ku τ -|Kv τ | s-2 Kv τ |Kw τ | dσ ≤2ν v -u 1 w 1 + c K L ∞ (Γ) Γ |Ku τ -Kv τ | s-1 |w τ | dσ ≤2ν v -u 1 w 1 + c K s L ∞ (Γ) u -v s-1 1 w 1 , hence A(v) -A(u) V ≤ 2ν v -u 1 + c K s L ∞ (Γ) u -v s-1 1
.

Existence of solutions

Here we construct a solution of (2.4) by using the theory of monotone operator. We first start by indicating the following result Proposition 2.1 The problem (2.5) is equivalent to the problem (2.4) .

This is a standard result for the Stokes problem and we refer the reader to [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF][START_REF] Lions | Quelques Méthodes de Résolution des Problèmes aux Limites Non Linéaires[END_REF] where similar proofs are given.

Owing to the equivalence in Proposition 2.1, we can restrict the analysis to problem (2.5). We shall prove that A is strictly monotone, coercive and hemi-continuous in V. Then the existence and uniqueness of (2.5) will follow by standard arguments, see [START_REF] Lions | Quelques Méthodes de Résolution des Problèmes aux Limites Non Linéaires[END_REF]Chapter 2]. We first start with the monotonicity property of A.

Lemma 2.3

The mapping v → Av defined by (2.9) is strictly monotone from

V div into V div : For all v, w ∈ V div , (Av -Aw)(v -w) > 0 .
Proof. Consider the functional J : V div → R defined by

J(v) = 1 2 a(v, v) + 1 s Γ |Kv τ | s dσ .
The first derivative of J is

J (u) • v = a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ , (2.10)
and the second derivative is

J (2) (u) • (v, w) =a(v, w) + (s -2) Γ |Ku τ | s-4 (Ku τ • Kw τ ) (Ku τ • Kv τ ) dσ + Γ |Ku τ | s-2 Kw τ • Kv τ dσ .
(2.11)

Note that J (2) is symmetric and positive definite. Indeed

J (2) (u) • (v, w) = J (2) (u) • (w, v) and
J (2) (u) • (v, v) =a(v, v) + (s -2) Γ |Ku τ | s-4 (Ku τ • Kv τ ) (Ku τ • Kv τ ) dσ + Γ |Ku τ | s-2 Kv τ • Kv τ dσ .
From the assumption on K, one deduces that K 2 is symmetric and uniformly positive definite. Let α min > 0 be the smallest eigenvalue of K 2 over Γ, that is

∀x ∈ Γ, ∀w ∈ R d , K 2 (x)w • w = |K(x)w| 2 ≥ α min |w| 2 . Then ∀x ∈ Γ, ∀w ∈ R d , (K 2 (x)w • w) s/2 = |K(x)w| s ≥ α s/2 min |w| s .
(2.12) Since 1 ≤ s < 2, one has

J (2) (u) • (v, v) ≥ a(v, v) - Γ |Ku τ | s-2 |Kv τ | 2 dσ + Γ |Ku τ | s-2 |Kv τ | 2 dσ = a(v, v) ≥ 0 .
One notes from (2.10) that

(A(v) -A(w))(v -w) =J (v) • (v -w) -J (w) • (v -w) = 1 0 J (2) (w + θ(v -w)) • (v -w, v -w)dθ ≥ 0
because the second derivative is positive definite. We then conclude that A is monotone (in fact strictly monotone).

Secondly, we have the following hemi-continuity result for A.

Lemma 2.4 A is hemi-continuous in V div , that is for u, v in V div , the mapping t -→ A(u + tv)v is continuous from R into R.
Proof. Indeed for any t 1 , t 2 in R, we can write from (2.10) that

(A(u + t 1 v) -A(u + t 2 v)) • v = (J (u + t 1 v) -J (u + t 2 v)) • v =(t 1 -t 2 ) 1 0 J (2) (u + t 2 v + θ(t 1 -t 2 )v) • (v, v)dθ .
Since u and v are fixed, it follows from (2.11) that the right hand side term on the second equality tends to zero with t 1 -t 2 .

Finally we show that A is coercive in the following way

Lemma 2.5 A is coercive, that is lim v 1→∞ 1 v 1 A(v)v = ∞ Proof.
From (2.12), (2.3) and (2.2) one has

1 v 1 A(v)v ≥ min(2ν, α s/2 min ) Dv 2 + v τ s L s (Γ) v 1 ≥ min(2ν, α s/2 min ) ∇v 2 + v τ s L s (Γ)
∇v .

From which we deduce that lim

∇v →∞ A(v)v ∇v = ∞ .
From Lemma 2.3, Lemma 2.5, and Lemma 2.4, we deduce the following existence result.

Theorem 2.1 If 1 < s < 2, then the variational problem (2.4) has a unique solution (u, p) ∈ V × M , which moreover satisfy the bounds

νc ∇u 2 + S |Ku τ | s dσ ≤ c f 2 , p ≤ 2ν ∇u + K s L ∞ (Γ) u τ s-1 L s (Γ) + c f . Remark 2.2
The contribution of the nonstandard boundary condition (1.7) to the a priori estimate in Theorem 2.1 is manifest by the presence of the boundary term Ku τ .

Variational analysis of the Navier-Stokes system

This part is devoted to the analysis of the problem related to (1.2)-(1.4), (1.7) and (1.8). The weak solution is constructed using a combination of approaches such as; Faedo-Galerkin's approximation, Brouwer's fixed point and compactness properties. The reader will find more details for example in [31, Chap 2].

Existence and uniqueness of solution

Turning to the Navier-Stokes system, the weak formulation associated with (1.2), (1.3), (1.4), (1.7) and (1.8) reads

         Find (u, p) ∈ V × M such that for all (v, q) ∈ V × M, a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ + d(u, u, v) -b(v, p) = (v), b(u, q) = 0 , (2.13) with d(u, v, w) = Ω (u • ∇)v • wdx . We recall that d(u, v, w) = -d(u, w, v) for all u, v, w ∈ V div × H 1 (Ω) × H 1 (Ω) d(u, v, v) = 0 for all u, v ∈ V div × H 1 (Ω) .
(2.14)

We discuss next, the existence and uniqueness of solution of (2.13). We state the following result.

Proposition 2.2 For any data f ∈ L 2 (Ω), the variational problem (2.13) has at least one solution

(u, p) ∈ V × M if 1 < s < 2, moreover there is a constant c such that νc u 2 1 + Γ |Ku τ | s dσ ≤ c ν f 2 , p ≤ 1 β 2ν u 1 + K s L ∞ (Γ) u τ s-1 L s-1 (Γ) + c u 2 1 + f . If moreover ν 2 > c f , (2.15)
then the solution is unique.

Proof. It is performed in several steps. We recall (2.13) is equivalent to

     Find u ∈ V div such that for all v ∈ V div , a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ + d(u, u, v) = (v) . (2.16)
We first construct the solution u ∈ V div solution of (2.16).

Step 1: Existence of approximate solutions. Since V div is separable, then there exist an increasing sequence a countable basis {ψ i } i≥1 . Let K m be the space spanned by the first m basis functions, {ψ i } 1≤i≤m . The problem (2.16) is discretized in K m by the square system of nonlinear equations:

Find u m ∈ K m such that for all v ∈ K m , a(u m , v) + Γ |Ku τ ,m | s-2 Ku τ ,m • Kv τ dσ + d(u m , u m , v) = (v) .
(2.17) to prove the existence of u m , we use fixed point by considering the mapping F : K m -→ K m defined as follows

F (u)v = a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ + d(u, u, v) -(v) .
We first show that F is continuous. Let u, v, w in K m , we have

|(F (u) -F (v))w| ≤2ν u -v 1 w 1 + |d(u -v, u, w)| + |d(v, u -v, w)| + Γ |Ku τ | s-2 Ku τ -|Kv τ | s-2 Kv τ |Kw τ | dσ ≤2ν u -v 1 w 1 + u -v L 4 (Ω) ∇u w L 4 (Ω) + v L 4 (Ω) ∇(u -v) w L 4 (Ω) + Γ |Ku τ | s-2 Ku τ -|Kv τ | s-2 Kv τ |Kw τ | dσ ,
together with the inequalities (2.8), implies that F is continuous. On the other hand, for v ∈ K m , one has

F (v)v =a(v, v) + Γ |Kv τ | s dσ -(v) ≥2νc v 2 1 + Γ |Kv τ | s dσ -f v 1 ≥ (2νc v 1 -f ) v 1 + Γ |Kv τ | s dσ .
Choosing r such that r > c ν f , it follows that for any v ∈ K m , with v 1 = r, we have F (v)v > 0. We deduce from a corollary of the Brouwer's fixed point (see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF]), that there exists u m ∈ K m solution of (2.17).

Step 2: Estimates and passage to the limit. We take v = u m in (2.17), using standard inequalities, this gives

νc u m 2 1 + Γ |Ku τ ,m | s dσ ≤ c ν f 2 . (2.18)
Hence there exists a subsequence, still denoted by u m such that

u m → u weakly in H 1 (Ω), Bu m → W weakly in V , (2.19)
because Bv ≤ c v s-1

1

, with

B(u)v = Γ |Ku τ | s-2 Ku τ • Kv τ dσ.
Owing to the compactness of the imbedding of H 1 (Ω) into L 4 (Ω), there exists a subsequence, still denoted by u m , which converges weakly to u ∈ H 1 (Ω) and strongly in L 4 (Ω) d . We easily check that div u = 0 in Ω, u • n| Γ = 0 and u| S = 0. Passing to the limit on m is direct for the linear term because one only needs weak convergence, but for the term (u m • ∇)u m , one needs the strong convergence of u m on L 4 (Ω) d (see [29, Chap 5, Thm 1.4]). Indeed

|d(u m , u m , v) -d(u, u, v)| = |d(u m -u, u m , v) + d(u, u m -u, v)| ≤ u m -u L 4 (Ω) ∇u m v L 4 (Ω) + |d(u, u m -u, v)| , (2.20)
and the right hand side of (2.20) tends to zero with m. Hence d(u m , u m , v) tends to d(u, u, v).

Finally for the term on Γ, we need to show that W = Bu. First, we deduce from the weak convergence on H 1 (Ω), that

u m → u strongly in L 2 (Ω)) .
The mapping u -→ u τ is continuous from H 1 (Ω) to L 2 (Γ). Then

u m | Γ → u τ | Γ weakly in L 2 (Γ).
Now one can use the fact that the mapping v → Bv is monotone, and hemi-continuous to show that (see [START_REF] Lions | Quelques Méthodes de Résolution des Problèmes aux Limites Non Linéaires[END_REF]Chap 2]) W = Bu. Therefore, with the help of the density of K m in V div , one deduces that the limit u satisfies the equation

     for all v ∈ V div , a(u, v) + Γ |Ku τ | s-2 Ku τ • Kv τ dσ + d(u, u, v) = (v) . (2.21)
Taking v = u in (2.21), and using standard inequalities, one obtains the a priori estimate on the velocity announced.

Step 3: Construction of the pressure. For 1 < s < 2, the inf-sup condition (2.6) holds and owing to (2.21), there exists a p in M see [29, Chap I, Lemma 4.1] such that

           for all (v, q) ∈ V × M, a(u, v) + |Ku τ | s-2 Ku τ , Kv τ Γ + d(u, u, v) -(v) = Ω div vpdx Ω q div udx = 0 . (2.22)
The a priori estimate on the pressure is obtained by using the inf-sup condition and the velocity equation.

Step 4: Uniqueness.

Let (u 1 , p 1 ) and (u 2 , p 2 ) be two solutions of (2.22). Setting for a while u = u 1 -u 2 and p = p 1 -p 2 , we deduce that

a(u, v) + Γ |Ku 1,τ | s-2 Ku 1,τ -|Ku 2,τ | s-2 Ku 2,τ • Kv τ dσ + d(u 1 , u 1 , v) -d(u 2 , u 2 , v) = Ω div vpdx .
We take successively v = u 1 , and v = u 2 , and take the difference of the resulting equations. This gives a(u, u)

+ Γ |Ku 1,τ | s-2 Ku 1,τ -|Ku 2,τ | s-2 Ku 2,τ • (Ku 1,τ -Ku 2,τ ) ≥0 dσ = -d(u, u 2 , u) , (2.23) which implies that 2νc u 2 1 ≤ u 2 L 4 (Ω) ∇u 2 ≤ c u 2 1 u 2 1 ≤ c f ν u 2 1 ,
which gives us the condition for the unique solvability. Next, the function p = p 1 -p 2 satisfies ∀v ∈ V , Ω div vpdx = 0 which together with (2.6), implies that p 1 = p 2 .

Remark 2.3

The smallness condition for uniqueness is not surprising given that we have a nonlinear problem. On the other hand, that condition is similar to the uniqueness condition for the Navier-Stokes equations with Dirichlet boundary condition. The similarity in our view is due to the fact that the "extra" nonlinear term is monotone and can be removed from the equation when analysing (see (2.23)).

The next section is about the construction of finite element solutions of (2.4) and (2.13).

Finite element approximation

This section is twofold. First, we introduce the basics tools for the formulation of the finite element approximation of (2.4) and (2.13). Secondly, we discuss the solvability of the nonlinear finite elements problems.

Preliminaries

From now on, we assume that Ω is a polygon in two dimensions or polyhedron in three dimensions, so it can be entirely triangulated by triangles or tetrahedron, according to the dimensions d. A regular family of triangulations (T h ) h of Ω, is a set of closed non degenerate triangles or tetrahedra, called elements, satisfying (see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]) (i) For each h, Ω is the union of all elements of T h , that is

Ω = K∈T h K .
(ii) The ratio of the diameter of an element K ∈ T h to the diameter of its inscribed circle or ball is bounded by a constant independent of h. This is to say that there exists a constant σ, independent of h and K such that

∀ K ∈ T h , h K ρ K = σ K ≤ σ ,
where h K is the diameter of K, and ρ K is the diameter of the sphere (circle) inscribed in K.

(iii) The intersection of two distinct elements of T h is either empty, a common vertex, or an entire common edge or face.

As usual, h is the maximal diameter of all elements of T h . For every K ∈ T h , we denote by P k the space of restrictions to K of polynomials of degree k in d variables and total degree at most k.

In what follows, c stand for a generic constant which may vary from line to line but are always independent of h.

Remark 3.1 In this work, we do not address the situation with curved boundary, but we refer the interested readers to the work in [START_REF] Dione | Penalty finite element approximation of Stokes equations with slip boundary conditions[END_REF][START_REF] Zhou | The CrouzeixRaviart element for the Stokes equations with the slip boundary condition on a curved boundary[END_REF].

Discretization of the Stokes system

Finite element approximation.We discretize the velocity and the pressure respectively in the spaces V h and M h defined as follows

V h = v h ∈ V ∩ C 0 (Ω) d , for all T ∈ T h , v h | K ∈ (P 1 (T ) + bubble) d , L 2 h = q h |q h ∈ C 0 (Ω), q h | T ∈ P 1 , for all T ∈ T h , M h = q h |q h ∈ L 2 h , Ω q h dx = 0 , (3.1)
where P 1 (T )+bubble is the sum of a polynomial of P 1 (T ) and a bubble function b T (x) = α 1 (x)α 2 (x)α 3 (x), for any T ∈ T h , and denoting the vertices of T ∈ T h by a i , 1 ≤ i ≤ 3, and its corresponding barycentric coordinates by α i . Note that b T (x) = 0 on T and that b T (x) > 0 on T .

With the above spaces it is quite natural to approximate (2.4) (with obvious notation) by:

         Find (u h , p h ) ∈ V h × M h such that for all (v h , q h ) ∈ V h × M h a(u h , v h ) + Γ |Ku τ ,h | s-2 Ku τ ,h • Kv τ ,h dσ -b(v h , p h ) = (v h ) , b(u h , q h ) = 0 . (3.2)
For 1 < s < 2, the discrete inf-sup holds, that is, there exists β > 0 such that inf

q h ∈M h sup 0 =v h ∈V h b(v h , q h ) v h 1 q h ≥ β . (3.3)
We recall that on the space

V div,h = v h ∈ V h : Ω q h div v h dx = 0 for all q h ∈ M h , (3.2) is equivalent to:      Find u h ∈ V div,h such that for all v h ∈ V div,h a(u h , v h ) + Γ |Ku τ ,h | s-2 Ku τ ,h • Kv τ ,h dσ = (v h ) . (3.4)
Remark 3.2 It is noted that a function in V h is reduced to polynomial of degree one on each side e of the element K and also on each boundary side. So, the condition u h • n| Γ = 0 implies that u h vanishes at each corner of Γ.

The velocity and pressure are interpolated respectively by the operators Π h ∈ L(V; V h ) and r h ∈ L(M ; M h ) constructed in Girault and Hecht [START_REF] Girault | Numerical Methods for Grade-Two Fluid Models: Finite-Element Discretizations and Algorithms In Handbook of Numerical Analysis[END_REF]Chap 5]. These operators satisfy, for each real number α ∈ [0, 1] and for each number m ≥ 2, the following properties. (1) There exists a constant c, independent of h such that

for all v ∈ W α+1,m (Ω) 2 ∩ V , |Π h v -v| W 1,m (Ω) ≤ ch α |v| W α+1,m (Ω) . (3.5)
(2) There exists a constant c, independent of h such that

for all q ∈ W α,m (Ω) 2 ∩ M , r h q -q L m (Ω) ≤ ch α |v| W α,m (Ω) . (3.6)
A key issue when studying the mixed formulation (3.2) with the spaces V h and M h is the discrete inf-sup condition (3.3). Its proof consists in using the trick introduce by Fortin (see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF]) based on the operator

H(v) = Π h (v) + K∈T h a K (v)b K where a K (v) = 1 K b K dx K (v -Π h (v))dx .
Existence and uniqueness. The Existence of a solution of (3.4) (and (3.2)) can be obtained by repeating the steps of paragraph 2.2. One observes that the variational problem (3.4) is a finite dimensional, square system of nonlinear equations. Thus just like the continuous analysis, we introduce the following mapping; for all (v h ,

w h ) ∈ V div,h × V div,h A h (v h )w h = a(v h , w h ) + Γ |Kv τ ,h | s-2 Kv τ ,h • Kw τ ,h dσ .
We can show that for all 1 < s < 2, A h is strictly monotone, coercive and hemi-continuous. Thus, the analog of theorem 2.1 reads as follows.

Theorem 3.1 For 1 < s < 2, the problem (3.2) has exactly one solution (u h , p h ) in V h × M h , and there exists c independent of h such that the following a priori estimates hold

∇u h 2 + Γ |Ku τ ,h | s dσ ≤ c f 2 , p h ≤ 2ν ∇u h + K s L ∞ (Γ) u τ ,h s-1 L s (Γ) + c f .

Analysis of the Navier-Stokes system

Finite element approximation. The velocity and pressure are taken as before. We consider the following finite element scheme

         Find (u h , p h ) ∈ V h × M h such that for all (v h , q h ) ∈ V h × M h , a(u h , v h ) + Γ |Ku τ ,h | s-2 Ku τ ,h • Kv τ ,h dσ + d(u h , u h , v h ) -b(v h , p h ) = (v h ), b(u h , q h ) = 0 , (3.7)
with the trilinear form d(•, •, •) given by R.Temam [START_REF] Temam | Navier-Stokes equations: Theory and Numerical Analysis[END_REF]

d(v h , u h , w h ) = d(v h , u h , w h ) + 1 2 ((div v h )u h , w h ) = 1 2 (d(v h , u h , w h ) -d(v h , w h , u h )) . It is noted that d(•, •, •) is consistent with d(•, •, •) in the sense that for all (v, u, w) ∈ V div × H 1 (Ω) × H 1 (Ω) , d(v, u, w) = d(v, u, w) . Furthermore d(•, •, •) is anti-symmetry meaning that for all (v h , u h , w h ) ∈ H 1 (Ω) × V h × V h , d(v h , u h , w h ) = -d(v h , w h , u h ) .
Existence and uniqueness. The kernel of divergence in V h is V div,h , and for 1 ≤ s < 2, (3.7) is equivalent to

     Find u h ∈ V div,h such that for all v h ∈ V div,h , a(u h , v h ) + Γ |Ku τ ,h | s-2 Ku τ ,h • Kv τ ,h dσ + d(u h , u h , v h ) = (v h ) . (3.8) (3.8
) is a finite dimensional, square system of nonlinear equations. We address next the solvability of (3.8) by Brouwer's Fixed point argument. To this end, for fixed u h ∈ V div,h , one introduces

F (u h ) in V div,h by (F (u h ), v h ) =a(u h , v h ) + d(u h , u h , v h ) + Γ |Ku τ ,h | s-2 Ku τ ,h • Kv τ ,h dσ - Ω f • v h dx .
From the analysis in the continuous formulation, we deduce that the mapping

F : V div,h -→ V div,h is continuous. Next, using the anti-symmetry property of d(•, •, •), (F (v h ), v h ) =a(v h , v h ) + Γ |Kv τ ,h | s dσ - Ω f • v h dx ≥ v h 1 (2νc v h 1 -f ) + Γ |Kv τ ,h | s dσ . Thus (F (v h ), v h ) is non-negative for all v h ∈ V h , with v h 1 = r > c ν
f . So, this proves existence of at least one solution of (3.8). Using the equivalence between (3.7) and (3.8), we construct the pressure and we have the following result. Proposition 3.1 The finite element problem (3.7) admits at least one solution (u h , p h ) ∈ V h ×M h , and there exist positive constant c such that

νc u h 2 1 + Γ |Ku τ ,h | s dσ ≤ c ν f 2 , β p h ≤ 2ν u h 1 + K s L ∞ (Γ) u τ ,h s-1 L s-1 (Γ) + c u h 2 1 + f .
The estimate above is obtained in the standard way (see proposition 2.2). One can observe the small difference between the estimate on the pressure in proposition 2.2 and proposition 3.1. This difference being the fact that in the discrete setting div u h is non zero. We end this analysis with the result below. Proposition 3.2 Let (u h , p h ) ∈ V h × M h be the solution of (3.7). There exists a positive constant c depending only on Ω such that if for ν, f , the relation

ν 2 ≥ c f ,
is satisfied, then the solution of (3.7) is unique.

Proof. Let (u 1 , p 1 ) and (u 2 , p 2 ) be two solutions of (3.7). Setting for a while u = u 1 -u 2 and p = p 1 -p 2 . the equation (2.23) is

a(u, u) + Γ |Ku 1,τ | s-2 Ku 1,τ -|Ku 2,τ | s-2 Ku 2,τ • (Ku 1,τ -Ku 2,τ ) ≥0 dσ = -d(u, u 2 , u) , (3.9) this leads to 2νc u 2 1 ≤ u 2 L 4 (Ω) ∇u 2 ≤ c u 2 1 u 2 1 ≤ c f ν u 2
1 , which gives us the condition for the unique solvability. Secondly, the function p satisfies ∀v ∈ V h , Ω div vpdx = 0 which together with the discrete inf-sup condition (3.3), implies that p 1 = p 2 .

Remark 3.3

The condition obtained for the unique solvability of (3.7) in Proposition 3.2, is not the same condition obtained in (2.15) for the unique solvability of (2.13). This difference being that in the discrete setting div u h is non zero. Hence the constants are not the same.

At this juncture, knowing that solution of (3.7) and (3.2) are well defined and computable, we next study the convergence of the finite element solution (u h , p h ). This exercise is very important in numerical analysis because by doing it we check the "reliability" of the approximation.

A priori error analysis

We now study the convergence of the finite element solution (u h , p h ) by computing an upper bound of the quantity uu h 1 + p -p h . The problems (3.7) and (3.2) are well defined, hence study the convergence of the finite element solution make perfect sense. It is important to recall that without uniqueness no convergence result can be formulated.

A priori error analysis: the Stokes equations

We start with the convergence of the solution (u h , p h ) of (3.2). It holds that Theorem 4.1 For 1 < s < 2, let (u, p) be the solution of (2.4) and (u h , p h ) the solution of (3.2). Then there is a positive constant c independent of h such that

u -u h 1 + p -p h + u τ -u τ ,h Γ ≤ c inf v h ∈V h u τ -v τ ,h 1/2 Γ + c inf v h ∈V h ∇(u -v h ) + c inf q h ∈M h q h -p . Proof.
Step 1: Inf-sup condition. Let (u h , p h ) be the unique solution of (3.2). Let q h ∈ M h , then p h -q h ∈ M h and in particular, there exists β independent of h such that

β p h -q h ≤ sup 0 =v h ∈V h ∩{v h |Γ=0} b(v h , p h -q h ) v h 1 . (4.1)
We take the difference between (2.4) and (3.2) and obtain (after observing that V h ⊂ V and

M h ⊂ M )              for all (v h , q h ) ∈ V h × M h , a(u -u h , v h ) + |Ku τ | s-2 Ku τ -|Ku τ ,h | s-2 Ku τ ,h , Kv τ ,h Γ = b(v h , p -p h ) , b(u -u h , q h ) = 0 . (4.2)
We solve for b(v h , p h -q h ) in the first equation of (4.2), report the result to (4.1), and using Cauchy-Schwarz's, Hölder's and trace's inequality gives

β p h -q h ≤ sup 0 =v h ∈V h ∩{v h |Γ=0} |a(u -u h , v h ) -b(v h , q h -p)| v h 1 ≤ 2ν ∇(u -u h ) + q h -p . (4.3)
Step 2:

Monotonicity. Set Bu τ = |Ku τ | s-2 Ku τ . Let w h ∈ V h , take v h = w h -u h , (4.2) leads to          for all (w h , q h ) ∈ V h × M h , a(u -u h , w h -u h ) + Γ (Bu τ -Bu τ ,h ) • (Kw τ ,h -Ku τ ,h )dσ = b(w h -u h , p -p h ) , b(u -u h , q h -p h ) = 0 .
which is re -written as follows

             a(w h -u h , w h -u h ) + Γ (Bu τ -Bu τ ,h ) • (Kw τ ,h -Ku τ ,h )dσ = a(w h -u, w h -u h ) + b(w h -u h , p -p h ) b(w h -u h , p -p h ) = b(w h -u, q h -p h ) + b(u h -w h , q h -p) . (4.4)
We use the second equation of (4.4) in the first one and obtain that

a(w h -u h , w h -u h ) + Γ (Bu τ -Bu τ ,h ) • (Kw τ ,h -Ku τ ,h )dσ =a(w h -u, w h -u h ) + b(w h -u, q h -p h ) + b(u h -w h , q h -p) . (4.5)
From linearity we have

Γ (Bu τ -Bu τ ,h ) • (Kw τ ,h -Ku τ ,h )dσ = Γ (Bu τ -Bu τ ,h ) • (Kw τ ,h -Ku τ )dσ + Γ (Bu τ -Bu τ ,h ) • (Ku τ -Ku τ ,h )dσ . (4.6)
Thus (4.5) becomes

a(w h -u h , w h -u h ) + Γ (Bu τ -Bu τ ,h ) • (Ku τ -Ku τ ,h )dσ =a(w h -u, w h -u h ) + b(w h -u, q h -p h ) + b(u h -w h , q h -p) + Γ (Bu τ -Bu τ ,h ) • (Ku τ -Kw τ ,h )dσ (4.7)
From (2.7), it was shown that (see [START_REF] Sandri | Sur l approximation des ecoulements numeriques quasi-Newtoniens dont la viscosite obeit a la loi de puissance ou de Carreau[END_REF]): there exist α 1 , α 2 two positive constants such that

α 1 u τ -u τ ,h 2 Γ α 2 + u τ 2-s Γ + u τ ,h 2-s Γ ≤ Γ (Bu τ -Bu τ ,h ) • (Ku τ -Ku τ ,h ) dσ . (4.8)
Now since 1 < s < 2, we use (2.8), (4.8), triangle's inequality and (2.3) in (4.7) and get

2νc ∇(w h -u h ) 2 + α 1 u τ -u τ ,h ≤2ν ∇(w h -u) ∇(w h -u h ) + ∇(w h -u) q h -p h + ∇(u h -w h ) q h -p + c u τ -u τ ,h s-1 Γ u τ -w τ ,h Γ ≤c u τ -w τ ,h 2(s-1) Γ + c ∇(u -w h ) 2 + c q h -p 2 , (4.13)
which together with the triangle's inequality gives the asserted result.

Remark 4.1 It is manifest that the rate of convergence in Theorem 4.2 depends on the parameter s and the regularity of the weak solution on the slip zone Γ. This is different to the conclusion obtained in Theorem 4.1. Hence the relations (4.6), and (4.10) play a big role in the analysis of the convergence of the finite element solution (u h , p h ) of (3.2).

Remark 4.2 From (4.13), we deduce that

w τ ,h -u τ ,h Γ ≤c c + w τ ,h -u τ (2-s)/2 Γ + u τ (2-s)/2 Γ + u τ ,h (2-s)/2 Γ × u τ -w τ ,h s-1 Γ + c ∇(u -w h ) + c q h -p . Remark 4.3 (a)
Taking (w h , q h ) = (Π h u, r h p), the right hand sides of the error estimates obtained in theorem 4.1 or theorem 4.2 tend to zero as h tends to zero.

(b) Let (u, p) in H 2 (Ω) 2 × H 1 (Ω)
, we take (w h , q h ) = (Π h u, r h p) and we recall that (see [34, p.39]), there exists c such that 1) from theorem 4.2

u τ -w τ ,h Γ ≤ c u -w h 1/2 u -w h 1/2 1 . Then u -u h 1 + p -p h ≤ c    h min(3(s-1)/2,
h 3/4 from theorem 4.1 . (c) If u τ ∈ H 2 (Γ), then w τ ,h -u τ Γ ≤ ch 2 .
So assuming that the weak solution (u, p) is in H 2 (Ω) 2 × H 1 (Ω), we take (w h , q h ) = (Π h u, r h p), and we deduce that

u -u h 1 + p -p h ≤ c    h min(2(
s-1),1) from theorem 4.2 h from theorem 4.1 .

A priori error analysis: the Navier-Stokes system

We address the convergence of the finite element approximation (u h , p h ), solution of (3.7) by showing that the upper bound of the term uu h 1 + p -p h tend to zero when h tends to zero. This a priori error estimate is obtained with the assumption that the problems (3.7) and (2.13) are uniquely solvable. In this sub-section, we continue in the footsteps of Theorem 4.1, and Theorem 4.2. It holds that Theorem 4.3 Let (u, p) be the solution given by (2.13) and (u h , p h ) the solution of (3.7). Assume that

c f ≤ ν 2 .
Then there is a positive constant c independent of h such that for 1 < s < 2

u -u h 1 + u τ -u τ ,h Γ + p -p h ≤c inf w h ∈V h u τ -w τ ,h 1/2 Γ + c inf w h ∈V h ∇(u -w h ) + c inf q h ∈M h q h -p .
Proof.

Step 1: Inf-sup condition. We subtract (2.13) to (3.7), and taking into account the fact that V h is a subset of V, we arrived at

         ∀v h ∈ V h , b(v h , p h -q h ) = b(v h , p -q h ) + a(u -u h , v h ) + Γ |Ku τ ,h | s-2 Ku τ ,h -|Ku τ | s-2 Ku τ • Kv τ ,h dσ + d(u -u h , u, v h ) + d(u h , u -u h , v h ) .
an adequate stopping condition is satisfied,

(u n+1 h , p n+1 h ) solution of          for all (v h , q h ) ∈ V h × M h , a(u n+1 h , v h ) + Γ |Ku n τ ,h | s-2 Ku n+1 τ ,h • Kv τ ,h dσ -b(v h , p n+1 h ) = (v h ) , b(u n+1 h , q h ) = 0 .
(5.2)

After the convergence of the system (5.2), we get the numerical solution u h of the discrete system (3.2).

The variational problem (5.1) is linear, and from Babuska-Brezzi's theory for mixed problems (see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF][START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]), (u n+1 h , p n+1 h ) is well defined. Continuing on the qualitative analysis on (5.1), we claim that Lemma 5.1 (see [START_REF] Girault | Numerical discretization of a Darcy-Forchheimer model[END_REF]) The algorithm (5.1) is consistent with (3.7) in the sense that if (u 1 h , p 1 h ) = (u h , p h ) with (u h , p h ) solution of (3.7), then for all n ≥ 2, (u n h , p n h ) = (u h , p h ).

Proof. It is done by indiction on n. So, we assume that (u n h , p n h ) = (u h , p h ), and we will show that (u n+1 h , p n+1 h ) = (u h , p h ).

We use the induction hypothesis, take the difference between (5.1) and (3.7) and obtain that

             for all (v h , q h ) ∈ V h × M h , a(u n+1 h -u h , v h ) + |Ku τ ,h | s-2 Ku n+1 τ ,h -|Ku τ ,h | s-2 Ku τ ,h , Kv τ ,h Γ + d(u h , u n+1 h , v h ) -d(u h , u h , v h ) -b(v h , p n+1 h -p h ) = 0 , b(u n+1 h -u h , q h ) = 0 . (5.3)
We take v h = u h -u n+1 h and q h = p h -p n+1 h in (5.3), and obtain

a(u n+1 h -u h , u n+1 h -u h ) + Γ |Ku τ ,h | s-2 Ku n+1 τ ,h -Ku τ ,h 2 dσ = 0 from which we deduce that u n+1 h = u h . Equations (5.3) become for all v h ∈ V h , b(v h , p n+1 h -p h ) = 0 . (5.4)
We use the inf-sup condition on b(•, •) and one obtains that p n+1 h = p h .

In order to justify why we use the iterative scheme (5.1) (implemented in paragraph 5.2) to compute the solution of (3.7), we shall prove a convergence result. First we state that the sequence generated by Algorithm (5.1) is bounded. Proof. We take v h = u n+1 h in (5.1) and obtain

a(u n+1 h , u n+1 h ) + Γ |Ku n τ ,h | s-2 Ku n+1 τ ,h 2 dσ = (u n+1 h ) ≤ c f Du n+1 h which implies that ∇u n+1 h ≤ c ν f .
From the inf-sup condition, we obtain that

p n+1 h ≤ c f + c ν 2 f 2 .
We have the following abstract convergence result ) ∈ V h × M h the solution of (5.1). Let (u h , p h ) solution of (3.7).

Assume that there exists c positive constant independent of h such that f ≤ ν 2 c . Then

lim n→∞ u n h = u h weakly in V .
Proof. By virtue of the bound obtained in Lemma 5.2, there is a subsequence, still denoted by {u n h } n such that lim n→∞ u n h = u h weakly in H 1 (Ω) .

(5.5)

Clearly, we need to show that u h = u h , meaning that

   for all v h ∈ V div,h , a( u h , v h ) + |K u τ ,h | s-2 K u τ ,h , Kv τ ,h Γ + d( u h , u h , v h ) = (v h ) .
(5.6)

Actually, since V div,h is closed, it follows from (5.5) that u h ∈ V div,h . Owing to the compactness of the imbedding of H 1 (Ω) into L 4 (Ω), there exists a subsequence, still denoted by {u n h } n , which converges weakly to u h ∈ H 1 (Ω) and strongly in L 4 (Ω). For v h ∈ V div,h , passing to the limit on n in (5.1) is direct for the linear term because one only needs weak convergence, but for the term d(•, •, •), one needs the strong convergence of u n h on L 4 (Ω) (see [29, Chap 5, Thm 1.4]). In fact

d(u n k , u n+1 h , v h ) -d( u h , u h , v h ) = d(u n k -u h , u n+1 h , v h ) + d( u h , u n+1 h -u h , v h ) ≤ u n k -u h L 4 (Ω) ∇u n+1 h v h L 4 (Ω) + d( u h , u n+1 h -u h , v h ) -→ 0 .
Next, for the term on Γ, we deduce from the weak convergence of u n h on H 1 (Ω), that

u n h → u h strongly in L 2 (Ω) ,
which together with the continuity of the trace implies that

u n h | Γ → u τ ,h | Γ weakly in L 2 (Γ).
But from [31, Chap 2], we have that

lim n |Ku n τ ,h | s-2 Ku n+1 τ ,h , Kv τ ,h Γ = |K u τ ,h | s-2 K u τ ,h , Kv τ ,h Γ .
Hence, u h satisfy (5.6), which has a unique

u h if f ≤ ν 2 c . Then u h = u h
Since the solution of (5.6) is unique (under the condition f ≤ ν 2 c ), then the entire sequence converges weakly to u h .

Simulations

In this section we study the numerical behaviour of the algorithm (5.1) which is more general than the algorithm (5.2). The test problems used are designed to illustrate the numerical behaviour of the algorithm in two dimensions.

All the numerical simulations showed in this section are performed with FreeFem code [START_REF] Hecht | New development in FreeFem++[END_REF] by taking u 0 h = 0 and by using the following stopping criterion:

||u n+1 h -u n h || 1 ||u n+1 h || 1 ≤ 10 -5 . (5.7)
In all examples, V h and M h are given by (3.1).

Rate of convergence

In this first example, we consider the following geometry Ω = (0, 1) 2 where each edge is divided into N segments of equal length. Thus the corresponding mesh contains 2N 2 elements. We set ν = 1, K = I and f = (x 2 + y 2 , 1) T . We take the anisotropic slip boundary condition (1.7), on right and top boundaries of Ω, that is on Γ = {(x, y)| x = 1 or y = 1} = Γ 2 ∪ Γ 3 , and the homogeneous Dirichlet boundary condition u = 0 on the bottom and left boundaries S = Γ 1 ∪ Γ 4 .

We have s ∈ (1, 2), and we are interested in computing the rate of convergence of the finite element solution (u h , p h ). We recall that (see Theorem 4.3)

u -u h 1 + p -p h ≤ ch if (u, p) ∈ H 2 (Ω) × H 1 (Ω) , u τ ∈ H 2 (Γ) . (5.8)
As we do not have the exact solution u corresponding to this case, we will approximate it by computing the numerical solution of the iterative system (5.1) for N = 400 which will be designated as the reference solution say, (u ref , p ref ) and which depends on s ∈ (1, 2). Next, we compute the numerical solution for N = 10, 15, 20, . . . , 50. Fig. 1 show the graph of the relative error

E err = ||u ref -u h || 1 + ||p ref -p h || ||u ref || 1 + ||p ref ||
with respect to h = 1 N in logarithmic scale for s = 1.2, 1.5, 1.8. These lines have slopes 1.07, 1.05, 1.04 respectively. We can deduce for this particular case that the numerical slope is one and is independent of s. In all these cases, the numerical slopes satisfy the condition (5.8)

Flow in rectangle

In this second example, we consider the following geometry Ω = (0, l 1 ) × (0, l 2 ) with l 1 = 10, l 2 = 1. We divide the edges of the boundary ∂Ω into segments of equal lengths, and we define N = 40 to be the number of such mesh segments per unit length. We set ν = 1, f = 0, s = 3/2 or s = 4/3, and the tensor K one of this two cases:

K 1 = 1 0 0 2 and K 2 = 5 -1 -1 4 .
We consider the isotropic slip boundary condition (1.4), (1.7), on bottom (Γ 1 ) and top (Γ 3 ) boundaries of Ω, that is on Γ = {(x, y)| y = 0 or y = l 2 } = Γ 1 ∪ Γ 3 . A non-homogeneous Dirichlet boundary condition is imposed on the left (Γ 4 ) and right (Γ 2 ) boundaries of Ω, that is

S = {(x, y)| x = 0 or x = l 1 } = Γ 2 ∪ Γ 4 : u| Γ2 = u| Γ4 = y(1 -y) -y(1 -y) .
(5.9)

Figure 1: The a priori error estimates with respect to the mesh step h in logarithmic scale for s = 1.2, 1.5, 1.8.

Thus, the velocity on S has a parabolic profile with a critical point at y = 1/2 with u| S = (1/4, -1/4) T . Thus it is manifest that anisotropic slip is expected on the upper and lower boundaries. In Figure 2 we have shown the velocity profiles for s = 3/2 and K = K 1 and K = K 2 . On S, one notes the parabolic profiles of the velocity. The fluid along Γ and far away from S is unidirectional.

Figure 2: The velocity plots for s = 3/2, K = K 1 (up) and K = K 2 (down).

Figure 3 shows the graphs x → u 1 (x, 1/2) and x → u 2 (x, 1/2) for s = 3/2 and K = K 1 or K 2 . One notes that u 1 (resp. u 2 ) is equal to 1/4 (resp. -1/4) on x = 0 (resp. x = 10) which is coherent with the boundary condition (5.9) for y = 1/2. Furthermore, Figure 4 plots the tangential velocity on Γ 1 and Γ 3 (u τ (x, 0) and u τ (x, 1)) for s = 3/2 and K = K 1 .

Next, Figure 5 show the velocity plots on the line y = 1/2 for K = K 1 and s = 4/3 or 3/2. We note the dependence on s of the velocity on the line y = 1/2.

Driven cavity flow

In this paragraph, we consider the numerical results for the Lid Driven cavity with the Navier-Stokes flows under the power law slip boundary condition by the numerical scheme (5.1). This problem is very popular example that can be seen in [START_REF] Owens | Computational Rheology[END_REF][START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF][START_REF] Djoko | Numerical methods for the Stokes and NavierStokes equations driven by threshold slip boundary conditions[END_REF]. The fluid is confined in the domain Ω = (0, 1) 2 and the boundary condition for the velocity field on the top is given by u = (16 x 2 (x -1) 2 , 0) T on y = 1 0 < x < 1 .

Thus one has a smooth start of the flow, and the lid velocity attains its maximum, u = (1, 0) T , at the centre x = 1/2. This boundary condition has the effect of avoiding local singularities at the top-right and top-left corners as it ensures both the velocity and velocity gradient vanish at the corners. On the other parts of the boundary, we consider the isotropic slip boundary condition. In all the numerical results considered in this section, we set K = I, f = 0 and ν = 1/Re where Re designates the Reynolds number.

Here also, each edge is divided into N = 100 segments of equal length. Thus the corresponding mesh contains 2N 2 elements. We begin ours numerical results with Re = 1. Figures 6 and7 shows the numerical velocity for s = 1.1, 1.5 and 1.9. One notes the formation of a vortex which is typical for this problem. Furthermore, we show in figure 8 the first and second components of the velocity on the line y = 1/2 for s = 1.1, 1.5 and 1.9. Again, one observes the dependence on s of the velocity on the line y = 1/2. Next, we take s = 3/2 and Re = 1000 and 2000. Figure 9 shows the repartition of the velocity for the corresponding Reynolds numbers. Thus, our algorithm converges when the Reynolds number

Conclusions

The goal of the present work is to propose a novel formulation for the Stokes and Navier-Stokes equations with a direction dependent slip on a part of the boundary. The variational formulations can be regarded are perturbations of classical variational formulations encountered with standard boundary conditions. In this work we have -Constructed the weak solutions of the Stokes and Navier-Stokes equations with power law slip boundary condition by using the same functional framework as for the standard boundary conditions.

-Formulated and studied the finite element approximations associated to these problems. In particular we have derived optimal convergence when the tangential part of the velocity on the slip zone has enough regularity.

-Constructed and studied an iterative scheme associated to the equations of Navier-Stokes, and derived some qualitative properties.

-The numerical experiments reveal optimal convergence rate, which is uniform with respect to the power law parameter s. As expected, one can observe the dependence of the velocity with respect to the power law parameter s in all examples treated. Finally, the convergence of the iterative scheme is observed for big Reynolds number.

The current methodology should be further extended, in particular by (i) devising an adaptive mesh strategy based on a posteriori error estimate; (ii) devising an adapted discretization approach for the case of a curved boundary and application in computational hemodynamics; (iii) implementation of the current problems by combining alternative directions method of multiplier(s) and augmented Lagrangian formulations; (iv) extending the present analysis to generalised non-Newtonian models; (v) extending the present analysis to discontinuous Galerkin methods and virtual element approach; (vi) extending the present analysis to time dependent problem.

ΩT:

  Dudx = 2ν Ω |Du|2 dx being the dissipation mechanisms in bulk, while Γ

. 8 )

 8 and (1.2)-(1.4), (1.7) are unchanged, and (v • ∇)v = d i=1 v i ∂v ∂x i is the convection term. The solvability of the Stokes and Navier-Stokes equations when (1.3), (1.4) and (1.7) are replaced by u • n = 0 and (T n) τ + H(u τ )u τ = g on ∂Ω ,

  and there is a constant c such that |v(x)| ≤ c a.e on Ω , with the associated norm

Remark 2 . 1 Lemma 2 . 1

 2121 In some contributions, the conservation of the momentum is expressed in terms of the Laplacian of the velocity which gives rise to the bilinear form a(u, v) = 2ν (∇v, ∇u), instead of a(•, •) defined before. Although at a continuous level div u = 0 implies that 2 div Du = ∆u, from a modelling viewpoint it may be important to work with symmetric tensor. For instance, the problem (1.1)-(1.4) and (1.7) gives directly the natural boundary condition in term of the force (traction force) exerted by the fluid on its boundary. Now with standard arguments, we can show the following equivalence result. The variational formulation (2.4) is equivalent in the sense of distribution to (1.1)-(1.4), (1.7).

Lemma 5 . 2

 52 Let (u n+1 h , p n+1 h) the solution of the algorithm (5.1). Then there exists a constant c independent of h such that ∇u

Theorem 5 . 1

 51 Let (u n+1 h , p n+1 h

Figure 3 :Figure 4 :

 34 Figure 3: The numerical velocity (u 1 right, u 2 left) on the line y = 1/2 for s = 3/2, K = K 1 and K = K 2 .

Figure 5 :

 5 Figure 5: The numerical velocity (u 1 right, u 2 left) on the line y = 1/2 for K = K 1 , s = 4/3 or s = 3/2.

Figure 6 :

 6 Figure 6: The numerical velocity for s = 1.1 (left) and s = 1.5 (right).

Figure 7 :

 7 Figure 7: The numerical velocity for s = 1.9.

Figure 8 :

 8 Figure 8: The numerical velocity (u 1 , u 2 ) for s = 1.1, 1.5 and 1.9.

Figure 9 :

 9 Figure 9: The numerical velocity for s = 3/2, Re = 1000 (left) and Re = 2000 (right).

Γ α 2 + u τ 2-s Γ + u τ ,h 2-s Γ

Γ α 2 + w τ ,h 2-s Γ + u τ ,h 2-s Γ ≤ Γ (Bw τ ,h -Bu τ ,h ) • (Kw τ ,h -Ku τ ,h ) dσ .(4.12)Now since 1 < s < 2, we use (2.8), (4.12), triangle's inequality and (2.3) in (4.11) and get2νc ∇(w h -u h ) 2 + α 1 w τ ,h -u τ ,h 2 Γ α 2 + c w τ ,h -u τ 2-s Γ + c u τ 2-s Γ + u τ ,h 2-s Γ ≤c Γ |Ku τ -Kw τ ,h | s-1 |Kw τ ,h -Ku τ ,h | dσ + 2ν ∇(u -w h ) ∇(w h -u h ) + ∇(w h -u) q h -p h + ∇(u h -w h ) q h -p ≤c u τ -w τ ,h s-1 Γ w τ ,h -u τ ,h Γ + 2ν ∇(u -w h ) ∇(w h -u h ) + ∇(w h -u) q h -p h + ∇(u h -w h ) q h -p ,from which after utilization of (4.3), the trace's inequality, and young's inequality leads to∇(w h -u h ) 2 + c 1 w τ ,h -u τ ,h 2 Γ c + c w τ ,h -u τ 2-s Γ + c u τ 2-s Γ + c u τ ,h 2-s Γ

Acknowledgments. The authors thank the referees for their constructive remarks and comments.

conflict of interest statement. None

from which after utilization of (4.3), the estimates in Theorem 2.1, Theorem 3.1, and young's inequality leads to

(4.9)

Hence the result announced is obtained after utilization of triangle's inequality and (4.9).

From Theorem 4.1, one concludes that the a priori error is independent of s, dominated by the interpolation error on the slip zone Γ, and depends on the regularity of the weak solution on the slip zone Γ. In the next result, we demonstrate that the a priori error may depend on s if one uses the monotonicity of B differently. Theorem 4.2 For 1 < s < 2, let (u, p) be the solution of (2.4) and (u h , p h ) the solution of (3.2). Then there is a positive constant c independent of h such that

Proof. We follow the proof of Theorem 4.1.

We have

(4.10) Thus from (4.10), inserting the second equation of (4.4) into the first one yields From (2.7), it was shown that (see [START_REF] Sandri | Sur l approximation des ecoulements numeriques quasi-Newtoniens dont la viscosite obeit a la loi de puissance ou de Carreau[END_REF]): there exist α 1 , α 2 two positive constants such that

Using Cauchy-Schwarz's, Hölder's and trace's inequality, one gets

which after inserting it in (4.1) and making use of (2.8), and the estimates in proposition 2.2 and proposition 3.1, one has

Step 2: Monotonicity. We let Bu τ = |Ku τ | s-2 Ku τ . We subtract (2.13) to (3.7) and obtain [START_REF] Najeeb Alam Khan | Effect of anisotropic slip and magnetic field on the flow and heat transfer of eyring powell fluid over an inifinite rotating disk[END_REF], and after using linearity of operators involved, one gets

(4.17) Now since 1 < s < 2, then using (2.8) and (2.3) in (4.17), one obtains

from which after utilization of (4.14), (4.12), the trace's inequality, and young's inequality leads to

which together with the triangle's inequality gives the result announced.

Following the proof of Theorem 4.2, we can show that Theorem 4.4 Let (u, p) be the solution given by (2.13) and (u h , p h ) the solution of (3.7). Assume that

Then there is a positive constant c independent of h such that for 1 < s < 2

The next section is concerned with the solution strategy for (3.2) and (3.7), and their numerical simulations.

Iterative schemes and Simulations

The finite element problems (3.2) and (3.7) are nonlinear. Hence iterative/incremental schemes need to be formulated for their effective resolution. It should be noted that since s is in (1, 2), the expression (|Ku τ ,h | s-2 Ku τ ,h , Kv τ ,h ) Γ has a poor differentiability. So a numerical solution of (3.2) or (3.7) with Newton's method or conjugate gradient may not converge at all. A direct strategy to make the above systems less nonlinear and weaken the nonlinearities in the systems is to linearize the nonlinear terms by considering the approximations

In the lines that follows, we concentrate on the iterative method for solving (3.7) since the nonlinear system of equations (3.2) can be viewed as a special case of (3.7).

Iterative scheme for (3.7)

We solve (3.7) by adopting the following strategy. Let u 0 h ∈ V h be a given initial guess, we generate the sequence (u n , p n h ) for n ≥ 1, until an adequate stopping condition is satisfied, by computing

, q h ) = 0 .

(5.1) After the convergence of the system (5.1), we get the numerical solution u h of the discrete system (3.7). One of the principal interests of the scheme (5.1) (which by the way is a fixed point algorithm) is that it converges to the unique fixed point for any starting point, but that convergence is slow.

Remark 5.1 Similar to the iterative scheme (5.1), the iterative scheme we may consider for (3.2) is the following: Let u 0 h a given element of V h . For n ≥ 1 and having u n h in V h , we compute until


V Scesa

P Henaff

F B Ouezdou

F Namoun

Time Window Width Influence On Dynamic BPTT(h) Learning Algorithm Performances: Experimental Study

The purpose of the research addressed in this paper is to study the influence of the time window width in dynamic truncated BackPropagation Through Time BPTT(h) learning algorithms. Statistical experiments based on the identification of a real biped robot balancing mechanism are carried out to raise the link between the window width and the stability, the speed and the accuracy of the learning. The time window width choice is shown to be crucial for the convergence speed of the learning process and the generalization ability of the network. Although, a particular attention is brought to a divergence problem (gradient blow up) observed with the assumption where the net parameters are constant along the window. The limit of this assumption is demonstrated and parameters evolution storage, used as a solution for this problem, is detailed.

Introduction

Born from the collaboration between a robotic research laboratory (LISV) and an industrial company (BIA), the supporting project of this study aims to shape a smart architecture able to learn to control non linear multi actuators system under real time constraints. Through the design and the implementation of this controller, we want to evaluate the ability of neural architectures to achieve the non linear control needs in real time. The project experiments are carried out on two structures: the ROBIAN biped robot from the LISV [START_REF] Mohamed | Emulation of the dynamic effects of human torso during walking gait[END_REF], and the BIA road simulator (www.bia.fr). The architecture of the developed algorithm is based on continuous time recurrent neural networks (CTRNN). To shape this architecture, a truncated BPTT algorithm was chosen for its ability of integrating the learning error and for its simplicity to be implemented for real time online applications. This well known gradient algorithm is widely used for system modeling and control [START_REF] Tsung | Modeling Dynamical Systems with Recurrent Neural Networks[END_REF], optimization applications, speech recognition [START_REF] Nguyen | Tau Net: A neural network for modeling temporal variability[END_REF], or meta learning [START_REF] Hochreiter | Learning to learn using gradient descent[END_REF]. A detailed description of this algorithm is given in [START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF], [START_REF] Werbos | Backpropagation through time: what it does and how to do it[END_REF] and [START_REF] Rumelhart | Learning internal representations by error propagation Parallel distributed processing: explorations in the microstructure of cognition[END_REF]. Nevertheless, as far as we know, no study aimed to extract the important role of the truncation width on the success of the learning. Obviously, the truncation width influences deeply the learning speed and quality. Our real time experimentations on ROBIAN biped also show us that, when the net parameters are considered as constant along a large window, the learning stability could be hardly spoiled.

Thus, through a statistical analysis of learning results, the links between the width of the time window, the stability of the learning, the "constant parameters" assumption, the convergence speed and the generalization abilities of the learned networks, have to be raised. To focus on this influence, a first experiment dealing with direct model identification rather than a controller shaping will be considered. The next section details the neural model and the BPTT(h) learning algorithm used. The experimental plant is described in section 3. In section 4, the influence of the window width will be studied. The "constant parameters" assumption with its influence on the learning stability will be discussed. A stabilizing modification of the learning equation, to take into account the parameters evolution history, will be proposed. It will allow raising the dependencies between the window width and the learning results. Finally, discussions about the parameters evolution storage and the window choice to find a compromise between learning speed and nets accuracy will be given.

2 Dynamic truncated Backpropagation Through Time algorithm

Neural Model

Following the classical CTRNN equation ( 1), the input data are propagated in the network to generate the neurons outputs. The net activities belong to the intrinsic neurons and network parameters : weights (w ij ), biases (b j ) and time parameters (T j ).

[

]

      + ⋅ + - = ∂ ∂ ⋅ ∑ i j i ij j j j b y w f y t y T
Where y j corresponds to the j neuron activity, f is the activation function (tanh). Using a time scale parameter S j ( classically S j =∆t/T j , ∆t is the time step ), the discrete corresponding equation can be written as follows:

( ) [ ] ( ) ( ) ∆t t y S b ∆t t y w f S t y j j j i i ij j j - ⋅ - +       + - ⋅ ⋅ = ∑       1
As the network is fully recurrent, each neuron receives the outputs of all the other ones. The weights, biases and time constants shape the temporal response of the net.

Learning algorithm

The objective of the parameters modification consists in minimizing a desired criterion. For an identification process, the criterion would be the gap between the neural model and the taught system. To carry out the adaptation of network parameters, BackPropagation Through Time algorithm, detailed in [START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF], [START_REF] Werbos | Backpropagation through time: what it does and how to do it[END_REF] and [START_REF] Rumelhart | Learning internal representations by error propagation Parallel distributed processing: explorations in the microstructure of cognition[END_REF], can be used. This algorithm is computing an error function that corresponds to the criterion to be minimized. The error function (E) is defined as the integral of each net output errors. The parameters modification (∆param) is leaded by the error gradient inverse value. The following equations give the error function and the parameters modification law:

(1) (2) ( ) [ ] ∫ ⋅ = t t j d e E 0 τ τ (3) param E ∆param ∂ ∂ ⋅ - = η (4)
Where e j (τ) is the output error of the neuron j stored at the current time step in the time window (τ). During the learning, to minimize E, the algorithm modifies the net parameters following the gradient descent (4), where param is w ij , b j or S j and η is the corresponding learning rate. To compute these "delta" values for continuous time neurons, a dynamic BPTT learning algorithm [START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF] is needed. The gradient values giving the influence of each parameter are computed with equations ( 5), ( 6) and ( 7):

( ) ( ) ( ) ( ) [ ] ∫ ⋅ - ⋅ - ′ ⋅ ⋅ = ∂ ∂ t t j k k k jk d ∆t y ∆t x f z S w E 0 τ τ τ τ ( ) ( ) ( ) [ ] ∫ ⋅ - ′ ⋅ ⋅ = ∂ ∂ t t k j j j d ∆t x f z S b E 0 τ τ τ ( ) ( ) ( ) ( ) ( ) [ ] ∫ ⋅ - - - ⋅ = ∂ ∂ t t j j j j d ∆t y ∆t x f z S E 0 τ τ τ τ
The backpropagated errors ( z j (τ) ) for a j th neuron can be written as following:

( ) ( ) ( ) ( ) ( ) [ ] ( ) ( ) ( ) ∆t e S ∆t z w x f S ∆t z y E z j j j k jk k k k j j ⋅ + - ⋅ + + ⋅ ′ ⋅ ⋅ + = ∂ ∂ = ∑ τ τ τ τ τ τ 1
Where f' is the derivative of the activation function (tanh). For each neuron, the z j is computed by considering the network states at each time step as successive layers. Thus, the output errors are backpropagated in the network and in time. In BPTT learning, z j is computed for each neuron and each time step since the starting. To prevent a memory explosion induced by the storage of every network states, a truncated BPTT algorithm (BPTT(h)) is proposed in [START_REF] Williams | Gradient-based learning algorithms for recurrent connectionist networks[END_REF]. This algorithm is keeping in memory only the past states included in a window following the current instant (t 0 = t -h with h the window width). The z j values are computed along this sliding window (τ∈[t-h;t]). Hence, only a gradient approximation is computed. The previous [START_REF] Pearlmutter | Gradient calculation for dynamic recurrent neural networks: a survey[END_REF] to [START_REF] Williams | Gradient-based learning algorithms for recurrent connectionist networks[END_REF] equations are constructed on the assumption that the parameters are constant along the time window. This assumption presented in [START_REF] Williams | An efficient gradient--based algorithm for on--line training of recurrent network trajectories[END_REF] and [START_REF] Campolucci | On-Line Learning Algorithms for Locally Recurrent Neural Networks[END_REF] minimizes the process memory needs as the parameters history is not stored. In our study, these algorithms are experimented on the ROBIAN biped robot described in the next section.

ROBIAN identification description 3.1 Plant description and perturbation signals

To focus on the time window width influence, a model identification rather than a controller shaping will be carried out. Fig. 1 shows how the identification of ROBIAN's torso influence on the ZMP (Zero Moment Point) will be performed. For more details concerning the ROBIAN biped and its torso mechanism see [START_REF] Mohamed | Emulation of the dynamic effects of human torso during walking gait[END_REF].

(6) (7) (8) (5)
The ZMP is a key notion in the balance control of walking robot. It corresponds to the position of the center of pressure on the ground. The ZMP algorithm, introduced thirty five years ago [START_REF] Vukobratovic | Zero-moment point -thirty five years of its life[END_REF], consists of controlling the equilibrium of the biped robot by keeping the ZMP inside the polygon defined by the contact points with the ground. The motion of the X and Y masses of the torso (Mx and My) are perturbing the robot balancing, leading to variations of the X and Y ZMP positions. During the learning process, the neural network (NN) computes the gap between the measured positions and its own outputs and modifies its parameters to vanish this error, by minimizing a cost function. This function is defined as the normalized squared errors sum on the two axes and the output errors (e j (t)) can be expressed as following:

Where y x , y y are the output activities and maxZMP axis is the maximum amplitude on each axis (for normalization). Representative situations should be given during the learning. They must fully characterize the behavior of the studied system. This implies to give the network patterns that express the dynamic of the ZMP for the X and Y axis. frequencies (1Hz-3Hz, containing the resonance of the robot), for both X and Y axes. In the test pattern, a sine command with a varying frequency is applied on both axes.

Example of learning

The cost function evolution is depicted on Fig. 3. In this experiment, the network (composed of 2 inputs, 15 hidden neurons and 2 outputs) was taught with one hundred loops on the learning pattern (learning rates: η 1 = 0.025 for weights and biases, η 2 = 1.25 for time constants) and a time window width h=20. At the beginning of the learning, the network parameters are randomly chosen to generate stable outputs. The time scale parameters (S j ) are taken in the interval [0.5;1], to obtain neurons with speeds close to the dynamic of the system. Fig. 3 gives the evolution of the cost during the learning. In the zoom area, the four peaks are due to the different kinds of command (amplitude changes: peaks 1 and 4, or frequency variation: peaks 2 and 3) of the learning pattern.

During the experiment, two main drops happened. They are respectively linked to a quick decrease of peak 3 (learning of the X direction behavior) and peak 2 (learning on Y direction). The instant when the second drop happened (SDI : Second Drop Instant) is a good factor to quantify the convergence speed. It happens when the network identifies correctly the behavior of ROBIAN torso in the both two directions.

Here it takes about 64 loops to arise. The SDI value can be interpreted as a physical expression of the convergence for our experiments. It can be related to the choice of a threshold below which the cost is considered as satisfactory.

After the learning, a relevant value for estimating the generalization ability of the learned network is the costs sum along the whole test pattern (TES : Test Errors Sum). It expresses the difference between an ideal model of the robot and the learned one. With a perfect model, the TES value is equal to zero. For the considered example, before the learning, TES = 363, and at the end of the learning process (after the 100 th loop), the remaining error decreases to TES = 40. In the usual gradient calculus, the parameters are considered as constant along the time window [START_REF] Williams | An efficient gradient--based algorithm for on--line training of recurrent network trajectories[END_REF], [START_REF] Campolucci | On-Line Learning Algorithms for Locally Recurrent Neural Networks[END_REF]. This approximation can be seen in the gradient equations ( 5) to [START_REF] Williams | Gradient-based learning algorithms for recurrent connectionist networks[END_REF] where the parameters (w ij , b j or S j ) are not indexed by time. The time window width reduces the use of this assumption by destabilizing the learning. The identification of the ROBIAN's torso influence on the ZMP positions is carried out for different time window width values (TW) from 5 to 60 states stored. For each TW, ten learning courses are carried out for a statistical analysis of the results. Each network is evaluated on the test pattern after learning. The learning and test results are depicted on the following figure: Here, it happens for TW=35 where there was only 10% of success This graph can be decomposed in three parts. When TW is less than 15, the learning is so slow that it doesn't manage to converge (no results on the graph). For TW∈[15;35], the SDI, and TES seem to decrease proportionally to TW, thus the convergence speed and the generalization ability are better for larger TW. Finally, if TW is greater than 35 states stored, the algorithm is diverging with a gradient blow up. For a stronger learning rate (η=0.1) the gradient blow up occurs since TW=15.

Gradient blow up divergence and parameters evolution storage

We called this divergence problem gradient blow up since it's caused by a sudden explosion of the neurons gradient values. In Fig. 5, the evolution in time of the gradi-ent values for each neuron are depicted. In this experiment, the time window width is TW=100 states stored (with a sampling rate=45Hz, it corresponds to an interval of 2.22s). The three graphs represent the initiation of the gradient blow up. Just after the last graph, all the gradient values are exploded. The blow up begins with a divergence of the oldest past values. Through the recurrent links raised in the network, the explosion is spread in the entire network. Finally, the divergence is so important that the parameters take infinite values and the algorithm is stopped. Our guess is that the divergence of the farer gradients is due to the "constant parameter" assumption. This assumption means that the parameters are the same along the entire time window. However, if the window width is large compared to the parameters variation, it's obviously false. Actually, if the width is important, the parameter values that contribute to create the farer past states stored could be strongly different from the current ones. This happens when either the first or the second drop is initiated, i.e. when the parameters are strongly modified. That's why the divergence is initiated in the oldest part of the time window. Hence, the assumption could only be used with short window or small learning rates. A way to avoid this blow up problem is to store the parameters history along the time window. So, the gradient equations will take into account their evolutions as in (5'), (6') and (7'). The weights, biases and time scale parameters become a function of time τ. The storage of the parameter history along the sliding time window implies an increase of memory needs. The number of stored values is multiplied by TW. Equation ( 7) is not changed since no parameter is involved. This modification of the classical BPTT learning equations will be used for the analysis of TW influence on the speed and accuracy of the learning in next section.

( ) ( ) ( ) ( ) [ ] ∫ ⋅ - ′ ⋅ ⋅ - = ∂ ∂ t t k j j j d ∆t x f z S b E 0 1 τ τ τ τ (6') (8') ( ) ( ) ( ) ( ) [ ] ( ) ( ) ( ) ∆t e S ∆t z w x f S ∆t z z j j j k jk k k k j ⋅ + - ⋅ + + ⋅ ′ ⋅ ⋅ + =                   ∑ τ τ τ τ τ τ τ τ 1 Gradient at time T
Gradient at time T+0.02s Gradient at time T+0.04s The previous graph can be decomposed in three parts depending on the TW value. If TW<15 states stored, the algorithm is not able to converge to a correct solution. Next, for TW varying between from 15 to 35, the convergence speed average is increased while the average TES is decreased. Then, with TW>35, the convergence process is more oscillating. The standard deviations are thus larger, and not all the 10 initial random networks lead to a correct solution. But, for these TW values, in case of convergence, the algorithm finds better solutions (i.e. the TES values are smaller). This is due to the amount of data taken into account. The modified algorithm allows also increasing the learning rates. Figure 7 represents the results obtained with η 1 =0.1: Fig. 7. Learning and test results with η 1 =0.1 for weights and biases and η 2 =1.25 for time constants. All the learning succeed for TW≤35, for TW= 40, 45, 50 and 55, the percentages of success are respectively 70%, 40%, 40%, 30% and 0% for 55 Compared to the experiments done with η 1 =0.025, the convergence speed gets faster since a TW = 15 states stored. Nevertheless, the learning is more unstable with this higher learning rate, and the algorithm meets convergence difficulties earlier. For larger TW values, the number of success is less important as the learning is oscillating, due to a too important learning rate. But, when they converge, the networks learned are more accurate. Here again, the TES obtained are the best for the biggest TW. As far as our experiment is concerned, the best configuration ensuring maximum speed, best quality and high percentage of success, is TW=30 and η 1 =0.1.

( ) ( ) ( ) ( ) ( ) [ ] ∫ ⋅ - ⋅ - ′ ⋅ ⋅ - = ∂ ∂ t t j k k k jk d ∆t y ∆t x f z S w E 0 1 τ τ τ τ τ (5')

Discussion and conclusion 5.1 Storing the parameters evolution?

The comparison between the results of the "constant parameters" assumption and those presented for the modified algorithm, shows that the second one is better for our identification experiment. First, the classical one is only valid for small TW values whereas the second one allows larger ones. Next, for the same TW range, the classical method finds networks with worst generalization abilities. The only advantage of the classical one seems to be a faster convergence speed. Nevertheless, as a stronger learning rate is not prohibited with the modified method, this speed advantage can be overcome. In that case, the modified algorithm results are comparable to the ones achieved with a smaller learning rate and a larger TW value. With the modified algorithm, for the largest TW values, the convergence is not always met. The algorithm fails to find a correct solution. But no gradient blow up occurs. The TW limit that leads to a gradient blow up for the classical method is probably linked to the dynamic of the studied system and the convergence speed defined by the learning rates values. If the learning rates are strong, the parameters modification will be fast and they couldn't be approximated as constant along the time window. In the same way, if the system is fast, the learning will tend to bring the time constants and the network to a faster dynamic. Thus, the errors will be quickly backpropagated in the time window, and the algorithm will behave as if the time window was bigger and will diverge for smaller TW. So, the choice between using or not the "constant parameters" assumption will depend on the dynamic of the system.

Convergence speed VS generalization ability

The comparison of the convergence speeds and the quality of the learned networks demonstrates that the TW value choice must be a trade-off between these two results. Choosing the largest one, i.e. the best quality, the speed could be strongly decreased or the convergence not guaranteed. Choosing the fastest convergence, with a shorter TW, could lead to non suitable networks. There is no general rule for optimizing the learning. The choice must be done following the needs and criteria fixed for the learning process. Here again, the dynamic of the studied system will act upon the TW that leads to the fastest convergence. For the same criteria, the chosen TW could be different for a slow or fast system. But, allowing the use of larger TW is obviously useful.

Conclusions

In this paper, we studied the influence of the time window width parameter upon the stability, the learning speed and the quality of the networks obtained. Based on statistical experiments, aiming to identify the links between the balancing of a biped robot and its torso motion, we discussed the net parameters evolution storage and the choice of an optimal TW value. We found that it could be useful to enlarge the window to reach faster and more accurate learning. As the "constant parameter" assumption classically adopted for BPPTT(h) is not adapted for larger window, a modified algorithm taking into account the parameters history can be advantageous. In the future, we will first carry out similar experiments on faster or slower systems to find out the influence of the system dynamic on the learning. Next, we will perform a mathematical analysis of the modified algorithm. The learning algorithm will be also applied to carry out inverse model identification and control.

Fig. 1 .

 1 Fig. 1. Learning architecture for the identification (left) of ROBIAN robot (right)

  Fig 2 represents the learning pattern adopted.

Fig. 2 .

 2 Fig.2. Masses positions and ZMP positions on the two axis for the learning pattern (sampling rate = 45Hz). In periods 1 and 4, the X and then the Y mass, are submitted to successive steps with various amplitudes. In the 2 and 3 periods, steps with varying frequency are appliedThe learning pattern contains a succession of four different excitations applied to the robot: squared commands with varying amplitudes (0.01m to 0.20m) and with varying

Fig. 3 .

 3 Fig. 3. Cost evolution, and a zoom on the first loop (top right corner), 1 loop =10200 iterations

1

 1 Stability limit with the "constant parameters" assumption

Fig. 4 .

 4 Fig. 4. Learning and test results following TW with η 1 =0.025 for weights and biases, η 2 =1.25 for time constants. The average value of SDI and TES are plotted. The error bars correspond to the standard deviation measured. They express the iteration range of convergence. If not all the 10 trials lead to a convergence, a mark is added on the graph giving the percentage of success.Here, it happens for TW=35 where there was only 10% of success

Fig. 5 .

 5 Fig. 5. Gradient blow up initiation. The little waves correspond to the "normal" backpropagation of the error in the time window on the neurons. The strong drops and jumps in the third graph correspond to the "abnormal" divergence called gradient blow up

4. 3

 3 Influence on the convergence speed and on the accuracy of the learned nets The same identification process is carried out with the modified learning algorithm, for time window width values from 5 to 60 states stored, and with 10 initial random nets for each TW value. The learning and test results are depicted on the fig 6:

Fig. 6 .

 6 Fig. 6. Learning and test results following the TW values with η 1 =0.025 and η 2 =1.25. All the learning succeed for TW≤45, for TW= 50 and 55, only 40% of the nets reach a convergence