Hyam Abboud
email: hyamabboud@usek.edu.lb

Toni Sayah
email: tsayah@fs.usj.edu.lb

Upwind discretisation of a time-dependent two dimensional grade-two fluid model

Keywords: Grade-two fluid, non-linear problem, incompressible flow, time and space discretizations

In this paper, we propose a finite-element scheme for solving numerically the equations of a transient two dimensional grade-two fluid non-Newtonian Rivlin-Ericksen fluid model. This system of equations is considered an appropriate model for the motion of a water solution of polymers. As expected, the difficulties of this problem arise from the transport equation. As one of our aims is to derive unconditional a priori estimates from the discrete analogue of the transport equation, we stabilize our scheme by adding a consistent stabilizing term. We use the IP 2 -IP 1 Taylor-Hood finite-element scheme for the velocity v and the pressure p, and the discontinuous IP 1 finite element for an auxiliary variable z. The error is of the order of h 3/2 + k, considering that the discretization of the transport equation loses inevitably a factor h 1/2 .

Introduction

This article is devoted to the numerical solution of the equation of a grade two fluid non-Newtonian Rivlin-Ericksen fluid ( [START_REF] Rilvin | Stress-deformation relations for isotropic materials[END_REF]) :

∂ ∂t (u -α∆u) -ν∆u + curl(u -α∆u) × u + ∇p = f in ]0, T [×Ω, (1)
with the incompressibility condition :

div u = 0 in ]0, T [×Ω, (2)
where the velocity u = (u 1 , u 2 , 0),

div u = ∂u 1 ∂x 1 + ∂u 2 ∂x 2 , curl u = (0, 0, curlu), curl u = ∂u 2 ∂x 1 - ∂u 1 ∂x 2 ,
here f denotes an external force, ν the viscosity and α is a constant normal stress modulus.

This model is considered an appropriate one for the motion of water solutions of polymers ( [START_REF] Dunn | Rajagopal Fluid of differential type: Critical review and thermodynamic analysis[END_REF]).

The case α = 0 represents the transient Navier-Stokes problem. Here, p is not the pressure, but the formula which gives the pressure from u and p is complex. To simplify, we refer to p as the "pressure" in the sequel. According to Dunn and Fosdick's work [START_REF] Dunn | Fosdic Thermodinamics, stability, and boundness of fluids of complexity two and fluids of second grade[END_REF], in order to be consistent with thermodynamics, a grade-two fluid must satisfy α ≥ 0 and ν ≥ 0. The reader can refer to [START_REF] Dunn | Rajagopal Fluid of differential type: Critical review and thermodynamic analysis[END_REF] for a discussion on the sign of α.

The equations of a grade two fluid model have been studied by many authors (Videmann gives in [START_REF]Videmann Mathematical analisys of visco-elastic non Newtonian fluid[END_REF] a very extensive list of references), but the best construction of solutions for the problem, with homogeneous Dirichlet boundary conditions and mildly smooth data, is given by Ouazar [START_REF] Ouazar | 3ème cycle de l'Université Pierre et Marie Curie Paris VI[END_REF] and by Cioranescu and Ouazar [START_REF] Cioranescu | Existence et unicité pour les fluides de second grade[END_REF], [START_REF] Cioranescu | Existence and uniqueness for fluids of second grade, Nonlinear Partial Differential Equations[END_REF]. They prove existence of solutions, with H 3 regularity in space, by looking for a velocity u such that

z = curl(u -α∆u)
has L 2 regularity in space, introducing z as an auxiliary variable and discretizing the equations of motion by Galerkin's method in the basis of the eigenfunctions of the operator curl curl(u -α∆u). This choice allows one to recover estimates from the transport equation in two dimensions

α ∂z ∂t + νz + αu • ∇z = ν curl u + α curl f , (3)
whenever curl f belongs to L 2 (Ω) 3 . In this case, z = (0, 0, z) with z = curl(u -α∆u). Hence, z is necessarily orthogonal to u.

In this article, we propose finite-element schemes for solving numerically the equation of a two dimensional grade-two fluid model. Defining z as above, the equations of motion becomes :

∂ ∂t (u -α∆u) -ν∆u + z × u + ∇p = f , (4)
and α ∂z ∂t

+ νz + αu • ∇z = ν curlu + α curlf , (5)
the Dirichlet boundary condition :

u = 0 on ]0, T [×∂Ω, (6)
and the initial conditions : u(x, t) = 0, and z(x, t) = 0.

This problem is analyzed by Girault and Saadouni [START_REF] Girault | On a time-dependent grade-two fluid model in two dimensions, Computers and Mathematics with Applications[END_REF]. If we want to derive an unconditional a priori estimate for the discrete analogue of (3), we add to the left-hand side of this last equation a stabilizing, consistent term, so it becomes

α ∂z ∂t + νz + αu • ∇z + α 2 (div u)z = ν curl u + α curl f . (8)
In this work, we propose to discretize this last equation, as Girault and Scott did in [START_REF] Girault | Upwind discretizations of a steady grade-two fluid model in tow dimensions[END_REF], by an upwind scheme based on the discontinuous Galerkin method of degree one introduced by Lesaint and Raviart in [START_REF] Lesaint | On a finite element method for solving the neutron transport equation, Mathematical aspects of finite Elements in Partial Differential Equation[END_REF]. Let X h , M h and Z h be the discrete spaces for the velocity and the pressure. We approximate the velocity and the pressure by the standard IP 2 -IP 1 Taylor-Hood scheme, where IP k denotes the space of polynomials of degree k in two variables. Also, in each element of the triangulation, z n h is a polynomial of degree one, without continuity requirement on interelement boundaries. Our discrete system corresponding to (4) and ( 8) is :

Find u n+1 h ∈ X h , p n+1 h ∈ M h and z n+1 h ∈ Z h such that ∀v h ∈ X h , 1 k (u n+1 h -u n h , v h ) + α k (∇(u n+1 h -u n h ), ∇v h ) + ν(∇u n+1 h , ∇v h ) +(z n h × u n+1 h , v h ) -(p n+1 h , div v h ) = (f n+1 , v h ), (9)
∀θ h ∈ Z h , α k (z n+1 h -z n h , θ h ) + ν(z n+1 h , θ h ) + c(u n+1 h ; z n+1 h , θ h ) = ν(curlu n+1 h , θ h ) + α(curlf n+1 , θ h ), (10)
where c(u n+1 h ; z n+1 h , θ h ) is the discrete non-linear part of the transport equation and the functions of X h vanish on ∂Ω. This system is linearized in the sense that in [START_REF] Girault | Raviart Finite element methods for the Navier-Stokes Equations. Theory and Algorithms[END_REF], knowing z n h , we calculate u n+1 h and p n+1 h with a linear equation. Then, we calculate z n+1 h with the second linear equation [START_REF] Girault | On a time-dependent grade-two fluid model in two dimensions, Computers and Mathematics with Applications[END_REF]. For both the velocity and pressure discretizations, the error is of order h 3/2 and k. This is the best that can be achieved, considering that the discretization of the transport equation loses inevitably a factor h 1/2 . Other finite elements can be used, cf. Crouzeix and Raviart [START_REF] Crouzeix | Raviart Conforming and non-conforming finite element methods for solving the stationnary Stokes problem[END_REF], Brezzi and Fortin [START_REF] Brezzi | Mixed and Hybrid finite element methods[END_REF] and Girault and Raviart [START_REF] Girault | Raviart Finite element methods for the Navier-Stokes Equations. Theory and Algorithms[END_REF]. Now, we recall some notation and basic functional results. As usual, for handling time-dependent problems, it is convenient to consider functions defined on a time interval ]a, b[ with values in a functional space, say X (cf. Lions and Magenes [START_REF] Lions | Magenes Problèmes aux limites non homogènes et applications I[END_REF]). More precisely, let . X denote the norm of X; then for any r, 1 ≤ r ≤ ∞, we define

L r (a, b; X) = {f mesurable in ]a, b[; b a f (t) r X dt < ∞} equipped with the norm f L r (a,b;X) = ( b a f (t) r X dt) 1/r , with the usual modifications if r = ∞. It is a Banach space if X is a Banach space. Let (k 1 , k 2 ) denote a pair of non-negative integers, set |k| = k 1 + k 2 and define the partial derivative ∂ k by ∂ k v = ∂ |k| v ∂x k 1 1 ∂x k 2 2 .
We denote by :

W m,r (Ω) = {v ∈ L r (Ω); ∂ k v ∈ L r (Ω) ∀|k| ≤ m},
This space is equipped with the seminorm

|v| W m,r (Ω) = [ |k|=m Ω |∂ k v| r dx] 1/r ,
and is a Banach space for the norm

v W m,r (Ω) = [ 0≤|k|≤m |v| r W k,r (Ω) dx] 1/r .
When r = 2, this space is the Hilbert space H m (Ω). In particular, the scalar product of L 2 (Ω) is denoted by (. , .).

Similarly, L 2 (a, b; H m (Ω)) is a Hilbert space and in particular L 2 (a, b; L 2 (Ω)) coincides with L 2 (Ω×]a, b[).
The definitions of these spaces are extended straightforwardly to vectors, with the same notation, but with the following modification for the norms in the non-Hilbert case. Let u = (u 1 , u 2 ); then we set

u L r (Ω) = [ Ω u(x) r dx] 1/r ,
where . denotes the Euclidean vector norm.

For functions that vanish on the boundary, we define for any r ≥ 1

W 1,r 0 (Ω) = {v ∈ W 1,r (
Ω); v| ∂Ω = 0} and recall Sobolev's imbeddings in two dimensions: for each r ∈ [2, ∞[, there exits a constant S r such that ∀v ∈

H 1 0 (Ω) , v L r (Ω) ≤ S r |v| H 1 (Ω) , (11)
where

|v| H 1 (Ω) = ∇v L 2 (Ω) . (12)
When r = 2, (11) reduces to Poincaré's inequality and S 2 is Poincaré's constant.

The case r = ∞ is excluded and is replaced by: for any r > 2, there exists a constant M r such that

∀v ∈ W 1,r 0 (Ω), v L ∞ (Ω) ≤ M r |v| W 1,r (Ω) . (13)
We have also in dimension 2,

||g|| L 4 (Ω) ≤ 2 1/4 g 1/2 L 2 (Ω) ∇g 1/2 L 2 (Ω) . (14)
Owing to Poincaré's inequality, the seminorm |.| H 1 (Ω) is a norm on H 1 0 (Ω) and we use it to define the dual norm:

f H -1 (Ω) = sup v∈H 1 0 (Ω) f, v |v| H 1 (Ω) ,
where ., . denotes the duality pairing between H -1 (Ω) and H 1 0 (Ω). Also, we introduce the space:

L 2 0 (Ω) = {q ∈ L 2 (Ω); Ω q dx = 0},

The exact problem

Let Ω be a bounded polygon in two dimensions with boundary ∂Ω and let ]0, T [ be a given time-interval. We want to find a vector velocity u, a scalar pressure p and an auxiliary scalar function z solution of

∂ ∂t (u -α∆u) -ν∆u + z × u + ∇p = f in ]0, T [×Ω, (15)
α ∂z ∂t + νz + αu • ∇z = ν curlu + α curlf in ]0, T [×Ω, (16)
u = 0 on ]0, T [×∂Ω, (17)
u(x, t) = 0 and z(x, t) = 0,

where z × u = (-zu 2 , zu 1 ). Here ν > 0 and α > 0 are given constants.

A straightforward formulation of ( 15)-( 18) is :

Find (u(t), p(t), z(t)) ∈ L ∞ (0, T ; H 1 0 (Ω) 2 )×L 2 (0, T ; L 2 0 (Ω))×L ∞ (0, T ; L 2 (Ω)), u ′ ∈ L 2 (0, T ; H 1 0 (Ω) 2 ) such that ∀v ∈ H 1 0 (Ω), (u ′ (t), v) + α(∇u ′ (t), ∇v) + ν(∇u(t), ∇v) +(z(t) × u(t), v) -(p(t), div v) = (f (t), v) in Ω×]0, T ], (19)
∀q ∈ L 2 0 (Ω), (q(t), div u(t)) = 0, (20)
α ∂z ∂t + νz + αu • ∇z = ν curlu + α curlf in ]0, T [×Ω, (21)
u(0) = 0 and z(0) = 0 in Ω. ( 22
)
The following theorem is established in [START_REF] Girault | On a time-dependent grade-two fluid model in two dimensions, Computers and Mathematics with Applications[END_REF]:

Theorem 2.1.
Let Ω be a lipschitz polygon. For all ν > 0 and f ∈ L 2 (0, T ; L 2 (Ω) 2 ) such that curl f ∈ L 2 (0, T ; L 2 (Ω)), ( 15)-( 18) has at least one solution (u, z, p) that satisfies the following estimates :

||z|| L ∞ (0,T ;L 2 (Ω)) ≤ √ 2 S 2 ν ||f || L 2 (0,T ;L 2 (Ω) 2 ) + |α| ν || curl f || L 2 (0,T ;L 2 (Ω)) , ||u|| L ∞ (0,T ;H 1 (Ω) 2 ) ≤ S 2 ν ||f || L 2 (0,T ;L 2 (Ω) 2 ) , ||p|| L 2 (0,T ;L 2 (Ω)) ≤ 1 β (S 2 ||f || L 2 (0,T ;L 2 (Ω) 2 ) + S 2 4 ||u|| L ∞ (0,T ;H 1 (Ω) 2 ) ||z|| L ∞ (0,T ;L 2 (Ω)) ).

A discontinuous upwind scheme

Let h > 0 be a discretization parameter and let T h be a regular family of triangulation of Ω, consisting of triangles κ with maximum mesh size h: There exists a constant σ 0 , independent of h, such that ∀κ ∈ T h , h κ ρ κ ≤ σ 0 , where h κ is the diameter and ρ κ is the diameter of the ball inscribed in κ. We introduce ρ min = min κ ρ k . As usual, the triangulation is such that any two triangles are either disjoint or share a vertex or a complete side. We first recall how upwinding can be achieved by the discontinuous Galerkin approximation introduced in [START_REF] Lesaint | On a finite element method for solving the neutron transport equation, Mathematical aspects of finite Elements in Partial Differential Equation[END_REF]. Let Z h be the discontinuous finite-element space :

Z h = {θ h ∈ L 2 (Ω); ∀κ ∈ T h , θ h | κ ∈ IP 1 }.
There exists an approximation operator, [START_REF] Clément | Approximation by finite element functions using local regularization[END_REF], R h ∈ L(W 1,p (Ω); Z h ∩ C 0 (Ω)) such that for any p ≥ 1, for m = 0, 1 and 0

≤ l ≤ 1 ∀z ∈ W l+1,p (Ω), |R h (z) -z| W m,p (Ω) ≤ Ch l+1-m |z| W l+1,p (Ω) .
Let u h be a discrete velocity in H 1 0 (Ω) 2 , and for each triangle κ, let

∂κ -= {x ∈ ∂κ; αu h • n < 0},
where n denotes the unit exterior normal to ∂k. Note that, for all triangles κ of T h , ∂κ -only involves interior segments of T h because u h = 0 on ∂Ω. Then, the non-linear term α[(u

• ∇z, θ) + 1 2 (div u z, θ)] is approximated by c(u n+1 h ; z n+1 h , θ n+1 h ) = α 2 Ω div u n+1 h z n+1 h θ n+1 h dx + κ∈T h κ α(u n+1 h • ∇z n+1 h )θ n+1 h dx + ∂κ - |αu n+1 h • n|(z n+1 h,int -z n+1 h,ext )θ n+1 h,int ds .
The subscript int (resp. ext) refers to the trace on the segment ∂κ of the function taken inside (resp. outside) κ. Note that in the above sum, the boundary integrations act in fact over complete interior segments.

On the other hand, let us recall the standard Taylor-Hood discretization of the velocity and pressure. The discrete space of the pressure is :

M h = {q h ∈ H 1 (Ω) ∩ L 2 0 (Ω); ∀κ ∈ T h , q h ∈ IP 1 }. There exists an operator r h ∈ L(L 2 0 (Ω); M h ) such that for 0 ≤ l ≤ 2, ∀q ∈ H l (Ω) ∩ L 2 0 (Ω), ||r h (q) -q|| L 2 (Ω) ≤ Ch l ||q|| H l (Ω) .
The discrete velocity space is :

X h = {v h ∈ C 0 (Ω); ∀κ ∈ T h , v h | κ ∈ IP 2 , v h | ∂Ω = 0},
and let V h = {v h ∈ X h ; (q h , div v h ) = 0 ∀q ∈ M h }.
There exists an operator P h ∈ L(H 1 0 (Ω) 2 ; X h ), such that

                 ∀v ∈ H 1 0 (Ω) 2 , ∀κ ∈ T h , ∀q h ∈ M h , κ q h div(P h (v) -v)dx = 0, for all number p ≥ 2; ∀v ∈ H 1 0 (Ω) 2 , ||P h (v) -v|| L p (Ω) ≤ Ch 2/p |v| H 1 (Ω) , for all number p ≥ 2, 1 ≤ s ≤ 3, m = 0 or 1, ∀v ∈ [W s,p (Ω) ∩ H 1 0 (Ω)] 2 , ||P h (v) -v|| W m,p (Ω) ≤ Ch s-m |v| W s,p (Ω) . (23) We take f n+1 (x) = 1 k t n+1
t n f (t, x)dt. Then the discrete system corresponding to the formulation ( 19)-( 22) is :

Given (u 0 h , z 0 h ) = (0, 0) and z n h ∈ Z h , find (u n+1 h , p n+1 h ) ∈ X h × M h such that : ∀v h ∈ X h , 1 k (u n+1 h -u n h , v h ) + α k (∇(u n+1 h -u n h ), ∇v h ) + ν(∇u n+1 h , ∇v h ) + (z n h × u n+1 h , v h ) -(p n+1 h , div v h ) = (f n+1 , v h ), ( 24
)
∀q h ∈ M h , (q h , div u n+1 h ) = 0. ( 25
)
Once we have u n+1 h , we compute z n+1 h by solving the system :

∀θ h ∈ Z h , 1 k (z n+1 h -z n h , θ h ) + ν(z n+1 h , θ h ) + c(u n+1 h ; z n+1 h , θ h ) = ν(curlu n+1 h , θ h ) + α(curlf n+1 , θ h ). (26)
In order to prove the existence of solutions of ( 24)-( 26), let us recall the following identity established by Lesaint and Raviart [START_REF] Lesaint | On a finite element method for solving the neutron transport equation, Mathematical aspects of finite Elements in Partial Differential Equation[END_REF] :

Lemma 3.1. For all v n h ∈ X h , z n h and θ n h in Z h , we have c(v n h ; z n h , θ n h ) = κ∈T h - κ α(v n h • ∇θ n h )z n h dx + ∂κ - α|v n h • n|z n h,ext (θ n h,ext -θ n h,int )ds - α 2 Ω (div v n h )θ n h z n h dx. For θ n h ∈ H 1 (Ω), we have c(v n h ; z n h , θ n h ) = - Ω α(v n h • ∇θ n h )z n h dx - α 2 Ω (div v n h )θ n h z n h dx. For θ n h = z n h ∈ Z h we have c(v n h ; z n h , z n h ) = 1 2 κ∈T h ∂κ - |αv n h • n|(z n h,ext -z n h,int ) 2 ds. Theorem 3.2. Given f n+1 ∈ L 2 (Ω) 2 with curl f n+1 ∈ L 2 (Ω), for all (u n h , z n h ) ∈ X h × Z h , there exists a unique solution (u n+1 h , p n+1 h , z n+1 h ) of problem (24)-(26) that belongs to X h × M h × Z h .
Proof. On the one hand, for z n h ∈ Z h , it is clear that problem (24)-( 25) has a unique solution (u n+1 h , p n+1 h ) as a consequence of the coerciveness of the corresponding bilinear form on X h ×X h . On the other hand, the last lemma proves that the bilinear form corresponding to the equation ( 26) is also coercive on Z h × Z h . Then (26) has a unique solution.

Theorem 3.3. We assume that f ∈ L 2 (0, T ; L 2 (Ω) 2 ) with curl f ∈ L 2 (0, T ; L 2 (Ω)).
The solution of the problem (24)-( 26) satisfies :

||u h || L ∞ (0,T ;H 1 (Ω) 2 ) ≤ C 1 ||f || L 2 (0,T ;L 2 (Ω) 2 ) , ||z h || 2 L ∞ (0,T ;L 2 (Ω)) ≤ C 2 ||f || 2 L 2 (0,T ;L 2 (Ω) 2 ) + C 3 || curl f || 2 L 2 (0,T ;L 2 (Ω)) , ||p h || 2 L 2 (0,T ;L 2 (Ω)) ≤ C 4 ||f || 2 L 2 (0,T ;L 2 (Ω) 2 ) + C 5 ||u h || 2 L ∞ (0,T ;H 1 (Ω) 2 ) + C 6 ||z h || 2 L ∞ (0,T ;L 2 (Ω)) ||u h || 2 L ∞ (0,T ;H 1 (Ω) 2 ) ,
where C i , i = 1, . . . , 6 are positive constants that depend on Ω and T .

Proof. On the one hand, we take v h = u n+1 h in (24) and we obtain :

1 2 ||u n+1 h || 2 L 2 (Ω) - 1 2 ||u n h || 2 L 2 (Ω) + α 2 |u n+1 h | 2 H 1 (Ω) - α 2 |u n h | 2 H 1 (Ω) + νk|u n+1 h | 2 H 1 (Ω) ≤ kε 2 ||f n+1 || 2 L 2 (Ω) + kS 2 2 2ε |u n+1 h | 2 H 1 (Ω) . We choose ε = S 2 2 2ν
and sum over n = 0, . . . i. We obtain :

1 2 ||u i h || 2 L 2 (Ω) + α 2 |u i h | 2 H 1 (Ω) ≤ i n=1 kS 2 2 4ν ||f n || 2 L 2 (Ω) ≤ N n=1 kS 2 2 4ν ||f n || 2 L 2 (Ω) .
This implies the first estimate :

||u h || 2 L ∞ (0,T ;H 1 (Ω) 2 ) = sup 0≤i≤N |u i h | 2 H 1 (Ω) ≤ S 2 2 2να ||f || 2 L 2 (0,T ;L 2 (Ω) 2 ) .
On the other hand, we choose θ h = z n+1 h in (26), use the third relation in Lemma 3.1 and we obtain :

α 2 ||z n+1 h || 2 L 2 (Ω) - α 2 ||z n h || 2 L 2 (Ω) + νk||z n+1 h || 2 L 2 (Ω) ≤ kν 2ε 1 || curl u n+1 h || 2 L 2 (Ω) + kνε 1 2 ||z n+1 h || 2 L 2 (Ω) + kα 2ε 2 || curl f n+1 || 2 L 2 (Ω) + kαε 2 2 ||z n+1 h || 2 L 2 (Ω) , taking ε 1 = 1, ε 2 = ν α
and summing over n = 0, . . . , i, this becomes :

α||z i h || 2 L 2 (Ω) ≤ νT ||u h || 2 L ∞ (0,T ;H 1 (Ω) 2 ) + α 2 ν || curl f || 2 L 2 (0,T ;L 2 (Ω)) .
Then we obtain the second estimate :

||z h || 2 L ∞ (0,T ;L 2 (Ω)) ≤ νT α ||u h || 2 L ∞ (0,T ;H 1 (Ω) 2 ) + α ν || curl f || 2 L 2 (0,T ;L 2 (Ω)) .
The third estimate is obtained in two steps: First, we take the function test

v n+1 h = u n+1 h -u n h k in (24). We obtain : 1 k 2 ||u n+1 h -u n h || 2 L 2 (Ω) + α k 2 |u n+1 h -u n h | 2 H 1 (Ω) + ν 2 |u n+1 h | 2 H 1 (Ω) + ν 2k 2 |u n+1 h -u n h | 2 H 1 (Ω) ≤ 1 2ε 1 ||f n+1 || 2 L 2 (Ω) + ε 1 2k 2 ||u n+1 h -u n h || 2 L 2 (Ω) + 1 2ε 2 ||z n h || 2 L 2 (Ω) ||u n+1 h || 2 L 4 (Ω) + S 2 4 ε 2 2k 2 |u n+1 h -u n h | 2 H 1 (Ω) .
Then by choosing ε 1 = 2 and

ε 2 = α 2S 2 4
we obtain :

| u n+1 h -u n h k | 2 H 1 (Ω) ≤ 1 2α ||f n+1 || 2 L 2 (Ω) + 2S 2 4 α 2 ||z n h || 2 L 2 (Ω) ||u n+1 h || 2 L 4 (Ω) + 2ν 2 α 2 |u n+1 h | 2 H 1 (Ω) .
Next, owing that the pair (X h , M h ) satisfies a uniform discrete inf-sup condition, we associate with

p n+1 h ∈ M h the function v h ∈ X h defined by        ∀w h ∈ V h , (∇v h , ∇w h ) = 0, ∀q h ∈ M h , (div v h , q h ) = (p n+1 h , q h ), |v h | H 1 (Ω) ≤ 1 β ||p n+1 h || L 2 (Ω) , (27)
we substitute this v h into (24) and we obtain :

||p n+1 h || 2 L 2 (Ω) ≤ P + α 2ε 1 | u n+1 h -u n h k | 2 H 1 (Ω) + P + α 2β 2 ε 1 ||p n+1 h || 2 L 2 (Ω)) + ν 2ε 2 |u n+1 h | 2 H 1 (Ω) + νε 2 2β 2 ||p n+1 h || 2 L 2 (Ω)) + 1 2ε 3 ||z n h || 2 L 2 (Ω) ||u n+1 h || 2 L 4 (Ω) + S 2 4 ε 3 2β 2 ||p n+1 h || 2 L 2 (Ω)) + 1 2ε 4 ||f n+1 || 2 L 2 (Ω) + P 2 ε 4 2β 2 ||p n+1 h || 2 L 2 (Ω)
.

By choosing ε 1 = β 2 4(P + α) , ε 2 = β 2 4ν , ε 3 = β 2 4S 2 4
and ε 4 = β 2 4P 2 and summing over n from 0 to N -1, we obtain the third estimate.

Error estimates

Theorem 4.1. Under the assumptions u ∈ L ∞ (0, T ; W 1,4 (Ω) 2 )∩L 2 (0, T ; H 3 (Ω) 2 ), u ′ ∈ L 2 (0, T ; H 3 (Ω) 2 ), p ∈ L 2 (0, T ; H 2 (Ω)), z ∈ L ∞ (0, T ; L 2 (Ω)) and z ′ ∈ L 2 (0, T ; L 2 (Ω))
, there exist positive constants C and C ′ that depend on u, z, Ω and T such that :

1 2 ||u N h -u(t N )|| 2 L 2 (Ω) + α 2 ||∇u N h -∇u(t N )|| 2 L 2 (Ω) + ν 2 N -1 n=0 k|u n+1 h -u(t n+1 )| 2 H 1 (Ω) ≤ C(h 4 + k 2 ) + C ′ N -1 n=0 k||z n+1 h -z(t n+1 )|| 2 L 2 (Ω) . (28)
Proof. We consider (19), choose the function test

v n+1 h = u n+1 h -P h u(t n+1
), integrate from t n to t n+1 and take the difference between this and (24) multiplied by k. We obtain :

(u n+1 h -u(t n+1 )) -(u n h -u(t n )), v n+1 h + α ∇(u n+1 h -u(t n+1 )) -∇(u n h -u(t n )), ∇v n+1 h +ν k∇u n+1 h - t n+1 t n ∇u(t)dt, ∇v n+1 h -kp n+1 h - t n+1 t n p(t)dt, div v n+1 h + kz n h × u n+1 h - t n+1 t n z(t) × u(t)dt, v n+1 h = 0.
Let us treat the terms in the left-hand side of this equation that we denote (a i ), i = 1, ..., 5.

For the first term, we insert P h u(t n+1 ) and P h u(t n ) and we split (a 1 ) into two terms that we treat separately. The first part is as follows :

(a 1,1 ) = 1 2 ||v n+1 h || 2 L 2 (Ω) - 1 2 ||v n h || 2 L 2 (Ω) + 1 2 ||v n+1 h -v n h || 2 L 2 (Ω) ,
and the second part is as follows :

|(a 1,2 )| = t n+1 t n (P h u ′ (τ ) -u ′ (τ ))dτ, v n+1 h ≤ 1 2ε 1 Ch 4 ||u ′ || 2 L 2 (t n ,t n+1 ;H 2 (Ω) 2 ) + S 2 ε 1 2 k|v n+1 h | 2 H 1 (Ω) .
We treat the second term (a 2 ) as the first one and we obtain :

(a 2,1 ) = α 2 |v n+1 h | 2 H 1 (Ω) - α 2 |v n h | 2 H 1 (Ω) + α 2 |v n+1 h -v n h | 2 H 1 (Ω) , and
|(a 2,2 )| ≤ Cα 2ε 2 h 4 ||u ′ || 2 L 2 (t n ,t n+1 ;H 3 (Ω) 2 ) + ε 2 α 2 k|v n+1 h | 2 H 1 (Ω) .
For the third term (a 3 ), we insert ∇P h u(t n+1 ) and ∇P h u(t) and we split it into three parts that are treated successively as follows :

(a 3,1 ) = νk|u n+1 h -P h u(t n+1 )| 2 H 1 (Ω) , |(a 3,2 )| = ν t n+1 t n ∇P h (u(t n+1 ) -u(t))dt, ∇v n+1 h = ν t n+1 t n ∇P h u ′ (τ )(τ -t n )dτ, ∇v n+1 h ≤ νε 3 2 √ 3 k|v n+1 h | 2 H 1 (Ω) + νC 2 k 2 2ε 3 √ 3 ||u ′ || 2 L 2 (t n ,t n+1 ;H 2 (Ω) 2 ) , and
|(a 3,3 )| = ν ∇ t n+1 t n (P h u(t) -u(t))dt, ∇v n+1 h ≤ νC 2 2ε 4 h 4 ||u|| 2 L 2 (t n ,t n+1 ,H 3 (Ω) 2 ) + νε 4 2 k|v n+1 h | 2 H 1 (Ω) .
To study the fourth term, we use the fact that

Ω p n+1 h div(P h u(t n+1 ) -u(t n+1 )) = 0, div u(t n+1 ) = 0 and Ω p n+1 h div u n+1 h
= 0 and we obtain :

|(a 4 )| = t n+1 t n (r h p(t) -p(t))dt, div v n+1 h ≤ C 1 2ε 5 h 4 ||p|| 2 L 2 (t n ,t n+1 ;H 2 (Ω)) + ε 5 2 k|v n+1 h | 2 H 1 (Ω) .
Finally, for the last term (a 5 ), we have (

z n h × u n+1 h , v n+1 h ) = (z n h × P h u(t n+1 ), v n+1 h ), because (a × b, b) = 0. But z n h × P h u(t n+1 ) -z(t) × u(t) = (z n h -z(t n )) × P h u(t n+1 ) + z(t n ) × P h (u(t n+1 ) -u(t)) +z(t n ) × (P h u(t) -u(t)) + (z(t n ) -z(t)) × (u(t) -u(t n )) + (z(t n ) -z(t)) × u(t n ),
than (a 5 ) is split into five parts that we treat successively.

The first part is as follows :

|(a 5,1 )| = t n+1 t n ((z n h -z(t n )) × P h u(t n+1 ), v n+1 h )dt ≤ S 2 4 ε 6 2 ||P h u|| 2 L ∞ (0,T ;H 1 (Ω) 2 ) k||z n h -z(t n )|| 2 L 2 + S 2 4 2ε 6 k|v n+1 h | 2 H 1 (Ω) . ≤ S 2 4 ε 6 2 ||u|| 2 L ∞ (0,T ;W 1,4 (Ω) 2 ) k||z n h -z(t n )|| 2 L 2 + S 2 4 2ε 6 k|v n+1 h | 2 H 1 (Ω) .
The second part is as follows :

|(a 5,2 )| = (z(t n ) × t n+1 t n (P h u(t n+1 ) -P h u(t))dt, v n+1 h ) = (z(t n ) × t n+1 t n P h u ′ (τ )(τ -t n )dτ, v n+1 h ) ≤ S 2 4 2 √ 3 C ′′′ ε 7 ||z|| 2 L ∞ (0,T ;L 2 (Ω)) ||u ′ || 2 L 2 (t n ,t n+1 ;H 2 (Ω) 2 ) k 2 + kε 7 |v n+1 h | 2 H 1 (Ω) .
For the third part, we have

|(a 5,3 )| = t n+1 t n (z(t n ) × (P h u(t) -u(t)), v n+1 h )dt ≤ S 2 4 ||z(t n )|| L 2 (Ω) |v n+1 h | H 1 (Ω) t n+1 t n |P h u(t) -u(t)| H 1 (Ω) dt ≤ C 2 S 2 4 2ε 8 ||z|| 2 L ∞ (0,T ;L 2 (Ω)) ||u|| 2 L 2 (t n ,t n+1 ;H 3 (Ω) 2 ) h 4 + C 2 S 2 4 ε 8 2 k|v n+1 h | 2 H 1 (Ω) .
The fourth part is treated as follows :

|(a 5,4 )| = t n+1 t n (z(t n ) -z(t)) × (u(t) -u(t n )), v n+1 h dt = t n+1 t n t t n z ′ (τ )dτ × t t n u ′ (τ )dτ , v n+1 h dt ≤ S 2 4 ε 9 2 √ 2 k |v n+1 h | 2 H 1 (Ω) + S 2 4 k 3 2 √ 2ε 9 ||u ′ || 2 L 2 (0,T ;H 1 (Ω) 2 ) ||z ′ || 2 L 2 (0,T ;L 2 (Ω)) .
Finally, for the last part, we have

|(a 5,5 )| = t n+1 t n (z(t n ) -z(t)) × u(t n ), v n+1 h dt = t n+1 t n z ′ (t)(t -t n+1 )dt × u(t n ), v n+1 h ≤ S 2 4 ε 10 2 √ 3 k|v n+1 h | 2 H 1 (Ω) + S 2 4 2 √ 3ε 10 k 2 ||z ′ || 2 L 2 (t n ,t n+1 ;L 2 (Ω)) ||u|| 2 L ∞ (0,T ;H 1 (Ω) 2 )
At the end, (28) follows easily after the decomposition

(a 1,1 ) + (a 2,1 ) + (a 3,1 ) ≤ |(a 1,2 )| + |(a 2,2 )| + |(a 3,2 )| + |(a 3,3 )| + |(a 4 )| + |(a 5 )|,
the sum over n = 1, . . . , N -1, a suitable choice of ε i , i = 1, . . . , 10 and by using the properties of P h in :

|u n+1 h -u(t n+1 )| H 1 (Ω) ≤ |u n+1 h -P h u(t n+1 )| H 1 (Ω) + |P h u(t n+1 ) -u(t n+1 )| H 1 (Ω) .
We define ρ h as the L 2 projection of z onto IP 1 in each triangle κ : for z ∈ L 2 (Ω),

∀q ∈ IP 1 , κ (ρ h (z) -z)qdx = 0.
This operator has locally the same accuracy as R h .

Theorem 4.2. We suppose that there exists a constant γ > 0 such that k ≤ γh. In addition to the assumptions of Theorem 4.1, we assume that u ∈ L ∞ (0,

T ; W 1,∞ (Ω) 2 ), z ∈ L ∞ (0, T ; W 1,4 (Ω)) and z ′ ∈ L ∞ (0, T ; L 4 (Ω)) ∩ L 2 (0, T ; W 1,4 (Ω)). Denoting θ n+1 h = z n+1 h - ρ h z(t n+1 ) we have : N -1 n=0 t n+1 t n c(u n+1 h ; ρ h z(t n+1 ) -z(t), θ n+1 h )dt ≤ L 1 (h 3 + k 2 ) + L 3 N -1 n=0 k||θ n+1 h || 2 L 2 (Ω) +L 2 N -1 n=0 k|u n+1 h -u(t n+1 )| 2 H 1 (Ω) + α 2 κ∈T h N -1 n=0 t n+1 t n ∂κ - |u n+1 h • n|(θ n+1 h,ext -θ n+1 h,int ) 2 dsdt
(29) where L i are constants that only depend on u, z, Ω, T and arbitrary coefficients ε i (i = 1, . . . , 8)

Proof. Owing to Lemma 3.1 and denoting ξ h = ρ h (z(t n+1 )), we have :

t n+1 t n c(u n+1 h ; ρ h z(t n+1 ) -z(t), θ n+1 h )dt = α κ∈T h - t n+1 t n κ (u n+1 h • ∇θ n+1 h )(ξ h -z(t))dxdt + t n+1 t n ∂κ - |u n+1 h • n|(ξ h -z(t)) ext (θ n+1 h,ext -θ n+1 h,int )dsdt - α 2 t n+1 t n Ω div(u n+1 h -u(t))(ξ h -z(t))θ n+1 h dxdt
In the left-hand side, we denote the terms by (d i ), i = 1 . . . , 3. as ∇(z hρ h (z)) is a constant vector, the first term (d 1 ), for any constant vector c, can be treated as :

t n+1 t n κ (u n+1 h • ∇θ n+1 h )(ξ h -z(t))dxdt ≤ α t n+1 t n κ (u n+1 h -u(t n+1 ))∇θ n+1 h (ξ h -z(t))dxdt +α t n+1 t n κ (u(t n+1 ) -c)∇θ n+1 h (ξ h -z(t n+1 ))dxdt + t n+1 t n κ u(t n+1 )∇θ n+1 h (z(t n+1 ) -z(t))dxdt . With α t n+1 t n κ (u n+1 h -u(t n+1 ))∇θ n+1 h (ξ h -z(t))dxdt ≤ αS 4 ρ k |u n+1 h -u(t n+1 )| H 1 (κ) ||θ n+1 h || L 2 (κ) c 1 k||ξ h -z(t n+1 )|| L 4 (κ) + c 2 k 2 √ 2 ||z ′ || L ∞ (0,T ;L 4 (κ)) ≤ αc 3 S 4 σ 0 ||z|| L ∞ (0,T ;W 1,4 (κ)) + γ||z ′ || L ∞ (0,T ;L 4 (κ)) k 2ε 1 ||θ n+1 h || 2 L 2 (κ) + ε 1 k 2 |u n+1 h -u(t n+1 )| 2 H 1 (κ) , and
α t n+1 t n κ (u(t n+1 ) -c)∇θ n+1 h (ξ h -z(t n+1 ))dxdt + t n+1 t n κ u(t n+1 )∇θ n+1 h (z(t n+1 ) -z(t))dxdt = α t n+1 t n κ (u(t n+1 ) -c)∇θ n+1 h (ξ h -z(t n+1 ))dxdt - t n+1 t n κ θ n+1 h u(t n+1 )∇(z(t n+1 ) -z(t))dxdt ≤ ||θ n+1 h || L 2 (κ) αc 3 ρ k kh 3/2 ||z(t n+1 )|| W 1,4 (κ) ||u(t n+1 ) -c|| L ∞ (κ) + ||u|| L ∞ (0,T ;κ) αk 3/2 √ 3 ||z ′ || L 2 (t n ,t n+1 ;H 1 (κ)) ≤ αc 4 σ 0 kh 3/2 ||θ n+1 h || L 2 (κ) ||u(t n+1 )|| W 1,∞ (κ) ||z|| L ∞ (0,T ;W 1,4 (κ)) + k 3/2 √ 3 ||θ n+1 h || L 2 (κ) ||z ′ || L 2 (t n ,t n+1 ;H 1 (κ)) ||u|| L ∞ (0,T ;κ) ≤ αc 5 ||u|| L ∞ (0,T ;W 1,∞ (κ) 2 ) k 2ε 2 ||θ n+1 h || 2 L 2 (κ) + ε 2 2 (k 2 ||z ′ || 2 L 2 (t n ,t n+1 ;H 1 (κ)) + kh 3 ||z|| 2 L ∞ (t n ,t n+1 ;W 1,4 (κ)) ) .
For the second part (d 2 ), we write :

α t n+1 t n ∂κ - |u n+1 h • n|(ξ h -z(t)) ext (θ n+1 h,ext -θ n+1 h,int )dsdt ≤ α 2 t n+1 t n ∂κ - |u n+1 h • n|(θ n+1 h,ext -θ n+1 h,int ) 2 dsdt + α 2 t n+1 t n ∂κ - |u n+1 h • n| (ξ h -z(t)) ext 2 dsdt
We keep the first term in the right-hand side of this inequality. The second term can be written as follows :

α 2 t n+1 t n ∂κ - |u n+1 h • n| (ξ h -z(t)) ext 2 dsdt ≤ α t n+1 t n ∂κ - |u n+1 h • n| (ξ h -z(t n+1 )) ext 2 dsdt + α t n+1 t n ∂κ - |u n+1 h • n| (z(t n+1 ) -z(t)) ext 2 dsdt with α t n+1 t n ∂κ - |u n+1 h • n| (ξ h -z(t n+1 )) ext 2 dsdt ≤ αc 6 k ||ξ h -z(t n+1 )|| 2 L 4 (∂κ -) ||u n+1 h -u(t n+1 )|| L 2 (∂κ -) + ||ξ h -z(t n+1 )|| 2 L 2 (∂κ -) ||u(t n+1 )|| L ∞ (ω k ) ≤ αc 7 k h 3/2 ||z(t n+1 )|| 2 W 1,4 (ω k ) ||u n+1 h -u(t n+1 )|| H 1 (ω k ) + h 3 ||z(t n+1 )|| 2 W 1,4 (ω k ) ||u(t n+1 )|| L ∞ (ω k ) ≤ αc 7 h 3 k ||z|| 2 L ∞ (0,T ;W 1,4 (ω k )) ||u|| L ∞ (0,T ;ω k ) + 1 2ε 3 ||z|| 2 L ∞ (0,T ;W 1,4 (ω k )) + ε 3 2 k||u n+1 h -u(t n+1 )|| 2 H 1 (ω k ) ,
where ω k denotes the union of triangles adjacent to κ and

α t n+1 t n ∂κ - |u n+1 h • n| (z(t n+1 ) -z(t)) ext 2 dsdt ≤ α ∂κ - |u n+1 h • n| t n+1 t n (t n+1 -t) t n+1 t |z ′ (τ )| 2 dτ ext dtds ≤ α ∂κ - |u n+1 h • n| t n+1 t n |z ′ (τ )| 2 dτ ext t n+1 t n (t n+1 -t)dtds ≤ αk 2 2 ∂κ - |(u n+1 h -u(t n+1 )) • n| t n+1 t n |z ′ (τ )| 2 dτ ext ds + ∂κ - |u(t n+1 ) • n| t n+1 t n |z ′ (τ )| 2 dτ ext ds ≤ αc 8 k 2ε 4 |u n+1 h -u(t n+1 )| 2 H 1 (ω k )) + ε 4 k 3 2 ||z ′ || 2 L 2 (t n ,t n+1 ;W 1,4 (ω k )) +k 2 ||u|| L ∞ (0,T ;H 1 (ω k ) 2 ) ||z ′ || 2 L 2 (t n ,t n+1 ;W 1,4 (ω k ))
For the third part (d 3 ), we insert respectively in the divergence term and the term in z, ±u(t n+1 ) and ±z(t n+1 ). We get four parts that we treat as follow :

|(d 3,1 )| = α 2 t n+1 t n Ω div(u n+1 h -u(t n+1 ))(ξ h -z(t n+1 ))θ n+1 h dxdt ≤ α 2 σ 0 C 9 ||z|| L ∞ (0,T ;W 1,4 (Ω)) k 2ε 5 |u n+1 h -u(t n+1 )| 2 H 1 (Ω) + ε 5 2 k||θ n+1 h || 2 L 2 (Ω) , |(d 3,2 )| = α 2 t n+1 t n Ω div(u(t n+1 ) -u(t))(ξ h -z(t n+1 ))θ n+1 h dxdt ≤ α 2 C 10 σ 0 k 3/2 ||z|| L ∞ (0,T ;W 1,4 (Ω)) ||θ n+1 h || L 2 (Ω) ||u ′ || L 2 (t n ,t n+1 ;H 1 (Ω) 2 ) ≤ α 2 C 10 σ 0 ||z|| L ∞ (0,T ;W 1,4 (Ω)) k 2 2ε 6 ||u ′ || 2 L 2 (t n ,t n+1 ;H 1 (Ω) 2 ) + ε 6 2 k||θ n+1 h || 2 L 2 (Ω) , |(d 3,3 )| = α 2 t n+1 t n Ω div(u(t n+1 ) -u(t))(z(t n+1 ) -z(t))θ n+1 h dxdt ≤ α 2 C 11 k 2 ρ min ||z ′ || L 2 (t n ,t n+1 ;L 4 (Ω)) ||θ n+1 h || L 2 (Ω) ||u ′ || L 2 (t n ,t n+1 ,H 1 (Ω) 2 ) ≤ α 2 C 11 σ 0 γ||u ′ || L ∞ (0,T ;H 1 (Ω) 2 ) k 2 2ε 7 ||z ′ || 2 L 2 (t n ,t n+1 ;L 4 (Ω)) + ε 7 2 k||θ n+1 h || 2 L 2 (Ω) , and
|(d 3,4 )| = α 2 t n+1 t n Ω div(u n+1 h -u(t n+1 ))(z(t n+1 ) -z(t))θ n+1 h dxdt ≤ α 2 C 13 k 3/2 ρ min ||θ n+1 h || L 2 (Ω) |u n+1 h -u(t n+1 )| H 1 (Ω) ||z ′ || L 2 (t n ,t n+1 ;L 4 (Ω)) ≤ α 2 C 13 γ||z ′ || L ∞ (0,T ;L 4 (Ω)) k 2ε 8 ||θ n+1 h || 2 L 2 (Ω) + ε 8 2 k|u n+1 h -u(t n+1 )| 2 H 1 (Ω) .
We deduce the result after summing over n = 0, . . . , N -1 Theorem 4.3. With the same assumptions of Theorem 4.2, we have :

N -1 n=1 k||z n+1 h -z(t n+1 )|| 2 L 2 (Ω) ≤ F 1 (h 3 + k 2 ) + F 2 N -1 n=0 k|u n+1 h -u(t n+1 )| 2 H 1 (Ω) , (30)
where F i are constants that only depend on u, z, Ω and T .

Proof. We consider (21), take the test function

θ h = θ n+1 h = z n+1 h -ρ h z(t n+1
), integrate from t n to t n+1 and subtract (26) multiplied by k. We obtain :

α (z n+1 h -z(t n+1 )) -(z n h -z(t n )), θ n+1 h + ν t n+1 t n (z n+1 h -z(t))dt, θ n+1 h + k c(u n+1 h ; z n+1 h , θ n+1 h ) -α t n+1 t n u(t)∇z(t) + 1 2 div u(t)z(t), θ n+1 h dt = ν t n+1 t n (curlu n+1 h -curlu(t))dt, θ n+1 h (31)
Let us treat each term of this equation that we denote by (b i ), i = 1, ..., 4.

For the first term, we follow the same steps as for the term (a 1 ) in the Theorem 4.1. We obtain :

(b 1,1 ) = α 2 ||θ n+1 h || 2 L 2 (Ω) - α 2 ||θ n h || 2 L 2 (Ω) + α 2 ||θ n+1 h -θ n h || 2 L 2 (Ω) , and
|(b 1,2 )| = α t n+1 t n (ρ h z ′ (τ ) -z ′ (τ ))dτ, θ n+1 h ≤ Cαh 3 2ε 9 ||z ′ || 2 L 2 (t n ,t n+1 ;W 1,4 (Ω)) + αε 9 2 k||θ n+1 h || 2 L 2 (Ω) .
For the second term (b 2 ), we write :

z n+1 h -z(t) = z n+1 h -ρ h z(t n+1 ) + ρ h z(t n+1 ) -ρ h z(t) + ρ h z(t) -z(t),
and we obtain three parts that we treat successively.

The first one gives : (b 2,1 ) = νk||θ n+1 h || 2 L 2 (Ω) . The second part is bounded as follows :

|(b 2,2 )| = ν t n+1 t n (ρ h z(t n+1 ) -ρ h z(t))dt, θ h = ν t n+1 t n ρ h z ′ (t)(t -t n )dt, θ n+1 h ≤ νk 2 2 √ 3ε 10 ||z ′ || 2 L 2 (t n ,t n+1 ;H 1 (Ω)) + ε 10 2 √ 3 k||θ n+1 h || 2 L 2 (Ω) ,
and the last part is bounded as follows :

|(b 2,3 )| = ν t n+1 t n (ρ h z(t) -z(t))dt, θ n+1 h ≤ νh 3 2ε 11 ||z|| 2 L 2 (t n ,t n+1 ;W 1,4 (Ω)) + νkε 11 2 ||θ n+1 h || 2 L 2 (Ω) .
The third term can be written as follows :

(b 3 ) = t n+1 t n c(u n+1 h ; z n+1 h , θ n+1 h )dt -α t n+1 t n u(t)∇z(t) + 1 2 div u(t)z(t), θ n+1 h dt = t n+1 t n c(u n+1 h ; θ n+1 h , θ n+1 h )dt + t n+1 t n c(u n+1 h ; ρ h z(t n+1 ) -z(t), θ n+1 h )dt + t n+1 t n c(u n+1 h ; z(t), θ n+1 h )dt -α t n+1 t n u(t)∇z(t) + 1 2 div u(t)z(t), θ n+1 h dt.
Owing to Lemma 3.1 and denoting ξ h = ρ h (z(t n+1 )), (b 3 ) becomes :

(b 3 ) = α 2 κ∈T h t n+1 t n ∂κ - |u n+1 h • n|(θ n+1 h,ext -θ n+1 h,int ) 2 dsdt + t n+1 t n c(u n+1 h ; ρ h z(t n+1 ) -z(t), θ n+1 h )dt +α t n+1 t n Ω (u n+1 h -u(t))∇z(t)θ n+1 h dxdt + α 2 t n+1 t n Ω div(u n+1 h -u(t))z(t)θ n+1 h dxdt.
We divide (b 3 ) into four terms (b 3,i ), i = 1, . . . , 4. We keep the term (b 3,1 ) in the left-hand side of (31). The second term (b 3,2 ) is bounded as in the previous theorem.

For the third part (b 3,3 ), we have :

|(b 3,3 )| = α t n+1 t n Ω (u n+1 h -u(t))∇z(t)θ n+1 h dxdt ≤ αC 14 ||z|| L ∞ (0,T ;W 1,4 (Ω)) k 2ε 12 |u n+1 h -u(t n+1 )| 2 H 1 (Ω) + ε 12 2 k||θ n+1 h || 2 L 2 (Ω) +αC 15 ||z|| L ∞ (0,T ;W 1,4 (Ω)) k 2 2ε 13 ||u ′ || 2 L 2 (t n ,t n+1 ;H 1 (Ω) 2 ) + ε 13 2 k||θ n+1 h || 2 L 2 (Ω) ,
and the last part of the third term (b 3 ) is :

|(b 3,4 )| = α 2 t n+1 t n Ω div(u n+1 h -u(t))z(t)θ n+1 h dxdt ≤ α 2 C 16 ||z|| L ∞ (0,T ;L ∞ (Ω)) k 2ε 14 |u n+1 h -u(t n+1 )| 2 H 1 (Ω) + ε 14 2 k||θ n+1 h || 2 L 2 (Ω) + α 2 C 17 ||z|| L ∞ (0,T ;L ∞ (Ω)) k 2 2ε 15 ||u ′ || L 2 (t n ,t n+1 ;H 1 (Ω) 2 ) + ε 15 2 k||θ n+1 h || 2 L 2 (Ω) .
For the last term (b 4 ), we split it into two parts, as follows :

Using || curl u h || 2 L 2 (Ω) ≤ 2|u h | 2 H 1 (Ω) , we have |(b 4,1 )| = ν t n+1 t n ( curl u n+1 h -curl u(t n+1 ))dt, θ n+1 h ≤ νk ε 16 |u n+1 h -u(t n+1 )| 2 H 1 (Ω) + νε 16 k 2 ||z n+1 h -ρ h z(t n+1 )|| 2 L 2 (Ω) , and
|(b 4,2 )| = ν t n+1 t n ( curl u(t n+1 ) -curl u(t))dt, θ n+1 h ≤ νk 2 √ 3ε 17 ||u ′ || 2 L 2 (t n ,t n+1 ;H 1 (Ω) 2 ) + νε 17 k √ 3 ||z n+1 h -ρ h z(t n+1 )|| 2 L 2 (Ω) .
Collecting all these results, we obtain :

(b 1,1 ) + (b 2,1 ) + (b 3,1 ) ≤ |(b 1,2 )| + |(b 2,2 )| + |(b 2,3 )| + |
Then (30) follows easily after the sum over n = 1, . . . , N -1, a suitable choice of ε i , i = 1, . . . , 17 and by applying a triangular inequality to ||z n+1 h z(t n+1 )|| L 2 (Ω) :

||z n+1 h -z(t n+1 )|| L 2 (Ω) ≤ ||z n+1 h -ρ h z(t n+1 )|| L 2 (Ω) + ||ρ h z(t n+1 ) -z n+1 h || L 2 (Ω)
and properties of P h .

Corollary 4.4. Under the assumptions of Theorem 4.1 and Theorem 4.3, and for k sufficiently small, there exist constants C 1 , C 2 and C 3 independent of h and k such that :

N -1 n=0 k|u n+1 h -u(t n+1 )| 2 H 1 (Ω) ≤ C 1 (h 3 + k 2 ), (32)
N -1

n=0 k||z n+1 h -z(t n+1 )|| 2 L 2 (Ω) ≤ C 2 (h 3 + k 2 ), ( 33
)
and sup

n |u n h -u(t n )| 2 H 1 (Ω) ≤ C 3 (h 3 + k 2 ). ( 34
)
Proof. On one hand, we consider (30). On the other hand, the only difference between this proof and that of Theorem 4.1 is the upper bound of the term (a 5,1 ). Here, using the inequality

||u|| 2 L 4 (Ω) ≤ C|u| H 1 (Ω) ||u|| L 2 (Ω) ,
we have :

|(a 5,1 )| ≤ S 4 ε 6 2 ||P h u|| 2 L ∞ (0,T ;H 1 (Ω) 2 ) k||z n h -z(t n )|| 2 L 2 (Ω) + S 4 4ε 6 k||u n+1 h -P h u(t n+1 )|| 2 L 2 (Ω) + S 4 4ε 6 k|u n+1 h -P h u(t n+1 )| 2 H 1 (Ω) .
Then, using this result with (30) and after a suitable choice of ε i , i = 1, . . . , 10 and ε6 , we obtain :

||u N h -P h u(t N )|| 2 L 2 (Ω) + N -1 n=0 ||(u n+1 h -P h (u(t n+1 ))) -(u n h -P h (u(t n )))|| 2 L 2 (Ω) +α|u N h -P h u(t N )| 2 H 1 (Ω) + α N -1 n=0 |(u n+1 h -P h u(t n+1 )) -(u n h -P h u(t n ))| 2 H 1 (Ω) +ν N -1 n=0 k|u n+1 h -P h u(t n+1 )| 2 H 1 (Ω) ≤ C(h 3 + k 2 ) + C ′ N -1 n=0 k||u n+1 h -P h u(t n+1 )|| 2 L 2 (Ω) .
Then by applying the discrete Gronwall lemma, we obtain, for k sufficiently small :

||u N h -u(t N )|| 2 L 2 (Ω) + ν N -1 n=0 k|u n+1 h -u(t n+1 )| 2 H 1 (Ω) ≤ Ce C ′ kN (h 3 + k 2 ),
and the results follow easily.

Remark 4.5. If we have, for example,

v m h 2 L 2 (Ω) + m-1 n=0 v n+1 h -v n h 2 L 2 (Ω) +ν m-1 n=0 k|v n+1 h | 2 H 1 (Ω) ≤ C 1 + C 2 m-1 n=0 k v n+1 h 2 L 2 (Ω) , by writing, v m h L 2 (Ω) ≤ v m h -v m-1 h L 2 (Ω) + v m-1 h L 2 (Ω) , we obtain C 2 k v m h 2 L 2 (Ω) ≤ 2C 2 k v m h -v m-1 h 2 L 2 (Ω) +2C 2 k v m-1 h 2 L 2 (Ω) .
By assuming k sufficiently small such that 2C 2 k ≤ 1, we obtain :

v m h 2 L 2 (Ω) + m-2 n=0 v n+1 h -v n h 2 L 2 (Ω) +ν m-1 n=0 k|v n+1 h | 2 H 1 (Ω) ≤ C 1 + 3C 2 m-1 n=1 k v n h 2 L 2 (Ω) ,
then we can apply the discrete classic Gronwall lemma.

Theorem 4.6. Under the assumptions of Corollary 4.4, we suppose that p ′ ∈ L 2 (0, T ; L 2 (Ω)).

Then the pressure satisfies the following estimate :

N -1 n=0 k||p n+1 h -p(t n+1 )|| 2 L 2 (Ω) ≤ 1 β ⋆ C(h 3 + k 2 ) +(α + S 2 2 ) N -1 n=0 k (u n+1 h -u(t n+1 )) -(u n h -u(t n )) k 2 H 1 (Ω) , (35)
where the coefficients C and S 2 are respectively the inf-sup constant and Poincaré's constant and are independent of h and k.

Proof. We consider again the first equation of the proof of Theorem 4.1, insert ±kr h p(t n+1 ) in the terms of the pression and we get :

t n+1 t n p n+1 h -r h p(t n+1 ), div v n+1 h dt = (u n+1 h -u(t n+1 )) -(u n h -u(t n )), v n+1 h +α ∇(u n+1 h -u(t n+1 )) -∇(u n h -u(t n )), ∇v n+1 h + ν t n+1 t n ∇(u n+1 h -u(t))dt, ∇v n+1 h + t n+1 t n (z n h ∧ u n+1 h -z(t) ∧ u(t))dt, v n+1 h - t n+1 t n r h p(t n+1 ) -p(t), div v n+1 h dt.
Owing to the inf-sup condition, ∀q h ∈ M h ,

∃w h ∈ V ⊥ h ; (div w h , q h ) = ||q h || 2 L 2 (Ω) and ||∇w h || L 2 (Ω) ≤ ||q h || L 2 (Ω) ,
and summing over n = 0, . . . , N -1, the left-hand side of this equation becomes

N -1 n=0 k||p n+1 h -r h p(t n+1 )|| 2 L 2 (Ω) .
Let us treat the terms in the right-hand side of the equation.

For the first term, we have

N -1 n=0 (u n+1 h -u(t n+1 )) -(u n h -u(t n )), v n+1 h ≤ S 2 N -1 n=0 k (u n+1 h -u n h ) -(u(t n+1 ) -u(t n )) k 2 L 2 (Ω) 1/2 N -1 n=0 k|v n+1 h | 2 H 1 (Ω) 1/2 ≤ S 2 2 N -1 n=0 k (u n+1 h -u n h ) -(u(t n+1 ) -u(t n )) k 2 H 1 (Ω) 1/2 N -1 n=0 k|v n+1 h | 2 H 1 (Ω) 1/2
, and for the second,

α N -1 n=0 ∇((u n+1 h -u(t n+1 )) -(u n h -u(t n ))), ∇v n+1 h ≤ α N -1 n=0 k (u n+1 h -u n h ) -(u(t n+1 ) -u(t n )) k 2 H 1 (Ω) 1/2 N -1 n=0 k|v n+1 h | 2 H 1 (Ω) 1/2
.

For the third term, we have :

ν N -1 n=0 t n+1 t n ∇(u n+1 h -u(t)), ∇v n+1 h dt ≤ ν N -1 n=0 t n+1 t n ∇(u n+1 h -u(t n+1 )), ∇v n+1 h dt + ν N -1 n=0 t n+1 t n t t n+1 ∇u ′ (τ )dτ, ∇v n+1 h dt ≤ ν C 1 (h 3 + k 2 ) 1/2 + C 2 k||u ′ || L 2 (0,T ;H 1 (Ω) 2 ) N -1 n=0 k|v n+1 h | 2 H 1 (Ω) 1/2
.

The fourth term is treated as the fifth term in the proof of Theorem 4.1 and by using the result of Theorem 4.3, the result is the following :

N -1 n=0 t n+1 t n z n h ∧ u n+1 h -z(t) ∧ u(t) dt, ∇v n+1 h ≤ C(h 3 + k 2 ) 1/2 N -1 n=0 k|v n+1 h | 2 H 1 (Ω) 1/2 .
Finally, the last term is treated as follows :

t n+1 t n r h p(t n+1 ) -p(t), div v n+1 h dt ≤ N -1 n=0 t n+1 t n (r h p(t n+1 ) -r h p(t))dt, div v n+1 h + N -1 n=0 t n+1 t n (r h p(t) -p(t))dt, div v n+1 h ≤ C N -1 n=0 t n+1 t n (τ -t n )||p ′ (τ )|| L 2 (Ω) |v n+1 h | H 1 (Ω) dτ + C 2 N -1 n=0 h 2 k 1/2 ||p|| L 2 (t n ,t n+1 ;H 2 (Ω)) |v n+1 h | H 1 (Ω) ≤ C 1 k √ 3 ||p ′ || L 2 (0,T ;L 2 (Ω)) + C 2 h 2 ||p|| L 2 (0,T ;H 2 (Ω)) N -1 n=0 k|v n+1 h | 2 H 1 (Ω) 1/2 .
Then (35) follows easily from these inequalities.

We still have to estimate

N -1 n=0 k (u n+1 h -u(t n+1 )) -(u n h -u(t n )) k 2 H 1 (Ω) 1/2 .
We introduce a variant of Stokes projection as follows :

∀(u, p) ∈ V × L 2 0 (Ω), S h u ∈ V h is defined by ∀v h ∈ V h , ν(∇(S h u -u), ∇v h ) = -(p, div v h ), (36)
Lemma 4.7. Let (u, p) ∈ V × L 2 0 (Ω). Then S h u defined by (36) satisfies :

|S h u -u| H 1 (Ω) ≤ 2|P h u -u| H 1 (Ω) + 1 ν ||r h p -p|| L 2 (Ω) . (37)
If, in addition, Ω is convex, there exists a constant C independent of h such that:

||S h u -u|| L 2 (Ω) ≤ Ch |S h u -u| H 1 (Ω) + ||r h p -p|| L 2 (Ω) . (38)
Theorem 4.8. Under the assumptions of Theorem 4.6 and assuming p ′ ∈ L 2 (0, T ; H 2 (Ω)), we have :

N -1 n=0 k (u n+1 h -u(t n+1 )) -(u n h -u(t n )) k 2 H 1 (Ω) ≤ C(h 3 + k 2 ). (39)
Proof. We consider, once more, the first equation in the proof of Theorem 4.6, choose v n+1 h ∈ V h , insert r h p(s) and S h u ′ = (S h u) ′ and we set e n h = u n h -S h u(t n ). We obtain :

(u n+1 h -S h u(t n+1 )) -(u n h -S h u(t n )), v n+1 h - t n+1 t n u ′ (s) -S h u ′ (s), v n+1 h ds +α ∇(u n+1 h -S h u(t n+1 )) -∇(u n h -S h u(t n )), ∇v n+1 h -α t n+1 t n ∇(u ′ (s) -S h u ′ (s)),

  We sum this above equation over n = 0, . . . , N -1 and we treat the terms denoted again(a i ), i = 1, . . . , 6 in the left-hand side. We take v n+1 -S h u ′ || L 2 (t n ,t n+1 ;L 2 (Ω) 2 ) k|| -S h u ′ || 2 L 2 (0,T ;L 2 (Ω) 2 ) + ||u ′ || 2 L 2 (0,T ;H 2 (Ω) 2 ) + ||p ′ || 2 L 2 (0,T ;H 1 (Ω)) +The second term (a 2 ) is treated as the first one. We have :(a 2 ) = α Cα 2ε 2 h 4 ||u ′ || 2 L 2 (0,T ;H 3 (Ω) 2 ) + ||p ′ || 2 L 2 (0,T ;H 2 (Ω)) + || 2 L 2 (Ω)and by defining S h u 0 with p = 0, we obtain :(tt n )|S h u ′ (t)| H 1 (Ω) dt |v n+1 h | H 1 (Ω) ||S h u ′ || L 2 (t n ,t n+1 ;H 1 (Ω) 2 ) ||S h u ′ || 2 L 2 (0,T ;H 1 (Ω) 2 ) .Using the definition of S h , we have(a 4 ) + (a 6 ) = νFinally, the last term (a 5 ) is bounded as previously in Theorem 4.6 :Theorem 4.9. Under the assumptions of Theorem 4.6, there exists a constant C that does not depend on h and k such that

	∇v n+1 h L 2 (Ω) . || 2 ds = 0 -e n ∇(S h u(s) -u(s))ds, ∇v n+1 h t n+1 t n+1 t n t n p(s), div v n+1 h ds = 0. h = e n+1 h -e n h k . The first term is composed of two parts as follows : +ν t n+1 t n ∇(u n+1 h -S h u(s))ds, ∇v n+1 h + ν + t n+1 t n (z n h ∧ u n+1 h -z(s) ∧ u(s))ds, v n+1 h + where (a 1,1 ) = N -1 n=0 k|| e n+1 h -e n h k L 2 (Ω) , and |(a 3,2 )| = ν N -1 n=0 t n+1 t n (∇S h (u(t n+1 ) -u(s)), ∇v n+1 h )ds || 2 |(a 1,2 )| ≤ N -1 n=0 ||u ′ e n+1 h -e n h k || 2 L 2 (Ω) 1/2 ≤ 1 2ε 1 ||u ′ ε 1 2 N -1 k|| e n+1 h -e n h k || 2 = ν N -1 t n+1 t n h )dt (∇S h u ′ (t)(t -t n ), ∇v n+1 n=0 ≤ ν N -1 t n+1 |v n+1 h | H 1 (Ω) t n n=0 L 2 (Ω) N -1 n=0 ≤ C 2ε 1 h 4 ε 1 2 N -1 n=0 k|| e n+1 h h k N -1 n=0 (∇e n+1 h -∇e n h , ∇v n+1 h )-α N -1 n=0 t n+1 t n ≤ ν n=0 k 3/2 ≤ ν 2ε 3 N -1 n=0 k| e n+1 h -e n h k | 2 H 1 (Ω) + νε 3 2 k 2 t n+1 t n ∇(S h u(s) -u(s))ds, ∇v n+1 h + t n+1 t n p(s), div v n+1 h ∇u N -1 n=0 k| e n+1 h -e n h k | 2 H αε 2 2 N -1 n=0 k|| e n+1 h -e n h k || 2 H 1 (Ω) . The third term is treated as : |(a 3 )| = ν N -1 n=0 t n+1 t n (∇(u n+1 h -S h u(s)), ∇v n+1 h )ds = ν N -1 n=0 (∇e n+1 h , ∇v n+1 h ) + ν N -1 n=0 t n+1 t n (∇S h (u(t n+1 ) -u(s)), ∇v n+1 h )ds = |(a 3,1 ) + (a 3,2 )|. Using the relation (a n+1 -a n , a n+1 ) = 1 2 ||a n+1 || 2 L 2 (Ω) -1 2 ||a n || 2 L 2 (Ω) + 1 2 |(a 5 )| ≤ Cε 4 2 (h 3 + k 2 ) + C 2ε 4 N -1 n=0 k| e n+1 h -e n h k | 2 H N -1 ||a n+1 -a n |(a 3,1 )| = ν 2 N -1 n=0 k| e n+1 h -e n h k | 2 H 1 (Ω) + ν 2 N -1 n=0 k | e n+1 h k | 2 H 1 (Ω) -| h k H 1 (Ω) | 2 e n n=0 k||p n+1 h -p(t n+1 )|| 2 L 2 (Ω) ≤ C(h 3 + k 2 ). (40)	ds
	(a 1 ) =	N -1 n=0 (e n+1 h =	-e n h , v n+1 h ) -ν 2 N -1 n=0 k| e n+1 N -1 n=0 h -e n h k	t n+1 t n | 2 H 1 (Ω) + u ′ (s) -S h u ′ (s), v n+1 h ν 2 k| e N h k | 2 H 1 (Ω) .	ds = (a 1,1 ) + (a 1,2 ),

′ (s)-∇S h u ′ (s), ∇v n+1 h ds = (a 2,1 )+(a 2,2 ), with (a 2,1 ) = α 1 (Ω) , |(a 2,2 )| ≤ 1 (Ω) .

Collecting these results, writing (a 1,1 ) + (a 2,1 ) + (a 3,1 ) ≤ |(a 1,2 )| + |(a 2,2 )| + |(a 3,2 )| + |(a 5 )|, choosing suitably ε i , i = 1, . . . , 4 and by applying the following triangular inequality

|u n+1 h u(t n+1 )| H 1 (Ω) ≤ |u n+1 h -S h u(t n+1 )| H 1 (Ω) + |S h u(t n+1 )u(t n+1 )| H 1 (Ω) ,

(39) follows easily.

i=2 (b 3,i )| + |(b 4 )|.


Sebastian Seeber
email: sebastian.seeber@unibw.de

Gabi Dreo
email: gabi.dreo@unibw.de

Towards an adaptive and effective IDS using OpenFlow

Processing huge amounts of traffic from core network components with respect to security remains a challenging task, since the amounts of data increase continuously. Therefore, new approaches need to be investigated to detect and handle attacks already in high-speed environments. In this PhD research, we will develop a new approach for detecting network attacks by processing data from core network components taking advantage of properties of OpenFlow in an SDN environment. Using this, we can collect metadata about forwarded traffic in an immediate and effective way. In addition, our solution will enable dynamic and adaptive redirection of traffic to various IDSs including cloud-based IDS solutions.

Introduction & Motivation

Security at network level states an important research area since consumers and companies push their data continuously into cloud environments [START_REF] Morris | Infographic: SMB Cloud Adoption Trends in 2014[END_REF]. A reason for this evolution is the growing popularity of cloud services as well as simplicity and fast dynamic expandability of resources on demand. In addition, distributed denial-of-service (DDoS) attacks increased dramatically during the last years [2]. A recent report from Akamai [START_REF]State of the Internet -Security Report[END_REF] shows 90 % growth of DDoS attacks in the last 12 month. The advent of software defined networking (SDN) promises a large variety of possibilities to improve network monitoring and traffic steering. Nevertheless, widespread deployments using SDN principles are still rare and mostly static (proactive) configured.

Network security solutions are primarily based on inspecting traffic. The most well-known approaches are on the one hand analyzing whole packets using deep packet inspection (DPI) and on the other hand using rather statistical data provided by NetFlow/IPFIX techniques [START_REF] Fry | Security Monitoring[END_REF]. The operational environment affects the chosen traffic inspection method. Environments carrying fewer amounts of data (average bandwidth from few Mb/s to 1 Gb/s) are suitable for DPI [START_REF] Abuhmed | A survey on deep packet inspection for intrusion detection systems[END_REF] whereas high-throughput (up to 100 Gb/s) environments like Internet Exchange Points (IXP) have only a chance inspecting traffic using NetFlow/IPFIX [START_REF] Sperotto | An Overview of IP Flow-Based Intrusion Detection[END_REF] and even sampled. A drawback of the latter approaches lays in their design, because statistics are calculated after a flow is terminated (by flow aging, TCP session termination or fixed interval). Therefore, these introduce a delay in a subsequent detection process.

In contrast, OpenFlow is able to raise an event or update a flow counter at arrival time of a packet depending on a match or mismatch with respect to an existing or non-existing flow. Deployments of commercial existing intrusion detection systems are mostly implemented in a static manner. If multiple Intrusion Detection Systems (IDSs) exist, traffic redirection is mainly based on subnets or IP addresses.

Our proposed effective IDS uses OpenFlow as a key-enabler and is adaptive since it involves multiple IDSs on-demand by taking into account immediate detection results. Keeping these thoughts in mind, the following research questions arise and need to be investigated: RQ1 How to use OpenFlow to gather statistics to detect attacks/anomalies? RQ2 How to find an optimal trade-off between performance and detection rate? RQ3 How to classify attack and reason suitable IDS and their sequence? RQ4 How to compare cloud security solutions and how to validate them? RQ5 How to verify the chosen path (Service Chain Verification)?

The remainder of the paper is organized as follows. Section 2 discusses related work. We introduce our approach in Section 3. Concluding remarks and future work are described in Section 4.

Related Work

The use of SDN concepts towards network security is not new. Kreutz et al. [START_REF] Kreutz | Towards secure and dependable softwaredefined networks[END_REF] argue for building dependable and secure SDN applications. Therefore, they identified and described current threat vectors in SDN environments that could be exploited and propose a general design to overcome these threats. Besides this Scott et al. [START_REF] Scott-Hayward | SDN security: A survey[END_REF] investigated possible new security issues introduced through SDN and identified affected layers. François et al. [START_REF] François | Network Security through Software Defined Networking: a Survey[END_REF] reviewed SDN security approaches according to their scope, practicability and advantages.

Research also has been done focusing on more specific attacks and their mitigation. Using self organizing maps the authors of [START_REF] Braga | Lightweight DDoS flooding attack detection using NOX/OpenFlow[END_REF] propose a method to detect DDoS attacks based on flow analysis. Feamster et al. [START_REF] Feamster | Outsourcing home network security[END_REF] and Schehlmann et al. [START_REF] Schehlmann | COFFEE: a Concept based on OpenFlow to Filter and Erase Events of botnet activity at high-speed nodes[END_REF] investigated possibilities to detect botnet traffic by using distributed monitoring approaches. Combining traditional network features (sFlow) and Open-Flow, Giotis et al. [START_REF] Giotis | Combining OpenFlow and sFlow for an effective and scalable anomaly detection and mitigation mechanism on SDN environments[END_REF] proposed a mechanism to detect anomalies and mitigate attacks by modifying flow tables.

Approach

The research in this PhD presents a new approach in providing security on the network level. To overcome existing static security function deployments, our approach provides security functions in a dynamic and free composable manner. Traditional approaches follow a static sequence of different security functions based on parameters like state of protective system, anticipated attack probability or expected risk-level. Instead, the mechanism concatenates security functions (e.g., various specialized IDSs) based on immediate detection results. OpenFlow is a well-established protocol in SDN deployments and therefore the most suitable enabler for our approach.

Overall Architecture Our proposed overall architecture is depicted in Figure 1. Starting from left, the input traffic reaches a network core OpenFlow enabled switch (OF-Switch). It is important that this switch is as near as possible to an IXP to cope with large amounts of attack traffic. On top of this switch a light IDS (OF-IDS) is implemented using OpenFlow counters (Meter-, Group-, Flow-Table -OpenFlow 1.3) as input data. RQ1 and RQ2 are related to the development of an effective IDS based on monitoring data from an OpenFlow enabled switch.

Triggered by events of OF-IDS the connected SDN-Controller (SDN-C) processes the event provided by the IDS and information from the associated flow in order to decide about the subsequent detection step. RQ3 will answer this question by using traffic and attack characterization techniques to introduce adaptiveness. After a decision is taken by the SDN-C, OpenFlow rules are applied to the underlying OF-Switch meeting the decision requirements. The procedure between incoming traffic at the OF-Switch and SDN-C decision (including applying OpenFlow rules) will be called detection cycle (DC). The described behavior repeats at every point an IDS (IDS-A, IDS-B, IDS-C) is present. Various DCs are conceivable, depending on the number of available IDSs and detection results. In cases where not only on-site IDSs, like Suricata [START_REF]Suricata IDS/IPS[END_REF], SNORT [START_REF]Snort[END_REF] or Bro [3] are considered, the architecture is flexible in including cloud-based IDS solutions (Cloud-IDS), e.g., Cloudflare [START_REF] Cloudflare | [END_REF]. For example, the red-dotted line shows a case where only on-site IDSs are involved, whereas the green path involves a cloudbased IDS. Therefore, RQ4 investigates techniques to evaluate cloud security solutions (cloud-based IDS) and comparison criteria. It's important to mention, that only attack traffic is redirected to cloud-based IDS solutions. Therefore, privacy aspects in using a cloud-based IDS are minimized compared to solutions that redirect the whole traffic to the cloud.

In the first step each SDN-C will decide autonomously, without including previous IDS results or decisions from earlier SDN-Cs. Later a connection between them will be established to enhance detection capabilities and design the SDN-Cs as a distributed controller. Use-Cases Meeting requirements of various deployment scenarios enhances applicability of our proposed architecture and helps to improve the continuous development process. However, we will start implementing a very narrow-focused deployment and enhance it based on our findings and additional use-cases.

-Software as a Service (SaaS) products like web-hosting, e-mail or financial management provider are usually accessible from all over the world. In addition, these services seem to be easily configurable, also by people not familiar with security principles. -Platform as a Service (PaaS) providers in most cases access to a preconfigured machine regardless of the underlying infrastructure (virtual/ physical). Users of PaaS require a rare understanding of configuring these systems and knowledge in securing them, but gain significant freedom in choosing applications running on top. -Company's Internet Access (CIA) Companies providing services to customers and make Internet access available to employees are responsible for (amongst others) securing their network against attacks or attackers trying to gather customer data. Furthermore, they have to satisfy their customers by meeting SLAs as contracted (e.g., availability, response time).

Service Chaining Composing a sequence of successive security functions raises immediately concerns about (a) trustworthiness of each security function and (b) the right and entire processing of a composed chain. Therefore, our approach will utilize the benefits of SDN and distributed controllers to establish a verifiable chain of security functions. RQ5 is the core of this aspect and enables adaptability.

Concluding Remarks & Perspectives

Since the importance of security raises if more and more companies push their data and processing capabilities into cloud services, our approach tries to provide a solution to inspect traffic by various IDSs including cloud-based IDSs. Finally, it will redirect or block attack traffic at the first point of occurrence in the monitored network. Our proposed solution is also relevant since DDoS attacks are getting more and more popular. As a next step, we need to investigate statistics derived from OpenFlow in order to detect attacks. At this point a trade-off between detection rate and performance of the OpenFlow device is essential to provide an effective solution. Furthermore, a mechanism needs to be investigated that decides about the next IDS in the detection process, to be adaptive. Later an evaluation needs to be done regarding the involvement of cloud-based IDSs. Therefore, we will investigate mechanisms to prove cloud-based IDSs performance and reliability. An interesting part of research is the verification of the path the traffic took. In that context, an investigation of service chaining and verification techniques needs to take place.

We envision, as extension of the work in this PhD, a thorough study of SLA restrictions and verification of legal requirements in introducing cloud-based IDS solutions needs to be investigated.

Fig. 1 :

 1 Fig. 1: Architecture Overview including Example Paths

Acknowledgment

This work was partly funded by FLAMINGO, a Network of Excellence project (ICT-318488) supported by the European Commission under its Seventh Framework Programme.


Jean Mailhé

Jean-Marc Linares
email: linares@iut.univ-aix.fr

Jean-Michel Sprauel
email: sprauel@iut.univ-aix.fr

VERIFICATION BY VIRTUAL GAUGE USING A STATISTICAL CRITERION

Keywords: measurement, verification, virtual gauge, interference probability map

For the most part, metrology software is currently based on the measurement of distances or angles between geometrical elements. If this method of verification is well adapted to geometrical specification without a virtual state modifier, this is not appropriate for specification based on envelop zone such as in maximal matter condition, for example, for the ISO2692 standard. Usually the least squares best fit method is used to estimate derived surfaces, but the statistical information contained in the acquired coordinates remains under-used. The aim of this paper is to present a new approach for the verification of a part, based on a virtual gauge and using a statistical criterion.

INTRODUCTION

The guiding principle of most verification processes is an assembly test between a set of points or its derived surface and a tolerance zone built from geometrical specifications. Hence the geometrical verification is articulated around three important points: the measurement process, part characterization and geometrical specification interpretation.

The measurement process has been the source of numerous research projects for the last decade. However, it will be assumed that the verification process will start with one or a set of real surfaces measurements taken on a classical coordinate measuring machine (CMM). Current metrology software is based on the measurement of distances or angles between geometrical elements. Usually, the least squares best fit method is used to associate a perfect feature to the acquired coordinates. However, the statistical information contained in the set of points is generally not put to use.

In this paper, a geometrical verification by virtual gauge using a statistical criterion is presented. The statistical information contained in a set of points representing the real surface, permits through uncertainty propagation to define the probability density of the position of the matter around the mean associated surface. Through uncertainty propagation method, the probability density function of the position of the matter around the mean associated surface can be defined by the statistical information contained in a set of points. On the other hand, a gauge model offers a closed integration domain which is perfectly adapted for the calculation of interference probability.

THE ISSUE OF GEOMETRICAL VERIFICATION

As mentioned in the introduction, the geometrical verifications issue could be summarized in three points (figure 1). First, the measurement of the real surfaces (here done on with CMM) will provide a set of digitized points. Secondly, the part characterization is achieved through a best fitting to obtain derived features or integral associated surfaces. Finally, the geometrical specifications and the building of tolerance zones (TZ) must be analyzed in order to perform a conformance test between the best fitted surfaces or the set of measured points, and the tolerance zones.

Part measurement process

Geometrical specification

Part characterization

Geometrical verification

Fig.1. Geometrical verification issue

Part measurement process and part characterization

Several research projects have been carried out on measurement processes. First kind of work has been performed on the improvement of the use of measurement instruments [START_REF] Weckenmann | Monitoring coordinate measuring machines by user-defined calibrated parts[END_REF]. This improvement could be made with more efficient methods of calibration [START_REF] Antunes | Estimation of precision and uncertainty of a calibration artefact for CMMS[END_REF]; other work proposes methods to decrease measurement errors [START_REF] Zhang | Error compensation of coordinate measuring machines[END_REF]. Then work has been done sampling strategy [START_REF] Edgeworth | Uncertainty Management for CMM Probe Sampling of Complex Surfaces[END_REF].

Nowadays measurements by CMM are widely used in manufacturing industries. The data supplied by CMM a set of acquired coordinates. To estimate an associated surface representing the real part surface, a best fitting must be done. This association is made by optimizing a criterion like least squares, infinite norm or likelihood functions [START_REF] Henke | Methods for evaluation of systematic geometric deviations in machined parts and their relationships to process variables[END_REF] [START_REF] Takamasu | Basic concept of feature-based metrology[END_REF]. Part characterization can be sorted into two categories: a best fit of a complex surface including classical shape error, or best fit of a simple feature and an estimation of statistical data. A method based on a likelihood function proposes a non-linear optimization association [START_REF] Sprauel | Contribution of nonlinear optimization to the determination of measurement uncertainties[END_REF].

Nevertheless, it is well known in measurement that the sampling of the real surface and hence the best fitting, imply that the estimation of the surface is provided with uncertainties.

Geometrical Dimensioning & Tolerancing

A geometrical specification describes a set of geometrical conditions which must be met by a set of real surfaces composing the workpiece. The geometrical product specification system (GPS) provides a set of efficient tools for univocal specification descriptions and verification processes [START_REF]Geometrical product specifications (GPS) -model for geometric specification and verification[END_REF]. These tools are based on the concept of tolerance zone, i.e. a domain of the 3D space related or not to a datum reference surface and where the real surface must be contained. However if international standards ISO provide accurate definitions for geometrical specifications, it is necessary to perform a mathematical translation of this requirement adapted to the measurement data. Aiming toward coherence between geometrical specifying and verifying, several mathematical-based models of specifications have been proposed. One approach involves a classification of symmetry groups and its impact on functional feature taxonomy, datum definition and parameterization [START_REF] Srinivasan | Geometrical Product Specification Language Based on a Classification of Symmetry Groups[END_REF]. As tolerance zone can be seen as geometrical boundaries, it is possible to describe the permissible geometrical variations of the parts by virtual gauges with internal mobility [START_REF] Dantan | Functional and product specification by Gauge with Internal Mobilities[END_REF]. Still based on geometrical limit, the concept of Virtual Boundary Requirement (VBR), generalized with Maximum Material Condition (MMC) and Least Material Condition (LMC), is now widely used.

Conformance test

The conformance test is a test which will validate that the specified surface of the part meets the geometrical requirement. The pre-condition of a conformance test is the characterization of the part and the interpretation of the geometrical specifications. In practical terms, the conformance test will determine if the estimation of the derived feature or the set of points can be contained in the tolerance zone. Most measurement software, in agreement with normalized specifications (ISO standard), is based on the measurement of distances or angles between geometrical elements. These geometrical quantities can be described as the following types:

-point/point distance, -point/plane distance, -point/line distance, -line/line angle, -line/plane angle, -plane/plane angle.

This will result in a set of inequations representing the tolerance zone, which must be satisfied. With the current industrial landscape, the competition in products manufacturing is becoming stiffer and stiffer due to the globalization and subcontractors management. This implies the fact that product manufacturing must always must continue to accelerate with a higher level of quality. Hence verification should be fast and reliable. As said earlier, tolerance zone can be seen as a space boundary, especially when modifiers of the envelope condition, maximum or least material condition, are specified. Consequently, gauges are actually a good solution to this issue. They are one of the most reliable tools for geometrical verifications. They can directly validate the largest part of assembly and functional requirements. There are two kinds of gauges: virtual gauges and physical gauges. Physical gauges have been used for a century. They permit extremely fast and reliable validation of functional requirements, such as assembly requirements. They are involved in specification verifying [START_REF] Bin | Polygon-based large diameter measurement with modular gauges[END_REF] [12] but also in machine tools and CMM [START_REF] Antunes | Estimation of precision and uncertainty of a calibration artefact for CMMS[END_REF] calibration. A virtual gauge is a set of ideal virtual features. It permits a comparison between the real part and a model of the nominal design. The goal of virtual gauge is to be defined as soon as possible in the conception process of a product, the product's verifier gauge and then to compare it to the obtained real part. With the improvement of CAD software and CMM, virtual gauges are showing all their potential.

CLASSICAL & STATISTICAL POINTS OF VIEW

If the concept of physical and virtual gauge is still the same, i.e., testing if the matter is completely inside/outside a spatial boundary, there is, however, different ways to perform it. The first method is to work directly on the real part with a material gauge. The second works on the set of digitized points without any kind of best fitting or on a derived feature element given by the best fit method.

Verification by virtual gauge without best fit

The aim of this method is to find the position of the virtual gauge where it will include all the set of points (figure 2.). Of course, this position shall be in agreement with the degrees of freedom given by the geometrical specifications.

Digitized points Virtual gauge

Datum reference di Mi

Fig.2. Cloud of point verification

From a mathematical point of view, an algebraic distance d i will be calculated between each digitized point M i and the surface of the virtual gauge. The conformance test will be to check if the whole distances d i are positive or negative, i.e. if the points are inside or outside the virtual gauge, according to geometrical specifications. The main advantage of this kind of conformance test is to avoid geometrical construction and hence to avoid uncertainties propagation.

Verification with associated feature

The principle behind this method is to perform an assembly test between an estimated and/or constructed feature and the virtual gauge.

Integral associated surface Virtual gauge

Datum reference

Feature virtual state di

Fig.3. Verification with virtual state of the matter

In order to perform the conformance test, a surface representing the virtual state of the feature is constructed. The estimation of the virtual state results from the addition of the perpendicularity error, and the estimate dimension of the extreme fit. Then, to be accepted, the virtual state dimension must be inferior or superior to the virtual gauge dimension, i.e. it must be inside or outside the matter according to the geometrical specifications.

The virtual gauge could also have a fitter behavior [START_REF] Pairel | The Gauge model: a new approach for coordinate measurement[END_REF]. In this case, an extreme fitting will be done with the whole set of digitized points. This extreme fitting of the set of features is done in a particular order defined by the geometrical specification. Then the same test as above will be applied between the fitter gauge and the limiter gauge obtained from the geometrical specifications.

Statistical point of view

However, all the previous kinds of verification do not take into account measurement uncertainties. In this case, each point is considered as being a point from the real surface. Actually it is well known fact that every measurement is tainted with uncertainties, and hence the acquired coordinates are a random sample of the real points of the surface.

In the next section, a method taking into account this statistical aspect is used.

STATISTICAL PART CHARACTERIZATION

As shown in the previous section, statistical parameters are unused in most of classical geometrical verification processes. With the method to be described in this section, a three-dimensional scalar field representing the probabilities of being inside the matter can be expressed.

The real surface to be verified is composed of an infinity of points which could be seen as the statistical population to be characterized. The main hypothesis of this method is that the set of digitized points is considered as a statistical sample of the real surface population. The first step of the statistical characterization process is the best fitting of a derived feature and a set of intrinsic parameters with the set of acquired coordinates. In fact, the derived feature is expressed by a random vector representing the mathematical expectation of the surface parameters. To complete the best fitting, the maximization of the likelihood method will be used.

Two sets of values are expressed:

-The first moment order representing the mean value of the derived feature and intrinsic parameter values.

-The second moment representing the variance and covariance of the random vector parameters and intrinsic parameters. These values will be provided on the form of a covariance matrix.

The second step is the uncertainties propagation to the complete associated feature. This feature can be constructed from the first moment order values, i.e. the mean values obtained at the best fitting. The aim of this step is to calculate the uncertainties directly on the complete associated surface. This is done via the propagation of a covariance matrix which is achieved with the following formula (1):

( ) ( ) ( ) ( ) ⇔ ⋅ ⋅ = ∑ T OM J p Σ OM J OM ( ) ( ) ( ) ( ) (residues) σ n v , OC n 2λ n v n λ n OC n OM σ 2 T T 2 T 2 n + ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ + ⋅ ⋅ = ∑ ∑ ∑ (1)
Where J represents the Jacobean operator, Σ the covariance matrix and λ the distance to the barycenter.

Next, the variance of the best fit residues should be added. In practice, the variance at the current point M belonging to the complete associated feature along the feature normal has been calculated in the equation ( 1).

This will bring about two important results: -For a fixed risk α, it is possible to express the statistical envelop containing the whole real surface of the part (figure 4). -It is possible to calculate the probability at any point to be inside the matter. This can be done by an integration of the normal law along the normal of the complete associated feature. The utilization of this result in a verification process will be developed in the next section.

STUDYING A SIMPLE CASE

In order to highlight this statistical approach, a simple case will be studied. In the next two subsections, the case of a geometrical verification problem of a bore in a parallelepiped will be considered. Two different kinds of specifications will be seen: firstly, envelop requirement specifications, and then, orientation specification using the maximum material condition.

Envelop zone specification:

As shown in figure 5, the case of a bore geometrically specified by an envelop requirement is taken into consideration. According to standards (ISO 8015 -1985), the principle of envelop requirement implies the nominally cylindrical real surface of the bore must be outside a perfect cylinder at the state of maximal material (figure 6.). In the case under consideration the envelop diameter is 30.93 mm (figure 6). Moreover, every distance between two points in opposition must be inferior to the maximum allowed distance 30.97 mm, in this case (ISO 14 660-1). However, only the constraint of the envelop zone will be examined in this sample. Basically, the virtual gauge will be a perfect cylinder with 6 degrees of freedom. Thus, the optimization of the gauge position will be gained on the 3 rotations and the 3 translations.

The first step of the verification process is the measurement of the real part. For this case, it has been carried out on a classic CMM. In order to show the influence of the number of digitized points on uncertainties, two sets of coordinate points have been acquired on the same bore: first with 16 points and then with 32 points. Acquired coordinates are provided by the CMM programs in a classic ASCII file of points.

The second part is the best fitting which is done with the maximization of the likelihood function. This step will provide the mean associated surface parameters (position and intrinsic dimension) and their covariance matrix.

Thirdly, the virtual gauge is made according to the constraints of the geometrical specification. Generally, the initial position and the dimension of an elementary virtual gauge should be raised from 3D part's design. An offset operation from nominal geometry is performed. For complex geometrical specification, a global virtual gauge could be made with an assembly of a set of elementary virtual gauges. According to the virtual gauge freedom degrees, a matrix representing the set of allowed solid displacements is generated for each elementary gauge. It will represent degrees of freedom defined in the geometrical specification. This matrix will be used in the optimization problem. As demonstrated in section 4, it is possible to express, for each point of the 3D space, the probability of being inside the matter. This will be used for the conformance test. Assuming the fact that the virtual gauge is a closed surface easily formulated by a parametric equation, a meshing of the gauge should be carried out. Thus, a probability of interference can be calculated for each node of the mesh. The result is given in a matrix permitting via linear interpolation the calculation of interference probability according to the curvilinear coordinate of the gauge.

Fourthly, the optimization of the virtual gauge position is performed aiming to minimize the highest interference probability.

Finally, the drawing of an interference probability map (IPM) is generated for each elementary virtual gauge for the optimal position. Using a grey scale, an IPM represents the interference probability according to the curvilinear coordinate of the gauge. This should be seen as an unfolded texture of the gauge.

Before seeing the statistical geometrical verification result, a classical verification process based on least square and shape default has been made. In the case of the 16 points, a diameter of 30.987 mm and a shape default of 0.021mm have been found. Hence the estimation of the diameter of the smallest envelop tangent to the matter is 30.987-0.021=30.966 mm, which is large enough to accept the part. In the case of the 32 points, a diameter of 30.988 mm and a shape default of 0.018mm have been found. Thus the part should be declared good with an estimation of the diameter of the smallest envelop tangent to the matter of 30.970 mm. It must be noted that these results do not take into account the measurement uncertainties.

In the case where the part is characterized from the set of 16 points, the following IPM is obtained after optimization (figure 7.):

Fig.7.Interference probability map for 16 points

This IPM shows three different zones. The centered one is the well-known zone of the part; here the probability of interference between real surface and virtual gauge is nearly null (below 2%). This zone is located around the barycenter of the digitized point. Next, there are the two extremity zones where the knowledge of the part is the worst (from 2% to 18.2% at edges); here uncertainty on the mean surface increases proportionally with the distance to the barycenter of acquired coordinates.

An obvious fact is the symmetry around the plane perpendicular to the cylinder axis. This shows that the result of gauge position optimization is that the gauge axis is collinear to the mean associated surface axis. This alignment seems foreseeable due to the envelop requirement verification leaving all degrees of freedom on gauge position.

Revolution symmetry is notable. This is a particular case. Indeed the covariance matrix is symmetric and defined positively. With propagation calculation using the Jacobian operator, the theoretical form of a statistical envelop is a hyperboloid with an elliptic base. At the moment only first and second moment order are used for part characterization. Moreover, for the best fit, the hypothesis that every point has the same standard deviation is made forward. Therefore, this symmetry is actually showing that the position of digitized points is symmetric.

In the case where the part is characterized from the set of 32 points, the following IPM is obtained after optimization (figure 8.): The maximal interference probability is 2.6% at the edge. This great decrease of uncertainties is due to the increase in point numbers. However, shape default also decreases from 0.021 mm with 16 points to 0.018 mm.

An important remark is the highlighting of the influence of orientation parameters on uncertainty due to the lever arm effect depending on the distance to barycenter of the set of points. If it is not so important in envelop verification, this will take its entire place in perpendicularity verification.

Perpendicularity specification:

In this subsection, the case of a perpendicularity requirement between a bore and a plane is considered (figure 9.). The plane A is taken as datum reference surface. Once again the case of 16 and 32 points will be studied.

According to geometrical specification standards, the datum reference surface is defined as the perfect plane extreme fitting to the real surface, nominally plane A, and minimizing the highest gap. To be accepted the nominally cylindrical real surface of the bore must be outside a perfect cylinder perpendicular to plane A with a diameter of 30.93 mm (ISO 2692 -1988) (figure 10.). The difference with the first case studied is the constraint put on the virtual gauge. In this case, degrees of freedom are the two translations leaving the reference plane A globally invariant. The surfaces characterization is achieved with two previous sets of digitized points. The datum reference surface is built on a cloud of 16 points. Although it is wrong, this surface will be admitted without uncertainty. Using the classical least square and shape default method, the following results are obtained:

-In the 16 points case the measured diameter is 30.987 mm, the shape default is 0.021 mm and the perpendicular default is 0.029 mm. Hence the diameter of the perfect cylinder representing the state of the maximum of matter is of 30.937 mm. So the part should be accepted and has a clearance of 7 µm.

-In the 32 points case the measured diameter is 30.988 mm, the shape default is 0.018 mm and the perpendicular default is 0.028 mm. Hence the diameter of the perfect cylinder representing the state of the maximum of matter is of 30.942 mm. So the part should be accepted and has a clearance of 12 µm.

In the case where the part is characterized from the set of 16 points, the following IPM is obtained after optimization (figure 11.): With this IPM, the two critical of the verified cylinder can be spotted. In these two zones, the maximal probability for the part to be nonconform is of 32.8%. This high risk is mainly due on the one hand to the uncertainties on cylinder direction and to shape error and on the other hand to the mean value showing a high orientation error.

In the case where the part is characterized from the set of 32 points, the following IPM is obtained after optimization (figure 12.): With a high number of digitized points, uncertainties on direction are highly decreased. Here the probability of interference decreases to 15.8% at critical zone.

CONCLUSION

This paper introduced a new method of part characterization based on statistical estimation of the real surface of the part, on the one hand, and, on the other hand, a new point of view in geometrical specification verification with the interference probability map. For characterization, the use of first and second order statistical moments permits the expression of the probability at any point of the space to be inside the matter. Thus, coupled with the concept of virtual gauge, it is possible to have an accurate statistical point of view of risk of interference between real surface and tolerance zone. The representation of the interference probability map permits a fast estimation of the efficiency of the measurement process planning and hence to correct it for critical cases.

However, the verification of complex parts composed of numerous elementary virtual gauges and a set of freedom degrees implies complex geometrical formulation. Moreover, due to the type of functions to optimize, which are continuous but not derivable (non-smooth optimization), the optimization problem must be solved with unconventional algorithms. A solution could be found in metric tensors and the perturbation matrix. This will be the object of future works.

Fig. 4 .

 4 Fig.4. Statistical limit of the matter

Fig. 5 .

 5 Fig.5. Nominal design and geometrical tolerance

Fig. 8 .

 8 Fig.8. Interference probability map for 32 points

Fig. 9 .

 9 Fig.9. Nominal design and geometrical tolerance

Fig. 11 .

 11 Fig.11. Interference probability map for 16 points

Fig. 12 .

 12 Fig.12. Interference probability map for 32 points