Mathieu Fabre
email: mathieu.fabre@insa-lyon.fr

Jérôme Pousin
email: jerome.pousin@insa-lyon.fr

Yves Renard
email: yves.renard@insa-lyon.fr

A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method

In this paper, we develop and analyze a finite element fictitious domain approach based on Nitsche's method for the approximation of frictionless contact problems of two deformable elastic bodies. In the proposed method, the geometry of the bodies and the boundary conditions, including the contact condition between the two bodies, are described independently of the mesh of the fictitious domain. We prove that the optimal convergence is preserved. Numerical experiments are provided which confirm the correct behavior of the proposed method.

Introduction

In the vast majority of finite element software, the contact conditions between deformable solids are taken into account through the introduction of Lagrange multipliers and/or penalization terms. The multipliers, which generally approximate the contact stresses, represent some additional unknowns. The approximated problem is then solved in a coupled way or iteratively on the multiplier using Uzawa's algorithm (see e.g. [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]). Recently in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF], it has been proposed an extension to the contact conditions of Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lsung von Dirichlet-Problemen bei Verwendung von Teilrumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF][START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF] which was originally dedicated to Dirichlet's condition. This method combines the advantages of both the penalty and Lagrange multiplier methods since it remains consistent, optimal and avoid the use of multipliers.

In a fictitious domain framework, this paper aims to adapt Nitsche's method to the case of frictionless contact of two elastic solids with the small deformations hypothesis. Frictionless contact is considered to keep the presentation as simpler as possible. However, the analysis extends without additional difficulties to the case of Tresca friction, in a similar way as in [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF]. One of the advantages of the fictitious domain approach comes from the possibility to work with structured meshes regardless of the complexity of the geometry of the bodies and of the potential contact zone. This approach is particularly advantageous in the case of free boundary problems such as shape optimization and fluid-structure interaction. In that case, it prevents the consecutive remeshing which can be very costly, in particular for three-dimensional problems, and which may also generates some instabilities. More generally, a fictitious domain method may be used in the presence of complex or moving geometries to avoid meshing them.

The fictitious domain approach we consider in this work is the one using "cut elements" which is currently a subject of growing interest and is closely related to XFem approach introduced in [START_REF] Möes | A finite element method for cracked growth without remeshing[END_REF] and widely studied since then (see for instance [START_REF] Möes | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF][START_REF] Pierres | A two-scale extended finite element method for modeling 3D crack growth with interfacial contact[END_REF][START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF][START_REF] Nicaise | Optimal convergence analysis for the eXtended Finite Element[END_REF]). The case of a body with a Dirichlet (or transmission) condition with the use of cut-elements is studied in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF] when Lagrange multipliers and a Barbosa-Hughes stabilization are used, and in [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF][START_REF] Annavarapu | A robust Nitsche's formulation for interface problems[END_REF] when Nitsche's method and an additional interior penalty stabilization are considered. This fictitious domain method is to be compared with more classical strategies (see [START_REF] Marchuk | Methods of numerical mathematics[END_REF][START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF][START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrange multiplier method[END_REF][START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Bertoluzza | The fat boundary method: Semi-discrete scheme and some numerical experiments[END_REF] and the references therein) where the elements are not cut. These more classical strategies offer the possibility to leave unchanged the stiffness matrix of the problem. The boundary conditions are then prescribed via additional penalty and Lagrange multiplier terms. However, in classical strategies, it is often quite difficult to obtain an optimal method regarding the convergence order which easily takes into account both Dirichlet and Neumann conditions. The Fictitious domain method with cut elements allows to consider both Dirichlet and Neumann conditions in a rather standard way. The main price to pay is the adaptation of integration methods on cut elements.

In that context of cut elements, our study is focused on the case of two bodies with Nitsche's method for both the Dirichlet condition and the frictionless contact condition.

The outline of the paper is the following. In Section 2, we introduce the contact problem and the fictitious domain situation. Then, in Section 3, the finite element approximation with the use of Nitsche's method is built. In particular, a specific, parameter free stabilization technique is introduced which is necessary to guarantee the optimal rate of convergence. The properties of the approximated problem are described in Section 4 including the existence and uniqueness of a solution to the discrete problem, the consistency and the a priori error analysis. Finally, in Section 5, some two and three-dimensional Hertz-type numerical experiments are presented which illustrate the optimality regarding the convergence of the method.

The unilateral contact problem in a fictitious domain framework

An example of fictitious domain situation is illustrated in Figure 1. Let Ω i , 1 i 2, be two possibly overlapping domains with piecewise C 1 boundaries included in R d , d = 2 or 3, representing the reference configurations of two elastics bodies. Let Ω be a simple shaped polygonal fictitious domain (typically allowing the use of a structured mesh) containing both Ω 1 and Ω 2 . The boundary Γ 1 of Ω 1 (respectively Γ 2 of Ω 2 ) is divided into three non overlapping parts: Γ 1,C the slave potential zone of contact with meas(Γ 1,C ) > 0 (respectively Γ 2,C with meas(Γ 2,C ) > 0); Γ 1,N the Neumann part (respectively Γ 2,N ) and Γ 1,D the Dirichlet part with meas(Γ 1,D ) > 0 (respectively Γ 2,D with meas(Γ 2,D ) > 0).

The two elastic bodies are subjected to volume forces

f = (f 1 , f 2 ) on Ω 1 × Ω 2 , to surface loads = ( 1 , 2 ) on Γ 1,N × Γ 2,N
and satisfy non homogeneous boundary Dirichlet conditions on Γ 1,D × Γ 2,D , the displacement being prescribed to the given value u D = (u 1,D , u 2,D ). We assume small elastic deformation for the two bodies. The linearized strain tensor field is given by ε(v) = 1 2 (∇v + ∇v T ) and the stress tensor field σ = (σ ij ) 1 i,j 2 is given by σ(v) = Aε(v) where A is the fourth order symmetric elasticity tensor satisfying the usual uniform ellipticity and boundedness properties. Consequently, the displacement (u 1 , u 2 ) on Ω 1 × Ω 2 has to satisfy the following set of equations, apart for the contact condition which will be described later:

           Find u = (u 1 , u 2 ) satisfying -divσ(u i ) = f i in Ω i , σ(u i ) = Aε(u i ) in Ω i , u i = u i,D on Γ i,D , σ(u i )n i = i on Γ i,N . (1)
Γ 1 Ω 1 n g Γ 2 Ω 2 Γ 2,D Ω Γ 2,C Γ 1,C Γ 1,N Γ 1,D Γ 2,N
Π : Γ 1,C → Γ 2,C x → Π(x). (2)
In order to simplify the mathematical analysis, the operator Π is assumed to be a C 1 one to one correspondence on Π(Γ 1,C ) (this hypothesis is satisfied, for instance, when Γ i,C are convex and C 1 for i ∈ {1, 2}). The outward unit normal vector n for the contact condition is chosen to be the one of Γ 2,C :

n : Γ 1,C → R d x → n 2 (Π(x)).
The initial gap g between Γ 1,C and Γ 2,C is defined to be the following distance function:

g : Γ 1,C → R x → (x -Π(x)) • n.
For (v 1 , v 2 ) a displacement field defined on Ω 1 × Ω 2 , the normal jump is defined on the slave boundary Γ 1 for the normal displacement as follows:

[[v • n]] = (v 2 • Π -v 1 ) • n.
Concerning the normal stress, we define

σ(v 1 )n 1 = -σ n (v 1 )n + σ t (v 1 ) with σ n (v 1 ) = -σ(v 1 )n 1 • n and σ(v 2 • Π)n 2 • Π = σ n (v 2 • Π)n + σ t (v 2 • Π) with σ n (v 2 • Π) = σ(v 2 • Π)n 2 • Π • n.
This allows to define the normal stress jump as

[[σ(u)n]] = σ(u 1 )n 1 + σ(u 2 • Π)n 2 • Π |det(J Π )| ,
with J Π the Jacobian matrix of Π. This latter expression is derived accordingly with Newton's second law (action-reaction principle) which is expressed on arbitrary elementary surfaces (see Figure 2):

∀ω ⊂ Γ 1,C , ω σ(u 1 )n 1 dΓ = - Π(ω) σ(u 2 )n 2 dΓ = - ω σ(u 2 • Π)n 2 • Π |det(J Π )| dΓ.
These jumps being defined, the unilateral frictionless contact conditions can be expressed on the slave boundary Γ 1,C as follows:

n 2 Ω 2 Ω 1 Γ 2,C σ 2 (u 2 ) σ 1 (u 1 ) ω Π(ω) Γ 1,C
           [[u • n]] g (i), σ n (u 1 ) 0 (ii), σ n (u 1 )([[u • n]] -g) = 0 (iii), [[σ(u)n]] = 0 (iv), σ t (u 1 ) = 0 (v). (3)
Now, let us introduce the Hilbert space V and the convex cone K of admissible displacements:

V = H 1 (Ω 1 ) d × H 1 (Ω 2 ) d , K = {v = (v 1 , v 2 ) ∈ V | v 1 = u 1,D on Γ 1,D and v 2 = u 2,D on Γ 2,D | [[v • n]] -g 0 on Γ 1,C }.
We assume that f belongs to

L 2 (Ω 1 ) d × L 2 (Ω 2 ) d , belongs to L 2 (Γ 1,N ) d × L 2 (Γ 2,N ) d and u D belongs to H 3 2 (Γ 1,D ) d × H 3 2 (Γ 2,D ) d .
We define the bilinear and the linear forms a(., .) and L(.) by

a(u, v) = i=1,2 Ω i σ(u i ) : ε(v i ) dΩ, L(v) = i=1,2 Ω i f i v i dΩ + i=1,2 Γ i,N i v i dΓ.
The weak formulation of Problem ( 1)-(3) as a variational inequality (see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational in-equalities and finite element methods[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]), reads:

Find u ∈ K such that a(u, v -u) L(v -u) ∀v ∈ K. (4)
Stampacchia's Theorem ensures that Problem (4) admits a unique solution.

3 A Nitsche-based finite element approximation

Nitsche's formulation

In this section, we assume that both the solution u and the test functions v are sufficiently regular (for instance, (u, v)

∈ (H 3/2+ν (Ω 1 ) d × H 3/2+ν (Ω 2 ) d ) 2 for ν > 0)
. From the equilibrium equations and Green's formula, we obtain:

a(u, v) - i=1,2 Γ i,D σ(u i )n i • v i dΓ - Γ 1,C σ n (u 1 )[[v • n]] dΓ = L(v).
In order to build Nitsche's formulations for the contact and Dirichlet conditions, the contact conditions are expressed in an equivalent way by extending to our case the formulation given in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. Denoting z + = max(z, 0) and for an arbitrary γ > 0, the contact conditions (3) on Γ 1,C can be equivalently rewritten:

σ n (u 1 ) = - 1 γ [[[u • n]] -g -γσ n (u 1 )] + (5)
Let θ ∈ R be a fixed parameter. This additional parameter for Nitsche's method determines the symmetry properties (see remarks (3.2) and [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]). Then by using [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF] and

[[v • n]] = ([[v • n]] - θγσ n (v)) + θγσ n (v))
, we obtain:

a(u, v) - Γ 1,C θγσ n (u 1 )σ n (v 1 ) dΓ - i=1,2 Γ i,D σ(u i )n i • v i dΓ + Γ 1,C 1 γ [[[u • n]] -g -γσ n (u 1 )] + ([[v • n]] -θγσ n (v 1 )) dΓ = L(v).
Using contact conditions (3), it holds σ n (u 1 ) = σ n (u 2 • Π) |det(J Π )|. In order to ensure the stability, we introduce a stabilized formulation for elements having a small contribution [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF][START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]. We replace σ n (u 1 ) by a convex combination of σ n (u 1 ) and σ n (u 2 •Π) |det(J Π )|. Namely, we define

σ n (u) = tσ n (u 2 • Π) |det(J Π )| + (1 -t)σ n (u 1 ), (6)
for a parameter t ∈ [0, 1] which may be different for an element to an other for the finite element approximation. Note that a similar approach has been developed in [START_REF] Annavarapu | A robust Nitsche's formulation for interface problems[END_REF] where an optimal choice of the fixed parameter t ∈ [0, 1] is proposed. We obtain:

a(u, v) - Γ 1,C θγσ n (u)σ n (v) dΓ - i=1,2 Γ i,D σ n (u)n i • v i dΓ + Γ 1,C 1 γ [[[u • n]] -g -γσ n (u)] + ([[v • n]] -θγσ n (v)) dΓ = L(v).
We did not treat yet the Dirichlet conditions. In order to be coherent with the fictitious domain approach, we also describe the Dirichlet conditions thanks to Nitsche's method [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF][START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF]. Then, writing v i = (v iθγσ(v i )n i ) + θγσ(v i )n i as in the formulation for the contact conditions, we deduce:

- Γ i,D σ(u i )n i • v i dΓ = Γ i,D 1 γ (u i -u i,D -γσ(u i )n i ) • (v i -γθσ(v i )n i ) dΓ - Γ i,D θγσ(u i )n i • σ(v i )n i dΓ. (7)
We obtain the following weak formulation:

a(u, v) + Γ 1,C 1 γ [[[u • n]] -g -γσ n (u)] + ([[v • n]] -θγσ n (v)) dΓ + i=1,2 Γ i,D 1 γ (u i -u i,D -γσ(u i )n i ) • (v i -γθσ(v i )n i ) dΓ - Γ 1,C θγσ n (u)σ n (v) dΓ - i=1,2 Γ i,D θγσ(u i )n i • σ(v i )n i dΓ = L(v) ∀v ∈ V. (8)
Finally, defining the bilinear form

A θγ (u, v) = a(u, v) - Γ 1,C θγσ n (u)σ n (v) dΓ - i=1,2 Γ i,D θγσ(u i )n i • σ(v i )n i dΓ,
our Nitsche-based method reads:

A θγ (u, v) + Γ 1,C 1 γ [[[u • n]] -g -γσ n (u)] + ([[v • n]] -θγσ n (v)) dΓ + i=1,2 Γ i,D 1 γ (u i -u i,D -γσ(u i )n i ) • (v i -γθσ(v i )n i ) dΓ = L(v) ∀v ∈ V. (9)

Discrete Nitsche's formulation

In what follows, Ciarlet's notations [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] are used. Let T h be a family of triangulations of the fictitious domain Ω such that Ω = K∈T h K. Let h K be the diameter of K ∈ T h and h = max K∈T h h K . The family of triangulations is assumed to be regular, i.e. it exists C > 0 such that h K ρ K C where ρ K denotes the radius of the ball inscribed in K. We suppose that the mesh is quasi uniform in the sense that it exists ζ > 0 a constant such that ∀K ∈ T h , h K ζ h. Let K be the fixed reference element (a triangle for d = 2, a tetrahedron for d = 3) and let T K be the geometric transformation which satisfies T K ( K) = K. The family of triangulations is supposed affine, i.e. T K reads as

∀K ∈ T h , T K (x) = J K x + b K , x ∈ K,
where J K ∈ R d,d is the Jacobian matrix of T K being invertible and b K ∈ R d . Thus, we have:

|det(J K )| = mes(K) mes( K) , J K h K /ρ K , J -1 K h K /ρ K .
Remark 3.1. The family of triangulations is regular and affine, so it holds:

|det(J K )| Ch d K , J K Ch K , J -1 K Ch -1 K .
We introduce U h ⊂ H 1 (Ω) a family of finite element spaces indexed by h coming from some order k 1 finite element method defined on T h . Consequently, we suppose the existence of a global interpolation operator π h : C 0 (Ω) → U h and a local one π h K on each element K ∈ T h such that:

∀u ∈ C 0 (Ω), π h (u) K = π h K (u K ) and ∀p ∈ P k (K), π h K (p) = p.
We assume that the finite element method satisfies the following classical local interpolation error estimate for k l 0, u ∈ H l+1 (Ω):

u -π h K u m,K Ch l+1-m |u| l+1,K , with 0 m l k.
Note that, in particular, the classical P k Lagrange finite element method [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] satisfies this estimate.

The approximation spaces for our problem are defined by

V h 1 = (U h ) d Ω 1 , V h 2 = (U h ) d Ω 2 and V h = (V h 1 × V h 2 ).
In the same way, we define the global operators

Π h i : H k+1 (Ω) d → V h i , i = {1, 2} and Π h : H k+1 (Ω) d × H k+1 (Ω) d → V h .
In order to write a discrete approximation of formulation [START_REF] Duvaut | Les inquations en mcanique et en physique[END_REF], let us introduce the following discrete linear operators:

P h τ : V h 1 × V h 2 → L 2 (Γ 1,C ) v → [[v • n]] -τ σ n (v), P h i,τ : V h i → L 2 (Γ i,D ) d v i → v i -τ σ(v i )n i .
Then, a finite element approximation of our Nitsche-based method reads as:

             Find u h ∈ V h such that A θγ (u h , v h ) + Γ 1,C 1 γ [P h γ (u h ) -g] + P h θγ (v h ) dΓ + i=1,2 Γ i,D 1 γ (P h i,γ (u h i ) -u i,D ) • P h i,γθ (v h i ) dΓ = L(v h ) ∀v h ∈ V h . (10)
In the following, we define γ = γ 0 h K .

Remark 3.2. The additional parameter θ is aimed to be chosen in [-1, 1]. The following values of θ are of particular interest: for θ = 1, we recover the symmetric method proposed and analyzed in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]; for θ = 0, we recover a non-symmetric version presented in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] and for θ = -1, we obtain a skew-symmetric version which has the remarkable property that convergence occurs for any value of γ 0 (see [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]).

Remark 3.3. Note that, concerning the Dirichlet conditions, we obtain Nitsche's classical reformulation since the terms on Γ i,D in (10) read

Γ i,D 1 γ (u i -u i,D ) • v i dΓ -θ Γ i,D (u i -u i,D ) • σ(v i )n i dΓ - Γ i,D σ(u i )n i • v i dΓ.
Indeed, the first term is a kind of penalty term for the Dirichlet condition, the second one ensure the symmetry when θ = 1 and the third one ensure the consistency.

Consistency

The advantage of Nitsche's method, compared to penalization, is the consistency of the approximation in the following sense.

Theorem 3.4. Let u be the solution to Problem (1)-( 3). Assume u is sufficiently regular (typically, (u

1 , u 2 ) ∈ H 2+ν (Ω 1 ) d × H 2+ν (Ω 2 ) d , for ν > 0)
, then u is also a solution to the discrete problem [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] replacing u h by u.

Proof. Let u be the solution to (1)-( 3) and take v h ∈ V h . We assume u sufficiently regular such that

σ n (u) ∈ L 2 (Γ 1,C ) and for i = 1, 2, σ n (u i ) ∈ L 2 (Γ i,D ). As a result, P h θγ (u) ∈ L 2 (Γ 1,C ), for i = 1, 2, P h i,θγ (u i ) ∈ L 2 (Γ i,D
) and A θγ (u, v h ) makes sense. On the one hand, we use the definition of P h θγ , P h i,θγ , the reformulations ( 5) and ( 7) to obtain:

A θγ (u, v h ) + Γ 1,C 1 γ [P h γ (u) -g] + P h θγ (v h ) dΓ + i=1,2 Γ i,D 1 γ (P h i,γ (u i ) -u i,D ) • P h i,γθ (v h i ) dΓ = a(u, v h ) - Γ 1,C θγσ n (u)σ n (v h ) dΓ - i=1,2 Γ i,D θγσ(u i )n i • σ(v h i )n i dΓ + Γ 1,C 1 γ [[[u • n]] -g -γσ n (u)] + ([[v h • n]] -θγσ n (v h )) dΓ + i=1,2 Γ i,D 1 γ (u i -u i,D -γσ(u i )n i ) • (v h i -γθσ(v h i )n i ) dΓ = a(u, v h ) - Γ 1,C θγσ n (u)σ n (v h ) dΓ + Γ 1,C 1 γ (-γσ n (u))([[v h • n]] -θγσ n (v h )) dΓ + i=1,2 Γ i,D 1 γ (-γσ(u i )n i ) • (v h i -γθσ(v h i )n i ) dΓ - i=1,2 Γ i,D θγσ(u i )n i • σ(v h i )n i dΓ = a(u, v h ) - Γ 1,C σ n (u)[[v h • n]] dΓ - i=1,2 Γ i,D σ(u i )n i • v h i dΓ = a(u, v h ) - Γ 1,C σ n (u 1 )[[v h • n]] dΓ - i=1,2 Γ i,D σ(u i )n i • v h i dΓ
On the other hand, multiplying by v h i and integrating (1), it holds:

- i=1,2 Ω i divσ(u i )v h i dΩ = i=1,2 Ω i f i v h i dΩ.
Using Green's formula, we have:

-

Ω i divσ(u i )v h i dΩ = Ω i σ(u i ) : ε(v h i ) dΩ - Γ i σ(u i )n i • v h i dΓ i = 1, 2,
with

-
Using the one to one correspondence of the projection, it holds:

Γ 2,C σ(u 2 )n 2 • v h 2 dΓ = Γ 1,C σ(u 2 • Π)n 2 • Π • v h 2 • Π |det(J Π )| dΓ. Hence - i=1,2 Ω i divσ(u i )v h i dΩ = Ω 1 σ(u 1 ) : ε(v h 1 ) dΩ + Ω 2 σ(u 2 ) : ε(v h 2 ) dΩ - Γ 1,C σ(u 1 )n 1 • v h 1 dΓ - Γ 1,C σ(u 2 • Π)n 2 • Π • v h 2 • Π |det(J Π )| dΓ - i=1,2 Γ i,D σ(u i )n i • v h i dΓ - Γ 1,N 1 v h 1 dΓ - Γ 2,N 2 v h 2 dΓ.
Using (3), it holds:

- i=1,2 Ω i divσ(u i )v h i dΩ = a(u, v h ) - i=1,2 Γ i,D σ(u i )n i • v h i dΓ - Γ 1,N 1 v h 1 dΓ - Γ 2,N 2 v h 2 dΓ - Γ 1,C σ n (u 1 )[[v h • n]] dΓ. So a(u, v h )- Γ 1,C σ n (u 1 )[[v h • n]] dΓ - i=1,2 Γ i,D σ(u i )n i • v h i dΓ = L(v h ).
Which ends the proof.

Moreover, formulation ( 9) is formally equivalent to (1) and (3) in the following sense.

Theorem 3.5. Let u ∈ H 2 (Ω 1 ) d × H 2 (Ω 2
) d be a solution to equation ( 9) then u is a solution to (1) and (3).

Proof. For u ∈ H 2 (Ω 1 ) d × H 2 (Ω 2 ) d a solution to (9) and whatever v ∈ H 2 (Ω 1 ) d × H 2 (Ω 2 ) d , it satisfies: Hence Γ 1,C 1 γ [[[u • n]] -g -γσ n (u)] + + σ n (u 1 ) v 1 • n dΓ = 0 ∀v 1 ∈ H 2 (Ω 1 )
and

Γ 1,C 1 γ [[[u • n]] -g -γσ n (u)] + -σ n (u 2 • Π) |det(J Π )| v 2 • Π • n dΓ = 0 ∀v 2 ∈ H 2 (Ω 2 ). Hence 1 γ [[[u • n]] -g -γσ n (u)] + = -σ n (u 1 ) a. e. on Ω 1 ,
which is a formulation equivalent to [START_REF] Brézis | quations et inquations non linaires dans les espaces vectoriels en dualit[END_REF]. Arguing in the same way as above the Neumann and Dirichlet conditions are recovered.

Stabilization method

A stabilization technique is necessary to control the possible bad quality of σ n (u h ) on elements having very small intersection with the real domains. The stabilization used is the one proposed in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF] which consists in using extension of the normal stress on a neighbor element having a sufficiently large intersection with the real domain. The advantage of this stabilization technique is the absence of parameter to fit, except the threshold under which an intersection is considered to be too small. Note that other stabilization techniques are available, such as the so-called ghost penalty stabilization considered in [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF]. For a given small radius 1 > ρ > 0, let R ρ (respectively R ρ) be an operator of approximation of the normal stress of displacements σ n (u h ) (respectively σ(u h i )) which we define thereafter. For K ∈ T h such that K ∩ Γ 1,C , we note S K = {K ∈ T h | K ∩ Π(K) = ∅}. We note also E K , the polynomial extrapolation of an element v h ∈ V h define from K to Ω.

We distinguish three cases to define the stabilized operator R ρ. Let K ∈ T h and K ∩ Γ 1,C = ∅ then:

• if the intersection between K and Ω 1 is sufficiently large i.e. it exists ŷK > 0 such that

B(ŷ K , ρ) ⊂ T -1 K (K ∩ Ω 1 ) (see Figure 3 a)), then R ρ(v h ) K = σ n (v h 1 K ), • otherwise, if it exists K ∈ S K intersecting Ω 2 such that it exists ŷ K > 0 with B(ŷ K , ρ) ⊂ T -1 K ( K ∩ Ω 2 ) (see Figure 3 b)), then R ρ(v h ) K = σ n (E K (v h 2 ) • Π) |det(J Π )|,
• otherwise, we suppose that it exits a neighbor element

K of K such that it exists ŷK > 0 with B(ŷ K , ρ) ⊂ T -1 K (K ∩ Ω 1 ) (see Figure 3 c)), then R ρ(v h ) K = σ n (E K (v h 1 )).
In the same way, we define the operator R ρ on Γ i,D for i = 1, 2: R ρ :

V h i → L 2 (Γ i,D ) d v i → R ρ(v h i ) = σ(v h i )n i ∃ ŷK > 0 such that B(ŷ K , ρ) ⊂ T -1 K (K ∩ Ω 1 ) σ(E K (v h i ))n i otherwise.
Let us introduce the stabilized discrete linear operators:

P h,ρ τ : V h 1 × V h 2 → L 2 (Γ 1,C ) v → [[v • n]] -τ R ρ(v), K K Ω 1 Γ 1,C Ω 2 Ω 2 Γ 1,C Ω 1 Ω 2 c) Otherwise a) If Ω 1 ∩ K is sufficiently large b) If ∃ K ∈ S K such that Ω 2 ∩ K is sufficiently large K K Γ 2,C Ω 1 Γ 2,C K' Γ 1,C Γ 2,C
V h i → L 2 (Γ i,D ) d v i → v i -τ R ρ(v i ).
We define the discrete form of A θγ (., .) as follows:

A h θγ (u h , v h ) = a(u h , v h ) - Γ 1,C θγR ρ(u h )R ρ(v h ) dΓ - i=1,2 Γ i,D θγR ρ(u h i )R ρ(v h i ) dΓ.
The stabilized version of our approximation [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] reads:

             Find u h ∈ V h such that A h θγ (u h , v h ) + Γ 1,C 1 γ [P h,ρ γ (u h ) -g] + P h,ρ θγ (v h ) dΓ + i=1,2 Γ i,D 1 γ (P h,ρ i,γ (u h i ) -u i,D ) • P h,ρ i,γθ (v h i ) dΓ = L(v h ) ∀v h ∈ V h . (11)
Note that strict consistency of this stabilized discrete problem do not occur. However, we have the following result. Theorem 3.6. Let u be the solution to Problem (1)-(3). Assume u is sufficiently regular (typically, (u 1 , u 2 ) ∈ H 2+ν (Ω 1 ) d × H 2+ν (Ω 2 ) d for ν > 0), then u is also a solution to the following problem:

                 a(u, v h ) - Γ 1,C θγσ n (u)R ρ(v h ) dΓ - i=1,2 Γ i,D θγσ(u i )n i • R ρ(v h i ) dΓ + Γ 1,C 1 γ [P h γ (u) -g] + P h,ρ θγ (v h ) dΓ + i=1,2 Γ i,D 1 γ (P h i,γ (u i ) -u i,D ) • P h,ρ i,γθ (v h i ) dΓ = L(v h ) ∀v h ∈ V h . ( 12
)
Proof. The proof can be straightforwardly deduced from the one of Theorem 3.4.

4 Analysis of the Nitsche-based approximation

Existence and uniqueness results

Theorem 4.1. Let γ = γ 0 h K . It exists a unique solution v h ∈ V h to the discrete problem [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF], for all γ 0 > 0 if θ = -1 and for γ 0 > 0 sufficiently small if θ = -1.

Proof. The proof is adapted from [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. The main adaptations concern the fictitious domain framework and in particular the stabilization term, the consideration of two elastic solids and the semi-coercivity of the bilinear form due to the fact that Dirichlet conditions are taken into account with Nitsche's method. We begin by providing some stability and approximation property for operators R ρ and R ρ in lemmas 4.2, 4.5 and 4.6. Then a coercivity property is proved in Lemma 4.7. Finally, the existence and uniqueness result is deduce from the hemi-continuity of the nonlinear operator which corresponds to [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF].

Lemma 4.2. Let v h ∈ V h , there exists a constant C > 0 independent of h such that R ρ(v h ) 2 0,Γ 1,C Ch -1 ( v h 1 2 1,Ω 1 + v h 2 2 1,Ω 2 ) ∀v h ∈ V h . ( 13
)
The proof of this lemma is detailed in the appendix.

Remark 4.3. The following more general operator R ρ could be considered:

R ρ(u h ) K = (1 -t)σ n (E K (u h 2 • Π)) |det(J Π )| + tσ n (E K (u h 1 )),
with t ∈ [0, 1], the element K being K itself or a neighbor element such as the intersection between K and Ω 2 is large enough and the element K being K itself or a neighbor element such as the intersection between K and Ω 1 is large enough. Lemma 4.2 can be easily extended to this operator. When the elastic coefficients in Ω 1 and Ω 2 are equal, a proposed optimum choice is given by (see [START_REF] Annavarapu | A robust Nitsche's formulation for interface problems[END_REF]):

t K = mes(Ω 1 ∩ K) mes(Ω 1 ∩ K) + mes(Ω 2 ∩ K) .
Remark 4.4. When the initial gap between the two bodies vanishes, for ρ sufficiently small either K ∩Ω 1 or K ∩Ω 2 is sufficiently large and thus it is not necessary to consider any neighbor element.

Lemma 4.5. Let u h ∈ V h , Γ i,D be Lipschitz continuous then it exits a constant C > 0 independent of h such that

R ρ(u h 1 ) 2 0,Γ 1,D Ch -1 u h 1 2 1,Ω 1 , and R ρ(u h 2 ) 2 0,Γ 2,D Ch -1 u h 2 2 1,Ω 2 .
The proof of this lemma can be straightforwardly deduced from the one of Lemma 4.2. Now, Let u h , v h ∈ V h and γ = h K γ 0 and using lemma 4.2, it holds:

γ 1 2 R ρ(u h -v h ) 2 0,Γ 1,C Cγ 0 i=1,2 u h i -v h i 2 1,Ω i , γ 1 2 R ρ(u h i -v h i ) 2 0,Γ i,D Cγ 0 u h i -v h i 2 1,Ω i .
Due to the know approximation properties of the stabilized operators on regular and quasi-uniform families of meshes (see [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]), one obtains the following lemma:

Lemma 4.6. For any v ∈ H k+1 (Ω 1 ) d × H k+1 (Ω 2 ) d R ρ(Π h (v)) -σ n (v) 2 0,Γ 1,C Ch 2k-1 i=1,2 v 2 k+1,Ω i , and
i=1,2 R ρ(Π h i (v i )) -σ(v i )n i 2 0,Γ i,D Ch 2k-1 i=1,2 v 2 k+1,Ω i .
The following coercivity property can be stated :

Lemma 4.7. For M > 0 fixed, it exists α > 0 such that for all γ with M γ > 0 the following coercivity property holds:

a(v, v) + 1 2 i=1,2 Γ i,D γ -1 v 2 i dΓ α i=1,2 v i 2 1,Ω i ∀v ∈ V. ( 14
)
The proof of this lemma is detailed in the appendix. Now, by defining the following operator

B h from V h to V h : (B h u h , v h ) 1,Ω = A h θγ (u h , v h ) + i=1,2 Γ i,D 1 γ (P h,ρ i,γ (u h i ) -u i,D ) • P h,ρ i,γθ (v h i ) dΓ + Γ 1,C 1 γ [P h,ρ γ (u h ) -g] + P h,ρ θγ (v h ) dΓ ∀u h , v h ∈ V h , (15)
it is sufficient to prove that B h is hemi-continuous (see the Corollary 15 p. 126 of [START_REF] Brézis | quations et inquations non linaires dans les espaces vectoriels en dualit[END_REF]) to end the proof of Theorem 4.1. The proof of hemi-continuity of B h is detailed in the appendix.

A priori Error analysis

In this section some optimal a priori error estimates are proved for the problem under consideration. The rate of convergence is the same as for standard finite element methods.

Theorem 4.8. Let u be a solution of the stabilized problem (4) belonging to H

3 2 +ν (Ω 1 ) d × H 3 2 +ν (Ω 2 ) d with ν > 0.
1. If θ = -1, we suppose γ 0 > 0 is sufficiently small. The solution u h of the stabilized problem [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF] satisfies the following error estimate:

i=1,2 u i -u h i 2 1,Ω i + γ 1 2 (σ n (u) + 1 γ [P h,ρ γ (u h ) -g] + ) 2 0,Γ 1,C + i=1,2 γ -1 2 (P h i,γ (u i ) -P h,ρ i,γ (u i )) 2 0,Γ i,D C inf v h ∈V h i=1,2 u i -v h i 2 1,Ω i + γ -1 2 (u -v h ) 2 0,Γ 1,C + γ 1 2 (σ n (u) -R ρ(v h )) 2 0,Γ 1,C + i=1,2 γ -1 2 (u i -v h i ) 2 0,Γ i,D + i=1,2 γ 1 2 (σ(u i ) -R ρ(v h i )) 2 0,Γ i,D (16)
with C > 0 a constant independent of h, u and γ 0 .

2. If θ = -1, then for all γ 0 > 0, the solution u h of the stabilized problem [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF] satisfies the error estimate [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF] with C > 0 a constant independent of h and u.

Proof. The proof is also an adaptation to our fictitious domain framework of the one in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. Let v h ∈ V h , using the coercivity inequality [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF] and continuity of the form a(., .) + 1 2

i=1,2 Γ i,D

γ -1 (.) 2 dΓ and Young's inequality, it holds:

α i=1,2 u i -u h i 2 1,Ω i a(u -u h , u -u h ) + 1 2 i=1,2 Γ i,D γ -1 (u i -u h i ) 2 dΓ = a(u -u h , (u -v h ) + (v h -u h )) + 1 2 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D C i=1,2 u i -u h i 1,Ω i u i -v h i 1,Ω i + a(u -u h , v h -u h ) + 1 2 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D α 2 i=1,2 u i -u h i 2 1,Ω i + C 2 2α i=1,2 u i -v h i 2 1,Ω i + 1 2 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + a(u, v h -u h ) -a(u h , v h -u h ).
Hence

α 2 i=1,2 u i -u h i 2 1,Ω i C 2 2α i=1,2 u i -v h i 2 1,Ω i + 1 2 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D +a(u, v h -u h ) -a(u h , v h -u h ). ( 17
)
Let u be the solution to (4), it verifies the stabilized formulation ( 12), thus we have:

a(u, v h -u h ) -a(u h , v h -u h ) = Γ 1,C θγ(σ n (u) -R ρ(u h ))R ρ(v h -u h ) dΓ - i=1,2 Γ i,D 1 γ (P h i,γ (u i ) -P h,ρ i,γ (u h i )) • P h,ρ i,γθ (v h i -u h i ) dΓ + i=1,2 Γ i,D θγ(σ(u i )n i -R ρ(u h i )) • R ρ(v h i -u h i ) dΓ + Γ 1,C 1 γ ([P h,ρ γ (u h ) -g] + -[P h γ (u) -g] + )P h,ρ θγ (v h -u h ) dΓ. (18)
First, using Cauchy-Schwarz and Young's inequalities for β 1 > 0, it holds:

Γ 1,C θγ(σ n (u) -R ρ(u h ))R ρ(v h -u h ) dΓ θ 2 2β 1 γ 1 2 (σ n (u) -R ρ(v h )) 2 0,Γ 1,C + (θ + β 1 2 ) γ 1 2 R ρ(v h -u h ) 2 0,Γ 1,C . (19)
For all a, b ∈ R, we have the following estimate:

([a] + -[b] + )(b -a) -([a] + -[b] + ) 2 . ( 20
)
Then, set:

τ 1 = Γ 1,C 1 γ ([P h,ρ γ (u h ) -g] + -[P h γ (u) -g] + )P h,ρ θγ (v h -u h ) dΓ = Γ 1,C (σ n (u) + 1 γ [P h,ρ γ (u h ) -g] + )P h,ρ θγ (v h -u h ) dΓ.
Using Cauchy-Schwarz and Young's inequalities for β 2 > 0 and β 3 > 0 and (20), it holds:

τ 1 -1 + 1 2β 2 + |1 -θ| 2β 3 γ 1 2 (σ n (u) + 1 γ [P h,ρ γ (u h ) -g] + ) 2 0,Γ 1,C + β 2 2 γ -1 2 (P h,ρ γ (v h ) -P h γ (u)) 2 0,Γ 1,C + |1 -θ| β 3 2 γ 1 2 R ρ(v h -u h ) 2 0,Γ 1,C . (21)
Moreover, set:

τ 2 = - i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D - i=1,2 Γ i,D γ -1 (v h i -u i ) • (P h i,γ (u i ) -P h,ρ i,γ (u h i )) dΓ + i=1,2 Γ i,D (u i -u h i ) • (σ(u i )n i -R ρ(v h i )) dΓ +(1 + θ) i=1,2 Γ i,D (u i -u h i ) • (R ρ(v h i -u h i )) dΓ (22)
Using Young's inequality for β 4 > 0, it holds:

- i=1,2 Γ i,D γ -1 (v h i -u i ) • (P h i,γ (u i ) -P h,ρ i,γ (u h i )) dΓ β 4 2 i=1,2 γ -1 2 (v h i -u i ) 2 0,Γ i,D + 1 2β 4 i=1,2 γ -1 2 (P h i,γ (u i ) -P h,ρ i,γ (u h i )) 2 0,Γ i,D β 4 2 i=1,2 γ -1 2 (v h i -u i ) 2 0,Γ i,D + 1 β 4 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + 1 β 4 i=1,2 γ 1 2 (σ(u i )n i -R ρ(u h i )) 2 0,Γ i,D , (23)
and for

β 5 > 0 i=1,2 Γ i,D (u i -u h i ) • (σ(u i )n i -R ρ(v h i )) dΓ 1 2β 5 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + β 5 2 i=1,2 γ 1 2 (σ(u i )n i -R ρ(v h i )) 2 0,Γ i,D , (24)
and for β 6 > 0

(1 + θ) i=1,2 Γ i,D (u i -u h i ) • (R ρ(v h i -u h i )) dΓ |1 + θ| 2β 6 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + β 6 |1 + θ| 2 i=1,2 γ 1 2 (R ρ(v h i -u h i )) 2 0,Γ i,D . (25)
Using inequalities ( 23), ( 24) and ( 25) in

τ 2 + 1 2 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D
, it holds:

τ 2 + 1 2 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D - 1 2 + 1 β 4 + 1 2β 5 + |1 + θ| 2β 6 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + β 4 2 i=1,2 γ -1 2 (v h i -u i ) 2 0,Γ i,D + 1 β 4 i=1,2 γ 1 2 (σ(u i )n i -R ρ(u h i )) 2 0,Γ i,D + β 5 2 i=1,2 γ 1 2 (σ(u i )n i -R ρ(v h i )) 2 0,Γ i,D + β 6 |1 + θ| 2 i=1,2 γ 1 2 (R ρ(v h i -u h i )) 2 0,Γ i,D . (26)
Gathering ( 18), ( 19), ( 21) and ( 26) in [START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF], it holds:

α 2 i=1,2 u i -u h i 2 1,Ω i C 2 2α i=1,2 u i -v h i 2 1,Ω i + θ 2 2β 1 γ 1 2 (σ n (u) -R ρ(v h )) 2 0,Γ 1,C + θ + β 1 2 |1 -θ| β 3 2 γ 1 2 R ρ(v h -u h ) 2 0,Γ 1,C + -1 + 1 2β 2 + |1 -θ| 2β 3 γ 1 2 (σ n (u) + 1 γ [P h,ρ γ (u h ) -g] + ) 2 0,Γ 1,C + β 2 2 γ -1 2 (P h,ρ γ (v h ) -P h γ (u)) 2 0,Γ 1,C + - 1 2 + 1 β 4 + 1 2β 5 + |1 + θ| 2β 6 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + β 4 2 i=1,2 γ -1 2 (v h i -u i ) 2 0,Γ i,D + 1 β 4 i=1,2 γ 1 2 (σ(u i )n i -R ρ(u h i )) 2 0,Γ i,D + β 5 2 i=1,2 γ 1 2 (σ(u i )n i -R ρ(v h i )) 2 0,Γ i,D + β 6 |1 + θ| 2 i=1,2 γ 1 2 (R ρ(v h i -u h i )) 2 0,Γ i,D . (27)
Using lemmas 4.2 and 4.5, we obtain:

γ 1 2 R ρ(v h -u h ) 2 0,Γ 1,C C 1 γ 0 i=1,2 v h i -u h i 2 1,Ω i 2C 1 γ 0 i=1,2 v h i -u i 2 1,Ω i + i=1,2 u i -u h i 2 1,Ω i , i=1,2 γ 1 2 R ρ(v h i -u h i ) 2 0,Γ i,D C 2 γ 0 i=1,2 v h i -u h i 2 1,Ω i 2C 2 γ 0 i=1,2 v h i -u i 2 1,Ω i + i=1,2 u i -u h i 2 1,Ω i ,
and we know

γ -1 2 P h,ρ γ (v h ) -P h γ (u) 2 0,Γ 1,C 2 γ -1 2 (v h -u) 2 0,Γ 1,C + 2 γ 1 2 (R ρ(v h ) -σ n (u)) 2 0,Γ 1,C
, and using Young's inequality for β > 0, it holds:

- i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D -(1 - 1 β ) i=1,2 γ 1 2 σ(u i )n i -R ρ(u h i ) 2 0,Γ i,D -(1 -β) i=1,2 γ -1 2 P h γ (u) -P h,ρ γ (u h ) 2 0,Γ i,D . (28)
Let θ ∈ R be fixed, if β 2 , β 3 , β 4 , β 5 and β 6 are chosen sufficiently large such that:

-1 + 1 2β 2 + |1 -θ| 2β 3 < - 1 2 , - 1 2 + 1 β 4 + 1 2β 5 + |1 + θ| 2β 6 < - 1 4 .
And if γ 0 is sufficiently small and β < 1, we get the inequality ( 16).

In the case θ = -1, thanks to [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF], it holds:

α 2 i=1,2 u i -u h i 2 1,Ω i C 2 2α i=1,2 u i -v h i 2 1,Ω i + 1 2β 1 γ 1 2 (σ n (u) -R ρ(v h )) 2 0,Γ 1,C +(-1 + β 1 2 + β 3 ) γ 1 2 R ρ(v h -u h ) 2 0,Γ 1,C + -1 + 1 2β 2 + 1 β 3 γ 1 2 (σ n (u) + 1 γ [P γ (u h ) -g] + ) 2 0,Γ 1,C + β 2 2 γ -1 2 P h,ρ γ (v h ) -P h γ (u) 2 0,Γ 1,C + - 1 2 + 1 β 4 + 1 2β 5 i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D + β 4 2 i=1,2 γ -1 2 (v h i -u i ) 2 0,Γ i,D + 1 β 4 i=1,2 γ 1 2 (σ(u i )n i -R ρ(u h i )) 2 0,Γ i,D + β 5 2 i=1,2 γ 1 2 (σ(u i )n i -R ρ(v h i )) 2 0,Γ i,D . (29)
Let η 1 > 0 and η 2 > 0, we take

β 1 = 2η 1 , β 2 = 1 + 1/η 1 , β 3 = 1 + η 1 , β 4 = 2(1 + η 2 ), β 5 = 2(1 + 1/η 2 ), then it holds: α 2 i=1,2 u i -u h i 2 1,Ω i C 2 2α i=1,2 u i -v h i 2 1,Ω i + 1 4η 1 γ 1 2 (σ n (u) -R ρ(v h )) 2 0,Γ 1,C + 2η 1 γ 1 2 R ρ(v h -u h ) 2 0,Γ 1,C - η 1 2(1 + η 1 ) γ 1 2 (σ n (u) + 1 γ [P γ (u h ) -g] + ) 2 0,Γ 1,C + 1 + η 1 2η 1 γ -1 2 P h,ρ γ (v h ) -P h γ (u) 2 0,Γ 1,C - η 2 4(1 + η 2 ) i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D +(1 + η 2 ) i=1,2 γ -1 2 (v h i -u i ) 2 0,Γ i,D + 1 2(1 + η 2 ) i=1,2 γ 1 2 (σ(u i )n i -R ρ(u h i )) 2 0,Γ i,D +(1 + 1/η 2 ) i=1,2 γ 1 2 (σ(u i )n i -R ρ(v h i )) 2 0,Γ i,D
.

Using [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF] and

β = η 2 2 + η 2
< 1, we have:

- η 2 4(1 + η 2 ) i=1,2 γ -1 2 (u i -u h i ) 2 0,Γ i,D 1 2(1 + η 2 ) i=1,2 γ 1 2 σ(u i )n i -R ρ(u h i ) 2 0,Γ i,D - η 2 2(1 + η 2 )(2 + η 2 ) i=1,2 γ -1 2 P h γ (u) -P h,ρ γ (u h ) 2 0,Γ i,D
, and i=1,2

γ 1 2 σ(u i )n i -R ρ(u h i ) 2 0,Γ i,D 2 i=1,2 γ 1 2 σ(u i )n i -R ρ(v h i ) 2 0,Γ i,D +2 i=1,2 γ 1 2 σ(R ρ(v h i -u h i ) 2 0,Γ i,D
.

Let γ 0 be positive. If we take η 1 = α/(32C 1 γ 0 ) and η 2 = C 2 γ 0 /(32α), then we get the inequality [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]. This ends the proof of Theorem 4.8

Theorem 4.9. Let u be a solution of the variational problem [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF]. Suppose that u belongs to (H

3 2 +ν (Ω 1 )) d × (H 3 2 +ν (Ω 2 )) d with 1/2 ν > 0 if k = 1 and with 1 > ν > 0 if k = 2.
Then, if additionally γ 0 > 0 is sufficiently small when θ = -1, the solution u h of the stabilize problem [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] satisfies the following a priori error estimate:

i=1,2 u i -u h i 2 1,Ω i + γ 1 2 (σ n (u) + 1 γ [P h,ρ γ (u h ) -g] + ) 2 0,Γ 1,C + i=1,2 γ -1 2 (P h,ρ i,γ (u h i ) -P h i,γ (u i )) 2 0,Γ i,D Ch 1+2ν i=1,2 u 2 3 2 +ν,Ω i (30)
with C > 0 a constant independent of h and u.

Proof. Now let us establish the inequality (30). Set v h i = Π h i (u i ), we have the following estimates: 16), γ = γ 0 h and we use the previous inequalities, we get (30). We can write:

u i -Π h i (u i ) m,Ω i Ch k+1-m u i k+1,Ω i , R ρ(Π h (u)) -σ n (u) 2 0,Γ 1,C Ch 2k-1 i=1,2 u 2 k+1,Ω i , and
i=1,2 R ρ(Π h i (u i )) -σ(u i )n i 2 0,Γ i,D Ch 2k-1 i=1,2 u 2 k+1,Ω i . If we replace v h i by Π h i (u i ) in (
i=1,2 u i -u h i 1,Ω i + γ 1 2 (σ n (u) + 1 γ [P h,ρ γ (u h ) -g] + ) 0,Γ 1,C + i=1,2 γ -1 2 (P h,ρ i,γ (u h i ) -P h i,γ (u i )) 0,Γ i,D Ch 1/2+ν i=1,2 u 3 2 +ν,Ω i (31)

Numerical study

This section is dedicated to some numerical experiments with isoparametric Lagrange P 1 or P 2 finite element methods. The accuracy of the method is discussed for the different cases with respect to the finite element used, the mesh size and the value of the parameter γ 0 . Note that the following results are obtained without the stabilization introduced in Section 3.4. From a numerical viewpoint, the stabilization seems not strictly necessary to obtain an optimal rate of convergence. This has already been observed in a linear case in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]. The numerical tests in two dimensions (resp. three dimensions) are performed on a fictitious domain Ω = ] -0.5, 0.5[ 2 (resp. Ω = ] -0.5, 0.5[ 3 which contains the first body Ω 1 , a circle of radius 0.25 and center (0, 0) (resp. a sphere of radius 0.25 and center (0, 0, 0)), and the second Ω 2 = ] -0.5, 0.5[×] -0.5, -0.25[ (resp. Ω 2 = ] -0.5, 0.5[ 2 ×] -0.5, -0.25[). A Dirichlet condition is prescribed on the bottom of the rectangle (resp. parallelepiped). The projector Π is defined from the lower part of the boundary of Ω 1 (i.e. for Γ 1,C = {x ∈ ∂Ω 1 : x d ≤ 0}) onto its projection on the top boundary of Ω 2 . All remaining parts of the boundaries of Ω 1 and Ω 2 are considered traction free.

Since no Dirichlet condition is applied on Ω 1 , the problem is only semi-coercive. In order to recover the uniqueness of the solution, it is needed to prescribe the horizontal rigid translation in 2D and two horizontal translations and one rotation in 3D. This is done by prescribing the displacement on some given convenient points.

We use a generalized Newton's method to solve the discrete problem [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] (see [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] for more details) and our finite element library GetFEM++ 1 . The tool for fictitious domain methods of GetFEM++ has been used which provides cut integration methods. The geometries are described with zero level sets of some signed distances to the domain boundaries. The distance functions are approximated by quadratic Lagrange finite elements. In order to build cut integration methods, each element of the mesh which crosses a domain boundary is cut into a set of sub-elements conforming to this boundary. Then, an integration method is produced on each sub-element lying on the interior of a domain and on each sub-element boundary lying on a domain boundary. In order to obtain a convenient order for the produced integration methods and for the approximation of the domain boundaries, curved sub-elements are used.

Moreover, no specific treatment have been considered for the fact that boundary terms for the contact condition approximated by Nitsche's method is non-regular (due to the positive part). We used an order four numerical integration method on each sub-element and we noted no improvement of the accuracy with higher order or refined numerical integration method. For simplicity, we consider a dimensionless situation with Lamé coefficients λ = 1 and µ = 1 and a vertical volume density of force -0.1.

The situation studied is not strictly speaking of Hertz type due to the fact that Ω 2 is bounded. The expression of the exact solution being unknown, the convergence is studied with respect to a reference solution computed with a P 2 isoparametric element on a very fine mesh (h = 1/200 in 2D and h = 1/30 in 3D) with the skew-symmetric method θ = -1 (see Figures 4 and5).

Numerical convergence in the two dimensional case

We perform a numerical convergence study on the three methods θ = 1, θ = 0 and θ = -1 for a fixed parameter γ 0 = 1/200 (chosen small in order to have the convergence for the three cases). On figures 6, 7 and 8, the relative error in percentage in L 2 and H 1 -norms on each bodies for P 1 Lagrange finite elements are plotted. As expected the optimal convergence is obtained in H 1 -norm for all methods in good accordance with Theorem 4.9. The rate of convergence in L 2 -norm is slightly sub-optimal on Ω 2 if one refers to Aubin-Nitsche lemma in the linear case. However, such a result is not available for the nonlinear contact problem. Moreover, this slight sub-optimal convergence may be caused by the Neumann-Dirichlet transition at the bottom of Ω 2 . On figures 9, 10 and 11, the same experiments are reported but for P 2 isoparametric Lagrange finite elements. The convergence rate for the three cases is close to 1.6 on Ω 1 and 1.3 on Ω 2 . This is also close to optimality if one takes into account that the expected maximal regularity of the displacement next to the transition between contact and non-contact should be H 5/2-η for any η > 0 (However, this result has only been proved in a scalar case in [START_REF] Moussaoui | Regularité des solutions d un problème mêlé Dirichlet-Signorini dans un domaine polygonal plan[END_REF]). Accordingly, one could expect that the convergence rate in the L 2 -norm would be close to 2.5. This is approximately the case with again some sub-optimal rates which may due to the nonlinear characteristic of the contact condition or to the presence of non-regularities on the transition between the Dirichlet and the Neumann condition.

Influence of the parameter γ 0

The influence of γ 0 on the H 1 -norm of the error is plotted in Figure 12 for P 1 elements and on Figure 13 for P 2 elements. The most affected method is the one for θ = 1. Indeed, it converges only for γ 0 very small. The large oscillation in the error norm comes from the fact that Newton's algorithm do not fully converge for all numerical experiments probably because there is no solution to the discrete problem in some cases. The method for θ = 0 gives a more regular error with respect to γ 0 . It is still important to have γ 0 small to keep a good solution but a larger value is allowed. Accordingly to the theoretical result of Theorem 4.9, the influence of γ 0 on the method θ = -1 is more limited. There is only a slight increase of the error for large values of γ 0 . Note that the nonlinear discrete system [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] becomes very stiff when γ 0 is very small. Thus, the possibility to have a large γ 0 is an advantage.

Numerical experiments in the 3D case

Due to the high number of degrees of freedom in 3D, it obviously has not been possible to produce convergence curves with a mesh size as small as in 2D. The convergence curves for 3D are shown in Figures 14, 15 and 16 only for P 1 elements. Although we also made some tests with P 2 elements and on the influence of γ 0 , we do not reproduce them for brevity of the paper. Indeed, the conclusions that can be drawn are were very similar to the 2D case.

Conclusion

In this paper, we developed a fictitious domain approach for the approximation in small deformations of the frictionless contact with nonzero initial gap of two elastic bodies. The main ingredients are the adaptation of Nitsche's method for the contact condition introduced in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] and the fictitious domain method (inspired by the X-fem) developed in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF] including the stabilization proposed for the elements having a small intersection with the real domains. Perspective of this works would be to weakened the conditions on the projection operator Π to include for instance non regular situations such as the one illustrated in Figure 17 where Π is only piecewise regular. Another possibility would be to consider a non-orthogonal projection. As already mentioned, the analysis can be easily adapted to Tresca friction similarly as it has been done in [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF] for the non-fictitious domain situation.

Γ 2,C Γ 1,C σ 1 (u 1 ) σ 2 (u 2 ) Π(ω) ω
From this study we conclude that the presented method allow an optimal approximation of unilateral contact problems for affine and quadratic finite element methods. The method for θ = 1 is symmetric which can be an advantage for the numerical solving but requires a very small parameter γ 0 which may lead to a very stiff discrete problem [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF]. The method for θ = 0 has the advantage of the simplicity and allows the use of a moderate γ 0 . Finally, the skew-symmetric method θ = -1 allows the use of larger value of γ 0 which can be a real advantage for the solving of the discrete problem.

We define Γ1 = T -1 K (Γ 1,C ∩ K) and σ(u 1 ) = σ(u h 1 ) • T K and n1 a unit normal vector on Γ1,C .

Γ 1,C ∩K σ n (u h 1 ) 2 dΓ = Γ1 σn (u 1 ) 2 |det(J K )| J -1 K n1 d Γ = Γ1 |σ(u 1 )n • n| 2 |det(J K )| J -1 K n1 d Γ Ch d-1 K | σ(u 1 ) | 2 ∞, K Γ1 . ( 32
) because |σ(u 1 )n • n| | σ(u 1 ) | F n 2 2 = | σ(u 1 ) | F . Moreover, Γ 1
is bounded, indeed the operator T K is a continuous one to one correspondence. Now using the equivalence of norms in P k ( K) d , we have:

| σ(u 1 ) | 2 ∞, K | σ(u 1 ) | 2 ∞,B(ŷ K ,2) = | σ(u 1 ) • t -ŷ K | 2 ∞,B(0,2) C | σ(u 1 ) • t -ŷ K | 2 2,B(0,ρ) = C | σ(u 1 ) | 2 2,B(ŷ K ,ρ) C | σ(u 1 ) | 2 2,T -1 K (Ω 1 ∩K)) = C T -1 K (Ω 1 ∩K) | σ(u 1 ) | 2 F dx. ( 33
)
Using the upper bound of Γ1 and the previous inequalities, it holds:

Γ 1,C ∩K R ρ(u h ) 2 dΓ C h d-1 K h d K T -1 K (Ω 1 ∩K) | σ(u 1 ) | 2 F |det(J K )| dx Ch -1 K Ω 1 ∩K σ(u h 1 ) 2 F dx Ch -1 K Ω 1 ∩K A∇u h 1 2 F dx Ch -1 K Ω 1 ∩K ∇u h 1 2 2 dx. • Otherwise, if ∃ K ∈ S K such as ∃ ŷ K > 0 such that B(ŷ K , ρ) ⊂ T -1 K ( K ∩ Ω 2 ), then R ρ(v h ) K = σ n (E K (v h 2 ) • Π) |det(∇P i)| and using the continuous of J Π i.e. |det(J Π )| C, it holds: R ρ(u h ) 2 0,Γ 1,C ∩K = Γ 1,C ∩K σ n (E K (u h 2 ) • Π) 2 |det(J Π )| 2 dΓ C Γ 1,C ∩K σ n (E K (u h 2 ) • Π) 2 |det(J Π )| dΓ C Π(Γ 1,C ∩K) σ n (E K (u h 2 )) 2 dΓ C K∈S K Γ 2,C ∩K σ n (E K (u h 2 )) 2 dΓ. We define Γ2 = T -1 K ( K∈S K Γ 2,C ∩ K) and σ(u 2 ) = σ(u h 2 )
• T K and n2 a unit normal vector on Γ2,C . As previously, we have Γ2 bounded. In the same way as in (32), we have:

K∈S K Γ 2,C ∩K σ n (E K (u h 2 )) 2 dΓ Ch d-1 K | σ(û 2 ) | 2 ∞, Γ2 Γ2 .
Now using the equivalence of norms in P k ( K) d and in the same way as in (33), we obtain:

| σ(û 2 ) | 2 ∞, Γ2 | σ(û 2 ) | 2 ∞,B(ŷ K ,2 ||Π|| ||T -1 K ||) C T -1 K (Ω 1 ∩ K) | σ(û 2 ) | 2 F dx.
Hence, using the previous inequalities, it holds:

Γ 1,C ∩K R ρ(u h ) 2 dΓ C h d-1 K h d K T -1 K (Ω 1 ∩ K) | σ(u 1 ) | 2 F det(J K ) dx Ch -1 K Ω 1 ∩ K ∇u h 2 2 2
dx.

• Otherwise, we suppose it exits an neighbor element K of K such that ∃ ŷK > 0 such that

B(ŷ K , ρ) ⊂ T -1 K (K ∩ Ω 1 ), then R ρ(v h ) K = σ n (E K (v h 1 )
). Then, it holds:

R ρ(u h ) 2 0,Γ 1,C ∩K = Γ 1,C ∩K σ n (E K (u h 1 )) 2 dΓ
We define by Γ

1 = T -1 K (Γ 1,C ∩ K) and σ (u 1 ) = σ(u h 1 )
• T K and by n1 a unit normal vector on Γ1,C . As previously, we have Γ 1 bounded. In the same way as in (32), we have:

Γ 1,C ∩K σ n (E K (u h 1 )) 2 dΓ = Γ 1 σ (u 1 ) 2 |det(J K )| J -1 K n1 d Γ Ch d-1 K | σ (u 1 ) | 2 ∞,T -1 K (K) Γ 1 .
Now using the equivalence of norms in P k ( K) d and in the same way as in (33), we have:

σ(û 1 ) 2 ∞,T -1 K (K) σ(û 1 ) 2 ∞,B(ŷ K ,4) C T -1 K (Ω 1 ∩K ) σ(û 1 ) 2 F dx.
Hence, using the previous inequalities, it holds:

Γ 1,C ∩K R ρ(u h ) 2 dΓ C h d-1 K h d K T -1 K (Ω 1 ∩K ) σ(û 1 ) 2 F |det(J K )| dx Ch -1 K Ω 1 ∩K ∇u h 1 2 2 dx.
Finally, by iterating on all the elements K intersecting Γ 1,C and using the quasi uniformity of the mesh, we obtain [START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF].

Proof of lemma 4.7. We argue by contradiction. It is sufficient to prove the result for γ = M . Suppose there exists (v n ) n∈N ⊂ V such that i=1,2 v i,n

Hence, it holds lim n→+∞ i=1,2 Γ i,D v 2 i,n dΓ = 0 and lim n→+∞ a(v n , v n ) = 0. From the weak sequential compactness of the unit ball of V , there exits v ∈ V and a subsequence still denoted by v n which weakly converges to v. The compact injection of H 1 into L 2 implies that up to a subsequence, v n converges to v strongly in L 2 (Ω 1 ) d × L 2 (Ω 2 ) d . First, we show that v = 0 and then that v n converges to v strongly in H 1 (Ω 1 ) d × H 1 (Ω 2 ) d . By using the lower semi-continuity of v → i=1,2 Γ i,D v 2 i,n dΓ, we have v i,n 0,Ω i = i=1,2 v i 0,Ω i . Similarly by using the weak lower semi-continuity of a(., .), we deduce a(v, v) = 0 and using the property of the fourth order tensor A, it holds:

0 = Ω i σ(v i ) : ε(v i ) dΩ = Ω i ε(v i ) : ε(v i ) dΩ = ε(v i ) 0,Ω i .
Let us finally show that v = 0. Since, the tensor A is uniformly elliptic, it holds:

a(v n , v n ) C i=1,2 Ω i ε(v i,n ) : ε(v i,n ) dΩ = C i=1,2 ε(v i,n ) 0,Ω i and lim n→+∞ a(v n , v n ) = 0. Hence lim n→+∞ i=1,2 ε(v i,n ) 0,Ω i = i=1,2 ε(v i ) 1,Ω i = 0.
Moreover, thanks to Korn's inequality (see [START_REF] Duvaut | Les inquations en mcanique et en physique[END_REF]), it holds: v i,n 0,Ω i .

We deduce: lim

n→+∞ i=1,2 v i,n 1,Ω i = i=1,2 v i 1,Ω i = 0 which contradicts i=1,2 v i,n 2 1
,Ω i = 1.

Proof for operator B h to be hemi-continuous (for the proof of Theorem 4.1). First, we need to prove B h is coercive which is a consequence of the previous lemmas. Then we establish an estimate which will imply the hemi-continuity. Let u h , v h ∈ V h , it holds:

(B h u h -B h v h , u h -v h ) 1,Ω = I + II + III (34)
with 1 γ ([P h,ρ γ (u h )g] + -[P h,ρ γ (v h )g] + )(P h,ρ θγ (u h ) -P h,ρ θγ (v h )) dΓ. Now, we need to bound I, II, III from below to prove the coercivity. Using Young's inequality for β > 0, it holds:

I = a(u h -v h , u h -v h ) - Γ 1,C θγR ρ(u h -v h )R ρ(u h -v h ) dΓ - i=1,2 Γ i,D θγR ρ(u h i -v h i ) • R ρ(u h i -v h i ) dΓ, II = i=1,2 Γ i,D 1
III (1 - |1 -θ| 2β ) γ -1 2 ([P h,ρ γ (u h ) -g] + -[P h,ρ γ (v h ) -g] + ) 2 0,Γ 1,C - |1 -θ| β 2 γ 1 2 R ρ(u h -v h ) 2 0,Γ 1,C (35)
Using Young's inequality for β > 0:

II (1 - |1 + θ| 2β ) i=1,2 γ -1 2 (u h i -v h i ) 2 Γ i,D +(θ - |1 + θ| β 2 ) i=1,2 γ 1 2 (R ρ(u h i ) -R ρ(v h i )) 2 Γ i,D . (36)
We deduced from the estimates of II and III that: .

(B h u h -B h v h , u h -v h ) 1,Ω a(u h -v h , u h -v h ) -θ γ 1 2 R ρ(u h -v h ) 2 0,Γ 1,C -θ i=1,2 γ 1 2 (R ρ(u h i ) -R ρ(v h i )) 2 Γ i,D +(1 - |1 -θ| 2β ) γ -1 2 ([P h,ρ γ (u h,ρ ) -g] + -[P h,ρ γ (v h ) -g] + ) 2 0,Γ 1,C -|1-θ|β 2 γ 1 2 R ρ(u h -v h ) 2 0,Γ 1,C +(1 - |1 + θ| 2β ) i=1,2 γ -1 2 (u h i -v h i ) 2 Γ i,D + (θ - |1 + θ| β 2 ) i=1,2 γ 1 2 (R ρ(u h i ) -R ρ(v h i )) 2 Γ i,D a(u h -v h , u h -v h ) + (1 - |1 + θ| 2β ) i=1,2 γ -1 2 (u h i -v h i ) 2 Γ i,D - |1 + θ| β 2 i=1,2 γ 1 2 (R ρ(u h i ) -R ρ(v h i ))
If θ = 1 and β = 2, we have:

(B h u h -B h v h , u h -v h ) 1,Ω a(u h -v h , u h -v h ) + 1 2 i=1,2 γ -1 2 (u h i -v h i ) 2 Γ i,D -2 i=1,2 γ 1 2 (R ρ(u h i ) -R ρ(v h i )) 2 Γ i,D + γ -1 2 ([P h,ρ γ (u h ) -g] + -[P h,ρ γ (v h ) -g] + ) 2 0,Γ 1,C -γ 1 2 R ρ(u h -v h ) 2 0,Γ 1,C
.

Thus, if γ 0 is sufficiently small and using the coercivity ( 14) for I and the previous lemma 4.7:

(B h u h -B h v h , u h -v h ) 1,Ω C i=1,2 u h i -v h i 2 1,Ω i . If θ = -1, choose β = |1-θ| 2 , it holds: (B h u h -B h v h , u h -v h ) 1,Ω a(u h -v h , u h -v h ) + i=1,2 γ -1 2 (u h i -v h i ) 2 Γ i,D

Figure 1 :

 1 Figure 1: Example of fictitious domain situation for a contact problem between two elastics bodies with an example of structured mesh.

Figure 2 :

 2 Figure 2: An example illustrating the action-reaction principle between the two bodies.

Figure 3 :

 3 Figure 3: The different cases for the definition of R ρ.

Figure 4 :

 4 Figure 4: 2D numerical reference solution with contour plot of Von Mises stress. Parameters h = 1/400, γ 0 = 1/200, θ = -1 and P 2 elements.

Figure 5 :

 5 Figure 5: Cross-section of 3D numerical reference solution with contour plot of Von Mises stress. Parameters h = 1/30, γ 0 = 1/200, θ = -1 and P 2 elements.

Figure 6 :

 6 Figure 6: Convergence curves in 2D for the method θ = 1, with γ 0 = 1/200 and P 1 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 7 :

 7 Figure 7: Convergence curves in 2D for the method θ = 0, with γ 0 = 1/200 and P 1 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 8 :

 8 Figure 8: Convergence curves in 2D for the method θ = -1, with γ 0 = 1/200 and P 1 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 9 :

 9 Figure 9: Convergence curves in 2D for the method θ = 1, with γ 0 = 1/200 and P 2 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 10 :

 10 Figure 10: Convergence curves in 2D for the method θ = 0, with γ 0 = 1/200 and P 2 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 11 :

 11 Figure 11: Convergence curves in 2D for the method θ = -1, with γ 0 = 1/200 and P 2 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

1 Figure 12 :

 112 Figure 12: Influence γ 0 on the relative H 1 -norm of the error on Ω 1 (on the left) and on Ω 2 (on the right) in 2D for h = 1/90 and P 1 elements.

1 Figure 13 :

 113 Figure13: Influence of γ 0 on the relative H 1 -norm of the error on Ω 1 (on the left) and on Ω 2 (on the right) in 2D for h = 1/90 and P 2 elements.

Figure 14 :

 14 Figure 14: Convergence curves in 3D for the method θ = 1, with γ 0 = 1/100 and P 1 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 15 :

 15 Figure 15: Convergence curves in 3D for the method θ = 0, with γ 0 = 1/100 and P 1 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 16 :

 16 Figure 16: Convergence curves in 3D for the method θ = -1, with γ 0 = 1/100 and P 1 finite elements for the relative L 2 -norm of the error (on the left) and the relative H 1 -norm of the error (on the right).

Figure 17 :

 17 Figure 17: Example of non regular situation on Γ 1,C .

i=1, 2 Γ i,D v 2 i

 22 dΓ = 0 with mes(Γ i,D ) = 0. Furthermore, due to theL 2 -convergence, one has lim n→+∞ i=1,2

i=1, 2 ε 2 v

 22 (v i,n ) 0,Ω i + i=1,i,n 0,Γ i,D C i=1,2

  h,ρ γ (u h )g] + -[P h,ρ γ (v h )g] + )

Γ i σ(u i )n i •v h i dΓ = -Γ i,D σ(u i )n i •v h i dΓ-Γ i,N σ(u i )n i •v h i dΓ-Γ i,C σ(u i )n i •v h i dΓ i = 1, 2, -Γ i σ(u i )n i • v h i dΓ = -Γ i,D σ(u i )n i • v h i dΓ -Γ i,N i v h i dΓ -Γ i,C σ(u i )n i • v h i dΓ i = 1, 2.

Ω i (divσ(u i ) + f i )v i dΩ = 0 ∀v i ∈ H 2 (Ω i ) i.e. -divσ(u i ) = f i a.e. in Ω i , 1 i 2.We have, forall v ∈ H 2 (Ω 1 ) d × H 2 (Ω 2 ) d : Γ 1,C 1 γ [[[u • n]]gγσ n (u)] + [[v • n]] dΓ + Γ 1,C σ(u 1 )n 1 • v 1 dΓ + Γ 2,C σ(u 2 )n 2 • v 2 dΓ = 0, Γ 1,C 1 γ [[[u • n]]gγσ n (u)] + (v 1v 2 • Π • n) dΓ -Γ 1,C σ n (u 1 )n • v 1 dΓ + Γ 1,C σ n (u 2 • Π)n • v 2 • Π |det(J Π )| dΓ = 0.

see http://download.gna.org/getfem/html/homepage/

1,Ω i = 1, for n ∈ N, which satisfies a(v n , v n ) + 1 2 i=1,2 Γ i,D M -1 v 2 i,n dΓ 1 n .

Appendix

Proof of lemma 4.2. First, we define the following matrix norms:

where | . | F is Frobenius' norm. If v is a fixed vector, we define the translation of a vector u, by t v (u) = u+v. In the following, the constant C may vary from a line to another but is independent of h. In order to prove (13), we distinguish the three different cases from the definition of R ρ. First, by using the geometric transformation, the integral is expressed on the reference element.

Then by using the equivalence of the infinity norm with the 2-norm located on a ball, we are able to deal with the 2-norm located on the current element. Finally by using the definition of the stress tensor, we obtain the result.

) and it holds:

and from the coercivity ( 14) for I, we obtain:

and β = |1 + θ|, it holds:

.

So, using γ 0 sufficiently small and using the coercivity ( 14) for I and previous lemma 4.7, it holds:

. Now, we prove the hemi-continuity of B h . Let t, s ∈ [0, 1] and u h , v h ∈ V h , we have:

For all a, b in R, we have the following estimate:

So we deduce that

Hence

Hence B h is hemi-continuous.


A Brownian quasi-helix in IR 4 , built from an automatic sequence The history of Brownian motion was described a number of times. In 1905, Einstein established his celebrated formula

∆X 2 = RT N 1 3πµa ∆t
for spherical particles of radius a suspended in a liquid of viscosity µ at temperature T ; the first member, ∆X 2 , is the average of the squares of their displacements during an interval of time ∆t ; R is the constant of perfect gaz, and N the Avogadro number. In the following years Jean Perrin made a series of experiments leading to a new determination of the Avogadro number, and observed that the very irregular motion of particles resembled the nowhere differentiable functions of mathematicians. Norbert Wiener introduced what he called "the fundamental random function" as a mathematical model for the physical Brownian motion. It was called immediately "the Wiener process", and later on, following Paul Lévy, "the Brownian motion". Wiener gave several versions of the construction and derived a number of fundamental properties, Lévy developed the theory to a high point of sophistication, and it is now a mathematical object of common use as well as a mine of interesting problems.

Here is the theory as it appears from the last exposition made by Norbert Wiener. The problem is to construct a random process X(t, ω), also denoted by X(t) (= X(t, •)), (t the time, ω ∈ Ω the probability space) such that 1) for almost all ω X(t, ω) is a continuous function of t 2) X(t) is a Gaussian process, meaning that the distribution of any n-uple (X(t 1 ), X(t 2 ), . . . X(t n )) is Gaussian 3) this Gaussian process has stationary increments, meaning that the distribution of X(t) -X(s) depends on ts only 4) it satisfies a normalized Einstein equation, that is

||X(t) -X(s)|| 2 2 = |t -s|
where the norm is taken in L 2 (Ω).

Here is such a construction. Let H be an infinite-dimensional subspace of L 2 (Ω) consisting of Gaussian centered variables, and W an isometric linear mapping of L 2 (I) (I = IR, or IR + , or [0, 1]) into H. Let χ t be the indicator function

1 [0,t] . Then X(t) = W (χ t )
satisfies all conditions 2) to 4). Moreover, given an orthonormal basis of L 2 (I), (u n ), its image by W is a normal sequence (sequence of independent Gaussian normalized random variables (ξ n )), and expanding χ t in the form

χ t = Σ a n (t)u n (in L 2 (I))
results in an expansion of X(t) as a random series of functions :

χ(t) = Σ a n (t)ξ n (in L 2 (Ω))
or, more explicitly,

X(t, ω) = Σ a n (t)ξ n (ω) .
To prove condition 1), it is enough to establish that the series in the second member converges uniformly in t for almost all ω, and this is done rather easily when the u n are classical orthonormal bases. By definition, an helix is a curve in a Hilbert space, parametrized by IR, such that the distance between two points depends only on the distance of the parameters :

||X(t) -X(s)|| 2 2 = ψ(ts) , and ψ(•) is called the helix function. A translation of the parameter results in a isometric motion of the curve onto itself. It is the abstract model for all Gaussian processes with independent increments.

When ψ(t) = |t| we say that the curve is a Brownian helix. In contrast with the realizations of the Brownian motions (the functions t -→ X(t, ω) when ω is fixed), the Brownian helix is a very regular curve. However some basic properties of the Brownian motion can be read on the Brownian helix : its Hausdorff dimension is 2, its 2-dimensional Hausdorff measure is nothing but dt, and any three points on the curve are the vertices of a rightangle triangle : the increments starting from a point are orthogonal to the past (therefore, independent from the past).

Simple examples of helices are : 1) the line (ψ(t) = a 2 t 2 ) 2) the circle (ψ(t) = r 2 sin 2 ωt) 3) the three-dimensional helices (ψ(t) = a 2 t 2 + r 2 sin 2 ωt) 4) generalizations of those, with

ψ(t) = a 2 t 2 + sin 2 ωt µ(dω) ,
where µ is a positive measure on IR + such that the integral is finite. Actually this is the general form of an helix function.

Except when µ is carried on a finite set, the helix cannot be imbedded in a finite dimensional Euclidean space.

At the end of the 70s, Patrice Assouad developed a theory of Lipschitz embeddings of a metric space into another [START_REF] Assouad | Espaces métriques, plongements, facteurs, Thèse d'Etat Orsay[END_REF]. He introduced and built quasi-helices in Euclidean spaces, meaning that

0 < a < ||X(t) -X(s)|| 2 2 ψ(t -s) < b < ∞
for some a and b and all t and s. When ψ(t) = |t| we call them Brownian quasihelices. Assouad constructed Brownian quasi-helices in Euclidean IR n for n ≥ 3, and this gives a new way to prove that the realizations of Brownian motion are continuous a.s.. He asked me the question whether a and b can be taken near 1 when n is large, that is, whether the Brownian helix can be approximated (in this sense) by Brownian quasi-helices. I gave a positive answer with an explicit construction, and it was published in my paper on Helices and quasi-helices [START_REF] Kahane | Hélices et quasi-hélices[END_REF].

A construction of Brownian quasi-helices by means of Walsh matrices

Let us consider IR 2 n (n ≥ 1) as a Euclidean space. Let N = 2 n . If we want to construct a function X : IN -→ IR N such that X(0) = 0 and ||X(t) -X(s)|| 2 = |t -s| when |t -s| ≤ N we have to choose an orthonormal basis u 0 , u 1 , . . . u N -1 , define u N +j = u j , and write

X(t) = 0≤j≤t-1 ±u j .
At this stage there is no restriction on the signs ±, and we may choose + when 0 ≤ j ≤ N -1. If we try to obtain ||X(2t) -X(2s)|| 2 = 2|t -s|, ||X(4t) -X(4s)|| 2 = 4|t -s| etc when |t -s| ≤ N , whe have more and more conditions on the ± and we are led to the following construction. We define the Walsh matrix of order N as the N × N matrix obtained as the n th tensor power of the matrix

1 1 1 -1 , that is 1 1 1 -1 ⊗ 1 1 1 -1 ⊗ • • • ⊗ 1 1 1 -1 (n times) .
For example, the Walsh matrix of order 4 is

M = M 2 =     1 1 1 1 1 -1 1 -1 1 +1 -1 -1 1 -1 -1 1    
and the matrix M n+1 of order 2 n+1 is obtained from M n as

M n+1 = M n M n M n -M n
The N 2 first signs ± are those of the entries of the Walsh matrix, read line by line. In order to obtain the following signs, we extend the Walsh matrix by a series of vertical translations and change of signs of some lines according to the following rule : the first row is nothing but the whole sequence of entries, written from line to line and from left to right.

With this procedure we define X(t) when t is an integer and we can extend the construction to all t > 0, then to all real t. It is proved in [START_REF] Kahane | Hélices et quasi-hélices[END_REF] that we obtain a quasi-helix with a and b close to 1 when n is large enough : it is the answer to the question of Assouad.

However, it was not proved that the construction provides a quasi-helix when n = 2 (it was remarked that it gives a Peano curve in the plane when n = 1). The aim of the present paper is to give a detailed exposition of the case n = 2 (most of it could be copied for n > 2) and to prove that we obtain a quasi-helix. Instead of t ∈ IR we shall consider only t ∈ IR + and a curve starting from 0 (X(0) = 0). We shall investigate the geometric properties of the curve, some of them leading to open questions of a combinatorial or arithmetical nature.

The sequences that we construct are automatic in the sense of [START_REF] Allouche | Automatic sequences[END_REF] 3 Description of the sequence 3.1 It is a sequence of +1 and -1 as described before, in case N = 4. We write it as a succession of + and -:

+ + + + + -+ -+ + --+ --+ + + + + • • •
The gaps between the blocks of four letters have no meaning, except a help to understand the construction. The construction proceeds as follows : given the initial word of length 4 j , we divide it into four words A, B, C, D of equal length 4 j-1 and write it ABCD ; then

A B C D A (-B) C (-D) A B (-C)(-D) A(-B)(-C) D
is the initial word of length 4 j+1 . We shall give several equivalent definitions, using substitutions, explicit expressions, or generating functions.

Beforehand let us write the sequence in a tabular form as in the previous section :

+ + + + a 0 a 1 a 2 a 3 + -+ - A a
n = n 0 + 4n 1 + • • • + 4 ν n ν (n ν = 1, 2, 3; n j = 0, 1, 2, 3 if j < ν)
the construction shows that

a n = a n0+4n1 a m , m = n 1 + 4n 2 + • • • + 4 ν-1 n ν , that is a n = a n0+4n1 a n1+4n2 • • • a nν-1+4nν .
In the second member we find a j s with j ≤ 15. Their value is -1 when j = 5, 7, 10, 11, 13 and 14 and +1 otherwise. Now let us express a n as a function of n written in the 4-adic system of numeration. We obtain the formula

a n = (-1) An
where A n is the number of 11, 13, 22, 23, 31, 32 in the 4-adic expansion of n. For example, if n = 1 3 2 0 0 1 1 1 0 2 3 1 1 1 2 2 the significant links are

1 -3 -2 0 0 1 -1 -1 0 2 -3 -1 -1 -1 2 -2
A n is nine and a n = -1.

3.3

Let me describe the sequence by means of a substitution rule. We start from an alphabet made of eight letters : +a, +b, +c, +d, -a, -b, -c, -d. The substitution rule is

(S 0 ) +a -→ +a + b + c + d +b +a -b + c -d +c +a + b -c -d +d +a -b -c + d -a -a -b -c -d -b -a + b -c + d -c -a -b + c + d -d -a + b + c -d
The infinite word beginning with +a and invariant under the substitution is

W = +a + b + c + d + a -b + c -d + a + b -c -d + a -b -c + d + • • •
Replacing a, b, c, d by 1 (or, in a graphic way, in suppressing them), we obtain our sequence of ±1 (or ±).

3.4

Actually there is a simpler substitution rule leading to the same result, namely

(S 1 ) +a -→ +a + b +b +c + d +c +a -b +d +c -d -a -a -b -b -c -d -c -a + b -d -c + d
It can be checked immediately that (S 1 )(S 1 ) = (S 0 ).

3.5

The generating function of the sequence (a n ) is

f (z) = a 0 + a 1 z + a 2 z 2 + • • •
It can be defined using partial sums of order 4 n . Let us introduce the matrix

M (z) =     1 z z 2 z 3 1 -z z 2 -z 3 1 z -z 2 -z 3 1 -z -z 2 z 3    
and define four sequences of polynomials by the formulas

    P 0 Q 0 R 0 T 0     =     1 1 1 1         P n+1 Q n+1 R n+1 T n+1     = M (z 4 n )     P n Q n R n T n     (n = 0, 1, . . .). Then     P n Q n R n T n     = M (z 4 n-1 )M (z 4 n-2 ) • • • M (z 4 )M (z)     1 1 1 1     .
When |z| = 1, we have M (z)M (z) = 4I, therefore the matrix 1 2 M (z) is unitary, hence

|P n | 2 + |Q n | 2 + |R n | 2 + |T n | 2 = 4(|P n-1 | 2 + |Q n-1 | 2 + |R n-1 | 2 + |T n-1 | 2 ) = 4 n (1 + 1 + 1 + 1) = 4 n+1
We obtain the generating function as

f (z) = lim n→∞ P n (z)
We can write the generating function in a more interesting form :

f (z) = f 0 (z 4 ) + zf 1 (z 4 ) + z 2 f 2 (z 4 ) + z 3 f 3 (z 4 ) ,
where the coefficients of the power series f 0 , f 1 , f 2 , f 3 are the columns of the table in 3.0. In order to obtain these coefficients, we can start from W in 3.2 and replace a, b, c, d by 1, 1, 1, 1 (for

f 0 ), 1, -1, 1, -1 (for f 1 ), 1, 1, -1, -1 (for f 2 ) and 1, -1, -1, 1 (for f 3 ). Then f 0 = f . Writing F (z) =     f 0 (z) f 1 (z) f 2 (z) f 3 (z)    
the functional equation of the generating functions of the columns is

F (z) = M (z)F (z 4 ) .
4 Description of the curve 4.1 Let u 0 , u 1 , u 2 , u 3 be an orthonormal basis of the Euclidean space IR 4 , and define u j+4 = u j (j = 0, 1, . . .). The partial sums of the series

a 0 u 0 + a 1 u 1 + a 2 u 2 + • • •
(that can be obtained from W in 3.2 by replacing a, b, c, d by u 0 , u 1 , u 2 , u 3 ) will be denoted by S(n). Then

S(n) = a 0 u 0 + a 1 u 1 + • • • + a n-1 u n-1 ∈ Z Z 4 .
It is easy to check on the table in 3.0 that

S(16n) = 4S(n) (n ∈ IN) .
Moreover it is not difficult to see (we shall be more specific later) that

||S(n) -S(m)|| 2 ≤ b |n -m|
for some b < ∞ and all n and m. This allows, first, to define S(t) when t is a binary number via a formula S(16 ν t) = 4 ν S(t), then to check that

S(16t) = 4S(t) ||S(t) -S(s)|| 2 ≤ b|t -s|
for such numbers, then to extend S(•) by continuity on IR + , and check the above formulas for all t ≥ 0 and s ≥ 0.

The curve we consider is the image of IR + by S(•).

Clearly (changing t into 16t) the curve is invariant under an homothety of center 0 and ratio 4. Our main purpose is to prove that it is a Brownian quasihelix. We shall point out some geometric properties first.

The matrix

M transforms u 0 , u 1 , u 2 , u 3 into u ′ 0 , u ′ 1 , u ′ 2 , u ′ 3 : (u ′ 0 , u ′ 1 , u ′ 2 , u ′ 3 ) = M (u 0 , u 1 , u 2 , u 3
) and the partial sums of order n of the series Σa j u ′ j are the partial sums of order 4n of the series Σa j u j . Therefore the equation

S(4t) = M S(t)
holds true when t = n ∈ IN and by extension for all t ∈ IR + .

It is easy to check that the eigenvalues of M are 2 and -2, and that

M     1 1 1 0 1 0 -1 1     =     2 2 2 0 2 0 -2 2     , M     1 0 -1 1 -1 -1 -1 0     =     -2 0 2 -2 2 2 2 0     . Let v 0 = 1 2 (u 0 + u 1 + u 2 -u 3 ) , v 1 = √ 2 2 (u 0 + u 3 ) v 2 = 1 2 (u 0 -u 1 -u 2 -u 3 ) , v 3 = √ 2 2 (u 1 -u 2 )
. They constitute an orthonormal system. The vectors v 0 and v 1 generate a plane, P , which is the eigenspace of the eigenvalue 2, and v 2 and v 3 a plane, Q, corresponding to the eigenvalue -2. Expressed via the orthonormal basis v 0 , v 1 , v 2 , v 3 , the operator M takes the form

M ′ = U M U -1 = 2     1 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 -1     ,
where U is the unitary matrix carrying u 0 , u 1 , u 2 , u 3 onto v 0 , v 1 , v 2 , v 3 . It means that the transformation S(t) -→ S(4t) is the product of a homothety of centre 0 and ratio 2 and an orthogonal symmetry with respect to the plane Q.

4.3

Clearly M 2 = 4I (I being the identity matrix), In turn, M is the square of another simple matrix, as it can be guessed from 3.2 and 3.3. Let us define

T =     1 0 1 0 1 0 -1 0 0 1 0 1 0 1 0 -1     . Then M = T 2 .
The eigenvalues of

T are √ 2, - √ 2, i √ 2 and -i √ 2.
The vectors

w 0 = √ 2 2 (v 0 + v 1 ) , w 1 = √ 2 2 (v 0 -v 1 )
are eigenvectors corresponding to √ 2 and -√ 2. Defining w 2 and w 3 in such a way that w 0 , w 1 , w 2 , w 3 is a direct orthonormal basis, and W being the unitary matrix carrying u 0 , u 1 , u 2 , u 3 onto w 0 , w 1 , w 2 , w 3 , we can write

W T W -1 = T ′ = √ 2     1 0 0 0 0 -1 0 0 0 0 0 -1 0 0 1 0     .
It means that T ′ is decomposed into : 1) an homothety of centre 0 and ratio √ 2 2) a rotation of π 2 of the orthogonal projection on Q 3) a symmetry with respect to w 1 of the orthogonal projection on P .

In the same way as we obtained the equation S(4t) = M S(t), we now have

S(2t) = T S(t)
and we have just given the interpretation of the transformation S(t) -→ S(2t) as a product of simple transformations.

4.4

We have investigated the properties of the transformations S(t) -→ S(16t), S(t) -→ S(4t), S(t) -→ S(2t) as products of homotheties and isometries. Now we shall look at the effect of a translation of t by an integer. We are interested in differences S(t) -S(s).

Let us begin with integers m < n < 16 k . Let us divide the series a 0 u 0 + a 1 u 1 + • • • into consecutive blocks of length 16 k , so that the series reads

+A + B + C + D + A -B + C -D + • • • If j = j 0 + 4j 1 (j 0 = 0, 1, 2, 3, j 1 ∈ IN), the j-th term is of type A, B, C, D
according to the value of j 0 and its sign is a j . Therefore

S(n + j • 16 k ) -S(m + j • 16 k ) = a j (S(n + j 0 16 k ) -S(m + j 0 16 k )) .
If 0 < s < t < 1, we can approximate s and t by m 16 -k and n 16 -k and we obtain S(t + j) -S(s + j) = a j (S(t + j 0 ) -S(s + j 0 )) (j 0 = 1, 2, 3, 4, j 0 = j modulo 4). This expresses that all arcs A j = S([j, j+1]) are isometric (actually translates or symmetric according to the value of a j ) of one of the arcs A 0 , A 1 , A 2 , A 3 (according to the value of j 0 ). Using 4.2, this holds true when we replace the A j s by A ν j = S([j2 ν , (j + 1)2 ν ]) whatever ν ∈ Z Z.

5 It is a Brownian helix 5.1 What we have to prove is that, writing

a = inf 0<s<t ||S(t) -S(s)|| |t -s| ≤ sup 0<s<t ||S(t) -S(s)|| |t -s| = b , we have a > 0 , b < ∞ .
We can write as well

a = inf m<n ||S(n) -S(m)|| |n -m| , b = sup m<n ||S(n) -S(m)|| |n -m| 5.2 The easy part is b < ∞. Let us first assume [m, n] = [j2 k , (j + 1)2 k ]. Then, according to 4.3, ||S(n) -S(m)|| = 2 k/2 .
In the general case, let us decompose [m, n] into such intervals in a minimal way, so that there are at most two intervals of the same length in the decomposition.

If the largest length is 2 k , we obtain

||S(n) -S(m)|| ≤ 2(2 k/2 + 2 (k-1)/2 + • • •) ≤ 2(1 + √ 2)2 k/2 therefore ||S(n) -S(m)|| ≤ 2(1 + √ 2)|n -m| 1/2 . This gives b ≤ 2(1 + √ 2)
5.3 To prove a > 0 is more tricky. We shall use two lemmas.

Lemma 1. There exists α > 0 such that

||S(n + h) -S(n)|| ≤ 1 -α for all n ∈ IN and h ∈ [-1 2 , 1 2 ] (with n + h ≥ 0).
Lemma 2. There exists an integer A such that

||S(n) -S(m)|| ≥ 2
for all integers m and n such that nm ≥ A.

Assuming that this is correct, the result is at hand : given t and s such that ts ≥ A + 1, we can write s = m + h and t

= n + h ′ with h, h ′ ∈ [-1 2 , 1 2 ] and m -n ≥ A, therefore ||S(t) -S(s)|| ≥ 2α . Whenever (A + 1)2 k < t -s ≤ (A + 1)2 k+1 (k ∈ Z Z), we have ||S(t) -S(s)|| ≥ 2α2 k/2 ≥ 2α 2(A + 1) |t -s| 1/2 therefore a ≥ 2α
2(A + 1) .

Proof of Lemma 1.

From now on it may be useful to represent S(n) on the table of 3.0, and also the differences S(n) -S(m), as figures consisting of consecutive lines plus or minus part of a line above and below, in such a way that each column in the figure has a sum equal to the corresponding coordinate of S(n) or S(n) -S(m).

+ + + + + -+ - + + -- + --+ + S (

Proof of Lemma 2

Here again we look at the table. We can compute ||S(n) -S(m)|| when nm is in a given interval, and moreover give the couples (m, n) for each the infimum is attained, and the expression of S(n) -S(m) (that is, the coordinates with respect to u 0 , u 1 , u 2 , u 3 ).

S(42) -S(22) = u 1 + u 2 -u 3 -u 4 , S(234) -S(214) = -u 1 -u 2 -u 3 + u 4 .
This proves Lemma 2 with A = 16.

Actually the proof can be given in a more concentrated form. It is enough to show that ||S(n) -S(m)|| 2 ≤ 3 is impossible when nm ≥ 16. Let us assume ab absurdo that ||S(n) -S(m)|| 2 ≤ 3. Let us add or remove the minimal number of terms in order to transform S(n) -S(m) into a difference of the form S(4n ′ ) -S(4m ′ ) (that is, to transform the figure S(n)-S(m) into a rectangle). In general, this minimal number is ≤ 4 and has the same parity as ||(S(n) -S(m)|| 2 ; here it is ≤ 3 and the resulting S(4n ′ ) -S(4m ′ ) has its squared norm ≤ 2, therefore ||S(n ′ ) -S(m ′ )|| 2 ≤ 3 and the process goes on until we reach S(n × ) -S(m × ) with n × ≤ 64. Then we know the possible pairs (m × , n × ), namely (5,11), (23,29), (35,41) and (53,59), and the reverse process never gives a squared norm ≤ 3.

Remarks and questions

The estimates we gave for b and a are quite rough. We can ask for better estimates and conjectures. The actual problem, of a combinatorial or arithmetical nature, is to compute these numbers exactly.

We were interested in estimating b from above and a from below. Examples provide estimates in the opposite direction : 6 Projections of the curve 6.1 The direction of u 0 is special : all first coordinates of the S(n) are ≥ 0. That means that the partial sums S 0 (n) of the original series described in 3.0 are positive.

A simple way to see it is to use Lemma 1 (of 5.2). Since S 0 (n) ≥ 1 for n = 1, 2, 3, 4, 5, 6, 7, 8, we have S 0 (t) > 0 for 1 2 ≤ t ≤ 8, therefore (changing t into 16 k t), S 0 (t) > 0 for all t > 0.

6.2 We are mainly interested in the three-dimensional projections of the curve. It seems likely that all parallel projections of the curve on a three-dimensional subspace of IR 4 have an infinity of double points. The question can be formulated in the equivalent forms : 1) is every direction in IR 4 the direction of some S(t) -S(s) ?

2) are the S(n)-S(m) √ n-m (n > m) dense on the sphere S 3 ? 6.3 Let us project the curve from 0 on the sphere S 3 , that is consider

S(t) = S(t) ||S(t)|| (t > 0) .
We obtain a closed C, image of any interval [a, 16a] by S(•).

C is invariant under the isometries of IR 4 defined by with respect to the orthonormal basis (w 0 , w 1 , w 2 , w 3 ). The vectors v 0 and v 1 (w 0 and w 1 as well) generate a plane P such that the mapping S(t) -→ S(4t) is an orthogonal symmetry with respect to P . For the projection of C on P , the change of t into 2t means a symmetry with respect to the line generated by w 0 .

C has a double point at t = 1 3 , t ′ = 4 3 : T ( 4 3 ) = T ( 1 3 ). In order to prove it, we expand t and t ′ in base 4 (we underline the expansion)

t = 0.1111 • • • t ′ = 1.1111 • • • .

  Séminaire sur l'interface entre l'analyse harmonique et la théorie des nombres Luminy 21 octobre 2005 J.-P. Kahane Université Paris-Sud, Laboratoire de Mathématique, Bâtiment 425, F-91405 Orsay Cedex, France 1 History and definitions

1 ) 4 1 -

 141 fig.1

  seems not impossible that the estimate for α is precise, that is a = 1 √ 5 . A careful investigation of the table would confirm or disprove this conjecture. It would lead also to a better estimate for b.

1 2 M and 1 √ 2 T

 12 (see 4.1 and 4.2). The first takes the form the orthonormal basis (v 0 , v 1 , v 2 , v 3 ) and the second 1

  4 a 5 ...........

	+ + --+ --+ + + + +		.................... ............. a 15 a 16 .............
	-+ -+ + + ---+ + -+ + + +	B	.................... .................... ............. a 31 a 32 .............
	+ -+ ---+ + -+ + -+ + + +	C	.................... .................... ............ a 47 a

48 ............. -+ -+ D .................... --+ + .................... + --+ ............ a 63 a 64 ............. A .................... .................... ............ a 79 a 80 ............. -B .................... .................... .................... 3.2 Let me give an explicit expression for a n . Writing

Using base 4 again, we easily obtain the figures and the values of S(1), S(11), S(111) and so on : S(1) = u 0 = (1 0 0 0)

S(11)

= S(10) + (S(11) -S(10)) = ( 1 The ratio between two consecutive vectors tend to 2 (meaning that the ratios of coordinates tend to 2), hence

.

By isometry we also have

These double points are contained in the plane P , and they are symmetric with respect to the line generated by w 0 .

I believe, but did not prove, that these are the only multiple points of the curve C. In that case, C is a Brownian quasi-helix (actually, a Brownian quasicircle) on some 4-covering of the sphere S 3 .

One can see the curve C in two other ways

First, taking into account that the first coordinate S 0 (t) is always positive, we can consider

and the curve C ′ described by R(•), projection of the original curve with a source at 0 and a screen at the hyperplane x 0 = 1. Symmetries and double points can be studied on this model as well.

Secondly, we obtain a projective model of C, say, C ′′ , on choosing four points A 0 , A 1 , A 2 , A 3 in IR 4 , defining A j+4 = A j (j = 0, 1, . . .), starting with a point M 0 = A 0 and defining the sequence of points

Some real figures would help. If a reader is willing to draw figures of the above curves, I'll appreciate to see them.