Jean-François Semblat
email: semblat@lcpc.fr

Roberto Paolucci

Anne-Marie Duval

Simplified vibratory characterization of alluvial basins

Keywords: earthquake, site effects, eigenmodes, vibration, resonance, amplification, seismic wave séisme, effets de site, modes propres, vibration, résonance, amplification, onde sismique

For the analysis of seismic wave amplification, modal methods are interesting tools to study the modal properties of geological structures. Modal approaches mainly lead to information on such parameters as fundamental frequencies and eigenmodes of alluvial basins. For a specific alluvial deposit in Nice (France), a simplified modal approach involving the Rayleigh method is considered. This approach assumes a set of admissible shape functions for the eigenmodes and allows a fast estimation of the fundamental frequency of the basin. The agreement between modal numerical results and experimental ones is satisfactory. The simplified modal method then appears as an efficient mean for the global vibratory characterization of geological structures towards resonance.

Caractérisation vibratoire simplifiée de remplissages sédimentaires

Résumé.

L'analyse de l'amplification du mouvement sismique peut être réalisée en étudiant les propriétés vibratoires des structures géologiques. Les approches modales peuvent donner des informations telles que les fréquences propres et les modes propres de bassins sédimentaires. Pour étudier la résonance d'un bassin sédimentaire du centre de

Nice, une approche modale simplifiée utilisant la méthode de Rayleigh est proposée. Cette approche suppose la connaissance d'une famille de fonctions modales admissibles et permet une détermination rapide de la fréquence fondamentale du bassin. La comparaison avec d'autres résultats numériques (obtenus par équations intégrales de frontière) et des mesures réalisées sur site est très favorable. Cette méthode conduit donc aisément et rapidement à une caractérisation vibratoire globale de structures géologiques vis-à-vis de la résonance.

Version française abrégée

L'amplification locale du mouvement sismique dans les bassins sédimentaires peut parfois être très importante [START_REF] Bard | The two dimensional resonance of sediment filled valleys[END_REF][START_REF] Bard | Effects of surface geology on ground motion : recent results and remaining issues[END_REF][START_REF] Duval | Détermination de la réponse sismique d'un site à l'aide du bruit de fond[END_REF][START_REF] Pitilakis | Site effects : recent considerations and design provisions[END_REF][START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF]. La propagation des ondes sismiques dans les structures géologiques de surface peut en effet conduire à une forte amplification du mouvement due au contraste de propriétés mécaniques entre un bassin sédimentaire et le substratum rocheux. L'amplification du mouvement sismique peut également s'interpréter comme la résonance des couches de surface à certaines fréquences. Les effets de site peuvent alors être analysés en termes de résonance vibratoire des structures géologiques de surface [START_REF] Guéguen | From soil-structure to site-city interaction[END_REF][START_REF] Semblat | Quantitative estimation of seismic site effects through modal methods[END_REF]. Dans cette note, nous considérons une approche modale simplifiée basée sur la méthode de Rayleigh [START_REF] Paolucci | Fundamental vibration frequencies of 2D geological structures[END_REF] pour estimer directement la fréquence fondamentale d'une structure géologique.

La méthode modale considérée ici vise à estimer la fréquence « fondamentale » d'une structure géologique. Le processus de propagation d'ondes conduit généralement à une forte amplification du mouvement à des fréquences variées [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF] mais, dans un but pratique, il est intéressant de déterminer une fréquence fondamentale unique par des approches simplifiées. La méthode de Rayleigh peut en effet permettre une estimation rapide et fiable de la fréquence de résonance d'une structure géologique [START_REF] Paolucci | Fundamental vibration frequencies of 2D geological structures[END_REF].

Simplified vibratory characterization of alluvial basins

J-F Semblat, R. Paolucci, A-M Duval
Cette méthode considère que le déplacement d'un système élastique dans un de ses modes propres peut être approché par celui d'un système à un degré de liberté. Nous étudions ici le premier mode propre caractérisé par la fréquence :   =2f 0 . En notant V l'énergie élastique du système et T son énergie cinétique, la conservation de l'énergie totale d'un système élastique implique que V max =T max . On peut alors écrire le déplacement s k ( x ,t) correspondant aux vibrations harmoniques à la fréquence   . La valeur exacte de la fréquence fondamentale du système peut ainsi être obtenue lorsque la forme modale réelle  k x ( ) est connue. Cependant, comme la solution exacte n'est généralement pas accessible, la valeur de   peut être estimée correctement en considérant une approximation réaliste de  k x ( ) . Cette approximation doit satisfaire à la fois les compatibilités géométriques et les conditions aux limites. Il a néanmoins été démontré [START_REF] Paolucci | Fundamental vibration frequencies of 2D geological structures[END_REF] que la seconde condition peut ne pas être satisfaite complètement et que la forme modale peut être choisie dans une large gamme de fonctions satisfaisant uniquement les compatibilités géométriques. On peut alors déterminer la fréquence fondamentale à l'aide de la relation suivante :

       d x x d x x x jl jl o k ) ( ) ( ) ( ) ( min 2 2       (1) où ij
 est le tenseur de contrainte, ij  le tenseur de déformation,  la masse volumique et

) ( x k  la forme modale.
En vue de l'analyse vibratoire simplifiée, nous considérons un bassin sédimentaire situé dans le centre de Nice. Différents types de mesures (bruit de fond, mouvements sismiques faibles) ont permis d'analyser l'amplification des ondes sismiques dans le bassin principal [START_REF] Duval | Détermination de la réponse sismique d'un site à l'aide du bruit de fond[END_REF]. Pour évaluer la robustesse de la méthode, le bassin est supposé homogène et sa résonance vibratoire est analysée pour une sollicitation sismique antiplane en cisaillement caractérisée par la célérité V S    / . Les caractéristiques mécaniques des deux milieux (bassin et substratum) sont les suivantes :

 bassin :  1 =2000 kg/m 3 ,  1 =180 MPa, =0.2 c'est-à-dire V S =300 m/s ;  substratum :  1 =2300 kg/m 3 ,  1 =4500 MPa, =0.2 d'où V S =1400 m/s.
où  est la masse volumique,  le module de cisailement,  le coefficient de Poisson et V S la célérité des ondes de cisaillement.

Comme le montre la figure 1, l'interface entre le bassin et le substratum est décrite à l'aide de deux fonctions cosinus (parties ouest et est de la vallée). A partir de formes modales admissibles choisies judicieusement, on estime les fréquences fondamentales de ces deux parties du bassin sédimentaire. Les résultats sont obtenus pour le premier mode en considérant différentes valeurs du module de cisaillement :

 1 =180 MPa (i.e. V S =300 m/s),  2 =120 MPa (i.e. V S =245 m/s) et  3 =90 MPa (i.e. V S =212 m/s).
Les résultats numériques obtenus par la méthode de Rayleigh sont comparés aux fréquences d'amplification maximale obtenues expérimentalement et numériquement (méthode des éléments de frontière). D'après ces résultats, nous pouvons constater que les fréquences fondamentales obtenues par la méthode de Rayleigh sont en bon accord avec les fréquences de forte amplification obtenues expérimentalement et par approche propagative. Les valeurs données dans le tableau I sont comparées pour la résonance de la partie profonde du bassin (ouest) et la partie peu profonde (est) et ce pour différentes valeurs du module de cisaillement. Pour un modèle de bassin homogène, la méthode de Rayleigh semble donc donner des résultats très intéressants qui sont non seulement en accord avec les valeurs obtenues à l'aide d'autres modèles mais également avec les résultats expérimentaux (figure 2).

Introduction

The local amplification of seismic motion in alluvial basins could sometimes be large [START_REF] Bard | The two dimensional resonance of sediment filled valleys[END_REF][START_REF] Bard | Effects of surface geology on ground motion : recent results and remaining issues[END_REF][START_REF] Duval | Détermination de la réponse sismique d'un site à l'aide du bruit de fond[END_REF][START_REF] Pitilakis | Site effects : recent considerations and design provisions[END_REF][START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF]. Seismic wave propagation in surface geological structures can indeed lead to a strong motion amplification due to the impedance constrast between the alluvial deposit and the bedrock. However, the amplification of seismic motion can also be considered as a result of the resonance of surface layers at peculiar frequencies. Site effects can then be analysed in terms of vibratory resonance of surface geological structures [START_REF] Guéguen | From soil-structure to site-city interaction[END_REF][START_REF] Semblat | Quantitative estimation of seismic site effects through modal methods[END_REF]. In this note, we consider a simplified modal approach based on the Rayleigh's method [START_REF] Paolucci | Fundamental vibration frequencies of 2D geological structures[END_REF] for the direct estimation of the fundamental frequency of a geological structure.

The modal method considered herein aims at estimating the « fundamental » frequency of an alluvial basin. The wave propagation process generally leads to strong motion amplifications at various frequencies [START_REF] Semblat | Numerical analysis of seismic wave amplification in Nice (France) and comparisons with experiments[END_REF] but, from a practical point of view, it is interesting to estimate a unique fundamental frequency by a simplified efficient method. We will show herein that Rayleigh's method allows a fast and reliable estimation of the resonance frequency of a geological structure [START_REF] Paolucci | Fundamental vibration frequencies of 2D geological structures[END_REF].

Simplified modal method

Theoretical point of view

The modal method considered herein vise à estimer la fréquence « fondamentale » d'une structure géologique. Le processus de propagation d'ondes conduit généralement à une forte amplification du mouvement à des fréquences variées (Semblat et al 00b) mais, dans un but pratique, il est intéressant de déterminer une fréquence fondamentale unique par des approches simplifiées. La méthode de Rayleigh permet par exemple une estimation rapide et fiable de la fréquence de résonance d'une structure géologique (Paolucci 99).

Cette méthode considère que le déplacement d'un système élastique dans un de ses modes propres peut être approché par celui d'un système à un degré de liberté. Nous étudions ici le premier mode propre caractérisé par la fréquence :   =2f 0 . En notant V l'énergie élastique du système et T son énergie cinétique, la conservation de l'énergie totale d'un système élastique implique que V max =T max .

Le déplacement s k ( x ,t) correspondant aux vibrations harmoniques à la fréquence   peut s'écrire :

s x t x e k k i t o ( , ) ( )    (2)
où x représente la coordonnée d'espace, i le nombre imaginaire, t le temps et  k x ( ) la forme modale suivant la direction k. On calcule alors l'énergie cinétique du système comme suit :

T t x s t d k ( ) ( )         1 2 2             0 2 2 2 0 1 2 e x x d i t k ( ) ( )   (3) d'où       d x x t T T k o t ) ( ) ( 2
1 ) ( max 2 2 max    (4) L'énergie élastique V s'écrit alors :     d x x t V jl jl ) ( ) ( 2

Il vient ainsi :

Simplified vibratory characterization of alluvial basins

J-F Semblat, R. Paolucci, A-M Duval        d x x d x x k jl jl o ) ( ) ( ) ( ) ( 2 2      (6)
L'équation 6 donne la valeur exacte de la fréquence fondamentale du système lorsque la forme modale réelle  k x ( ) est connue. Cependant, comme la solution exacte n'est généralement pas accessible, la valeur de   peut être estimée correctement en considérant une approximation réaliste de  k x ( ) . Cette approximation doit satisfaire à la fois les compatibilités géométriques et les conditions aux limites. Il a néanmoins été démontré (Paolucci 99) que la seconde condition peut ne pas être satisfaite complètement et que la forme modale peut être choisie dans une large gamme de fonctions satisfaisant uniquement les compatibilités géométriques. On peut alors déterminer la fréquence fondamentale à l'aide de la relation suivante :

       d x x d x x x jl jl o k ) ( ) ( ) ( ) ( min 2 2       (7)

Estimation of the fundamental frequency of a geological structure

En vue de l'analyse vibratoire simplifiée, nous considérons un bassin sédimentaire situé dans le centre de Nice. Ce site est bien connu puisque de nombreuses expérimentations ont été réalisées ainsi que plusieurs modélisations numériques (Semblat et al 00b). Dans ce texte, le bassin est supposé homogène et sa résonance vibratoire est analysée pour une sollicitation sismique antiplane en cisaillement caractérisée par la célérité V S    / . Les caractéristiques mécaniques des deux milieux (bassin et substratum) sont les suivantes :

 bassin :  1 =2000 kg/m 3 ,  1 =180 MPa, =0.2 c'est-à-dire V S =300 m/s ;  substratum :  1 =2300 kg/m 3 ,  1 =4500 MPa, =0.2 d'où V S =1400 m/s.
où  est la masse volumique,  le module de cisailement,  le coefficient de Poisson et V S la célérité des ondes de cisaillement.

Comme le montre la figure 1, l'interface entre le bassin et le substratum est décrite à l'aide de deux fonctions cosinus :

 partie ouest de la vallée : f(x,z) = (h 1 +1).cos(2,7.10 -3 x+1,55)  partie est de la vallée : g(x,z) = (h 2 +2).cos(2,8.10 -3 x -1,3)

Parmi les formes modales admissibles, on choisit les suivantes (Paolucci 99) :

      2 2 1 2 1 1 2 1 2 1 2 ( , ) cos ( , ) sin ( ) cos (
) x z f x z n x a m z h r s t                                 (8)
où f(x,z) est la fonction donnée précédemment. Les résultats sont obtenus pour le premier mode en considérant différentes valeurs du module de cisaillement :  1 =180 MPa (i.e. V S =300 m/s),  2 =120 MPa (i.e. V S =245 m/s) et  3 =90 MPa (i.e.

V S =212 m/s). Les résultats numériques obtenus par la méthode de Rayleigh sont comparés aux fréquences d'amplification maximale données par la méthode des éléments de frontière à travers une analyse explicite de la propagation des ondes (Semblat et

Comparison between simplified modal method and amplification curves

Conclusion

L'analyse des effets de site par des approches modales est possible tant du point de vue qualitatif que quantitatif. Le site considéré dans cette note est un remplissage sédimentaire du centre de Nice pour lequel les mesures expérimentales prévoient des facteurs d'amplification variant entre 10 et 30 pour des fréquences comprises entre 1 et 2 Hz.

Les approches modales permettent de déterminer différentes fréquences propres d'un remplissage sédimentaire. Pour une application concrète de ces approches modales, quelques points devront encore être approfondis : la distribution de vitesse dans le remplissage, l'effet de l'amortissement...

Figure 1 .

 1 Figure 1. Homogeneous geological profile for the simplified modal analysis.Description du profil géologique homogène pour l'approche modale simplifiée.

Figure 2 .

 2 Figure 2. Fundamental frequencies from the simplified modal approach (vertical lines) compared with experimental (top) and propagative BEM results (bottom). Fréquences fondamentales estimée par méthode modale simplifiée et comparaison avec les mesures (haut) et les simulations par la méthode des équations intégrales de frontière (bas).

Modal estimation of the fundamental frequency 3.1 Dependency of fundamental frequencies on shear modulus D

  'après ces résultats, nous pouvons constater que les fréquences fondamentales obtenues par la méthode de Rayleigh sont en bon accord avec les fréquences de forte amplification déterminées par approche propagative. Les valeurs données dans le tableau I sont comparées pour la résonance de la partie profonde du bassin (ouest) et la partie peu profonde (est) et ce pour différentes valeurs du module de cisaillement. Pour un modèle de bassin homogène, la méthode de Rayleigh semble donc Duval donner des résultats très intéressants qui sont non seulement en accord avec les valeurs obtenues à l'aide d'autres modèles mais également avec les résultats expérimentaux présentés en détail dans (Semblat et al 00b).

			al 00b).
	TABLEAU I : Comparaisons entre fréquences de référence obtenues par méthode propagative et
	fréquences fondamentales issues de l'approche modale simplifiée.
		partie ouest		partie est
	fréquence	fréquence	fréquence	fréquence
	référence	fondamentale	référence	fondamentale
	1,35	1,50	2,42	2,86
	1,30	1,23	2,13	2,34
	1,13	1,07	1,75	2,02
	3.

Simplified vibratory characterization of alluvial basins

J-F Semblat, R. Paolucci, A-M


Ferruh Özbudak
email: ozbudak@metu.edu.tr

Zülfükar Saygı
email: zsaygi@etu.edu.tr

Counting quadratic forms of codimension 2 in characteristic 2 and relations to maximal curves

Keywords:

Let Fq be an arbitrary finite field of characteristic 2 and k ≥ 2 be an arbitrary integer. We count the number of quadratic forms of codimension 2 on F q k over Fq and we give some related results, including some relations to certain maximal curves over finite fields.

Introduction

Let q = 2 t with t ≥ 1 (a power of 2), k ≥ 2 be an integer and m = ⌊k/2⌋, the integer part of m/2. For a positive integer ℓ, let F 2 ℓ denote a finite field with 2 ℓ elements (of characteristic 2). Let Q : F q k → F q be an arbitrary quadratic form. We assume without loss of generality (see Section 2 below) that there exist uniquely determined ǫ 0 , ǫ 1 , . . . , ǫ m-1 ∈ F q k and some ǫ m ∈ F q k such that for x ∈ F q k we have

Q(x) = Tr x ǫ 0 x + ǫ 1 x q + • • • + ǫ m x q m . (1)
Here ǫ m is uniquely determined if k is odd and ǫ m is determined uniquely only modulo F q m if k is even. Throughout the paper Tr denotes the trace map from F q k onto F q , while Tr Fq/F2 denotes the trace map from F q onto F 2 .

In this paper we count the number of quadratic forms on F q k over F q of codimension 2. We also count the subset of them with the invariant Λ(Q) = 1. Our results depend on the function S(ǫ) we introduce in Definition 1 below. We also give some related results on the maximal curves of the form

χ : y q + y = x ǫ 0 x + ǫ 1 x q + • • • + ǫ m x q m .
Our results generalize some of the results of Fitzgerald in [START_REF] Fitzgerald | Highly Degenerate Quadratic Forms over F2[END_REF].

Background

In this section we recall basic definitions and some facts that we use in this paper. A quadratic form on F q k over F q is a map such that -Q(αx) = α 2 Q(x) for all α ∈ F q and x ∈ F q k , and -The related map B(x, y) on F q k × F q k defined by

B(x, y) = Q(x + y) + Q(x) + Q(y) is a bilinear map over F q .
The radical W of Q is an F q -linear subspace of F q k given by

W = x ∈ F q k : B(x, y) = 0 for all y ∈ F q k .
Let w be the F q -dimension of W . Its codimension k -w is always even. By codimension of Q we mean the difference k -w.

There is another invariant Λ(Q) of Q, which is an integer in the set {-1, 0, 1}.

In fact if w = 0, then Q is called non-degenerate and Λ(Q) is the Arf invariant of Q [1], which is in the subset {-1, 1}. Let N (Q) denote the cardinality N (Q) = x ∈ F q k : Q(x) = 0 .
The invariant Λ(Q) maybe defined as the integer satisfying

N (Q) = q k-1 + Λ(Q)(q -1)q k+w 2 -1 (2)
(see, for example, [5, Section 6.2]).

The following is a useful characterization of quadratic forms, which allows us to make the assumption on Q in the beginning of Section 1.

Theorem 1 (see [START_REF] Fitzgerald | Highly Degenerate Quadratic Forms over finite fields of characteristic 2[END_REF]). Let Q : F q k → F q be a quadratic form and m = ⌊k/2⌋. There exist ǫ 0 , ǫ 1 , . . . , ǫ m ∈ F q k such that

Q(x) = Tr x ǫ 0 x + ǫ 1 x q + • • • + ǫ m x q m . ( 3
)
Moreover ǫ 0 , ǫ 1 , . . . , ǫ m are uniquely determined, except when k is even in which case ǫ m is only unique modulo

F q m .
If the codimension of the radical is 2, then we have further information on Q.

Theorem 2 (see [START_REF] Fitzgerald | Highly Degenerate Quadratic Forms over finite fields of characteristic 2[END_REF]). Let ǫ 0 , ǫ 1 , . . . , ǫ m be the coefficients corresponding to Q as in (3). Then we have w = k -2 if and only if there exist a, b ∈ F q k such that the set {a, b} is linearly independent over F q and for 0 ≤ i ≤ ⌊(k -1)/2⌋ we have

ǫ i = a q i b + ab q i , (4)
and if k is even, then furthermore

ǫ m -ab q m ∈ F q m . ( 5
)
Moreover we have the following:

-If Λ(Q) = 1, then ǫ 0 = ab. -If Λ(Q) = -1
, then there exists s ∈ F q such that Tr Fq/F2 (s) = 1 and

ǫ 0 = a 2 + ab + sb 2 .
-If Λ(Q) = 0, then there exists c ∈ F q k such that the set {a, b, c} is linearly independent over F q and

ǫ 0 = c 2 + ab.
The following result will be used in a proof below.

Lemma 1. Let F be a finite field of characteristic 2. Let x, y ∈ F and t be a positive integer. We put

u = x + y and v = xy.
Then the following holds:

x 2 t +1 +y 2 t +1 = u 2 t +1 + vu 2 t +1-2 +v 2 u 2 t +1-2 2 +v 2 2 u 2 t +1-2 3 +• • •+v 2 t-1 u 2 t +1-2 t .

Main Results

In this section we present our main results. First we give an "if and only if" characterization of quadratic forms of codimension 2. The following theorem will be useful in our counting results.

Theorem 3. Q has a codimension 2 radical if and only if there exists a nonzero v ∈ F * q k such that each of the following holds: i) ǫ 1 = 0, and we have

v = ab and ǫ 1 v = a q-1 + b q-1 , for some a, b ∈ F q k . ii) v q ǫ 2 = ǫ q+1 1 1 + t-1 i=0 v q+1 ǫ -2 1 2 i . iii) For all 3 ≤ i ≤ k-1 2
, we have

ǫ i = ǫ q-1 1 q i-2 ǫ i-2 + ǫ 2 ǫ -1 1 q i-2 ǫ i-1 . iv) If k = 2m then ǫ m -ab q m ∈ F q m ,
where a, b are the elements in F q k satisfying i).

Proof. (⇒) Assume that Q has a codimension 2 radical, then using Theorem 2 we know that there exist a, b ∈ F q k such that the set {a, b} is linearly independent over F q and (4) holds. Set u = a q-1 + b q-1 and v = ab. As a, b are linearly independent over F q , we have u = 0. Then we obtain that ǫ 1 = a q b + ab q = uv = 0, and hence item i) holds. We know that

ǫ 2 = a q 2 b + ab q 2 = ab a q-1 q+1 + b q-1 q+1 . ( 6
)
Then using Lemma 1 and ( 6) we get that

v q ǫ 2 = ǫ q+1 1 1 + t-1 i=0 v q+1 ǫ -2 1 2 i , (7)
which shows item ii). Now consider the equation

0 = y + a q b ǫ 1 y + ab q ǫ 1 = y 2 + y + v q+1 ǫ 2 1 . (8)
Using (8), its square, its 4-th power, . . ., and its q/2-th power we obtain that

y 2 + y = v q+1 ǫ 2 1 , y 4 + y 2 = v q+1 ǫ 2 1 2 ,
. . .

y q + y q/2 = v q+1 ǫ 2 1 q/2 . ( 9
)
Summing the equations in (9) we get that

y q + y = t-1 i=0 v q+1 ǫ 2 1 2 i . (10)
Combining ( 7) and (10), we get that

y q + y = 1 + v q ǫ 2 ǫ q+1 1 . (11)
We know that the set of solutions of ( 8) is a q b ǫ 1 , ab q ǫ 1 . Now we first put a q b ǫ 1 in (11). Then we have

a q 2 b q ǫ q 1 + a q b ǫ 1 = 1 + v q ǫ 2 ǫ q+1 1 , (12)
which implies that

a q 2 = ǫ q-1 1 a + a q ǫ 2 ǫ -1 1 . (13)
Similarly putting ab q ǫ 1 in (11) we obtain that

b q 2 = ǫ q-1 1 b + b q ǫ 2 ǫ -1 1 . (14)
Now using ( 13) and ( 14) together with (4) we see that item iii) holds. Finally item iv) follows from ( 5).

(⇐) Now suppose that items i),ii),iii) and iv) hold. Let a, b ∈ F q k satisfying v = ab and

ǫ 1 v = a q-1 + b q-1 . ( 15
)
Then we see that a and b are linearly independent over F q as v ∈ F * q k shows a = 0 and b = 0, while ǫ 1 = 0 shows a q-1 = b q-1 , that is, a/b / ∈ F * q . Using (15) we obtain that

ǫ 1 = ab a q-1 + b q-1 = a q b + ab q . ( 16
)
Now putting v = ab in item iii) and using ( 16), we obtain that

ǫ 2 = a q 2 b + ab q 2 . ( 17
)
Using ( 16) and ( 17) together with item iii), for 3

≤ i ≤ k-1 2
we get that

ǫ i = a q i b + ab q i . ( 18
)
Finally using Theorem 2, the linear independence of a and b, (18) and item iv) we get that Q has a codimension 2 radical.

The following definition is crucial.

Definition 1. For ǫ ∈ F * q k , let S(ǫ) be the number of v ∈ F * q k such that there exist a, b ∈ F q k satisfying v = ab and ǫ v = a q-1 + b q-1 . If q = 2, then S(ǫ) reduces to the number of v ∈ F 2 k such that there exist a, b ∈ F 2 k satisfying v = ab and ǫ v = a + b.
This condition is equivalent to the condition that a and b are the roots of the equation

y 2 + ǫ v y + v = 0.
Furthermore using (19) we have

(a ′ ) q 2 b ′ + a ′ (b ′ ) q 2 = a q 2 b + ab q 2 = ǫ 2 (v),
which shows that all the choices in (19) give the same ǫ 2 .

Here we note that using item ii) of Theorem 3 there are at most q(q +1)/2 v's giving the same ǫ 2 . Now we will show that there are exactly q(q + 1)/2 different v's giving the same ǫ 2 in (19). Let α 1 , β 1 , α 2 , β 2 ∈ F q and for i ∈ {1, 2}

a i = a + β i (α i a + b) b i = α i a + b. (20)
Then if

a 1 b 1 = a 2 b 2 we have (a + β 1 (α 1 a + b)) (α 1 a + b) = (a + β 2 (α 2 a + b)) (α 2 a + b)
which implies that

β 1 = β 2 and α 1 = α 2 or α 1 = α 2 + β -1 1 . (21)
Therefore using (21) we see that there are q + (q -1)q/2 = (q + 1)q/2 different v's giving the same ǫ 2 .

We count also the number of quadratic forms

Q of codimension 2 having invariant Λ(Q) = 1. Corollary 1. The number of quadratic forms Q of codimension 2 on F q k with Λ(Q) = 1 is ǫ∈F * q k S(ǫ).
Proof. We fix ǫ = ǫ 1 ∈ F * q k . By the proof of Theorem 4, the number of possible choices for ǫ 2 is S(ǫ 1 ) q(q + 1)/2

. Moreover there is a q(q + 1)/2-to-1 map between the set of possible values of v (whose cardinality is S(ǫ 1 )) and the set of possible choices of ǫ 2 . Also Λ(Q) = 1 if and only if ǫ 0 = v by Theorem 2. Therefore the number we are counting is q(q + 1) 2

ǫ1∈F * q k S(ǫ 1 ) q(q + 1)/2 .

Relations to Certain Maximal Curves

In this section we give some related results to maximal curves of the form

χ : y q + y = x ǫ 0 x + ǫ 1 x q + • • • + ǫ m x q m . (22)
Assume that Q is of codimension 2 and k ≥ 6. Let ℓ be the integer

ℓ = max{0 ≤ i ≤ m : ǫ i = 0}.
By Theorem 3 we have that ǫ 1 = 0 and hence ℓ ≥ 1, χ is algebraically closed over F q k and its genus is g = (q -1)q ℓ 2 .

Next we show that ℓ ∈ {m -1, m} using Hasse-Weil bound [6, Theorem 5.2.3].

Proposition 1. Assume that Q is a quadratic form on F q k over F q of codimension 2 and k ≥ 6. Let Q be as in (1) without loss of generality. Then ǫ m-1 = 0 or ǫ m = 0.

Proof. Assume the contrary and hence for the genus g of the curve χ in (22) we have g ≤ (q -1)q m-2 2 .

By choosing ǫ 0 appropriately we can assume that Λ(Q) = 1 without loss of generality. Using Hilbert Theorem 90 and (2), for the number #χ of F q k rational points of χ we obtain that #χ = 1 + qN (Q) = 1 + q k + (q -1)q k-1 .

Moreover by Hasse-Weil inequality we have #χ ≤ 1 + q k + 2gq k/2 = 1 + q k + (q -1)q m+k/2-2 .

(24)

Comparing ( 23) and (24) we get that

k -1 ≤ m + k/2 -2.
which is a contradiction.

Assume further that k ≥ 6 and Q is a quadratic form of codimension 2 on Fq k over F q . Using Theorem 3 we choose and fix ǫ 1 , a, b ∈ F * q k . Furthermore we determine the other coefficients ǫ 2 , . . . , ǫ m-1 ∈ F q k uniquely and ǫ m ∈ F q k uniquely if k is odd and uniquely modulo F q m if k is even following Theorem 3. If k is even and ab q m ∈ F q m we put ǫ m = 0 without loss of generality. Let χ be the curve in (22) with these coefficients. Proposition 2. Under the notation and assumptions as above, χ is maximal if and only if each of the following items hold: i) k is even,

Using Hilbert's Theorem 90 and a result of Klapper [START_REF] Klapper | Cross-correlations of geometric sequences in characteristic 2[END_REF], Fitzgerald computed S(ǫ) exactly in [START_REF] Fitzgerald | Highly Degenerate Quadratic Forms over F2[END_REF]Lemma 2.2].

For general F q of characteristic 2, it seems more difficult to compute S(ǫ) exactly. Using some computer experiments we currently can only conjecture that for ǫ ∈ F * q k it holds that

This conjecture is true when q = 2 (see [START_REF] Fitzgerald | Highly Degenerate Quadratic Forms over F2[END_REF]) and, for example,

that we verified by a computer. Now we determine the number of quadratic forms of codimension 2.

Theorem 4. The number of quadratic forms of codimension 2 on F q k over F q is q k ǫ∈F * q k S(ǫ) q(q + 1)/2 .

Proof. We know that ǫ 0 ∈ F q k . So there are q k choices for ǫ 0 . Now for any fixed ǫ = ǫ 1 ∈ F * q k , we will show that there are S(ǫ 1 ) q(q + 1)/2 choices for ǫ 2 . For each v ∈ F * q k satisfying v = ab and ǫ 1 v = a q-1 + b q-1 for some a, b ∈ F q k , we get an

Let v, a, b and ǫ 2 (v) be one choice. Now for any α, β ∈ F q we set

Using ( 19) we obtain that

which gives

Note that Corollary 1 is useful for counting such quadratic forms satisfying items i) and ii) of Proposition 2. However to determine how many of them also satisfy item iii) of Proposition 2 seems to be difficult (see also [START_REF] Fitzgerald | Highly Degenerate Quadratic Forms over F2[END_REF]Section 3]).

Acknowledgments. The authors would like to thank the anonymous reviewers for their valuable comments. The authors were partially supported by T ÜB İTAK under Grant No. TBAG-109T672. The work of Z. Saygı was also supported by T ÜB İTAK under Grant No. TBAG-109T344.


Risto M Hakala
email: risto.m.hakala@aalto.fi

A Lower Bound for the Nonlinearity of Exponential Welch Costas Functions

Keywords: Nonlinearity, Welch Costas functions, exponential function, Fourier transform

We study the nonlinearity of Exponential Welch Costas functions using the Fourier transform on Zn. Exponential Welch Costas functions are bijections from Zp-1 to Zp-1 defined using the exponential function of Zp, where p is an odd prime. Their linearity properties were recently studied by Drakakis, Requena, and McGuire, who conjectured that the absolute values of the Fourier coefficients of an Exponential Welch Costas function are bounded from above by O(p 0.5+ǫ ), where ǫ is a small constant. In this paper, we establish an upper bound of order O( √ p log p), which is asymptotically strictly less than the bound conjectured by Drakakis et al.

Introduction

Linear cryptanalysis [START_REF] Matsui | A new method for known plaintext attack of FEAL cipher[END_REF][START_REF] Matsui | Linear cryptanalysis method for DES cipher[END_REF] is a prominent cryptanalytic method that exploits linearity properties of S-boxes in the cipher. Traditional linear cryptanalysis is intended for ciphers containing S-boxes that can be viewed as vector-valued Boolean functions. It relies on biased linear combinations between the input and output bits of these functions. Fourier analysis is commonly used in studying the resistance of a function to linear cryptanalysis. It allows us to determine the most biased linear combination for a function by finding the largest absolute value of the Fourier coefficients of the function. This quantity is called the linearity of the function and it indicates the maximum correlation between the function and affine functions. The nonlinearity of the function is defined as the minimum distance of a function to the set of affine functions: higher distance means higher nonlinearity and better resistance to linear cryptanalysis. Optimal resistance against linear cryptanalysis is achieved by bent functions [START_REF] Rothaus | On "bent" functions[END_REF], which have a flat Fourier spectrum.

A Fourier transform on binary vector spaces is utilized when linearity properties of Boolean functions are studied using Fourier analysis. Since non-Boolean functions are also used in symmetric ciphers, e.g., in SAFER [START_REF] Massey | SAFER K-64: A byte-oriented block-ciphering algorithm[END_REF][START_REF] Massey | SAFER K-64: One year later[END_REF], linear cryptanalysis and the notion of nonlinearity have been extended to include these types of functions as well. For example, Carlet and Ding [START_REF] Carlet | Highly nonlinear mappings[END_REF] studied non-Boolean functions that have optimal nonlinearity; Baignères, Stern, and Vaudenay [START_REF] Baignères | Linear cryptanalysis of non binary ciphers[END_REF] studied how linear cryptanalysis can be used to analyze ciphers containing non-Boolean functions. The generalized notion of nonlinearity is defined using a Fourier transform on arbitrary finite groups, where the chosen groups depend on the definition of the studied function.

In this paper, we study the nonlinearity of Exponential Welch Costas (EWC) functions and their inverses, Logarithmic Welch Costas (LWC) functions. EWC and LWC functions are bijections from Z p-1 to Z p-1 , where p is an odd prime, and they are used as S-boxes, e.g., in SAFER with p = 257. The nonlinearity of EWC and LWC functions has been previously studied by Drakakis, Requena, and McGuire [START_REF] Drakakis | On the nonlinearity of Exponential Welch Costas functions[END_REF] using the Fourier transform on Z p-1 . Using a heuristic argument, they estimated that the absolute values of the Fourier coefficients of an EWC function are upper bounded by O(p 0.5+ǫ ), where ǫ is a small constant. We continue the research initiated by Drakakis et al. and derive an upper bound for the Fourier coefficients of an EWC function. The bound is of order O( √ p log p), which improves the estimate given before.

Studying nonlinearities using Fourier analysis reduces into determining the maximum absolute value of the exponential sum defined by the Fourier transform. The exponential sum studied in this paper is closely related to the exponential sum introduced and studied by Mordell [START_REF] Mordell | On the exponential sum P X x=1 exp(2πi(ax + bg x )/p)[END_REF]. The main difference between these sums is that Mordell's sum is incomplete and taken over pth roots of unity, while the sum in this paper is complete and taken over (p -1)th roots of unity. Our analysis is also different, but contains analogical elements since the results depend on similar trigonometric sums.

The paper is organized as follows. In Sect. 2, we introduce the definitions needed in the paper and set up our notation. In Sect. 3, we derive an upper bound for the linearity of an Exponential Welch Costas function. We conclude the paper in Sect. 4.

Preliminaries

In this section, we introduce some necessary definitions and set up the notation. Let n and m be positive integers. We use Z n to denote the ring of integers modulo n, p to denote an odd prime, and g to denote a generator of the multiplicative group Z * p . We also denote e(z) = exp(2πiz) and e n (z) = e(z/n) for a real number z.

Welch Costas Permutations

We recall the definition of an Exponential Welch Costas function from [START_REF] Drakakis | APN permutations on Zn and Costas arrays[END_REF]. We consider Z p-1 to be the set {0, 1, . . . , p -2} and Z * p to be the set Z p \ {0} = {1, 2, . . . , p -1}. Using this convention, we define f :

Z p-1 → Z p-1 as f (x) = (g x mod p) -1.
This function is called the Exponential Welch Costas function, and it defines a Costas permutation. Its inverse function f -1 (x) = log g (x + 1) is called the Logarithmic Welch Costas function.

Linearity and Nonlinearity

Nonlinearity is defined the same way as Drakakis et al. [START_REF] Drakakis | On the nonlinearity of Exponential Welch Costas functions[END_REF]. The Fourier transform of f :

Z n → Z m is defined by f (α, β) = x∈Zn e m (βf (x))e n (αx) for α ∈ Z n and β ∈ Z m . Definition 1. The linearity of f : Z n → Z m is defined by L(f ) = max α∈Zn max β∈Zm β =0 | f (α, β)|. Definition 2. The nonlinearity of f : Z n → Z m is defined by N L(f ) = n -L(f ) m .
It is not hard to show that a bijection and its inverse have the same nonlinearity [4, Theorem 2]. Also, the nonlinearity of an EWC function depends only on the prime p, not on the generator g of Z * p [4, Theorem 3]. It follows that all EWC and LWC functions over Z p-1 have the same nonlinearity.

Nonlinearity of Exponential Welch Costas Functions

In this section, we derive an upper bound for the linearity of an EWC function. The bound can be used to obtain a lower bound for the nonlinearity. Let p be an odd prime and f : Z p-1 → Z p-1 be an EWC function defined by f (x) = (g x mod p) -1.

We obtain an upper bound for

| f (α, β)| = p-2 x=0 e p-1 (βf (x) + αx) for α, β ∈ Z p-1 with β = 0.
We split the sum into separate parts and obtain estimates for each individual part. These estimates are given in Lemmas 1 and 2 and then used in the proof of Theorem 1 to obtain a bound for the linearity of f . Lemma 1 was shown by Drakakis et al. [START_REF] Drakakis | On the nonlinearity of Exponential Welch Costas functions[END_REF] and used to approximate the average linearity of f , but we also present the proof here for completeness. We get

|W | 2 = p-2
x,y=0

e p (r(g xg y ))e p-1 (α(xy))

= p-2 z=0 e p-1 (αz) p-2 y=0
e p (rg y (g z -1)).

Since

p-2 y=0 e p (rg y (g z -1)) = p -1 if z = 0 or r = 0, -1 otherwise,
we obtain

|W | 2 =          (p -1) 2 if α = r = 0, 0 if α = 0, r = 0, 1 if α = 0, r = 0, p if α, r = 0,
and the result follows.

⊓ ⊔

We use the ideas presented in the proof of Lemma 8.80 in [START_REF] Lidl | Finite fields[END_REF] to prove the following result.

sin π s ≥ π 6 -1 sin π 6 , so sin π s ≥ 3 s . It follows that 1 |sin πϕ(r)| = 1 sin π|ϕ(r)| ≤ 1 3|ϕ(r)| (1)
holds if |ϕ(r)| ≤ 1/6. Since 1 ≤ β ≤ p -2, we have |ϕ(r)| < 2/p for all r ∈ R ′ .
Therefore, (1) holds for p ≥ 12 and r ∈ R ′ . Since also |sin θ| ≤ |θ| for all θ, we get

|S(r)| = |sin πϕ(r)(p -1)| |sin πϕ(r)| ≤ |πϕ(r)(p -1)| 3|ϕ(r)| = π 3 (p -1) for p ≥ 12 and r ∈ R ′ . (2)
We then estimate the remaining part of the sum. Suppose that p ≥ We denote u = π(t/pβ/(p -1)), so t = pu/π + βp/(p -1) and dt = (p/π) du. Let

θ 1 = π β + 2 p - β p -1 = π 2p -2 -β p(p -1) and
θ 2 = π β + p -1 p - β p -1 = π 1 - p -1 + β p(p -1) . Since 1 ≤ β ≤ p -2,
csc u du = - 2p π ln tan π 2(p -1) = 2p π ln cot π 2(p -1) ≤ 2p π ln 2(p -1) π = 2p π ln 2 π + 2p π ln(p -1) for p ≥ 5. (3)
From ( 2) and ( 3), we get

p-1 r=0 |S(r)| ≤ 4π 3 (p -1) + 2p π ln 2 π + 2p π ln(p -1) = 4π 3 + 2 π ln 2 π p + 2p π ln(p -1) - 4π 3 < 2 π p ln p + 4p for p ≥ 12.
This inequality can be quickly checked for 3 ≤ p ≤ 11, so the result follows. ⊓ ⊔ Theorem 1. Let f be an EWC function and p be an odd prime. Then ⊓ ⊔

L(f ) < 2 π √ p ln p + 4 √ p.

Conclusion

We derived an upper bound of order O( √ p log p) for the linearity of an Exponential Welch Costas function. For reasonable values of p, the bound is not very tight as is easily shown by experiments. However, the bound shows that the nonlinearity of an Exponential Welch Costas function is asymptotically better than conjectured by Drakakis et al. [START_REF] Drakakis | On the nonlinearity of Exponential Welch Costas functions[END_REF]. These results apply for Logarithmic Welch Costas functions as well.

Lemma 1 .

 1 For any integers α ∈ Z p-1 and r ∈ Z p with r = 0, we have p-2 x=0 e p (rf (x))e p-1 (αx) ≤ √ p. rf (x))e p-1 (αx).

Lemma 2 .

 2 For any integers p ≥ 3 and β ∈ Z p-1 with β = 0, we have 1 (βy)e p (-ry) < 2 π p ln p + 4p. 1 (βy)e p (-ry). For an integer r ∈ Z p , we have |sin πϕ(r)(p -1)| |sin πϕ(r)| . Let R ′ denote the set {β -1, β, β + 1, β + 2} ⊆ Z p . We will find an upper bound for p-1 r=0 |S(r)| by comparing sums with integrals. For this purpose, we first find individual upper bounds for |S(r)|, r ∈ R ′ , since the divisor |sin πϕ(r)| in |S(r)| is close to zero when r ∈ R ′ . For s ≥ 6, we have π s -1

  5 and let R denote the set Z p \ R ′ . Since |sin θ| = |sin(-θ)| = |sin(πθ)| for all θ, we obtain r∈R

Proof.

  Let α, β ∈ Z p-1 be integers with β = 0. We have f (α, β) = p-2 x=0 e p-1 (βf (x) + αx) = p-2 x=0 e p-1 (βf (x))e p-1 (αx).For integers z, y, and n with n ≥ 2, we haven-1 r=0 e n (r(zy)) = n if y ≡ z mod n, 0 if y ≡ z mod n, 1(βy)e p (rf (x))e p (-ry).The value r = 0 can be omitted from the sum above since p-2 y=0 e p-1 (βy) = 0. By Lemmas 1 and 2, we obtain | f (α, β)| = 1 (βy)e p (rf (x))e p (-ry)e p-1 (αx)

  we have θ 1 ≥ π/(p -1) and πθ 2 ≥ π/(p -1). Therefore,

	r∈R	|S(r)| ≤	p π	θ2 θ1	csc u du ≤	2p π	π/2 π/(p-1)

Acknowledgements. The author would like to thank Gary McGuire for bringing [4] to his attention and Kaisa Nyberg for helpful suggestions. The research work has been supported by Helsinki Graduate School in Computer Science and Engineering, Academy of Finland (project #122736), Nokia Foundation, and KAUTE Foundation.


Bastien Durand

Franck Delvare

Patrice Bailly

Numerical solution of Cauchy problems in linear elasticity in axisymmetric situations

Keywords: Cauchy problem, Linear elasticity, Inverse problem, Regularization, Data completion, Boundary conditions, Friction coefficient identification

An iterative method for solving axisymmetric Cauchy problems in linear elasticity is presented. This kind of problem consists in recovering missing displacements and forces data on one part of a domain boundary from the knowledge of overspecified displacements and forces data on another part of this boundary. Numerical simulations using the finite element method highlight the algorithm's efficiency, accuracy and robustness to noisy data as well as its ability to deblur noisy data. An application of the inverse technique to the identification of a friction coefficient is also presented.

Introduction

Inverse problems can be defined by opposition to direct problems [START_REF] Kubo | Inverse problems related to the mechanics and fracture of solids and structures[END_REF] and characterized by the lack of knowledge of one of the following elements of information: the geometry of the domain, the equilibrium equations, the constitutive equations, the boundary conditions on the whole boundary of the domain and the initial conditions. According to this definition, many mechanical problems, for instance, identification of material parameters, identification of unknown boundaries (such as cracks or cavities), identification of initial boundary conditions, identification of inaccessible boundary conditions can be considered as inverse problems and more specific examples relating to elasticity problems can be found in [START_REF] Bonnet | Inverse problems in elasticity[END_REF].

In a mathematical sense, direct problems can be considered as well-posed problems. In linear cases, these problems have a unique solution which is stable (continuously dependent on the data). Conversely, inverse problems are generally ill-posed problems in the Hadamard sense [START_REF] Hadamard | Lectures on Cauchy's problem in linear partial differential equations[END_REF], since the existence or uniqueness or the continuous dependence on the data of their solutions may not be ensured. This paper examines, in axisymmetric situations, an inverse boundary value problem in linear elasticity, namely known as a Cauchy problem. It consists in recovering missing displacement and force data on some part of the boundary of a domain from overspecified displacement and force data on another part. In this case, the equilibrium equations, the constitutive equations, the domain and its boundary are known.

In order to solve Cauchy problems in linear elasticity, many regularization methods have been introduced which can be classified as Tikhonov type methods [START_REF] Bilotta | A numerical study on the solution of the Cauchy problem in elasticity[END_REF][START_REF] Koya | An inverse problem in elasticity with partially overspecified boundary conditions. II. Numerical details[END_REF][START_REF] Maniatty | Finite element analysis of some elasticity problems[END_REF]Marin and Lesnic, 2002, 2003, 2004;Marin, 2005;[START_REF] Schnur | Finite element solution of two-dimensional elastic problems using spatial smoothing[END_REF][START_REF] Tikhonov | Solution of ill-posed problems[END_REF][START_REF] Yeih | An inverse problem in elasticity with partially overspecified boundary conditions. I. Theoretical approach[END_REF][START_REF] Zabaras | Spatially regularized solution of inverse elasticity problems using the BEM[END_REF] or iter-ative methods [START_REF] Andrieux | An energy error-based method for the resolution of the Cauchy problem in 3D linear elasticity[END_REF]Delvare et al., 2010;Ellabib and Nachaoui, 2008;[START_REF] Marin | Boundary element method for the Cauchy problem in linear elasticity[END_REF]Marin et al, 2002b;Marin and Lesnic, 2005;Marin, 2009;Marin and Johansson, 2010a,b),... Tikhonov regularization methods present the advantage of leading to well-posed problems where the equilibrium equations have been modified. Some iterative methods are based on the use of a sequence of well-posed problems and others on the minimization of an energy-like functional. Numerical algorithms are implemented using different numerical methods, such as the finite element method (FEM) [START_REF] Andrieux | An energy error-based method for the resolution of the Cauchy problem in 3D linear elasticity[END_REF][START_REF] Bilotta | A numerical study on the solution of the Cauchy problem in elasticity[END_REF]Delvare et al., 2010;[START_REF] Maniatty | Finite element analysis of some elasticity problems[END_REF][START_REF] Schnur | Finite element solution of two-dimensional elastic problems using spatial smoothing[END_REF], the boundary element method (BEM) (Ellabib and Nachaoui, 2008;[START_REF] Koya | An inverse problem in elasticity with partially overspecified boundary conditions. II. Numerical details[END_REF][START_REF] Marin | Boundary element method for the Cauchy problem in linear elasticity[END_REF]Marin and Lesnic, 2002;Marin, 2002;Marin et al, 2002a,b;Marin andLesnic, 2003, 2005;Marin, 2009;Marin and Johansson, 2010a;[START_REF] Yeih | An inverse problem in elasticity with partially overspecified boundary conditions. I. Theoretical approach[END_REF][START_REF] Zabaras | Spatially regularized solution of inverse elasticity problems using the BEM[END_REF] or meshless methods (Marin and Lesnic, 2004;Marin, 2005;Marin and Johansson, 2010b). Some papers present comparisons between different numerical methods (Marin et al, 2002a;Marin, 2009).

A somewhat different resolution approach was introduced in [START_REF] Cimetière | Une methode inverse avec regularisation evanescente[END_REF](Cimetière et al., , 2001;;Delvare et al., 2002) to solve the Cauchy problem for the Laplace equation and was extended to solve the Cauchy problem in linear elasticity by Delvare et al. (2010). This approach reduced the resolution of the Cauchy problem to the resolution of a sequence of optimization problems under equality constraints. The functional is composed of two terms. At each step of the resolution, the first term (relaxation term) gives the gap between the optimal element and the overspecified boundary data, the second one (regularization term) gives the gap between the optimal element and the previous optimal element. The equality constraints are the equilibrium equations. So, at each step an optimal element is obtained which is an exact solution to the equilibrium equations and is nearer to the overspecified data than the previous optimal element calculated. In the case of compatible data, it was also proved that the sequence converges and its limit is the solution to the Cauchy problem. The additional regularization term tends to zero as its iterations continue.

In this paper, this inverse technique is extended to solve axisymmetric Cauchy problems in linear elasticity. Section 2 is devoted to the formulation of the Cauchy problem in linear elasticity in axisymmetric situations. Section 3 describes the iterative inverse method and Section 4 is devoted to its numerical implementation using the finite element method. In Section 5, several numerical simulations are presented and in Section 6, the application of the method to the identification of the friction coefficient is presented.

The Cauchy problem in linear elasticity

Let us consider an axisymmetric linear elastic material which occupies the domain Ω, with a smooth boundary Γ. We assume that the boundary is divided into three disjoint parts Γ d , Γ p and Γ u , where Γ d ∪ Γ p ∪ Γ u = Γ. We also assume that the loadings are axisymmetric.

With no body force, the equilibrium equations in cylindrical coordinates are given by:

∂σ rr ∂r + ∂σ rz ∂z + σ rr -σ θθ r = 0 ∂σ rz ∂r + ∂σ zz ∂z + σ rz r = 0 (1)
where r is the radial coordinate, θ is the angular one and z is the longitudinal one.

The Cauchy stress tensor components σ ij are related to the infinitesimal strain tensor components ε ij by the following constitutive equations:

σ ij = 2 µ ε ij + λ δ ij ε kk (2)
where λ and µ are the Lamé constants. These are related to Young's modulus E and Poisson's ratio ν as:

λ = νE (1 + ν) (1 -2 ν) µ = E 2 (1 + ν)
The strain tensor components ε ij , related to the displacement components and to those gradients, are given by:

ε rr = ∂u r ∂r ε θθ = u r r ε rθ = 0 ε θz = 0 2ε rz = ∂u r ∂z + ∂u z ∂r ε zz = ∂u z ∂z (3)
By substituting Hooke's law (2) into the governing equations (1), the Lamé equations are obtained:

(λ + 2µ) ∂ 2 u r ∂r 2 + 1 r ∂u r ∂r - u r r 2 + (λ + µ) ∂ 2 u z ∂r∂z + µ ∂ 2 u r ∂z 2 = 0 (λ + 2µ) ∂ 2 u z ∂z 2 + (λ + µ) ∂ 2 u r ∂r∂z + 1 r ∂u r ∂z + µ ∂ 2 u z ∂r 2 + 1 r ∂u z ∂r = 0 (4)
These equations can be divided by E and lead to the system of equations

L(u) = 0: 2(1 + ν) ∂ 2 u r ∂r 2 + 1 r ∂u r ∂r - u r r 2 + ∂ 2 u z ∂r∂z + (1 -2ν) ∂ 2 u r ∂z 2 = 0 2(1 + ν) ∂ 2 u z ∂z 2 + ∂ 2 u r ∂r∂z + 1 r ∂u r ∂z + (1 -2ν) ∂ 2 u z ∂r 2 + 1 r ∂u z ∂r = 0 (5) At a point x ∈ Γ, n(x)
is the outward unit normal vector as well as P(x) is the stress vector whose components are defined by:

P i (x) = P i u(x) = σ ij u(x) n j (x) x ∈ Γ
We define the adimensional stress vector by:

p i (x) = P i (x) E x ∈ Γ
It is assumed that both the displacement vector u = (u r , u z ) and the adimensional stress vector p = (p r , p z ) are given or known on the boundary part Γ d . It is also assumed that only the adimensional stress vector p is given or known on the boundary part Γ p but no condition is prescribed on the remaining part Γ u :

u(x) = φ d x ∈ Γ d p(x) = ψ d x ∈ Γ d ∪ Γ p (6)
where φ d and ψ d are prescribed vector functions. The Lamé (or Navier) system (5) and the boundary conditions (6) lead to the formulation of the Cauchy problem in linear elasticity:

         L(u) = 0 x ∈ Ω u(x) = φ d x ∈ Γ d p(x) = ψ d x ∈ Γ d ∪ Γ p (7)
This problem is difficult to solve, since it is ill-posed. When it admits a solution, its solution is unique [START_REF] Yeih | An inverse problem in elasticity with partially overspecified boundary conditions. I. Theoretical approach[END_REF], but it is known to be very sensitive [START_REF] Hadamard | Lectures on Cauchy's problem in linear partial differential equations[END_REF] to small perturbations on boundary conditions (6).

The evanescent regularization method

The iterative algorithm

Let us introduce the space H(Ω) of solutions of the equilibrium equations

(5):

H(Ω) = v ∈ H 1 (Ω) satisfying L(v) = 0 in Ω
Next, let us denote H(Γ) the space composed of couples of restrictions on Γ of elements v in H(Ω) and of their associated stress vector p(v).

H(Γ) = U = (u, p) ∈ H 1 2 (Γ) × H -1 2 (Γ) ; ∃v ∈ H(Ω) such that v| Γ = u , p(v)| Γ = p
An equivalent formulation of problem (7) reads:

         F ind U = (u, p) ∈ H(Γ) such as : u = φ d on Γ d p = ψ d on Γ d ∪ Γ p (8)
The problem ( 8) is also ill-posed even if it admits a unique solution. So then an iterative regularizing method is introduced to solve it. Given c > 0 and U 0 ∈ H(Γ) the iterative algorithm reads:

Find U k+1 ∈ H(Γ) such as :

J(U k+1 ) ≤ J(V) ∀V ∈ H(Γ) with J(V) = v -φ d 2 Γ d + p -ψ d 2 Γ d ∪Γp + c V -U k 2 Γ (9)
where the norms are defined by:

v 2 Γ d = Γ d v 2 ds p 2 Γ d ∪Γp = Γ d ∪Γp p 2 ds V 2 Γ = Γ v 2 ds + Γ p 2 ds
This iterative process is also used in Delvare et al. (2010) to solve the Cauchy problem for the linear elasticity and is a generalization of the inverse technique introduced by [START_REF] Cimetière | Une methode inverse avec regularisation evanescente[END_REF]Cimetière et al. ( , 2001) ) to solve the Cauchy problem for the Laplace equation. It can be considered as an iterative Tikhonov-type method.

In this iterative process, the equilibrium equations ( 5) are taken into account exactly since at each step the search for the optimal element is performed in space H(Γ). The functional is composed of three terms which play different roles. The first one (respectively the second one) acts only on Γ d (respectively only on Γ d ∪Γ p ). These terms represent the gap between the optimal element and the overspecified boundary data. They relax the overspecified data which can be possibly blurred by measurement noises (relaxation terms). The third term of the functional acts on the whole boundary Γ and not only on the boundary Γ u where the boundary conditions are to be completed. This term is a regularization term and controls the distance between the new optimal element and the previous one. This term tends to zero as the iterations continue.

So, at each step the optimal element obtained is an exact solution of the equilibrium equations (5) and is close to the overspecified data Φ d =(φ d , ψ d ).

Convergence results

The unique optimal element U k+1 = (u k+1 , p k+1 ) is characterized, for all V = (v, p) ∈ H(Γ), by:

u k+1 -φ d , v Γ d + p k+1 -ψ d , p Γ d ∪Γp + c U k+1 -U k , V Γ = 0 (

Properties of the iterative process

Some properties of the functional terms in the minimizing sequence can be easily established without the assumption that the data

Φ d = (φ d , ψ d )
are compatible:

• The sum of the relaxation terms S R U k = u k -φ d 2 Γ d + p k -ψ d 2 Γ d ∪Γp
is monotonically decreasing :

S R U k+1 ≤ S R U k (11)
• The regularization term J Γ (U k ) is monotonically decreasing as soon as c > 0:

U k+1 -U k 2 Γ ≤ U k -U k-1 2 Γ ( 12
)
• The sequence defined by the values of the functional J for each optimal element U k is also monotonically decreasing as soon as c > 0:

J(U k+1 ) ≤ J(U k ). ( 13
)
4. The implementation using the Finite Element Method

Discrete solutions space

The first issue in this section is to discretize space H(Γ). Our main concern has been to make use of any ordinary finite element code. Computations were run using Cast3M 1 (CAST3M, 1998) and piecewise linear finite elements, which means a piecewise constant approximation for the stress vector. Let us now discretize the domain Ω, h being the discretization parameter standing for the element size, leading to n nodes and n elements on the boundary, and m nodes inside the domain. Let V h be the space of continuous piecewise linear functions with respect to the mesh, and let us define V h (Γ) and W h (Γ) as the space of continuous piecewise linear functions and the space of piecewise constant functions on the boundary. Traces of functions belonging to V h (Ω) span the space V h (Γ), whereas the associated stress vectors belong to the space W h (Γ) of piecewise constant functions. Defining U and P as the 2n vectors standing respectively for the 2n nodal values of u and the 2n discrete values of p on the boundary, and U * the 2m-vector of internal nodal values of u, the discrete equilibrium equations read as:

  A ii A T ei A ei A ee     U * U   =   0 -BP   ( 14
)
A ii is the stiffness matrix corresponding to the Dirichlet problem and is thus invertible. Expressing the internal unknowns U * in terms of the boundary ones U, i.e performing a condensation, equation ( 14) reduces to:

( A ee -A ei A -1 ii A T ei ) U + B P = 0 (15)
The matrix form of (15) reads:

A B   U P   = 0 (16)
The finite element method leads to the definition of the following discrete compatible pairs space which reads:

H h (Γ) =    (U, P ) ∈ IR 2n × IR 2n such that E h (U, P ) = AU + BP = 0    ( 17
)
where E h denotes a linear operator mapping IR 2n × IR 2n onto IR 2n .

Discretization of the (k + 1) th iteration

Given now c > 0 and (U k , P k ) ∈ H h (Γ), iteration (k +1) of the discretized iterative algorithm reads as follows:

Find (U k+1 , P k+1 ) ∈ IR 2n × IR 2n such that :

J U k+1 , P k+1 ≤ J(V, Q) ∀(V, Q) ∈ IR 2n × IR 2n
under the 2n scalar equality constraints

E h (V, Q) = 0 (18)
Problem ( 18) is a minimization problem in R 2n × R 2n under the 2n equality constraints expressed by ( 16). Its solution is given by:

Initializing with (U 0 = 0 , P 0 = 0)

Find U k+1 , P k+1 , λ k+1 ∈ IR 2n × IR 2n × IR 2n such that :

∇J U k+1 , P k+1 + λ k+1 T ∇E h U k+1 , P k+1 = 0 E h U k+1 , P k+1 = 0 (19)
where λ k+1 is a 2n-vector of Lagrange multipliers introduced to take the equality constraints ( 16) into account. Each iteration in the iterative algorithm needs to solve a system of 6n linear equations with 6n unknowns. The matrice of this linear system is independent of the iterations and needs to be computed only once. For this reason, a direct algorithm (the Crout factorization) has been preferred to iterative methods. The factorization, which is obtained at the first step, is also used at each following step.

In the case of compatible data, the proof of the convergence of the discrete algorithm is similar to that established for the corresponding algorithm used to solve the Cauchy problem for the Laplace equation [START_REF] Cimetière | Une methode inverse avec regularisation evanescente[END_REF].

This proof is valid for all c > 0 and the c value only influences the convergence rate of the algorithm.

Numerical procedure

The procedure used during the numerical simulations is as follows:

i The meshing of the boundary is made using SEG2 elements. The SEG2 element is a finite element with two nodes which leads to a linear interpolation of the displacements. This induces a piecewise constant interpolation of the stress vector components.

ii The user specifies the meshing of the boundary specifying the number and the distribution of the finite elements on each part of the boundary.

iii The mesh of the entire domain is generated automatically by a routine included in the Cast3m software. This mesh is constituted by 4-node quadrilaterals.

iv The computation and the assembly of the stiffness matrix corresponding to the domain is performed thanks to the standard routines of the Cast3m software.

v This software generates a superelement based on the boundary and computes the corresponding stiffness matrix. This leads to the condensed stiffness matrix A.

vi The stiffness matrix A is then used by the specific code that implements the inverse method introduced.

Note that all numerical computations have been performed on a machine with a 2.20GHz Intel Core T M 2 Duo processor T7500.

The following control quantities are used to estimate the accuracy of the method:

-the L 2 (Γ) relative error made on u:

u error = Γ u -u an 2 ds Γ u an 2 ds

13

-the L 2 (Γ) relative error made on p:

p error = Γ p -p an 2 ds Γ p an 2 ds
where u an and p an denote the analytical solution.

Numerical results using analytical reference solutions

The purpose of this section is to present numerical results obtained with the method introduced. Firstly, different cases are considered for which an analytical solution is known.

Examples

A two-dimensional isotropic linear elastic medium in an axisymmetric stress state characterized by the material constants E = 200 GP a and ν = 0.34 corresponding to a steel is studied.

The domain Ω (Figure 1) is defined by:

Ω = (r, z) / R 1 < r < R 2 , - h 2 < z < h 2 with R 1 = 0.01 m , R 2 = 0.014 m and h = 0.04 m.
The boundary part Γ d is defined by:

Γ d = {(r, z) ∈ Γ/ r = R 2 }
Γ d is discretized using 160 finite elements SEG2. All finite elements have the same length and the nodes are uniformly distributed. The boundary part Γ u is defined by:

Γ u = {(r, z) ∈ Γ/ r = R 1 }
and discretized using a regular mesh with 160 finite elements SEG2. The boundary part Γ p is defined by:

Γ p = (r, z) ∈ Γ/ z = ± h 2
and each side is discretized using a regular mesh with 20 finite elements SEG2. It is used to find the displacement and the stress vectors on Γ u , from the knowledge of the displacement on Γ d and the stress vectors on Γ d ∪ Γ p .

Example 1

The data are built using the following analytical solution for the displacements:

u r (r, z) = KR 1 (1 -ν) r E∆ + KR 1 R 2 2 (1 + ν) rE∆ (20) u z (r, z) = 2νKR 1 (H -z) E∆ (21)
where

∆ = R 2 2 -R 2 1 R 1
and with K = 0.01 E. The corresponding components of the stress tensor are:

σ rr (r, z) = K R 1 ∆ 1 - R 2 2 r 2 (22) σ θθ (r, z) = K R 1 ∆ 1 + R 2 2 r 2 (23) σ rz (r, z) = 0 (24) σ zz (r, z) = 0 (25)
which correspond to a constant internal radial stress P r (z) = K (or its

dimensionless expression p r (z) = K E ) on Γ u and a free-force boundary Γ d ∪Γ p .
Example 2

The data are built using the following analytical solution for the displacements:

u r (r, z) = 2K (1 + ν) E∆    A 1 z 2 r + R 2 2 r -A 2 r 3 +A 3 R 2 2 r ln r R 2 + A 4 R 2 2 r    (26) u z (r, z) = 2K (1 + ν) E∆    z B 1 r 2 -B 2 R 2 2 ln r R 2 + B 3 R 2 2 -B 4 z 3    (27)
where

∆ = R 2 2 -R 2 1 R 1 , A 1 = 1 4 , A 2 = 1 16 , A 3 = 1 + ν 4 (1 -ν) , A 4 = 2ν 2 -2ν + 1 (1 -ν) (1 -2ν) A 3 - 3 + ν 16 (1 -ν) , B 1 = 1 4 , B 2 = 3 2 , B 3 = 3 (1 -2ν) 16ν - 1 -ν ν A 3 - 1 ν A 4 , B 4 = 1 6
and with K = 220 GP a.m -1 . The corresponding components of the stress tensor are:

σ rr (r, z) = K ∆     z 2 2 1 - R 2 2 r 2 + 3 8 R 2 2 -r 2 + 1 -3ν 2 (1 -ν) R 2 2 ln r R 2     (28) σ rz (r, z) = Kz ∆ r - R 2 2 r (29) σ zz (r, z) = K ∆ R 2 2 -r 2 2 - 3 -ν 1 -ν R 2 2 ln r R 2 -z 2 (30)
which correspond to the following stress vector P(z) on Γ u :

P r (z) = - K ∆     z 2 2 1 - R 2 2 R 2 1 + 3 8 R 2 2 -R 2 1 + 3ν -1 2 (1 -ν) R 2 2 ln R 2 R 1     P z (z) = Kz
and to its dimensionless expression p(z) = P(z) E .

Stopping criterion and influence of parameter c

In a first step, a reliable stopping criterion is needed to stop the iterative process. As in Delvare et al. ( 2010), the determination of the iteration to stop the iterative process is made using the quantity J(U k+1 ) where U k+1 is the optimal element obtained at the (k + 1) th step. The evolution of J versus the number of iterations k follows an L-curve [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the Lcurve[END_REF]. As expected from relation ( 13), the control quantity J decreases and then becomes almost constant. So the iterative process is stopped when J becomes almost constant. This stopping criterion is blind because when calculating the quantity J it does not need to know the analytical solution. Indeed, it only needs to know two successive optimal elements, the displacements data on Γ d and the pressure data on Γ d ∪ Γ p . As expected from relation ( 12), we may also notice that J Γ (the regularization term) decreases as the iterations continue. This term becomes negligible compared to the sum S R of the relaxation terms and tends to zero. This proves that the algorithm converges. As expected from relation (11), the sum S R decreases during the iterative process. After convergence, this sum remains constant and corresponds to the approximation error of the finite element method when the data are not noisy. When the displacements data φ d are noisy, after convergence, the residual value of the sum S R corresponds to the distance on the supports of data between the deblurred reconstructions and the noisy data.

It is also necessary to look at the influence of parameter c which defines the relative weight of the regularization term compared to the relaxation terms.

Different values of the parameter c are tested. Table 1 lists the results obtained for each value of parameter c by specifying the number of iterations necessary to achieve convergence, the u error , the p error and the CPU time.

The errors on u and p are quite identical for each value of the parameter c.

This confirms that the algorithm converges to the same solution whatever the value of c. However, the choice of the parameter c affects the number k of iterations needed to obtain convergence. The number of iterations necessary to achieve convergence evolves linearly with c. From the evolution of the CPU time with c (or with the number of iterations k) it can be deduced that the CPU time taken by each iteration for k > 1 is roughly 1.15 10 -2 s. This CPU time is less than the CPU time taken to achieve both the preliminary computations and the first iteration (roughly 4.23 s). For a small value of c, the convergence only takes a few iterations, however the p error increases a little. This may be explained by the fact that the regularization term becomes too weak and induces some instabilities on the reconstructions. Moreover, it is preferable to have a greater value of the parameter c which induces little additional CPU time in order to have more accurate reconstructions. Subsequently, the value of c and the number of iterations required to achieve convergence will be no longer specified.

c
It can also be observed that the errors in the numerical stresses obtained (p error ) using the iterative method are larger than those corresponding to the reconstructed displacements (u error ). This last remark is also valid for all the following numerical examples which will be analysed.

Reconstruction on Γ u with noisy displacement data φ d

It is necessary to see how the reconstructions are influenced when both components of the displacement data φ d are noisy. The noisy displacement data φ d are generated by:

φ d = φ d an + δ η φ d max (31)
where -1 ≤ η ≤ 1 is a random value, δ is the noise level in % and φ d max is the maximal value of the data of Γ d . All the reconstructions of the components of the displacement obtained can be seen to be very accurate. It can also be observed that the reconstructions of the components of the stress vector obtained using the iterative method are less accurate than those corresponding to the displacement vector. The functional is composed of terms which play different roles. As in most inverse methods, there is a regularization term which tends to zero as the iterations continue. But, in the present case, there are also relaxation terms that allow data blurred by noise to be taken into account. We therefore seek a solution which is close to the data but not a solution that exactly fits the data. The algorithm then recomputes, at each step, a solution on the whole boundary. Figure 14 represents, for Example 1, the evolution of the functional terms S R , J Γ and J versus the number of iterations k when the data φ d is noisy (δ = 5%). After convergence, the regularization term J Γ is negligible and the residual term J is equal to the sum of the relaxation terms. This residual term corresponds to the distance on Γ d between the deblurred reconstruction and the noisy data. It can also be noted that the relations (11-13) are also verified when the data φ d is noisy. This is not surprising because the inequalities were established without any assumptions on the data (for instance that Φ d is compatible).

Behaviour with respect to the mesh refinement

It is necessary to see how the reconstructions are influenced when the mesh refinement of the boundary Γ increases. In this section, the boundary part Γ d is discretized using 8N finite elements SEG2 where all finite elements have the same length and the nodes are uniformly distributed. The boundary part Γ u is discretized using a regular mesh with 8N finite elements SEG2 and both sides of the boundary part Γ p are discretized using a regular mesh with N finite elements SEG2. The level of noise δ added into the boundary displacement data φ d is fixed to 5%.

Figure 15 shows the reconstructions on Γ u of the u r -component of the displacement u obtained with N = 5, N = 10 and N = 20. Figure 16 gives the corresponding reconstructions of the p r -component of the stress vector.

These figures confirm that the inverse method is stable with respect to the mesh refinement. Similar results, which are not presented here, are obtained for the u z -component of the displacement and for the p z -component of the stress vector.

Behaviour with respect to the extension of the supports of data

In this section, we look at the influence of the type of prescribed data on the numerical solution. We investigate the numerically retrieved solutions corresponding to the following three cases associated with the given data:

case 1: u = φ d ∈ Γ d and p = ψ d ∈ Γ d ∪ Γ p case 2: u = φ d ∈ Γ d ∪ Γ p and p = ψ d ∈ Γ d case 3: u = φ d ∈ Γ d ∪ Γ p and p = ψ d ∈ Γ d ∪ Γ p
For the three cases, the same mesh is used and the level of noise δ added into the boundary displacement data φ d is fixed to 5%. The boundary part Γ d is discretized using 160 finite elements with a regular mesh. The boundary part Γ u is discretized using a regular mesh with 160 finite elements and both sides of the boundary part Γ p are discretized using a regular mesh with 20 finite elements.

Figure 17 shows the reconstructions on Γ u of the u r -component of the displacement u obtained for the three cases. Figure 18 gives the corresponding reconstructions of the p r -component of the stress vector. These figures confirm that the inverse method is accurate in each case but the reconstructions are more accurate when the extension of the supports of data is greater (case 3). Similar results, which are not presented here, are obtained for the u zcomponent of the displacement and for the p z -component of the stress vector.

6. Application of the method to the identification of a friction coefficient 6.1. Scope of the study Some materials, like concretes, rocks or geomaterials, have a material behavior which depends on the first invariant I 1 = 1 3 Trace σ of the stress tensor. In order to determine this dependence, some experimental devices have been developed to perform triaxial tests in quasi-static or dynamic situations. For one of them, the confining stress is obtained by placing the specimen in a metallic sleeve to achieve passive confinement (Figure 19a).

The compressive stress is applied directly on the specimen. The metallic ring is then subjected to internal pressure due to the lateral expansion of the specimen. This technique allows a higher level of I 1 to be obtained and has been used in quasi-static (Forquin et al., 2007) or dynamic situations [START_REF] Bailly | Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests[END_REF]Gary et al., 1998;Forquin et al., 2007).

Assuming a frictionless contact between the specimen and the metallic ring and assuming an elastic or a perfect elastoplastic behavior of the ring, the analytical solution of a hollow cylinder submitted to inner radial pressure, enables the radial stress on the inner surface of the ring to be deduced from the deformation recorded by a unique gauge glued at the outer surface [START_REF] Bailly | Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests[END_REF].

The objective is to prove that the inverse data completion method proposed could be very useful when friction at the inferface between the specimen and an elastic ring occurs. The specimen could then be replaced by the unknown loading it applies to the elastic ring (Figure 19b). This unknown loading must be identified by an inverse technique, leading to a Cauchy problem (7).

Determination of the numerical reference solution using FEM

In order to validate the procedure, we will use a numerical reference solution to a similar problem. The response of the elastic ring made of steel (Young's modulus E = 220 GP a and ν = 0.34) is studied. The loading is defined by the following dimensionless stress distributions along Γ u :

p r (z) =              a exp -β z + h 4 - h 2 ≤ z ≤ - h 4 b cosh(αz) - h 4 ≤ z ≤ h 4 a exp β z - h 4 h 4 ≤ z ≤ h 2 (32) p z (z) =            f p r (z) - h 2 ≤ z ≤ - h 4 c sinh(αz) - h 4 ≤ z ≤ h 4 -f p r (z) h 4 ≤ z ≤ h 2 (33)
where f = 0.13, a = 1.031, b = 0.963, c = -0.350, α = 37.435 m -1 and β = 13.394 m -1 . The values of these parameters are fixed assuming:

-the contact law is governed by a Coulomb-type criterion:

|p z (z)| ≤ f |p r (z)| ( 34
)
where f is the friction parameter.

-the distribution of the radial stress p r and of the axial stress p z are continuous along Γ u -the zone defined by -

h 4 ≤ z ≤ h 4 is an adherence zone -the zones defined by - h 2 ≤ z ≤ - h 4 and h 4 ≤ z ≤ h 2 are sliding zones.
The boundary part Γ d ∪ Γ p is also assumed to be a free-force boundary.

Then, a standard direct finite element simulation is performed determining a numerical reference solution. The FEM computing code used is always

Cast3m. The boundary and the domain meshes used are the same as those used in Section 5. Subsequently, a restriction of this numerical reference solution will be used to determine the input data for the Cauchy problem (7).

Example 3

Given that Γ d ∪ Γ p is a free-force boundary, ψ d is zero. For this example, only the u r -component of the displacement is noisy. This noisy data u d r is generated by:

u d r = u ref r + η δ (35)
where -1 ≤ η ≤ 1 is a random value, δ is the noise level due to uncertainties on the radial displacement measurements and u ref r

is the restriction to Γ d of the radial displacement obtained by direct simulation.

Reconstructions on Γ u

R(z) = p z (z) p r (z) (36)
The ratio R obtained is quite accurate and allows the extensions of the sliding zones where R = f , the extension of the adherence zone where R < f and the value of the friction parameter f to be identified a posteriori with relative precision. It may however be noted that some instabilities appear at the ends of the boundary part Γ u .

Conclusion

This paper introduces an iterative method for solving the Cauchy problem in linear elasticity in axisymmetric situations. This problem consists in recovering missing displacements and forces on some part of a boundary domain from the knowledge of overspecified data on another part of the boundary.

This approach reduces the resolution of the Cauchy problem to the resolution of a sequence of optimization problems under equality constraints. The algorithm reads as a least square fitting. The functional is composed of different terms. Some terms are relaxation terms which represent the gap between the optimal element and the overspecified boundary data. The other term is a regularization term which represents the gap between the optimal element and the previous optimal element. At each step, the optimal element obtained is an exact solution of the equilibrium equations and is close to the overspecified data. The regularization term vanishes when the iterations continue.

Numerical simulations using the finite element method have highlighted the accuracy and robustness of the inverse method to noisy data as well as its ability to deblur noisy data. For all the situations analysed, it can be observed that the errors in the force reconstructions obtained using the method are higher than those corresponding to the displacements.

On a numerical aspect, further developments will concern the improvement of the inverse method in order to obtain more accurate reconstructions at the ends of the boundary part Γ u . For the Cauchy problem for the Laplace equation, a first order method was introduced in Delvare and Cimetière (2008) to improve the reconstruction of the normal derivatives when the boundary part Γ u had corners. An extension of this inverse technique to our problem could be very useful.

On the experimental aspect, in further works, the inverse method will be combined with experimental techniques like digital image correlation in order to deal with experimental displacement data.

A.

  10) Theorem: Convergence of the sequence Let Φ d = (φ d , ψ d ) be compatible data associated with the compatible pair U e ∈ H(Γ). Then the sequence produced by the iterative scheme (9) strongly converges on Γ d and weakly converges to U e on Γ, where U e is the solution of the problem (8). The proof of the theorem is similar to that established for the algorithm introduced in Cimetière et al. (2000) to solve the Cauchy problem associated with the Laplace equation. This proof is valid for all c > 0. The c value only influences the convergence rate of the algorithm.

Figure 2 (

 2 Figure 2 (respectively Figure 3) shows the reconstructions on Γ u of the u rcomponent (respectively u z -component) of the displacement obtained with different noise levels (δ = 1% and δ = 5%) for Example 1. On the same figures, these reconstructions are compared with the reconstruction obtained with no noisy data. Figure 4 (respectively Figure 5) gives the corresponding reconstructions of the p r -component (respectivelyp z -component) of the stress vector. Figures 6-9 give the corresponding reconstructions for Example 2.

5. 4 .

 4 Reconstruction on Γ d : Deblurring the noisy displacement data φ d

Figure 10 (

 10 Figure 10 (respectively Figure 11) gives the reconstruction of the u r -component of the displacement (respectively the u z -component) on Γ d and the noisy data φ d used (δ = 5%) for Example 1. Figures 12 and 13 give the corresponding reconstructions for Example 2. It can be noted that all these reconstructions correspond to the analytical solution and that the noise in the data has been deleted by the algorithm.

Figure 20 (

 20 Figure20(respectively Figure21) shows the reconstructions on Γ u of the u r -component (respectively u z -component) of the displacement obtained with different uncertainty levels (δ = 5.10 -5 , δ = 1.10 -4 and δ = 5.10 -4 ).On the same figures, these reconstructions are compared on one side to the reconstruction obtained with no noisy data and on the other side to the reference solution obtained by the direct simulations. Each reconstruction of the components of the displacement can be seen to be very accurate. Figure22(respectively Figure23) gives the corresponding reconstructions of the component p r (respectivelyp z ) of the stress vector. The reconstruction of the components of the load is less accurate because at the end of the boundary part Γ u some instabilities appear in these reconstructions.

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 12 :Figure 14 :Figure 15 :Figure 16 :

 234567891012141516 Figure 2: The analytical solution and the numerical reconstructions of the u r -component of the displacement obtained on Γ u , for δ = 0%, δ = 1% and δ = 5%, for the Cauchy problem considered in Example 1.

  k u error in % p error in % CPU time in s

	1e-4 182404 0.0138653 0.543446	2103.65
	1e-5 18242 0.0138653 0.543447	214.10
	1e-6 1825	0.0138649 0.543439	25.21
	1e-7 184	0.0138617 0.543480	6.34
	1e-8 20	0.0138480 0.544349	4.45
	1e-9 3	0.0140402 0.549715	4.23

Table

1

: influence of c on the number of iterations k to achieve convergence -influence of c on the u error and on the p error

  Cimetière, F. Delvare, M. Jaoua, F. Pons, Solution of the Cauchy problem

	L. Marin, L. Elliott, D.B. Ingham, D. Lesnic, Boundary element regulariza-0.00035
	using iterated Tikhonov regularisation, Inverse Problems, 17, (2001), 553-570. tion methods for solving the Cauchy problem in linear elasticity, Inverse Problems in Engineering, 10, (2002), 335-357. 0.000346 0.000348 Analytic solution Reconstruction (0%) Reconstruction (1%) Reconstruction (5%)
	F. Delvare, A. Cimetière, F. Pons, An iterative boundary element method for L. Marin, D.N. Hao , D. Lesnic, Conjugate gradient-boundary element 0.000344
	Cauchy inverse problems. Computational Mechanics, 28, (2002), 291-302. method for the Cauchy problem in elasticity, Quarterly Journal of Me-0.000342
	F. Delvare, A. Cimetière, A first order method for the Cauchy problem for chanics and Applied Mathematics 55, (2002), 227-247. u 0.00034 r
	the Laplace equation using BEM, Computational Mechanics ,41 ,(2008), L. Marin, D. Lesnic, BEM first-order regularization method in linear elastic-0.000338
	789-796. ity for boundary identification, Computer Methods in Applied Mechanics 0.000336
	F. Delvare, A. Cimetière, J.L. Hanus, P. Bailly, An iterative method for and Engineering, 192, (2003), 2059-2071. 0.000334
	0.000332
	the Cauchy problem in linear elasticity with fading regularization effect, L. Marin, D. Lesnic. The method of fundamental solutions for the Cauchy 0.00033 Comput. Methods Appl. Mech. Engrg, 199, (2010), 3336-3344. -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 problem in two-dimensional linear elasticity, International Journal of Solids z
	A. Ellabib and A. Nachaoui, An iterative approach to the solution of an and Structures, 41, (2004), 3425-3438.
	inverse problem in linear elasticity", Mathematics and Computers in Sim-L. Marin, A meshless method for solving the Cauchy problem in three-
	ulation,77,(2008),189 -201. dimensional elastostatics, Computers and Mathematics with Applications
	P. Forquin, A. Arias, R. Zaera, An experimental method of measuring the 50, (2005), 73-92.
	confined compression strength of geomaterials, International Journal of L. Marin, D. Lesnic. Boundary element-Landweber method for the Cauchy
	Solids and Structures,44, (2007), 4291-4317. problem in linear elasticity, IMA Journal of Applied Mathematics 18,
	P. Forquin, G. Gary, F. Gatuingt, A testing technique for concrete under (2005), 817-825.
	confinement at high rates of strain, International Journal of Impact Engi-L. Marin, The minimal error method for the Cauchy problem in linear elastic-neering, 35,(2008), 425-446 . ity. Numerical implementation for two-dimensional homogeneous isotropic
	G. Gary, P. Bailly, F. Gatuingt, Testing concrete at high strains and high linear elasticity, International Journal of Solids and Structures, 46, (2009),
	rates of strain, 3th International Symposium of Impact Engineering, Sin-957-974.
	gapore (1998).

The FE code Cast3M is developed by the Department of Mechanics and Technology (DMT) of the French Atomic Energy Agency (CEA -DEN/DM2S/SEMT), http://www-cast3m.cea.fr

acknowledgements

The authors would like to thank the reviewers for their valuable comments.

000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000