Kadri Souleymane
email: souleymane.kadri-harouna@imag.fr

Harouna

Val Érie Perrier
email: valerie.perrier@imag.fr

DIVERGENCE-FREE AND CURL-FREE WAVELETS ON THE SQUARE FOR NUMERICAL SIMULATIONS

We present a construction of divergence-free and curl-free wavelets on the square, that could satisfy suitable boundary conditions. This construction is based on the existence of biorthogonal multiresolution analyses (BMRA) on [0, 1], linked by differentiation and integration. We introduce new BMRAs and wavelets for the spaces of divergence-free and curl-free vector functions on the square. The interest of such constructions is illustrated on examples including the Helmholtz-Hodge decomposition of vector flows and the Stokes problem.

Introduction

In many physical problems, like the numerical simulation of incompressible flows or in electromagnetism, the solution has to fulfill a divergence-free condition. For the numerical treatment of the relevant equations (Navier-Stokes equation in fluid mechanism or Maxwell's equation in electromagnetism) it is helpful to have at hand bases satisfying a divergence-free or a curl-free condition. In the context of solution schemes for Partial Differential Equations, wavelet bases provide very efficient algorithms, characterized by a reduced computational complexity, with respect to standard methods [START_REF] Cohen | Wavelet methods in numerical analysis[END_REF]. Divergence-free wavelet bases on R d , with compact support, were originally defined by Lemarié-Rieusset in 1992 [START_REF] Lemarié-Rieusset | Analyses multi-résolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF] and applied by Urban to the numerical solution of the Stokes-problem [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF]. In the periodic case, anisotropic divergence-free wavelets have been constructed in [START_REF] Deriaz | Divergence-free and curl-free wavelets in 2D and 3D, application to turbulent flows[END_REF], and firstly used to compute the numerical solution of the Navier-Stokes equation in velocity-pressure formulation [START_REF] Deriaz | Direct Numerical Simulation of Turbulence using divergencefree wavelets[END_REF]. Such numerical scheme requires, at each time-step, the Helmholtz-Hodge decomposition of the nonlinear term, which is no more divergence-free. In Fourier space, this decomposition writes explicitly, whereas in wavelet domain, it can be computed using divergence-free and curl-free wavelets [START_REF] Urban | Wavelet Bases in H(div) and H(curl)[END_REF][START_REF] Deriaz | Orthogonal Helmholtz decomposition in arbitrary dimension using divergence-free and curl-free wavelets[END_REF]. In the general case with physical boundary conditions, it is the key of Navier-Stokes numerical simulations to have at hand an explicit and efficient procedure to compute the Helmholtz-Hodge decomposition of the nonlinear term.

Precisely, the Helmholtz-Hodge decomposition of a vector field u on the square Ω = [0, 1] 2 consists in splitting u into a divergence-free part and a curl-free part [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. A first formulation leading to an orthogonal splitting is the following: there exist a stream function ψ and a potential field q such that:

u = curl ψ + ∇q in Ω curl ψ • ν = 0 on Γ = ∂Ω (1)
where ν is the outward normal to Γ. This decomposition corresponds to the orthogonal splitting of the space (L 2 (Ω)) 2 :

(L 2 (Ω)) 2 = H div,Γ (Ω) ⊕ H ⊥ div (Ω) (2)
where

H div,Γ (Ω) = {u ∈ (L 2 (Ω)) 2 : div(u) = 0, u • ν| Γ = 0}
is the divergence-free function space with velocity tangent to the boundary. It can also be seen as the "curl" space:

H div (Ω) = {u = curl ψ ; ψ ∈ H 1 0 (Ω)} On the other side the space of gradient functions H ⊥ div (Ω) = {∇q ; q ∈ H 1 (Ω)} corresponds to a curl-free function space [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF].

Other types of boundary conditions for the divergence-free space can be considered. For instance, the decomposition:

(L 2 (Ω)) 2 = H div (Ω) ⊕ {∇q ; q ∈ H 1 0 (Ω)} (3)
where now

H div (Ω) = {u ∈ (L 2 (Ω)) 2 : div(u) = 0} = {u = curl ψ ; ψ ∈ H 1 (Ω)}
does not incorporate boundary condition on Γ. In fluid mechanism, this type of boundary condition is less considered since in general Γ corresponds to a physical wall that cannot be crossed by fluid particles (except for porous media). A more useful condition corresponds to the homogeneous Dirichlet boundary condition on Γ, which leads to:

(H 1 0 (Ω)) 2 = H div,0 (Ω) ⊕ H ⊥ div,0 (Ω) (4)
where now H div,0 (Ω) = H div (Ω) ∩ (H 1 0 (Ω)) 2 while H ⊥ div,0 (Ω) is a subspace of H ⊥ div (Ω) ∩ (H 1 0 (Ω)) 2 , see [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] for details. For sake of simplicity, we will focus in this article to the divergence-free spaces being involved in decompositions [START_REF] Beylkin | On the representation of operator in bases of compactly supported wavelets[END_REF][START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF].

Accordingly, the objective of the present paper is to provide multiresolution analyses and wavelet bases of the spaces H div (Ω) and H ⊥ div (Ω). We present in the next section a new construction, based on wavelets on the interval [0, 1] that should satisfy homogeneous boundary conditions, like in [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF]. The key idea of our construction is based on a couple of wavelet bases on the interval, linked by differentiation, like in the theoretical approach of Jouini-Lemarié-Rieusset [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF]. The construction of divergence-free and curl-free approximation spaces and wavelets satisfying suitable boundary conditions are then straightforward. Moreover, our method extends readily to the cube [0.1] d by tensor product [START_REF] Souleymane | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF].

The outline of the paper is as follows. Section 2 details the principles of the construction of divergence-free and curl-free BMRAs and wavelets on the square. Section 3 is dedicated to the description of the divergence-free fast wavelet transform. Finally, section 4 presents two examples of use of these new wavelets in numerical simulations: the Helmholtz decompositions and the Stokes problem.

2 Divergence-free and curl-free wavelets on the square Divergence-free wavelets on the whole space R d have been firstly constructed by Battle-Federbush [START_REF] Battle | Divergence-free vector wavelets[END_REF] in the orthogonal case. Since these previous functions don't have compact support they where not implemented, contrarily to the biorthogonal bases arising from the construction proposed by Lemarié-Rieusset in [START_REF] Lemarié-Rieusset | Analyses multi-résolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF]. Urban was the first who used them in a practical problem, the Stokes problem [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF]. Later, Urban proposed to extend this construction to derive curl-free wavelets [START_REF] Urban | Wavelet Bases in H(div) and H(curl)[END_REF]. An alternative fast decomposition into divergence-free wavelets was proposed by Deriaz-Perrier, based on anisotropic (tensorproduct) wavelets in the periodic case [START_REF] Deriaz | Divergence-free and curl-free wavelets in 2D and 3D, application to turbulent flows[END_REF]. The objective below is to extend these constructions to the square [0, 1] 2 , following the construction principle already exposed in the theoretical work of Jouini-Lemarié-Rieusset [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF].

Construction principle

The construction of divergence-free wavelets on the cube [0, 1] d uses the same arguments as in the whole domain R d [START_REF] Lemarié-Rieusset | Analyses multi-résolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF][START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF]. The key ingredient is to have at hand two onedimensional multiresolution analyses (V 1 j ) and (V 0 j ) of L 2 (0, 1) linked by differentiation:

d dx V 1 j = V 0 j (5)
On the interval [0, 1], following [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF], the biorthogonal spaces should satisfy:

Ṽ 0 j = H 1 0 (0, 1) ∩ x 0 Ṽ 1 j = f : f ∈ Ṽ 1 j and f (0) = f (1) = 0 (6)
The existence of such spaces follows from the fundamental proposition of Lemarié-Rieusset, used at the begining to construct divergence-free wavelets on the whole space R d [20]:

Proposition 2.1 Let (V 1 j (R)
) be a multiresolution analysis (MRA) of L 2 (R), with differentiable and compactly supported scaling function ϕ 1 and associated wavelet ψ 1 . Then there exists a MRA (V 0 j (R)), with associated scaling function ϕ 0 and wavelet ψ 0 , such that:

(ϕ 1 ) (x) = ϕ 0 (x) -ϕ 0 (x -1) and (ψ 1 ) (x) = 4 ψ 0 (x) (7)
Similar relations hold for the dual functions φ1 , ψ1 and φ0 , ψ0 of the primal ones (ϕ 1 , ψ 1 ) and (ϕ 0 , ψ 0 ):

x+1 x φ1 (t) dt = φ0 (x) and ( ψ0 ) (x) = -4 ψ1 (x) (8)
In [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF], Jouini-Lemarié-Rieusset prove that from the scaling functions (ϕ 1 , φ1 ), (ϕ 0 , φ0 ) and wavelets (ψ 1 , ψ1 ), (ψ 0 , ψ0 ) of proposition 2.1, it is possible to exhibit multiresolution analyses on the interval satisfying the relations [START_REF] Cieselski | Spline bases in classical function spaces on compact manifolds, part II[END_REF][START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF].

Our objective in the present paper is first to provide an effective construction of such multiresolution analyses, which enable boundary conditions, fast wavelet algorithms, approximation results and practical computations. This is done in section 2.2.1.

Then the construction of biorthogonal MRAs and wavelets bases of H div (Ω) (with suitable boundary conditions) and H ⊥ div (Ω), are obtained by considering resp. the curl of (

V 1 j ⊗ V 1 j ), and the grad of (V 1 j ⊗ V 1 j
). This will be described in section 2.3.

Construction of biorthogonal MRAs on [0, 1] linked by differentiation

We detail now the construction of spaces (V 1 j , Ṽ 1 j ) and (V 0 j , Ṽ 0 j ) satisfying ( 5) and [START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF]. We proceed in two steps. The first step (construction of (V 1 j , Ṽ 1 j ), section 2.2.1) is classical and based on biorthogonal multiresolution analyses on the interval reproducing polynomials [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF][START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF][START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF][START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF], but it is needed to introduce the sekels. The second step (section 2.2.2) introduces a practical and new way to provide spaces (V 0 j , Ṽ 0 j ) and associated wavelets. It is based on proposition 2.1.

Construction of

(V 1 j , Ṽ 1 j ) on [0, 1] with polynomial reproduction (r, r)
We first recall the definition of a biorthogonal multiresolution analysis (BMRA) on [0, 1] [18, 8]:

Definition 2.1
The sequence (V 1 j , Ṽ 1 j ), j ≥ j min (j min ∈ N * ) is a biorthogonal MRA of approximation order (r, r) on the interval [0, 1] associated to the generators (ϕ 1 , φ1 ), if it satisfies:

(i) V 1 j ⊂ V 1 j+1 , Ṽ 1 j ⊂ Ṽ 1 j+1 and ∪ j≥j min V 1 j = ∪ j≥j min Ṽ 1 j = L 2 (0, 1). (ii) V 1,int j ⊂ V 1 j ⊂ V 1 j (R)| [0,1] , Ṽ 1,int j ⊂ Ṽ 1 j ⊂ Ṽ 1 j (R)| [0,1] . (iii) V 1
j and Ṽ 1 j are finite dimensional biorthogonal spaces spanned by biorthogonal bases

{ϕ 1 j,k : k ∈ ∆ j } and { φ1 j,k : k ∈ ∆ j }:< ϕ 1 j,k , φ1 j,k >= δ k,k , ∀ k, k ∈ ∆ j . (iv) V 1
j and Ṽ 1 j have respectively r and r polynomial exactness.

In point (ii), V 1 j (R)| [0,1] means the restriction of V 1 j (R)-functions to the interval [0, 1], whereas V 1,int j means interior functions of V 1 j (R)
, as introduced below (definitions 2.2,2.4), and same for the biorthogonal spaces. The dimension ∆ j ≈ 2 j of spaces V 1 j and Ṽ 1 j will be explicited later.

To construct such spaces (V 1 j , Ṽ 1 j ), as described in the numerous and now classical approaches [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF][START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF][START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF][START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF], we start with generators (ϕ 1 , φ1 ), that are biorthogonal scaling functions of a BMRA on R. We suppose that ϕ 1 is compactly supported on [n min , n max ] and reproduces polynomials up to degree r -1:

0 ≤ ≤ r -1, x ! = +∞ k=-∞ p1 (k) ϕ 1 (x -k) (9)
with p1 (k) = x ! , φ1 (x -k) . We suppose also that φ1 reproduces polynomials up to degree r -1:

0 ≤ ≤ r -1, x ! = +∞ k=-∞ p 1 (k) φ1 (x -k) (10)
with

p 1 (k) = x ! , ϕ 1 (x -k) .
To define a BMRA on [0, 1] we first define the set of interior scaling functions:

Definition 2.2
Let δ , δ ∈ N be two fixed parameters. For j ≥ 0, interior scaling functions of V 1 j are defined as scaling functions ϕ 1 j,k (x) = 2 j/2 ϕ 1 (2 j x -k) whose supports are included into

[ δ 2 j , 1 - δ 2 j ] ⊂ [0, 1]. If supp ϕ 1 = [n min , n max ], they correspond to indices k such that: δ -n min ≤ k ≤ 2 j -δ -n max
The space generated by interior scaling functions is then given by:

V 1,int j = span{ϕ 1 j,k ; k = k , 2 j -k } with k = δ -n min and k = δ + n max .
Similarly, let δ , δ ∈ N be two parameters. Interior scaling functions of Ṽ 1 j are defined as scaling functions φ1 j,k (x) = 2 j/2 φ1 (2 j x-k) whose supports are included into

[ δ 2 j , 1- δ 2 j ].
The space generated by interior scaling functions is then given by:

Ṽ 1,int j = span{ φ1 j,k ; k = k , 2 j -k } with k = δ -ñmin and k = δ + ñmax , if the support of φ1 is [ñ min , ñmax ].

Remark 1

The parameters (δ , δ , δ , δ ) are "free" parameters (chosen as small as possible), and chosen in practice to adjust the dimension of the spaces V 1 j and Ṽ 1 j . To preserve the polynomial reproduction (9, 10) on the interval [0, 1], we follow the approach of [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF][START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF] and define edge scaling functions at the edge 0:

Definition 2.3
The edge scaling functions at the edge 0 are defined by:

0 ≤ ≤ r -1, Φ 1, (x) = k -1 k=1-nmax p1 (k) ϕ 1 k (x) χ [0,+∞[
and the duals:

0 ≤ ≤ r -1, Φ1, (x) = k -1 k=1-ñmax p 1 (k) φ1 k (x) χ [0,+∞[
At the edge 1, the edge scaling functions Φ 1, j, are constructed on ] -∞, 1] by symmetry, using the transform T f (x) = f (1 -x).

As usual, one define the multiresolution spaces V 1 j on [0, 1], by the direct sum:

V 1 j = V 1, j ⊕ V 1,int j ⊕ V 1, j (11)
where:

V 1, j = span{2 j/2 Φ 1, (2 j x) ; = 0, • • • , r -1} V 1, j = span{2 j/2 Φ 1, (2 j (1 -x)) ; = 0, • • • , r -1}
are the edge spaces. In practice we have to choose j ≥ j min where j min is the smallest integer which verifies:

j min > log 2 [n max -n min + δ + δ ]
This condition ensures that the supports of edge scaling functions at 0 do not intersect the supports of edge scaling functions at 1. The polynomial reproduction in V 1 j is then satisfied since, for 0 ≤ ≤ r -1 and x ∈ [0, 1] we have:

2 j/2 (2 j x) ! = 2 j/2 Φ 1, (2 j x) + 2 j -k k=k p1 (k) ϕ 1 j,k (x) + 2 j/2 Φ 1, (2 j (1 -x)) (12)
Similarly, multiresolution spaces Ṽ 1 j are defined with the same structure, allowing the polynomial reproduction up to degree r -1:

Ṽ 1 j = span{ Φ1, j, } =0,r-1 ⊕ Ṽ 1,int j ⊕ span{ Φ1, j, } =0,r-1 (13)
In order to obtain the equality between dimensions of V 1 j and Ṽ 1 j , we have to adjust the parameters δ = k -ñmax and δ = k + ñmin in the definition 2.2 such that:

k -r = k -r and k -r = k -r (14)
We get:

dim(V 1 j ) = dim( Ṽ 1 j ) = 2 j -(δ + δ ) -(n max -n min ) + 2r + 1 where (δ , δ ) remain "free" parameters of the construction. Like for V 1 j we have to choose j ≥ jmin with jmin > log 2 [ñ max -ñmin + δ + δ ].
The last point of the construction lies in the biorthogonalization process of the new basis functions, since edge scaling functions of V 1 j and Ṽ 1 j are not biorthogonal. Several biorthogonalization methods exist [START_REF] Andersson | Wavelets on closed subsets of the real line[END_REF][START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF][START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF], here we apply on one hand the method proposed by Dahmen and al. [START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF] when using B-spline generators, and on the other hand a Gram-Schmidt process with Daubechies orthogonal generators [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF]. In both cases, it requires the inversion of the Gram matrix associated with the two systems, which for orthogonal and B-Spline generators is non singular [START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF]. Finally, the spaces (V 1 j , Ṽ 1 j ) j≥max{j min , jmin } , constitute a biorthogonal MRA of L 2 (0, 1) in the sense of definition 2.1.

Moreover, homogeneous boundary conditions can be simply imposed to

V 1 j , of the form f (λ) (α) = 0 at point α = 0 or 1, with λ = 0, • • • , r -1, by removing the edge scaling function Φ 1, λ if α = 0 or Φ 1, λ if α = 1 in the definition 2.
3 of edge spaces (see [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF] for more details). In such case, we also remove the function Φ1, λ or Φ1, λ from Ṽ 1 j prior to biorthogonalization.

2.2.2 Construction of (V 0 j , Ṽ 0 j ) on [0, 1] linked by differentiation /integration with (V 1 j , Ṽ 1 j )
We will now construct spaces (V 0 j , Ṽ 0 j ), related to the spaces (V 1 j , Ṽ 1 j ) of section 2.2.1 by the relations [START_REF] Cieselski | Spline bases in classical function spaces on compact manifolds, part II[END_REF][START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF] of differentiation/integration. Given (ϕ 1 , φ1 ) biorthogonal scaling functions with approximation orders (r, r) (r > 1), and compact supports [n min , n max ], [ñ min , ñmax ], we consider (ϕ 0 , φ0 ) arising from proposition 2.1. Then (ϕ 0 , φ0 ) satisfy some properties on R, that we recall below.

First, the scaling functions (ϕ 0 , φ0 ) are defined such that:

d dx ϕ 1 (x) = ϕ 0 (x) -ϕ 0 (x -1) and x+1 x φ1 (t)dt = φ0 (x) (15)
This implies that ϕ 0 has for compact support [n min , n max -1], and reproduces polynomials up to degree r -2:

0 ≤ ≤ r -2, x ! = +∞ k=-∞ p0 (k) ϕ 0 (x -k) (16)
with p0 (k) = x ! , φ0 (x -k) . Equation (8) implies: p0 (k) = p1 +1 (k) -p1 +1 (k -1) for = 0, • • • , r -2.
In the same way φ0 has for compact support [ñ min -1, ñmax ], and reproduces polynomials up to degree r:

0 ≤ ≤ r, x ! = +∞ k=-∞ p 0 (k) φ0 (x -k) (17)
with

p 0 (k) = x ! , ϕ 0 (x -k) . Equation (7) implies: p 1 (k) = p 0 +1 (k + 1) -p 0 +1 (k) for = 1, • • • , r.
Like for V 1 j , we first define the set of interior scaling functions of V 0 j : Definition 2.4 Let δ , δ ∈ N and k = δ -ñmin and k = δ + ñmax be the parameters introduced in definition 2.2. For j ≥ 0, the interior scaling functions of V 0 j are defined as scaling functions ϕ

0 j,k (x) = 2 j/2 ϕ 0 (2 j x -k) whose supports are included into [ δ 2 j , 1 - δ 2 j ] ⊂ [0, 1]. Since supp ϕ 0 = [n min , n max -1]
, the space generated by interior scaling functions is given by:

V 0,int j = span{ϕ 0 j,k ; k = k , 2 j -k + 1}
Similarly, let δ , δ ∈ N, and k = δ -ñmin , k = δ + ñmax be fixed by definition 2.2 and relation [START_REF] Deriaz | Direct Numerical Simulation of Turbulence using divergencefree wavelets[END_REF]. Interior scaling functions of Ṽ 0 j are defined as scaling functions

φ0 j,k (x) = 2 j/2 ϕ 0 (2 j x-k) whose supports are included into [ δ 2 j , 1- δ 2 j ].
The space generated by interior scaling functions is then given by: Ṽ 0,int [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF] on [0, 1] in (V 0 j , Ṽ 0 j ), we define edge scaling functions at the edge 0:

j = span{ φ0 j,k ; k = k + 1, 2 j -k } To preserve the polynomial reproduction
Definition 2.5
The edge scaling functions at the edge 0 are defined by:

0 ≤ ≤ r -2, Φ 0, (x) = k -1 k=2-nmax p0 (k) ϕ 0 k (x) χ [0,+∞[
and the biorthogonal ones, vanishing at 0:

1 ≤ ≤ r, Φ0, (x) = k k=1-ñmax p 0 (k) φ0 k (x) χ [0,+∞[
At the edge 1, the edge scaling functions Φ 0, j, and Φ0, j, are constructed by symmetry, using the transform

T f (x) = f (1 -x).
Remark 2 Following Jouini-Lemarié-Rieusset [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF], to preserve the commutation between the derivation and the multiscale projectors, the space Ṽ 0 j should satisfy ( 6): Ṽ 0 j ⊂ H 1 0 (0, 1) . Indeed we impose by construction homogeneous Dirichlet boundary conditions to Ṽ 0 j , since we do not consider Φ0, 0 and Φ0, 0 ( = 0) in definition 2.5. The multiresolution spaces V 0 j on [0, 1] are then defined by:

V 0 j = V 0, j ⊕ V 0,int j ⊕ V 0, j
where:

V 0, j = span{2 j/2 Φ 0, (2 j x) ; = 0, • • • , r -2} V 0, j = span{2 j/2 Φ 0, (2 j (1 -x)) ; = 0, • • • , r -2}
The polynomial reproduction up to degree r-2 in V 0 j is then ensured. In practice j > j min where the parameter j min is now adapted to both BMRA (V 1 j , Ṽ 1 j ) and (V 0 j , Ṽ 0 j ) by:

j min > max{log 2 [n max -n min + δ + δ + 1], log 2 [ñ max -ñmin + δ + δ + 1]}
Multiresolution spaces Ṽ 0 j are defined with similar structure, allowing a polynomial reproduction up to degree r, and satisfying vanishing boundary conditions at 0 and 1.

Ṽ

0 j = span{ Φ0 j, } =1,r ⊕ Ṽ 0,int j ⊕ span{ Φ1, j, } =1,r (18)
A simple calculation shows that:

dim(V 0 j ) = 2 j -k -k + 2r and dim( Ṽ 0 j ) = 2 j -k -k + 2r (19)
Since the parameters k , k , k and k are chosen to satisfy equation ( 14), we obtain:

dim(V 0 j ) = dim( Ṽ 0 j ) The following proposition proves that d dx V 1 j = V 0 j and d dx Ṽ 0 j ⊂ Ṽ 1 j . Proposition 2.2 (i)
The interior scaling functions of (V 1 j , V 0 j ) and ( Ṽ 1 j , Ṽ 0 j ) introduced in definitions 2.2, 2.4 satisfy:

d dx (ϕ 1 j,k ) = 2 j [ϕ 0 j,k -ϕ 0 j,k+1 ] d dx ( φ0 j,k ) = 2 j [ φ1 j,k-1 -φ1 j,k ] and dim(V 0,int j ) = dim(V 1,int j ) + 1, dim( Ṽ 0,int j ) = dim( Ṽ 1,int j ) -1. (ii) The edge scaling functions of (V 1 j , V 0 j ) introduced in definitions 2.3, 2.5 satisfy: for = 1, r -1, (Φ 1, 0 ) (x) = -ϕ 0 k , (Φ 1, ) (x) = Φ 0, -1 -p1 (k -1) ϕ 0 k (Φ 1, 0 ) (1 -x) = ϕ 0 2-k , (Φ 1, ) (1 -x) = -Φ 0, -1 (1 -x) + p1 (2 -k ) ϕ 0 2-k
whereas those of ( Ṽ 1 j , Ṽ 0 j ) are linked by: for = 1, r -1,

( Φ0, 0 ) = -φ1 k , ( Φ0, ) = Φ1, -1 -p 0 ( k ) φ1 k ( Φ0, 0 (1 -x)) = φ1 1-k , ( Φ0, (1 -x)) = -Φ1, -1 (1 -x) + p 0 (2 -k ) φ1 1-k Proof 2.1
The point (i) follows from [START_REF] Deriaz | Divergence-free and curl-free wavelets in 2D and 3D, application to turbulent flows[END_REF] and since interior functions are defined by ϕ 1 j,k (x) = 2 j/2 ϕ(2 j x -k) (and the same for ϕ 0 j,k , φ1 j,k , φ0 j,k ). The point (ii) comes directly from definitions 2.3, 2.5 of edge scaling functions. We focus on the first line of equalities: at edge 0 and for 0 ≤ ≤ r -1, the edge scaling functions of (V 1 j ) are defined by dilation of:

Φ 1, (x) = k -1 k=1-nmax p1 (k) ϕ 1 k (x)χ [0,+∞[ with ϕ 1 k (x) = ϕ 1 (x -k)
Differentiating in ]0, +∞[, one obtains for = 0 (p 1 0 (k) = 1, ∀k):

(Φ 1, 0 ) = k -1 k=1-nmax (ϕ 0 k -ϕ 0 k+1 ) χ ]0,+∞[ = ϕ 0 1-nmax χ ]0,+∞[ -ϕ 0 k = -ϕ 0 k since supp ϕ 0 1-nmax = [n min -n max + 1, 0].
In the same way, for = 1, r -1:

(Φ 1, ) = k -1 k=1-nmax p1 (k) (ϕ 0 k -ϕ 0 k+1 ) χ ]0,+∞[ = k -1 k=2-nmax [p 1 (k) -p1 (k -1)] ϕ 0 k χ [0,+∞[ -p 1 (k -1) ϕ 0 k From (8), since p0 -1 (k) = p1 (k) -p1 (k -1)
we get:

(Φ 1, ) = k -1 k=2-nmax p0 -1 (k) ϕ 0 k χ [0,+∞[ -p1 (k -1) ϕ 0 k = Φ 0, -1 -p1 (k -1) ϕ 0 k
This proves the relation between edge scaling functions at 0 of V 1 j and V 0 j . The proof for edge scaling functions at edge 1 and in the biorthogonal spaces Ṽ 1 j and Ṽ 0 j is obtained with similar arguments.

For easy reading, the two pairs of biorthogonal bases of (V 1 j , Ṽ 1 j ) and (V 0 j , Ṽ 0 j ) will be denoted by (ϕ 1 j,k , φ1 j,k ) and (ϕ 0 j,k , φ0 j,k ) respectively. The oblique projector on V 1 j parallel to ( Ṽ 1 j ) ⊥ will be denoted by P 1 j :

P 1 j : L 2 (0, 1) → V 1 j , f → P 1 j (f ) = k f, φ1 j,k ϕ 1 j,k (20)
while P1 j will denote its adjoint, and P 0 j , P0 j the biorthogonal projectors associated with (V 0 j , Ṽ 0 j ).

The following proposition proves that the constructed biorthogonal MRAs take place in the theoretical framework of Jouini-Lemarié-Rieusset in [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF].

Proposition 2.3

The two pairs of biorthogonal spaces (V 1 j , Ṽ 1 j ) and (V 0 j , Ṽ 0 j ) are related to:

d dx V 1 j = V 0 j and Ṽ 0 j = H 1 0 ∩ x 0 Ṽ 1 j , with Ṽ 0 j ⊂ H 1 0 (0, 1) Proof 2.2 The inclusions d dx V 1 j ⊂ V 0 j and d dx Ṽ 0 j ⊂ Ṽ 1 j
are straightforward according to proposition 2.2. Moreover the equality of dimensions between spaces ends the proof.

We then define the change of bases between the spaces

( d dx V 1 j , d dx Ṽ 0 j ) and (V 0 j , Ṽ 1 j ) as follows.
Definition 2.6 Let (L 1 j , L 0 j ) and ( L1 j , L0 j ) be the two pairs of sparse matrices defined by the change of bases between spaces involved in proposition 2.3 :

d dx ϕ 1 j,k = n (L 1 j ) k,n ϕ 0 j,n , d dx φ0 j,k = n ( L0 j ) k,n φ1 j,n (21)
and -

x 0 ϕ 0 j,k = m (L 0 j ) k,m ϕ 1 j,m , - x 0 φ1 j,k = m ( L1 j ) k,m φ0 j,m (22)
Remark 3

The matrices (L 1 j , L 0 j ) and ( L1 j , L0 j ) are rectangular and from the biorthogonality of spaces (V 0 j , Ṽ 0 j ), it comes:

L 0 j ( L0 j ) T = I dim(V 0 j )
where I denotes the matrix identity. Except for the first scaling function Φ 1, j,0 , which does not satisfy homogeneous Dirichlet boundary condition at 0, similarly we have:

L 1 j ( L1 j ) T = I dim(V 1 j )-1
In addition, the definition of L1 j [START_REF] Souleymane | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF] must include the scaling functions Φ0, j,0 and Φ0, j,0 , if not, this definition leads to 1 0 φ1 j,m = 0, which is not true. With this definition, we prove the commutation between multiscale projectors and differentiation.

Proposition 2.4

Let (P 1 j , P1 j ) and (P 0 j , P0 j ) be the multiscale projectors defined by (20):

(i) ∀ f ∈ H 1 (0, 1), d dx • P 1 j f = P 0 j • d dx f . (ii) ∀ f ∈ H 1 0 (0, 1), d dx • P0 j f = P1 j • d dx f . Proof 2.3 (i) The relation of commutation d dx • P 1 j f = P 0 j • d dx f
was demonstrated first in [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF] in the general setting. We will now prove the commutation of the projectors with derivation (ii). Let (ϕ 0 j,k , φ0 j,k ) be a pair of biorthogonal scaling functions of (V 0 j , Ṽ 0 j ) and as φ0 j,k ∈ H 1 0 (0, 1), from proposition 2.2 we obtain:

- d dx φ0 j,k , x 0 ϕ 0 j,k = φ0 j,k , ϕ 0 j,k = δ k,k = n ( L0 j ) k,n (L 0 j ) k ,n = [ L0 j L 0T j ] k,k
where L 0 j and L0 j are introduced in (21) and [START_REF] Souleymane | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF].

For f ∈ H 1 0 (0, 1), since f, ϕ 0 j,k = -d dx f,
x 0 ϕ 0 j,k we get:

d dx P0 j (f ) = k f, ϕ 0 j,k d dx φ0 j,k = n m k (L 0 j ) k,m ( L0 j ) k,n d dx f, ϕ 1 j,m φ1 j,n = n m δ n-m d dx f, ϕ 1 j,m φ1 j,n = n d dx f, ϕ 1 j,n φ1 j,n = P1 j ( d dx f )
This proves the relation (ii).

Wavelet spaces

We begin with the construction of wavelet bases of the biorthogonal MRA (V 1 j , Ṽ 1 j ). This point is classical, although different kinds of wavelets may be designed [START_REF] Andersson | Wavelets on closed subsets of the real line[END_REF][START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF][START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF][START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF][START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF]. For j ≥ j min , the biorthogonal wavelet spaces associated to V 1 j and Ṽ 1 j are in all cases defined by:

W 1 j = V 1 j+1 ∩ ( Ṽ 1 j ) ⊥ and W 1 j = Ṽ 1 j+1 ∩ (V 1 j ) ⊥
The wavelet space W 1 j has the following structure:

W 1 j = W 1, j ⊕ W 1,int j ⊕ W 1, j where    W 1, j = span{2 j/2 Ψ 1, (2 j x) ; = 0, p -1} W 1,int j = span{ψ 1 j,k ; k = p , 2 j -p -1} W 1, j = span{2 j/2 Ψ 1, (2 j (1 -x)) ; = 0, p -1} (23)
p and p introduced above are suitable integers to ensure that the support of each interior wavelet

ψ 1 j,k (x) = 2 j/2 ψ 1 (2 j x -k) of W 1,int j is included into [ δ 2 j , 1 - δ 2 J ]. Recall that the support of ψ 1 (wavelet on R) is [ n min -ñmax+1 2 , nmax-ñ min +1 2
], then we deduce:

p = ñmax + k -1 2 and p = k -ñmin + 1 2
To construct the edge wavelets Ψ 1, j, and Ψ 1, j, we have followed the work of Grivet-Talocia and Tabacco [START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF]. The biorthogonal spaces W 1 j have the same structure and the wavelet bases of the two spaces must to be biorthogonalized identically as the scaling functions. The resulting wavelet bases are denoted by {ψ 1 j,k } k=0,2 j -1 and { ψ1 j,k } k=0,2 j -1 without distinction.

Another advantage of this construction is the existence of a fast wavelet transforms because the scaling functions and wavelets satisfy both a two-scale relation. Indeed, there are sparse matrices H 1 j , H1 j , G 1 j and G1 j such that:

ϕ 1 j,k = n (H 1 j ) k,n ϕ 1 j+1,n and ψ 1 j,k = n (G 1 j ) k,n ϕ 1 j+1,n φ1 j,k = n ( H1 j ) k,n φ1 j+1,n and ψ1 j,k = n ( G1 j ) k,n φ1 j+1,n
The main objective of this section is now to construct biorthogonal wavelet bases of W 0 j and W 0 j , that will be linked to ψ 1 j,k and ψ1 j,k by differentiation/integration. A first result in this direction is given by the following proposition:

Proposition 2.5 Let (V 1 j , Ṽ 1 j
) and (V 0 j , Ṽ 0 j ) be two biorthogonal multiresolution analyses satisfying proposition 2.3. The spaces W 0 j and W 0 j defined by:

W 0 j = d dx W 1 j and W 0 j = x 0 W 1 j (24)
correspond to biorthogonal wavelet spaces associated to (V 0 j , Ṽ 0 j ) in the classical sense:

W 0 j = V 0 j+1 ∩ ( Ṽ 0 j ) ⊥ and W 0 j = Ṽ 0 j+1 ∩ (V 0 j ) ⊥
Proof 2.4 First we prove the relation:

d dx W 1 j = V 0 j+1 ∩ ( Ṽ 0 j ) ⊥ . Let w 1 j ∈ W 1 j
, we get from proposition 2.4 (i):

P 0 j ( d dx w 1 j ) = d dx P 1 j (w 1 j ) = 0 thus d dx w 1 j ∈ ( Ṽ 0 j ) ⊥ . Moreover, following the proposition 2.3 (i): d dx w 1 j ∈ d dx V 1 j+1 = V 0 j+1 This implies d dx W 1 j ⊂ V 0 j+1 and thus d dx W 1 j ∈ V 0 j+1 ∩ ( Ṽ 0 j ) ⊥ .
The equality is obtained by the equality between both space dimension. To prove the relation between W 0 j and W 1 j , as W 0 j ⊂ H 1 0 (0, 1), one can remark that:

W 0 j = x 0 W 1 j ⇔ d dx W 0 j = W 1 j
The relation d dx W 0 j = W 1 j is then proved as before.

As suggested in [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF], the wavelet bases of W 0 j and W 0 j are constructed directly by differentiating and integrating the wavelets of W 1 j and W 1 j :

Definition 2.7 Let {ψ 1 j,k } k=0,2 j -1 and { ψ1 j,k } k=0,2 j -1 be two biorthogonal wavelet bases of W 1 j and W 1 j respectively. The wavelets of W 0 j and W 0 j are defined by:

ψ 0 j,k = 2 -j (ψ 1 j,k ) and ψ0 j,k = -2 j x 0 ψ1 j,k
Interior wavelets ψ 0 j,k (x) = 2 j/2 ψ 0 (2 j x -k) in this definition correspond to the classical wavelets arising from previous constructions [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF][START_REF] Dahmen | Biorthogonal Spline-wavelets on the intervalstability and moment conditions[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF][START_REF] Grivet-Talocia | Wavelet on the interval with optimal localization[END_REF][START_REF] Chiavassa | On the Effective Construction of Compactly Supported Wavelets Satisfying Homogeneous Boundary Conditions on the Interval[END_REF], ψ 0 being a wavelet on R associated to the derivative ϕ 0 . On the other hand, in standard constructions, the edge wavelets do not verify the relations:

d dx Ψ 1, j,k = 2 j Ψ 0, j,k or Ψ0, j,k = -2 j x 0 Ψ1, j,k (25)
The following proposition guarantees [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF] and the new edge wavelets provided by definition 2.7 preserve fast algorithms since they satisfy a two-scale relationship.

Proposition 2.6

Let {ψ 1 j,k } k=0,2 j -1 and { ψ1 j,k } k=0,2 j -1 be two biorthogonal wavelet bases of W 1 j and W 1 j associated respectively to filters G 1 j and G1 j :

ψ 1 j,k = n (G 1 j ) k,n ϕ 1 j+1,n and ψ1 j,k = n ( G1 j ) k,n φ1 j+1,n
Then the following propositions hold: (i) The system {ψ 0 j,k = 2 -j (ψ 1 j,k ) } k=0,2 j -1 and { ψ0 j,k = -2 j x 0 ψ1 j,k } k=0,2 j -1 form biorthogonal wavelet bases of W 0 j and W 0 j respectively. (ii) There exist sparse matrices G 0 j and G0 j defined by:

G 0 j = 2 -j G 1 j L 1 j+1 and G0 j = -2 j G1 j L 0T j+1 ( 26
)
and the wavelets ψ 0 j,k and ψ0 j,k satisfy:

ψ 0 j,k = n (G 0 j ) k,n ϕ 0 j+1,n and ψ0 j,k = n ( G0 j ) k,n φ0
j+1,n Proof 2.5 (i) The construction of biorthogonal bases by differentiation /integration is an idea of Cieselski-Figiel [START_REF] Cieselski | Spline bases in classical function spaces on compact manifolds, part II[END_REF] used by Jouini-Lemarié-Rieusset [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF]. Moreover, from proposition 2.5 we can see that the two systems are independent and they generate the corresponding spaces.

(ii) To derive the filters, we only use the definition of wavelets. So, we have:

2 j ψ 0 j,k = n (G 1 j ) k,n (ϕ 1 j+1,n ) = n,m (G 1 j ) k,n (L 1 j+1 ) n,m ϕ 0 j+1,m = m [G 1 j L 1 j+1 ] k,m ϕ 0 j+1,m = 2 j m (G 0 j ) k,m ϕ 0 j+1,m
which gives the filter of ψ 0 j,k . Similarly, for ψ0 j,k we obtain:

d dx ψ0 j,k = m ( G0 j ) k,m ( φ0 j+1,m ) = m,n ( G0 j ) k,m ( L0 j+1 ) m,n φ1 j+1,n = n [ G0 j L0 j+1 ] k,n φ1 j+1,n = -2 j n ( G1 j ) k,n φ1 j+1,n
Taking into account the relation [ L0 j L 0T j ] k,k = δ k,k , there are summary:

G 0 j = 2 -j G 1 j L 1 j+1 and G0 j = -2 j G1 j L 0T j+1
This completes the proof.

Remark 4

The above construction of wavelets ψ 0 j,k and ψ0 j,k has two main interests: their filters are directly accessible from those of ψ 1 j,k and ψ1 j,k and there is no need for biorthogonalization as for classical constructions.

Example 1

To illustrate, we give the plot of edge scaling functions and wavelets at 0 in (V 1 j , Ṽ 1 j ) on Figure 1. The generators (ϕ 1 , φ1 ) used are biorthogonal B-Spline with r = r = 3. Then we have: n min = -1, n max = 2, ñmin = -3 and ñmax = 4. The "free" integer parameters are chosen as δ min = δ max = 2 and δmin = δmax = 0. On Figure 2, we plot the corresponding edge scaling functions and wavelets of (V 0 j , Ṽ 0 j ) at 0. The Figure 3 and Figure 4 show the non zeros elements of filters G 0 j and G0 j and matrices L 0 j and L 1 j respectively, for j = 6.

Biorthogonal MRA of H div (Ω)

Let Ω be the square [0, 1] 2 . The aim of the present section is to provide a divergence-free MRA and wavelet bases of the space H div (Ω) [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]:

H div (Ω) = {u ∈ (L 2 (Ω)) 2 : div(u) = 0 and u • ν = 0}
Since this space is equal to:

H div (Ω) = {u = curl χ ; χ ∈ H 1 0
(Ω)} our construction consists in taking the curl of a regular MRA of the two-dimensional scalar space H 1 0 (Ω).

Such MRA of H 1 0 (Ω) is usually defined as tensor-product of one-dimensional MRA of H 1 0 (0, 1). We now consider a regular one-dimensional MRA satisfying homogeneous boundary conditions:

V D j = V 1 j ∩ H 1 0 (0, 1
) as constructed in section 2.2, and which takes the form:

V D j = span{Φ 1, j, ; = 1, r -1} ⊕ V 1,int j ⊕ span{Φ 1, j, ; = 1, r -1}
To simplify, we denote by ϕ D j,k the scaling functions of V D j :

V D j = span{ϕ D j,k ; k = 0, 2 j -k -k + 2r -2}
and ψ D j,k the corresponding wavelets. With these notations, the divergence-free scaling function spaces are defined below.

Definition 2.8

For j ≥ j min , the divergence-free scaling function spaces V div j are defined by:

V div j = curl(V D j ⊗ V D j ) = span{Φ div j,k } (27)
where the divergence-free scaling functions are given by:

Φ div j,k := curl[ϕ D j,k 1 ⊗ ϕ D j,k 2 ], j ≥ j min (28)
The spaces V div j defined above constitute an increasing sequence of subspaces of (L 2 (Ω)) 2 :

V div j ⊂ V div j+1 of dimension: dim(V div j ) = dim(V D j ) 2 = (2 j -k -k + 2r -1) 2 = 2 j -(n max -n min ) -(δ + δ ) + 2r -1 2
We will also consider a more standard multiresolution analysis V j of (L 2 (Ω)) 2 defined as:

V j = (V 1 j ⊗ V 0 j ) × (V 0 j ⊗ V 1 j ) ( 29
)
V 0 j being the spaces defined in section 2.2.2. By proposition 2.3, this discrete space V j preserves the divergence-free condition, as stated by Jouini-Lemarié-Rieusset [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF]:

u ∈ (L 2 (Ω)) 2 , div(u) = 0 ⇒ div[ P j (u)] = 0 ( 30
)
where P j is the biorthogonal projector on V j :

P j = (P 1 j ⊗ P 0 j , P 0 j ⊗ P 1 j ) (31)
In the same way, we now introduce anisotropic divergence-free wavelets and wavelet spaces:

Definition 2.9

The anisotropic divergence-free wavelets and wavelet spaces are given by:

Ψ div,1 j,k := curl[ϕ D j min ,k 1 ⊗ ψ D j 2 ,k 2 ] and W div,1 j = span{Ψ div,1 j,k }, j 2 ≥ j min Ψ div,2 j,k := curl[ψ D j 1 ,k 1 ⊗ ϕ D j min ,k 2 ] and W div,2 j = span{Ψ div,2 j,k }, j 1 ≥ j min Ψ div,3 j,k := curl[ψ D j 1 ,k 1 ⊗ ψ D j 2 ,k 2 ] and W div,3 j = span{Ψ div,3 j,k }, j 1 , j 2 ≥ j min
The following proposition proves that (V div j ) j≥j min is a multiresolution analysis of H div (Ω). Proposition 2.7

The divergence-free scaling functions spaces V div j and wavelet spaces W div,ε j for ε = 1, 2, 3, satisfy:

(i) V div j min ⊂ • • • ⊂ V div j ⊂ V div j+1 ⊂ • • • ⊂ H div (Ω) and ∪V div j = H div (Ω). (ii) V div j = V div j min j min ≤j 1 ,j 2 ≤j-1 (⊕ ε=1,2,3 W div,ε j ).
(iii) For all j and ε = 1, 2, 3, {Ψ div,ε j,k } is a Riesz basis of W div,ε j .

Then each vector function u of H div (Ω) has an unique decomposition into the basis {Φ div j min ,k , Ψ div,ε j,k } j 1 ,j 2 ≥j min ;ε=1,2,3 :

u = k c div j min ,k Φ div j min ,k + j,k ε=1,2,3 d div,ε j,k Ψ div,ε j,k
with the norm-equivalence:

u L 2 ∼ k |c div j min ,k | 2 + j,k ε=1,2,3 |d div,ε j,k | 2
Proof 2.6 (i) Let V j be the spaces defined in (29). Since the spaces H div (Ω)∩ V j provide a multiresolution analysis of H div (Ω) [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF], point (i) is reduced to prove that:

V div j = H div (Ω) ∩ V j .
According to proposition 2.2, we have V div j ⊂ V j and V div j ⊂ H div (Ω) by construction. Conversely, let u ∈ H div (Ω) ∩ V j , and P j be the biorthogonal projector on V j defined in (31). We are going to prove that u ∈ V div j . On one side, as u ∈ V j we have u = P j (u), on the other hand due to u ∈ H div (Ω) we have u = curl (χ) with χ ∈ H 1 0 (Ω), and thus:

u = P j [curl (χ)] Since the spaces (V D j ⊗ V D j ) j≥j min form a MRA of H 1 0 (Ω)
, we can decompose χ as:

χ = P D j (χ) + j 1 ,j 2 ≥j Q D 1 (χ) + Q D 2 (χ) + Q D 3 (χ)
where

P D j (χ) = k c k ϕ D j,k 1 ⊗ ϕ D j,k 2 , Q D 2 (χ) = j 1 ≥j,k d 2 j 1 ,k ψ D j 1 ,k 1 ⊗ ϕ D j,k 2 Q D 1 (χ) = j 2 ≥j,k d 1 j 2 ,k ϕ D j,k 1 ⊗ ψ D j 2 ,k 2 , Q D 3 (χ) = j 1 ,j 2 ≥j,k d 3 j,k ψ D j 1 ,k 1 ⊗ ψ D j 2 ,k 2
are the biorthogonal projectors on respectively

V D j ⊗ V D j , W D j 1 ⊗ V D j , V D j ⊗ W D j 2 and W D j 1 ⊗ W D j 2 . Proposition 2.2 implies that: curl [ϕ D j,k 1 ⊗ ψ D j 2 ,k 2 ] ∈ (V D j ⊗ W 0 j 2 ) × (V 0 j ⊗ W D j 2 ) hence: P j (curl [ϕ D j,k 1 ⊗ ψ D j 2 ,k 2 ]) = 0 and same for P j (curl [ψ D j 1 ,k 1 ⊗ ϕ D j,k 2 ]) and P j (curl [ψ D j 1 ,k 1 ⊗ ψ D j 2 ,k 2 ]
). This leads to:

P j (curl (χ)) = P j (curl [P D j (χ)]) = curl [P D j (χ)] By construction we have curl [P D j (χ)] ∈ V div j , which implies u ∈ V div j
and then completes the proof:

V div j = H div (Ω) ∩ V j .
(ii) The spaces V D j are a multiresolution analysis of H 1 0 (0, 1), and we can write:

V D j ⊗ V D j = (V D j min j-1 j 1 =j min W D j 1 ) ⊗ (V D j min j-1 j 2 =j min W D j 2 )
By definition of V div j , we obtain:

V div j = curl   (V D j min ⊗ V D j min ) j min ≤j 1 ,j 2 ≤j-1 [(V D j min ⊗ W D j 2 ) ⊕ (W D j 1 ⊗ V D j min ) ⊕ (W D j 1 ⊗ W D j 2 )]   which is exactly V div j = V div j min j min ≤j 1 ,j 2 ≤j-1 ⊕ ε=1,2,3 W div,ε j . (iii)
Following [START_REF] Lemarié-Rieusset | Analyses multi-résolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF][START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF], this point is a consequence of the proposition 2.8 below.

We now introduce the biorthogonal divergence-free scaling functions and wavelets. Let:

Φdiv j,k := φD j,k 1 ⊗ γj,k 2 -γ j,k 1 ⊗ φD j,k 2 , Ψdiv,1 j,k := 2 j 2 φD j min ,k 1 ⊗ ψ0 j 2 ,k 2 -γ j min ,k 1 ⊗ ψD j 2 ,k 2 (32) Ψdiv,2 j,k :
= ψD j 1 ,k 1 ⊗ γj min ,k 2 -2 j 1 ψ0 j 1 ,k 1 ⊗ φD j min ,k 2 , Ψdiv,3 j,k := 2 j 2 ψD j 1 ,k 1 ⊗ ψ0 j 2 ,k 2 -2 j 1 ψ0 j 1 ,k 1 ⊗ ψD j 2 ,k 2 (33)
where: γj,k = -

x 0 φD j,k . Remark that γj,k (0) = γj,k (1) since φD j,k ∈ H 1 0 (0, 1). Proposition 2.8
For a fixed j ≥ j min , the normalized families

1 √ 2 Φ div j,k , 1 √ 4 j 2 +1 Ψ div,1 j,k , 1 √ 4 j 1 +1 Ψ div,2 j,k , 1 √ 4 j 1 +4 j 2 Ψ div,3 j,k ; j1, j2 ≥ j , k and 1 √ 2 Φdiv j,k , 1 √ 4 j 2 +1
Ψdiv,1 j,k ,

1 √ 4 j 1 +1 Ψdiv,2 j,k , 1 √ 4 j 1 +4 j 2 Ψdi j,k
are biorthogonal in (L 2 (Ω)) 2 , then they form Riesz sequences of (L 2 (Ω)) 2 .

Remark 5

Contrarily the usual definitions of MRAs, the (L 2 -normalized) divergence-free scaling functions (Φ div j,k ) k don't form a Riesz basis of the space V div j , since they do not verify:

k c div j,k Φ div j,k 2 (L 2 (Ω)) 2 ∼ k |c div j,k | 2
for all (c div j,k ) ∈ 2 and independently of j. A counterexample is given by:

u = k ≤k 1 ,k 2 ≤2 j -k Φ div j,k 1 ,k 2 ,
which satisfy:

k |c div j,k | 2 = (2 j -k -k + 1) 2 ∼ 2 2j
. On the other hand:

u = k ≤k 1 ,k 2 ≤2 j -k 2 -j √ 2 ϕ D j,k 1 ⊗ (ϕ D j,k 2 ) -(ϕ D j,k 1 ) ⊗ ϕ D j,k 2 k ≤k 2 ≤2 j -k (ϕ D j,k 2 ) = k ≤k 2 ≤2 j -k 2 j (ϕ 0 j,k 2 -ϕ 0 j,k 2 +1 ) = 2 j (ϕ 0 j,k -ϕ 0 j,2 j -k +1 ) Let h(x) = k ≤k≤2 j -k ϕ D j,k (x) ∼ 2 j/2 χ [δ ,δ ] , ∀x ∈]0, 1[ and u 2 (L 2 (Ω)) 2 = 2 -2j 1 0 h 2 1 0 2 2j (ϕ 0 j,k -ϕ 0 j,2 j -k +1 ) 2 ∼ 2 j
2.4 Biorthogonal MRA of H ⊥ div (Ω) In this section, the curl-free function space H curl (Ω) that we consider is the following:

H curl (Ω) = {u = ∇ q : q ∈ H 1 0 (Ω)}
This space is a proper subspace of H ⊥ div (Ω):

H ⊥ div (Ω) = H curl (Ω) ⊕ H ∆ (Ω), with H ∆ (Ω) = {∇ q : q ∈ H 1 (Ω) and ∆ q = 0}
To construct a multiresolution analysis of H curl (Ω), it suffices to consider the gradient of a MRA of H 1 0 (Ω). With the same notations as in previous section 2.3, curl-free scaling functions spaces are defined by: Definition 2.10 For j ≥ j min , a curl-free scaling function space V ∇ j is defined by:

V ∇ j = ∇(V D j ⊗ V D j ) = span{Φ ∇ j,k } (34)
where the curl-free scaling functions are given by:

Φ ∇ j,k := ∇[ϕ D j,k 1 ⊗ ϕ D j,k 2 ], j ≥ j min (35)
The spaces (V ∇ j ) constitute an increasing sequence of subspaces of (L 2 (Ω)) 2 , of dimension:

dim(V ∇ j ) = dim(V D j ⊗ V D j ) = (2 j -k -k + 2r -1) 2 . Let V *
j be the standard multiresolution analysis of (L 2 (Ω)) 2 defined by:

V * j = (V 0 j ⊗ V 1 j ) × (V 1 j ⊗ V 0 j ) (36)
By proposition 2.3, the spaces V ∇ j are contained in V * j . We now define the corresponding irrotational wavelets.

Definition 2.11

The anisotropic curl-free wavelets and wavelet spaces are defined by:

Ψ ∇,1 j,k := ∇[ϕ D j min ,k 1 ⊗ ψ D j 2 ,k 2 ] and W ∇,1 j = span{Ψ ∇,1 j,k }, j 2 ≥ j min Ψ ∇,2 j,k := ∇[ψ D j 1 ,k 1 ⊗ ϕ D j min ,k 2 ] and W ∇,2 j = span{Ψ ∇,2 j,k }, j 1 ≥ j min Ψ ∇,3 j,k := ∇[ψ D j 1 ,k 1 ⊗ ψ D j 2 ,k 2 ] and W ∇,3 j = span{Ψ ∇,3 j,k }, j 1 , j 2 ≥ j min
The following proposition holds:

Proposition 2.9

The spaces V ∇ j and W ∇,ε j for = 1, 2, 3 verify:

(i) V ∇ j min ⊂ • • • ⊂ V ∇ j ⊂ V ∇ j+1 ⊂ • • • ⊂ H curl (Ω) and ∪V ∇ j = H curl (Ω) (ii) V ∇ j = V ∇ j min j min ≤j 1 ,j 2 ≤j-1 (⊕ ε=1,2,3 W ∇,ε j )
(iii) For all j and ε = 1, 2, 3, {Ψ ∇,ε j,k } is a Riesz basis of W ∇,ε j .

The proof uses the same arguments as for proposition 2.7. In addition, these construction of curl-free scaling functions and wavelets can readily be extended to higher dimensions, for more details see [START_REF] Souleymane | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF][START_REF] Souleymane | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets[END_REF].

Example 2

3 Fast divergence-free wavelet transform

We describe in this section the practical computation of divergence-free scaling function and wavelet coefficients of a vector field u ∈ V div j . We use the same notations as in previous sections. The starting point is the decomposition of u in the MRA of (L 2 (Ω)) 2 provided by V j :

V j = (V 1 j ⊗ V 0 j ) × (V 0 j ⊗ V 1 j )
On the scaling functions basis of V j , the vector field u = (u 1 , u 2 ) can be written as:

u 1 = k c 1 j,k ϕ 1 j,k 1 ⊗ ϕ 0 j,k 2 and u 2 = k c 2 j,k ϕ 0 j,k 1 ⊗ ϕ 1 j,k 2 (37)
The computation of divergence-free coefficients will use a change of bases between ( d dx ϕ 1 j,k ) and (ϕ 0 j,k ). This needs to construct the matrices L 0 j and L 1 j introduced in ( 21) and ( 22), which we recall the definition:

d dx ϕ 1 j,k = dim(V 0 j ) n=1 (L 1 j ) k,n ϕ 0 j,n
and -

x 0 ϕ 0 j,n = dim(V 1 j ) k=1 (L 0 j ) n,k ϕ 1 j,k
So we first give a more precise result on the computation of the elements of these matrices.

Proposition 3.1

Let k * j denoted the dimension of V 1 j . The only nonzero elements of matrices L 0 j and L 1 j correspond to: (i) For edge scaling functions and for 2 ≤ k ≤ r, we have:

(L 1 j ) 1,r = -1, (L 1 j ) k,k-1 = 1, (L 1 j ) k,r = -p 1 k-1 (k -1) (L 1 j ) k * j ,k * j -r = 1, (L 1 j ) k * j -k+1,k * j -k+1 = -1, (L 1 j ) k * j -k+1,k * j -r = p1 k-1 (2 j -k + 1) and
(L 0 j ) r,1 = 1, (L 0 j ) k-1,k = -1, (L 0 j ) k-1,1 = -p 1 k-1 (k -1) (L 0 j ) k * j -r,k * j = -1, (L 0 j ) k * j -k+1,k * j -k+1 = 1, (L 0 j ) k * j -k+1,k * j = -p 1 k-1 (2 j -k + 1)
(ii) For interior scaling functions and for r + 1 ≤ k ≤ k * j -r, we have:

(L 1 j ) k,k-1 = 1, (L 1 j ) k,k = -1 and
(L 0 j ) k,m = -1, (L 0 j ) k,k * j = -1, k ≤ m ≤ k * j -r Proof 3.1
To obtain the relation on L 1 j , it suffices to use proposition 2.2 written for j ≥ j min . To obtain the relation on L 0 j , as it is still true by differentiation we get:

-Φ 0, j, = r-1 k=0 (L 0 j ) ,k (Φ 1, j,k ) and -Φ 0, j, = r-1 k=0 (L 0 j ) k * j + ,k * j +k (Φ 1, j,k )
Then using the proposition 2.2 again, we have:

Φ 0, j, = 2 -j [(Φ 1, +1 ) + p1 +1 (k -1) (Φ 1, j,0 ) ] and Φ 0, j, (1 -.) = 2 -j [(Φ 1, +1 (1 -.)) -p1 +1 (2 j -k + 1) (Φ 1, j,0 (1 -.)) ] because: ϕ 0 j,k = -2 -j (Φ 1, j,0
) and ϕ 0 j,2 j -k +1 = 2 -j (Φ 1, j,0 (1 -.)) Similarly for interior scaling functions, by proposition 2.2 we get:

ϕ 0 j,k = 2 -j (ϕ 1 j,k ) + ϕ 0 j,k+1
and recursively for k + 1 ≤ k ≤ 2 j -k we deduce that ϕ 0 j,k satisfy:

ϕ 0 j,k = 2 -j (ϕ 1 j,k ) + ϕ 0 j,k+1 = 2 -j [(ϕ 1 j,k ) + • • • + (ϕ 1 j,2 j -kmax ) ] + ϕ 0 j,2 j -kmax+1 = 2 -j [(ϕ 1 j,k ) + • • • + (ϕ 1 j,2 j -kmax ) + (Φ 1, j,0 (1 -.
)) ] This completes the proof. Now, using matrices L 0 j and L 1 j we can rewrite the components of u as follows:

u 1 = k c 1 j,k ϕ 1 j,k 1 ⊗ ϕ 0 j,k 2 = - k [(c 1 j,k )L 0 j ] k ϕ 1 j,k 1 ⊗ (ϕ 1 j,k 2 ) (38) and u 2 = k c 2 j,k ϕ 0 j,k 1 ⊗ ϕ 1 j,k 2 = - k [L 0 T j (c 2 j,k )] k (ϕ 1 j,k 1 ) ⊗ ϕ 1 j,k 2 (39)
If u ∈ V div j , it can be uniquely written as:

u = k c div k Φ div j,k (40)
Therefore, we have the following proposition.

Proposition 3.2

The matrices of coefficients [c 1 j,k ] and [c 2 j,k ] are linked to the matrix of coefficients [c div j,k ] by:

2 1/2 [c div j,k ] = L 0 T j [c 2 j,k ] -[c 1 j,k ]L 0 j (41)
and conversely:

[c 1 j,k ] = 2 -1/2 [c div j,k ]L 1 j and [c 1 j,k ] = -2 -1/2 L 1 T j [c div j,k ] (42) Proof 3.2
We assume that the basis of divergence-free scaling functions used is that of the lemma 2.8. With the help of (38) and (39), by computing directly the inner product we get:

u/ Φdiv j,k = 2 -1/2 L 0 T j [c 2 j,k ] -[c 1 j,k ]L 0 j k
Which proves (41). The second relation (42) is nothing but the change of basis described by the previous definition of L 1 j .

Now the objective is to compute the divergence-free wavelet coefficients of u:

u = k c div j,k Φ div j,k + j 2 ≥j,k d div,1 j,k Ψ div,1 j,k + j 1 ≥j,k d div,2 j,k Ψ div,2 j,k + j 1 ,j 2 ≥j,k d div,3 j,k Ψ div,3 j,k
We start with the standard wavelet decomposition of u = (u 1 , u 2 ) in V j :

u 1 = k c 1 j,k ϕ 1 j,k 1 ⊗ ϕ 0 j,k 2 + j 2 ≥j,k d 1,1 j 2 ,k ϕ 1 j,k 1 ⊗ ψ 0 j 2 ,k 2 + j 1 ≥j,k d 1,2 j 1 ,k ψ 1 j 1 ,k 1 ⊗ ϕ 0 j,k 2 + j 1 ,j 2 ≥j,k d 1,3 j,k ψ 1 j 1 ,k 1 ⊗ ψ 0 j 2 ,k 2 and u 2 = k c 2 j,k ϕ 0 j,k 1 ⊗ ϕ 1 j,k 2 + j 2 ≥j,k d 2,1 j 2 ,k ϕ 0 j,k 1 ⊗ ψ 1 j 2 ,k 2 + j 1 ≥j,k d 2,2 j 1 ,k ψ 0 j 1 ,k 1 ⊗ ϕ 1 j,k 2 + j 1 ,j 2 ≥j,k d 2,3 j,k ψ 0 j 1 ,k 1 ⊗ ψ 1 j 2 ,k 2
Which rewrites, using the matrices L 0 j and L 1 j :

u 1 = - k [(c 1 j,k )L 0 j ] k ϕ 1 j,k 1 ⊗ (ϕ 1 j,k 2 ) + j 2 ≥j,k d 1,1 j 2 ,k ϕ 1 j,k 1 ⊗ ψ 0 j 2 ,k 2 - j 1 ≥j,k [(d 1,2 j 1 ,k )L 0 j ] k ψ 1 j 1 ,k 1 ⊗ (ϕ 1 j,k 2 ) + j 1 ,j 2 ≥j,k d 1,3 j,k ψ 1 j 1 ,k 1 ⊗ ψ 0 j 2 ,k 2 and u 2 = - k [L 0 T j (c 2 j,k )] k (ϕ 1 j,k 1 ) ⊗ ϕ 1 j,k 2 - j 2 ≥j,k [L 0 T j (d 2,1 j 2 ,k )] k (ϕ 1 j,k 1 ) ⊗ ψ 1 j 2 ,k 2 + j 1 ≥j,k d 2,2 j 1 ,k ψ 0 j 1 ,k 1 ⊗ ϕ 1 j,k 2 + j 1 ,j 2 ≥j,k d 2,3 j,k ψ 0 j 1 ,k 1 ⊗ ψ 1 j 2 ,k 2
We now prove the following proposition.

analyses on the interval allowing the reproduction of polynomials, respects the theoretical framework of the previous work of Jouini and Lemarié-Rieusset [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF]. Moreover our construction can incorporate homogeneous boundary conditions in the basis functions, which allows the representation of more physical divergence-free vector functions. This ability is not present, for instance, in the attempt addressed by Stevenson [START_REF] Stevenson | Divergence-free wavelet bases on the hypercube[END_REF]. Associated fast wavelet transforms have been implemented satisfactory, opening new prospects for the realistic simulation of incompressible flows. First attempts have successfully been presented in this article with the Helmholtz decomposition of a vector flow, or with the computation of a Stokes problem solution. Work on more complex problems are underway, such as the direct simulation of turbulence, and this will the subject of a forthcoming paper.

Figure 1 :

 1 Figure 1: Scaling functions Φ 1, (three first left) and wavelets Ψ 1, (last three left), their duals scaling functions Φ1, (three first right) and wavelets Ψ1, (last three right).

Figure 2 :

 2 Figure 2: Scaling functions Φ 0 (three first left) and wavelets Ψ 0 (last three left), their duals scaling functions Φ0 (three first right) and wavelets Ψ0 (last three right) provided Dirichlet boundary condition.

Figure 3 :

 3 Figure 3: Wavelet filter G 0 j (left) and wavelet filter G0 j (right), for j = 6.

Figure 4 :Figure 5 :

 45 Figure 4: Matrices of change of bases L 0 j (left) and L 1 j (right), for j = 6.

Figure 8 :Figure 9 :

 89 Figure 8: Vector field of ∇[Ψ 1, 1 ⊗ Ψ 1, 1 ] and ∇[Ψ 1, 2 ⊗ Ψ 1, 2 ].

3 Figure 10 :

 310 Figure 10: Non linear error approximation: first component u 1 left and second component u 2 right.

Figure 11 :

 11 Figure 11: Error on u 1 (left) and u 2 (right) reconstructed from 20% of their divergence-free wavelet coefficients.

  pour f(x,y)=sin(2*pi*x)x 2 (1-x) 2 y 2 (1-y) 2 k=-2.54

Figure 12 :

 12 Figure 12: Relative l 2 norm of the projection on this base of curlf (x, y) with f (x, y) = sin(2πx)x 2 (1x) 2 y 2 (1 -y) 2 .

Figure 13 :

 13 Figure13: Relative L 2 errors on u and ∇u according to the resolution j, at t = 10 -2 and η = 10 -3 .

Inversely we have:

We assume also that the basis of divergence-free wavelet used is that of the lemma 2.8. The formula are obtained by considering the inner products:

The algorithm of reconstruction is still a consequence of proposition 2.2.

Example 3

We start with a vector field u arising from a numerical simulation of lid driven cavity flow. Then, we compute its divergence-free scaling function and wavelet coefficients using proposition 3.2 and proposition 3.3. The divergence-free wavelets are constructed from the B-Spline generators of Figure 1 with the same parameters. Figure 9 shows the plot of this vector field and corresponding coefficients.

Applications

In this section, we illustrate some practical uses of the divergence-free and curl-free wavelets constructed before. We first show on a numerical example their powerful properties of nonlinear approximation. Then, we present their application on two problems relevant for the numerical simulation of incompressible flows: the Helmholtz decomposition and the Stokes problem, with homogeneous boundary conditions.

Nonlinear approximation

Divergence-free and curl-free wavelet bases provide nonlinear approximation estimates, governed by the approximation orders of one-dimensional spaces involved in their construction. In this part, we investigate the convergence rate obtained from the N -best terms approximation.

For numerical test, we used the vector field and wavelets of example 3. On Figure 10, we plot the 2 error on each component through a solution recovered by non-linear approximation on the divergence-free wavelet basis. We get the same result as on the standard multiresolution analysis provided by V j . This is consistent with the theoretical result proved in [START_REF] Souleymane | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF]. Figure 11 show the classical boundary wavelet error phenomena. It is well know that this error does not prevent the convergence of the multi-scale projector of these bases.

Helmholtz decomposition

The Helmholtz decomposition of a vector field u of (L 2 (Ω)) 2 , is a unique decomposition of u of the form:

with χ ∈ H 1 0 (Ω) and q ∈ H 1 (Ω). The objective in this section is to compute an approximation u div j in V div j of the divergence-free part u div = curl (χ), using the divergence-free bases built in section 2.3. For simplicity, we use the scaling function basis (Φ div j,k ) , since the wavelet one can be deduced using one-dimensional fast wavelet transform along each direction. u div j is searched as its decomposition onto divergence-free scaling functions:

By orthogonality of the decomposition (49) in (L 2 (Ω)) 2 , one obtains:

where M the Gram matrix of the basis {Φ div j,k }. The computation of the coefficients (c div j,k ) is then reduced to the resolution of a linear system of matrix M. This system can be easily inverted, since M is no more than the stiffness matrix of a standard Laplacian onto the scalar scaling function basis {ϕ D j,k 1 ⊗ ϕ D j,k 2 }. Indeed:

More details on the implementation and resolution of (51) can be found in [START_REF] Souleymane | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets[END_REF]. Figure 12 presents the 2 -error of convergence of the algorithm, according to the space resolution j, using wavelets of example 3. The exact solution u corresponds to:

Stokes problem

The Stokes problem is a simple test case for the simulation of incompressible flows. In the non-stationary case, and for a velocity u vanishing at the boundary, it is described by the following equations:

where p is the corresponding pressure.

K. Urban was the first who uses interior divergence-free wavelets for the resolution of the stationary case [START_REF] Urban | Wavelets in Numerical Simulation[END_REF]. His method uses a variational method [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] in H div (Ω) and thus requires the inversion of the stiffness matrix in the divergence-free wavelet basis.

We propose here to use the Helmholtz decomposition to simplify a classical method of resolution, called the Chorin projection method [START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF]. Our algorithm to compute u n (x) ≈ u(x, nδt) is the following:

Starting with initial values u 0 = u(0, x), repeat for 1 ≤ n ≤ N

Step 1: Find a(x) solution of

a = 0 on Γ (55)

Step 2 Find u n+1 solution of

where P is the orthogonal projector from (L 2 (Ω)) 2 onto H div (Ω) computed in practice by the Helmholtz decomposition described in section 4.2.

This method has the advantage of decoupling the resolution of the diffusion term and the incompressibility constraint. Moreover, in more general boundary condition (u = 0 on Γ) there is no need of homogenization for the divergence-free basis like in [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF][START_REF] Urban | Wavelets in Numerical Simulation[END_REF], we can incorporate this boundary condition directly in the basis of H div (Ω) used in the computation of P(a).

On Figure 13, we plot the 2 error on u and on its gradient ∇u in V j . The exact solution u is taken from [START_REF] Krell | Stabilized DDFV schemes for Stokes problem[END_REF]:

We used the divergence-free wavelets of example 3 to compute the projector P.

Conclusion

In this article we have presented a practical construction of divergence-free and irrotational multiresolution analyses and wavelets. Our construction, based on one-dimensional


François Combet
email: francois.combet@lis.inpg.fr

Nadine Martin

Pierre Jaussaud

François Léonard

DETECTION OF A GEAR COUPLING MISALIGNMENT IN A GEAR TESTING DEVICE

On the gear testing device of the LIS laboratory (Grenoble, France), a total of 20 signals are synchronously recorded including shaft acceleration signals in several positions, torque, rotation speed, optical encoder signal of both shafts and currents and tensions of the non-synchronous training motor. After about 3500 hours of using this device, the driving gear shaft broke at the gear coupling position with the training motor shaft, which created an impressive helicoïdal crack . This failure could be due to a flexion strain caused by a misalignment of the two shafts, but this fault was not detected before. In this paper we compare results of different signal processing methods for the detection of this fault. Particularly, we use the phase spectrogram whose advantage is to be able to highlight a slight phase modulation in signals that is not detectable on the amplitude spectrogram. We also attempt to detect the fault using spectral analysis, a non-stationary modelling based on Prony's model and comparison of signals recorded at different times before the crack: the acceleration signals, the torque, and the torsional vibration signal obtained from the shaft encoder signal. The detection results are compared and discussed between the different methods.

INTRODUCTION

In rotating machinery rotor misalignment is known to generate vibration problems. It is measured that about 60 per cent of faults was caused by misalignment [Li 01]. The moment created by a misalignment in a gear coupling is quite complex but its fluctuations occur at twice the shaft speed [Dewell 84]. Thus, the presence of the 2 times shaft speed vibration is known to be an indicator of a misalignment. Furthermore, coupling misalignment does not significantly alter the 1x vibration response but the 2x vibration response is clearly a characteristic signature of misaligned shafts [Sekhar 95]. Misalignment initially affects the 1x response but in the case of severe misalignment, a 2x component appears in the response [Jordan 93]. So an increase of the 2x shaft speed vibration corresponds to an increase of the misalignment. The physical source of these effects may be modelled as a rotor bend for the 1x response and as a rotor asymmetry for the 2x response [Sinha 04].

In this paper, different signals measured on our gear testing device are analysed in order to attempt to find criteria for an early detection of faults that could appear in the gear coupling. Two signal processing methods are used. The first method takes advantage of the non-stationarity of the observed measured. It is based on a classical spectrogram followed by a time-frequency tracking of one harmonic modulation. Then a spectral analysis of the extracted modulation yields interesting parameters for supervising a running system. The second method takes advantage of a known model of the observation and proposes the tracking of local model parameters. In this paper, we aim at studying time-frequency local approaches.

The proposed methods actually do not need data means or synchronous measures. Results are presenting on the gear testing device of the lab (GOTIX bench of LIS-Grenoble France).

In this paper, we focus on the results and references are given for going insight the methods we do not describe. We first describe the experiment and we analyse the torsional vibration and the torque signals after synchronous averaging. We then use the frequency spectrogram for analysing slight frequency modulations in the acceleration signals. A time-frequency Prony analysis is finally applied in order to track the evolution of Prony's parameters. The gear testing device of the LIS laboratory (Grenoble, France) is used to study and develop new signal processing methods adapted to the early detection of a crack or the wear of the tested gears. Twenty signals are synchronously recorded on this device: accelerations at the bearing positions of the shafts, torque, rotation speeds, encoder signals and tachometer signals of both shafts and currents and tensions of the non-synchronous training motor. The sampling frequency is 25 kHz except for the encoder and tachometer signals (100 kHz). This device (fig. 1) has been running for 3500 hours when the gear coupling that links the motor shaft to the gear shaft suddenly broke, creating an impressive helicoidal crack of the gear shaft at its coupling position (fig. 1 right). This appears while we were testing new gears for 300 hours under normal load conditions. The gear coupling has never been changed for the 3500 hours. GOTIX testing device is described on our website (http://www.lis.inpg.fr/english/research/projets/aproj21.htm) where all data can be loaded. Only one personal identification is required.

DESCRIPTION OF THE EXPERIMENT

FIRST OBSERVATIONS

Torsional vibration signal

The torsional vibration signal of the shaft is obtained by phase demodulating the shaft encoder signal [Sweeney 96]. The encoder signal is first band-pass filtered around the carrying frequency (1024x12.3=12.6kHz). The band-pass filter is chosen of width +/-8kHz, which allows high modulation frequency analysis. The spectrum of the torsional vibration is rich in harmonics: it contains harmonics of the driving shaft speed (12.3Hz) of the driven shaft speed (46.7Hz, gear ratio is 57/15=3.8), and also harmonics of the 50Hz (current frequency). The time waveform of the vibration is so complex to analyse. In order to retain the driving shaft synchronous components only and to filter other components, we perform a synchronous averaging based on the tachometer signal of the driving shaft. Figure 2 shows the torsional vibration at time 3200h averaged on 125 rotations of the shaft (25 times the gear coincidence period). The main harmonics are the 1x shaft speed, 57x (gear frequency) and also 18x and 78x. These two last harmonics are likely to be due to the motor constitution.

In order to detect a fault that could appear in the gear coupling, we computed the averaged torsional vibration at different instants from 3200h (beginning of the new test) to 3465h (a few hours before the crack), but we did not find particular evolution of the time waveform and of the main harmonic amplitudes. Since a misalignment fault is characterised by an increase of the second harmonic, we plotted the evolution of the 2x frequency modulation amplitude and compared it to the fundamental one (figure 3). The 2x modulation slightly increases at 3385h and 3455h, which could be an indicator of a misalignment fault. We also plotted the 2x/1x amplitude ratio that may be more valuable for the detection.

Torque signal

As for the torsional vibration signal, the torque is perturbed by harmonics of the driven shaft frequency and of the 50Hz frequency. We so performed a synchronous averaging in order to retain the driving shaft harmonics only. The averaged and zeromean torque at time 3200h is shown on figure 4: the 18x shaft speed harmonic is here preponderant, which confirms that this component must be due to the training motor. The 1x and 57x are also present.

We computed the averaged torque at different instants from 3200h to 3465h, and we plotted the evolution of the 1x and 2x shaft speed modulation amplitudes (figure 5). The 2x modulation turns out to decrease at the end of the experiment, unlike the torsional vibration. The torque seems here to be difficult to interpret for the detection of a fault.

ACCELERATION ANALYSIS BY THE FREQUENCY SPECTROGRAM

The frequency spectrogram

On the well known power spectrogram, only the amplitude of the short-time Fourier transform (STFT) is shown. Thus, a slight frequency modulation may not be clearly visible on the power time-frequency plane. The frequency spectrogram [Léonard 00] was developed as a complement to the power spectrogram. It exploits the phase information of the STFT. The frequency spectrogram is simply computed by the phase difference of the STFT between two consecutive segments divided by the time shift between the segments. Furthermore, its sensibility can be increased by applying an integer gain to the result. The advantage is to be able to highlight a slight frequency or phase modulation that does not clearly appears on the amplitude spectrogram. Figure 6 shows the power and frequency spectrograms of a simulated frequency modulated sinusoid. The slight modulation is not evident on the power spectrogram itself but is better visible on the frequency spectrogram. A cut can then be applied on it at the modulated frequency in order to recover the modulation signal. This modulation can actually be recovered by tracking the maxima position of the power spectrogram, but this needs to highly interpolate the frequency axis in order to limit quantization [Léonard 00]. Thus, the frequency spectrogram permits to decrease the calculation time for recovering a frequency modulation. It also has the advantage to be more visual in a first analysis. We measured the frequency modulation of the second gear harmonic at 1400Hz which here stronger than the fundamental. The rotation period of the driving shaft corresponds to around 2032 samples. In order to track a modulation within the rotation period, a size of 512 samples was chosen for the spectrogram gliding window. Once the modulation signal is recovered by a cut of the frequency spectrogram at 1400Hz, a spectral analysis can be applied on it. Figure 7 shows the spectrum of the modulation signal of acceleration CH6 at the beginning of the new test (3200h) and before the crack (3465h). By comparison of these two spectrums, the frequency modulation amplitudes of the 1x shaft speed (12.3Hz) and of the 2x clearly evolved from 3200h to 3465h: the 2x becomes higher than the 1x at 3465h. Figure 8 shows evolution of these harmonics. A fault seems to appear at 3300-3350h due to an increase of the 2x/1x amplitude ratio. It also increases just before the crack (3465h) for signal CH6.

The originality of the approach is the extraction of the modulation signal in the time-frequency plane followed by a spectral analysis in order to estimate the evolution of the modulation amplitude versus time.

PRONY ANALYSIS

We finally performed a time-frequency Prony analysis [Martin 04] to the synchronous averaged acceleration signals CH5 and CH6. We observed the evolution of the Prony parameters (frequencies, dampings, amplitudes and phases) versus the time of analysis. The main evolution was observed for the damping parameter. Figure 9 shows evolution of the mean value and of the standard deviation of the estimated damping parameter. Damping is often higher for a cracked structure.

CONCLUSION

In a diagnostic and surveillance purpose, this paper presents two signal processing approaches in order to attempt to find criteria for early fault detection. The frequency spectrogram is interesting in order to demodulate frequency modulations when the modulation frequency is not known and/or apt to change before the occurring of a crack. The only a priori knowledge is the modulation bandwidth, whose choice is not critical. The algorithm is simple with a low computation cost. The Prony analysis is a local estimation of model parameters and is suitable for detecting the modification of shocks induced in a system. In our experimentation, the tracking of the damping factor seems to be decisive. Contrary to the first one, this method has a higher computation cost. The proposed methods have been tested on GOTIX bench of LIS-Grenoble. We observed from time 3200h (beginning of experiment) to 3465h (before the crack) an increase of the 2x shaft speed modulation amplitude in the torsional vibration signal and in the acceleration signals close to the gear coupling. This could be an indication that the misalignment increased in the gear coupling, and the friction between the teeth due to the created moment finally led to the crack of the shaft in the coupling.

Both approaches are working in time-frequency space and can be applied to all observations including harmonic modulations or shocks. They are adapted for tracking changes in a great amount of data. The next step now is to deduce a detection algorithm. Works are in progress in that way trying to define a detector, which will not be sensitive to the choice of a threshold, as in [Hory 02].

Fig. 1 :

 1 Fig.1: Left: gear testing device, middle: tested gears, right: cracked gear coupling. The gear coupling that links the training motor shaft to the gear shaft appears in the lower right corner of the middle figure.

Fig. 2 :

 2 Fig.2: Synchronous averaged torsional vibration signal at 3200h and its harmonics: main harmonics are 1x shaft speed, 57x (gear frequency), 18x and 78x.

Fig. 3 :Fig. 4 :

 34 Fig.3: Evolution of the 1x and 2x shaft speed frequency modulation amplitudes estimated from the torsional vibration signal. The 2x modulation slightly increases at 3385h and 3455h (crack is at 3485h).

Fig. 5 :

 5 Fig.5: Evolution of the 1x and 2x shaft speed modulations of the torque.

Fig. 6 :

 6 Fig.6: top: power spectrogram of a simulated frequency modulated sinusoid, middle: frequency spectrogram, bottom: modulation signal recovered by a cut of the frequency spectrogram. Signal to noise ratio is 10dB.

Fig. 7 :

 7 Fig.7: Spectrum of the modulation signal obtained from the frequency spectrogram of acceleration CH6 at time 3200h (left) and 3465h before the crack(right).

Fig. 8 :

 8 Fig.8: Frequency modulation amplitude at 1x and 2x shaft speed of the 2 nd gear frequency harmonic (1400Hz) for accelerations CH5 (left) and CH6 (right).

Fig. 9 :

 9 Fig.9: Evolution of the mean value and standard deviation of the damping parameter of Prony analysis applied to the averaged acceleration signals CH5 and CH6.


Nicolas Praly
email: nicolas.praly@sysnav.fr

Nicolas Petit
email: nicolas.petit@mines-paristech.fr

Julien Laurent-Varin
email: julien.laurent-varin@cnes.fr

Using distributed magnetometry in navigation of heavy launchers and space vehicles

Recently, a new technique has emerged to address the general problem of reconstructing the inertial velocity of a rigid body moving in a magnetically disturbed region. The contribution of this paper is to apply the developed method, in a prospective spirit, to a case of space navigation in view of estimating the performance improvement that could be obtained using state-of-the-art magnetometer technology onboard heavy launchers and other space vehicles. The main underlying idea of the approach is to estimate the inertial velocity by readings of the magnetic field at spatially distributed (known) locations on the rigid body. Mathematically, through a chain-rule differentiation involving variables commonly appearing in classic inertial navigation, an estimate of this velocity can be obtained. In this paper, we show the potential of this method in the field of navigation of heavy launchers passing through particular regions of the Earth magnetosphere as considered, e.g., for upcoming Galileo missions. Numerical results based on the specifications of candidate embedded magnetic sensors stress the relevance of the approach.

Introduction

The overall life expectation of a spacecraft is strongly related to the amount of fuel that is available on-board. Interestingly, the flight phase of injection into orbit of the payload is a critical phase with respect to this factor. The main reason for this is that the spacecraft may have to use a vast amount of its propellant to compensate for launchers position offsets. To minimize the undesired usage of propellants during this phase, a new generation of Ariane 5 launcher (Ariane 5 ME -Midlife Evolution) will be equipped with a re-ignitable engine: the Vinci rocket engine which can be used for the injection task. In this context, a new standard mission appears as very promising: the GTO+. It is a new version of the usual Geosynchronous Transfer Orbit (GTO) which is the historical mission of Ariane 5. The GTO apogee equals the altitude of the geosynchronous orbit, while its perigee lies between 200 km and 250 km of altitude. The GTO+ has the same apogee, but its perigee is significantly higher (approx. 6000 km). To reach this new orbit, the launcher trajectory includes a long ballistic phase lasting approximately 5.5 hours. This very long phase is inconsistent with the accuracy of current inertial navigation systems, whose drift are proportional to the total flight time. Importantly, navigation devices are critical for the success of the mission as they are needed to estimate the (generalized) position of the launcher with respect to an inertial frame of reference. To maintain a good level of accuracy during such long-lasting launch missions, hybridizing the inertial measures with other sources of positioning information seems necessary. Unfortunately, classic technologies such as Global Navigation Satellite Systems (GNSS) or GPS are discarded because they are easily blurred. Lately, CNES and SYSNAV have proposed a new paradigm of navigation exploiting the Earth magnetic field to address this problem. This technique, referred to as "magneto-inertial navigation" and its application to the field of space navigation is the subject of this article. The paper is organized as follows. In Section 2, we expose the principles of the magneto-inertial navigation technique for near-Earth application and report a state-of-the-art on the useable magnetometers. In Section 3, we give a brief overview of the geomagnetic field model used to evaluate the spatial derivatives (gradient) of the magnetic field, which is one of the key factors for magneto-inertial navigation. In Section 4, we evaluate the potential gain of the magneto-inertial navigation for a launcher during a GTO+ mission, and we sketch a method of hybridization.

Magnetic navigation

Notation and reference frames

We consider the following frames of reference:

• The inertial Earth frame (index IF ) which has its origin at the center of the Earth. Its axes are not rotating with the Earth. The IF Oz axis coincides with the Earth's polar axis (which is assumed to have a constant direction).

• The Earth rotating frame (index EF ) which has its origin at the center of the Earth. Its axes are attached to the Earth. The EF Oz axis coincides with the Earth's polar axis.

• The body frame (index BF ) which is fixed to the body under consideration (vehicle). The orthogonal axes are aligned with the roll, pitch and yaw axes of the vehicle, respectively.

Notations are reported in Table 1.

Basic equations of magneto-inertial navigation

Magneto-inertial navigation is a newly introduced technique to evaluate the velocity of a rigid body in a spatially disturbed and time-invariant magnetic field ([R1, R6]). It has been introduced for navigation applications on the ground (for pedestrians [R7] and automotive vehicles [R8], mostly). Interestingly, it can be generalized to space vehicles, at the expense of more sophisticated equations modelling the Earth magnetic field as rotating with the Earth. As will be discussed, magneto-inertial navigation uses inertial measurements to evaluate some parameters. Therefore, the body we study is assumed to be equipped with an Inertial Measurement Unit (IMU).

The main equation used is the differentiation of the magnetic field

= dB B B dX dt t X dt ∂ ∂   +   ∂ ∂   r r r (1)
Let us now detail each term in equation ( 1):

• dt B d
r can be measured with a 3-axes magnetometer attached to the body. The measure is

BF dt B d r . • t B ∂ ∂ r
can not be easily measured. We make the following assumption: at the vicinity of the Earth (<800km), the magnetic field is rotating with the Earth. Therefore,

0 = r r EF t B ∂ ∂ . • [ ] X B ∂ ∂
is the magnetic field gradient. This quantity can be estimated from measurements obtained with an array of 3-axes magnetometers located at known distinct locations in the body

Legendre polynomials

Considering the above measurement principles, we express the vector equation (1) in the Earth frame EF

EF EF EF dt X d X B t B dt B d r r r       ∂ ∂ + ∂ ∂ = (2)
This equation can be transformed under the form to make the velocity with respect to the inertial frame IF appear

        ⊗ Ω +       ∂ ∂ + ∂ ∂ → X dt X d X B t B dt B d EF IF IF EF EF r r r r r = (3)
As one can measure

        ⊗ Ω +       ∂ ∂ + ⊗ Ω - → → X dt X d X B B dt B d EF IF IF EF BF BF r r r r r r = (5)
Let us now discuss each term appearing in equation ( 5):

• EF BF → Ω r
is the sum of : Because all the magnetic measurements are performed in the body frame BF, we now express equation [START_REF]Ariane 5 : User's Manual[END_REF] in this reference frame

(

)

BF IF EF BF BF IF BF BF IF BF BF BF B M B dt B d r r r r r ⊗ Ω - ⊗ Ω - → → → = IF IF BF IF BF dt X d M X B r →       ∂ ∂ + (
)

IF IF EF IF BF IF BF X M X B r r ⊗ Ω       ∂ ∂ + → → (6)
To evaluate the potential accuracy improvement obtained by hybridization of the inertial navigation and the magneto-inertial technique, we detail the accuracy of the standard IMU and the accuracy of embeddable magnetic sensors. We use data from the Ariane 5 IMU and analyse their impact through equation [START_REF] Vissière | Using magnetic disturbances to improve IMU-based position estimation[END_REF]. Preliminary to this analysis, we present a state-of-the-art of magnetometers in the following section.

Magnetometers state-of-the-art

Various sensing technologies could be embedded onboard a heavy launcher such as Ariane 5. Because we need to evaluate the three coordinates of the magnetic field, only vector magnetometer are considered in the study. Table 2 reports a list of the main technologies along with their respective advantages and drawbacks. (*) already used in space applications

As will appear, SQUID sensors seem to be the only available technology measuring the magnetic field with enough accuracy to exploit the principles of magneto-inertial navigation in space applications.

To evaluate the magnetic field gradient, an array of 3-axes magnetometers is used. Two configurations can be considered:

• A tight network of magnetometers will ensure a thermal stability and an increased knowledge of the mechanical properties (distance between the magnetometers, orientation, limited mechanical deformation...). On the other hand, the sensitivity of this arrangement is greatly reduced.

• A spread-out network of magnetometers will ensure a greater sensitivity, but thermal and mechanical fluctuations may reveal much larger than in the compact network.

Model of Earth magnetic field and evaluation of its gradient

There exist numerous models of the Earth magnetic field (e.g. [START_REF]International Geomagnetic Reference Field: the eleventh generation[END_REF], [START_REF] Tsyganenko | A model of the near magnetosphere with a dawn-dusk asymmetry 1[END_REF] and [START_REF] Olson | Quantitative Model of the Magnetospheric Magnetic Field[END_REF]). Earlier, we have formulated the assumption of a time-invariant magnetic field in the Earth frame at low altitudes. An internal-source field model is consistent with the hypothesis. Among such models is the IGRF-2011 ( [START_REF]International Geomagnetic Reference Field: the eleventh generation[END_REF]). This model describes the geomagnetic internal field, and its secular variations, as being the gradient of a scalar potential V that is expressed through a spherical harmonic expansion with respect to spherical geocentric coordinates:

( )

)) ( cos ( ) ( sin ) ( cos = ) , , ( , , , 1 0 = 1 = θ φ φ φ θ m n m n m n n n m N n P m h m g R a a R V +       + ∑ ∑ (7)
Notations are given in Table 1.

The potential model ( 7) can be used to evaluate the spatial derivatives of the magnetic field at any point. In Figure 1, the variations of the components of the magnetic field gradient are reported for altitude variations typical of a launch vehicle performing a GTO mission (from liftoff to GEO altitude). A typical value of 10 -2 to 10 -3 nT/m is observed for altitudes below 800 km. Such values can be effectively measured with an embedded network of SQUID sensors located at 1m from each other, because their theoretical sensitivity is 10 -6 nT/m. Interestingly, thanks to this high sensitivity, the measure of the magnetic field gradient could also be achieved at GEO altitude. However, at this altitude, the assumption of a stationarity of the magnetic field in the Earth frame may prove untrue for multiple reasons, e.g. solar activity.

Performance of magneto-inertial navigation during a GTO mission

Accuracy

Equation ( 6) is the projection of equation ( 5) in the body frame axes set. We now wish to estimate the propagation of errors (uncertainties) in this equation. Using a current Ariane 5 IMU, the orientation error is small (≈10 -8 rad at 1500s). Developing equation ( 5) at first order gives:

B B B dt B d EF IF IF BF IF BF BF r r r r r r r ∆ ⊗ Ω - ∆ ⊗ Ω - ⊗ Ω ∆ - ∆ → → → =         ⊗ Ω +       ∂ ∂ ∆ + → X dt X d X B EF IF IF r r r         ∆ ⊗ Ω + ∆       ∂ ∂ + → X dt X d X B EF IF IF r r r (8)
An estimation of the accuracy for an injection at a GTO injection point (≈700 km from Earth) is given in Table 3, using equation [START_REF] Bristeau | Navigation system for ground vehicles using Temporally Interconnected Observers[END_REF] and the sensors accuracies. The current Ariane 5 injection accuracy [START_REF] Broquet | HiNAV Inertial / GNSS Hybrid Navigation System for Launchers and Re-entry Vehicles[END_REF] is 30 to 10 m/s. Considering the data reported in Table 3, one can conclude that magneto-inertial navigation techniques have the potential to improve the injection accuracy by a factor of 10 to 100.

Typical use of magneto-inertial navigation in launch vehicle and other spacecraft

By contrast to pure inertial systems, magneto-inertial navigation devices can provide an absolute estimate of velocity (which may be uncertain but does not drift as it is not the result of an integration over time). As long as the orders of magnitude of the various parameters remain unchanged, the accuracy on the estimated velocity is constant. Yet, during the launch mission, the magnetic field gradient decreases as discussed earlier and reported in Figure 1. Then, the obtained velocity estimate becomes less and less accurate. With respect to the preceding discussion, one can define a 3 phases navigation strategy for a launch mission. This strategy is illustrated in Figure 2. It combines traditional inertial navigation and magneto-inertial techniques:

• The first phase consists in using the inertial navigation estimates to calibrate the magneto-inertial navigation system. During this phase, the inertial navigation is accurate as the drift is small. • After some time, due to the inertial navigation drift, the magneto-inertial navigation technique provides a better estimate of the velocity than the inertial system. Hybridization of the inertial and magneto-inertial information yields an improvement of the flight parameter estimates. • As the altitude rises, the magneto-inertial navigation velocity estimates loses accuracy. Then, the inertial navigation technique, whose error has been kept small during the previous phase (second phase), is used for the rest of the mission. inertial navigation's velocity time magnetic navigation's velocity magnetic system calibration magnetic/inertial hybridization inertial navigation

Figure 2: Typical inertial and magnetic navigation on a space launch.

Conclusion

Information about absolute and Earth's relative velocity is critical for numerous space applications including putting payloads into orbit. The magneto-inertial navigation technique presented in this paper brings a new solution for this problem. Considering the theoretical performance of SQUID magnetometers, we have shown that this technology has the potential of improving the current launch vehicle navigation accuracy (and thus the injection accuracy) by a factor of 10 to 100. The potential gain is related to i) the magnitude of the magnetic field gradient observed along the trajectory ii) the accuracy of the employed magnetometers. This preliminary study and its promising results call for further investigations. Information on the real near-Earth magnetic gradient, and the impact of the space vehicle electric environment onto the magnetometers sensors array are critical to evaluate the actually obtainable performance of the system. Finally, it is interesting to note that other potential applications for magnetic navigation can be found in the domain of interplanetary space probe navigation. As magneto-inertial navigation takes advantage of the disturbances of the magnetic field, every celestial body with a sufficiently strong and stationary magnetic field and a known motion could allow us to compensate the drift of classic inertial navigation systems.

r

  Time derivative of the vector A in the body frame, expressed in the inertial frame axis set s A / | | A r ∆ First order error on the component of A

  equation (3), and we finally obtain 03.07 Navigation, Estimation, and Filtering Techniques

  angular velocity of the body with respect to IF (measured with a 3-axes gyrometer) and EF IF → Ω r which is the angular velocity of the Earth with respect to IF. • X r is the position of the body and can be determined by integrating the acceleration (measured by the IMU) or the velocity.

Figure 1 :

 1 Figure 1: Variations of the magnetic field gradient along a GTO launch mission.

  IMU parameter taken from[START_REF] Broquet | HiNAV Inertial / GNSS Hybrid Navigation System for Launchers and Re-entry Vehicles[END_REF] (**) injection accuracy taken from[START_REF]Ariane 5 : User's Manual[END_REF]

Table 1 :

 1 Nomenclature and notations

	Symbol	Quantity

. • dt X d r is the body velocity with respect to IF. To use this variable in the navigation algorithms one needs to consider IF dt X d r .

Table 2 :

 2 Magnetometers state-of-the-art

	Technology	Advantages	Drawbacks	Accuracy (nT)
		• Good accuracy
	Fluxgate*	• Most exploited technology in	• creates a magnetic field	10 -2 -10 7
		space
	GMR sensors	• Small size	• poor accuracy	10 1 -10 8
	Magnetoelectric sensors	• Good accuracy • Small size	• nonlinearity for measuring small magnetic field	10 -3 -10 8
	MTJ sensors	• Small size	• poor accuracy	10 1 -10 8
	Search-Coil*	• Very good accuracy	• time dependent magnetic field	10 -5 -10 8
	SQUID sensors	• Best accuracy • Low noise	• low working temperature	10 -6 -10 0

Table 3 :

 3 Accuracy of magnetic navigation at an actual injection point

	Parameter r	Order of Magnitude	Unit
	∆	B dt d B r ∆	BF	10 -6	nT/s

The presented methodology is patent pending and was financed by the French