Christel Pequet

Patrice Lasne

Makhlouf Hamide

Elisabeth Massoni

Michel Bellet
email: michel.bellet@ensmp.fr

A COUPLED APPROACH FOR THE MODELLING OF ARC WELDING PROCESSES

Keywords: welding, finite elements, thermal-metallurgical-mechanical coupling, material deposit

A 3D finite element model is presented, addressing some major phenomena arising in arc welding as well as their interaction: heat input, metal deposit, solidification, phase transformations in the solid state and material behavior. As a result of a simulation, the shape of the weld bead, the phase distribution, residual stresses and distortions are obtained. Some application and validation cases are discussed.

Introduction

The arc welding process allows assembling parts with or without material deposit. The high temperatures reached induce structural transformations in the metal such as solid-liquid or metallurgical solid state phase changes. Deformations in the assembly and residual stresses are also generated during the welding process. All these phenomena are complex and vary according to the process parameters (type of welding process, heat source power, welding velocity, clamping, number of pass, etc). The material behavior varies during welding due to the temperature evolution and phase transformations. Numerical simulation may help in the design and the optimization of welded assemblies. The first objectives of the simulation must be:

• To predict the geometry of the welded assembly, including weld bead shape.

• To follow the thermal and mechanical evolutions in the weld bead and in the part as well as the structural transformations and to evaluate the risks of defects during the process and in the final assembly. • To predict the final mechanical properties of the assembly, starting from the prediction of residual stresses and structural variables. Despite an abundant literature on welding numerical simulation (one can refer to the review of Runesson et al [START_REF] Runesson | Nonlinear finite element analysis and application to welded structures[END_REF]), there is still a real challenge in the numerical modeling of welding processes, because of the diversity of the phenomena involved at different time and length scales. This paper presents a contribution in that direction, coupling thermal, mechanical and metallurgical resolutions in a 3D finite element code. An original method to simulate the metal deposit will be also described.

Metallurgical Model

During the welding operation, very high heating rate and cooling rate are obtained in the heataffected zone of the parts and in the weld bead. This temperature evolution may induce transformations of the initial phases of the steel in austenite during heating, then in martensite, bainite, ferrite and pearlite during cooling. The austenizing in the part during heating is partial or total. The volume fraction of austenite g aust is a function of the maximum temperature reached and of the heating rate. The austenizing may be described with the following equation:

β α ) ( ) 1 ( 3 c aust aust A T g t g - - = ∂ ∂ (1)
With this simple model, the partial austenizing at temperatures lower than austenizing temperature A c3 is neglected and the austenizing rate is assumed to increase with temperature.

The simulation of the metallurgical transformation during cooling is based on TTT diagrams [START_REF] Inoue | Metallo-thermo-mechanical simulation of quenching process, theory and implementation of computer code 'hearths[END_REF][START_REF] Aliaga | Simulation numérique par éléments finis 3D du comportement thermomécanique au cours du traitement thermique d'acier : application à la trempe de pièces forgées ou coulées[END_REF], as illustrated in Figure 1. The incubation of ferrite, pearlite and bainite, is represented using Scheil parameter. For the growth, a Johnson-Mehl Avrami law [START_REF] Johnson | Reaction kinetics in processes of nucleation and growth[END_REF] is used to compute the fraction of ferrite, pearlite or bainite transformed. The use of the TTT diagram in the case of a non-isothermal transformation is done considering its decomposition into successive isothermal transformations. An additivity principle and the notion of fictitious time are used [START_REF] Fernandes | Prévision de l'évolution thermique et structurale des aciers au cours de leur refroidissement continu[END_REF].

Incubation and growth in the case of an isothermal transformation Principle of the simulation of a non-isothermal transformation Figure 1. TTT approach (temperature, time, transformations) for modeling solid state transformations.

The martensitic transformation is predicted using the Koistinen-Marburger law [START_REF] Koistinen | A new kinetic model for anisothermal metallurgical transformation in pure Fe-C alloys and plain carbon steels[END_REF]. The martensite volume fraction g mart is a function of the temperature when the temperature is lower than the martensite transformation temperature M s .

( ) ( )

T M A aust mart S M e g g - - - = 1 (2)
The metallurgical transformation is coupled with the thermal computation through the transformation enthalpy of each phase change (see next section). The variation of density between the different phases induces volume changes during allotropic transformations. The transformation plasticity is also taken into account to represent the abnormal deformation created during the metallurgical transformation under stress [START_REF] Gautier | Martensitic transformation under stress in ferrous alloys. Mechanical behaviour and resulting morphologies[END_REF]. In this case, the decomposition of the strain rate tensor is generalized with the following equation:

] [ ] [ tr tp th p e ε ε ε ε ε ε + + + + = (3

Coupling Between Metallurgical and Thermal Resolutions

As expressed above, metallurgical transformations highly depend on temperature history. Conversely, the impact of phase transformations (liquid-solid and solid-solid) on heat transfer should be considered. Applying the spatial averaging method on an elementary representative volume (REV) of the multiphase material, we get [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF]:

( ) ( ) k k k k k k k k k Q h g h g t = ⋅ ∇ + ⋅ ∇ + ∂ ∂ q v ρ ρ (4)
where k denotes either the liquid ( l k = ) or one of the possible solid phases (

N k , 1 =
), h the specific enthalpy, g the volume fraction and ρ the density. Assuming the temperature uniform in the REV, neglecting advection effects and phase density differences, it yields:

0 ) ( , 1 = ∇ ⋅ ∇ - ∂ ∂ +         ∂ ∂ + ∂ ∂ ∑ = T t T c h t g h t g p l l N i i i λ ρ ρ (5) Denoting 1 = k
the primary solid phase, and noticing that 1

, 2 1 = + + ∑ = N i i l g g g
,we get:

1 , 2 , 2 1 ) ( ) ( h t g h t g T h h t g t T c N i i N i i i l l p ∑ ∑ = = ∂ ∂ + ∂ ∂ - = ∇ ⋅ ∇ - - ∂ ∂ + ∂ ∂ ρ ρ λ ρ ρ (6) Denoting now 1 h h L l - =
the latent heat of solidification, and

L g c T h l T T p + = ∫ 0 d ) ( ) ( ~τ τ
, the final form of the conservation equation is very similar to the one that would be used for an enthalpy resolution with a unique solid phase, but it includes an additional right hand side:

1 , 2 , 2 ) ( ~h t g h t g T t h N i i N i i i ∑ ∑ = = ∂ ∂ + ∂ ∂ - = ∇ ⋅ ∇ + ∂ ∂ ρ ρ λ ρ (7)
This right hand side term is a source term composed of the enthalpy changes that are associated with each solid state phase change (i-j). It can be noticed that after solidification of the REV, we have 1

, 1 = ∑ = N i i g
. The previous equation becomes then:

∑ = ∂ ∂ - = ∇ ⋅ ∇ + ∂ ∂ N i i i h t g T t h , 1 ) ( ~ρ λ ρ (8)
Regarding implementation, a weak form of this equation is discretized with finite elements. The non linear resolution provides the nodal values of h ~. From these values, the unknown nodal values of variables T and g l can be calculated by solving a linear or non linear scalar equation (depending on the case: liquid, mushy, solid). This formulation permits then to take into account the energy changes associated with the solid state phase changes, while taking advantage of the stability and robustness of the enthalpy formulation for the liquid-solid phase change.

Coupling Between Metallurgical and Mechanical Resolutions

In the framework of a velocity-pressure formulation, the weak form of the momentum and mass conservation equations can be expressed as follows:

       = ∀ = ⋅ + ⋅ - ⋅ - ⋅ ∇ - ∀ ∫ ∫ ∫ ∫ ∫ ∫ Ω Ω Ω Ω ∂ Ω Ω 0 d tr 0 d d d d d d d : ) ( * * * * * * * * V p p V t V S V p V vp ε v v v g v T v ε v s v & & ρ ρ (9)
where s is the stress deviator, p the pressure, ρ the density, g the gravity, v the velocity field and T the external stress vector. The constitutive equations used in solidification problems have been presented in previous papers [START_REF] Bellet | Finite elements for a thermomechanical of analysis of solidification processes[END_REF][START_REF] Bellet | An ALE-FEM Approach to the Thermomechanics of Solidification Processes with Application to the Prediction of Pipe Shrinkage[END_REF]. Here it is just reminded that in the liquid and mushy states, steels are modeled as Newtonian and non Newtonian fluids, respectively, using a temperature dependent viscoplastic model. In the solid state they are modeled using a temperature dependent elastic-viscoplastic model. After spatial discretization using finite elements and time discretization, this results in a set of non linear equations 0 ) , ( = P V R mech , the unknowns of which are the velocity components and the pressure value at each node of the finite element mesh. While the influence of stress state on solid state phase changes has not been taken into account in this work, the reverse coupling has been treated as follows. In the context of a multiphasic solid, the variables s and p have been homogenized in a very simple way, as a first approach: the velocity field is supposed identical in the different phases. Then the deviatoric stress tensor s and the pressure p can be calculated for each phase, using constitutive parameters attached to this phase. This step is followed by an averaging of s and p, according to the linear mixing law:

k k k k p g p g = = s s ( 10
)

Modeling of Material Deposit

First of all, since the calculations are carried out in a Eulerian frame, the whole structure including the weld joint must be meshed. For that purpose, at the beginning of the simulation, a box, which will contain the final weld joint must be defined and meshed, as shown in Figure 2 (left). It is important to notice that the shape of the box will not drive the shape of the weld joint. Then, a geometrical object, C (a cone for instance) and its kinematics are defined. At each time step, elements of the mesh belonging to both the box boundary and the cone define the inlet surface for the filler metal. It is assumed here that during the material deposit, there is no formation of multiple droplets, but a continuous flow of filler material. The flow of the filler liquid metal is solved concurrently with the mechanical problem in the assembly. This is done through a VOF-like formulation in the reception box for the weld bead. In the elements of this box, a VOF (volume of fluid) function α is considered. Hence the deviatoric stress tensor s and the density ρ appearing in (9a) result from the following mixing law:

air metal air metal ρ α ρ α ρ µ α α ) 1 ( ) 1 ( 2 ) ( 2 - + = - + = ε v s s & (11)
where ) (v s metal is obtained by a viscoplastic or an elastic-viscoplastic constitutive equation, according to the local temperature. The air is considered as a compressible Newtonian medium with viscosity air µ . Accordingly, a mixture rule is applied to the mass conservation, extending equation (9b) to:

v v 0 v 0 ( ) 0 d ) 1 ( d d 1 * * =  ⌡ ⌠         + ⋅ ∇ - +         + ⋅ ∇ ∀ Ω V p t p p air metal metal χ α ρ ρ α v v (12)
where air χ is the (high) compressibility modulus in the air surrounding the filler material.

In order to track the interface between the flow of filler material and the air, an additional transport equation is solved for the VOF function α :

0 = ∇ ⋅ + ∂ ∂ α α v t ( 13
)
Figure 3 illustrates a material deposit on a plate without chamfer at different times during the process.

a) The source cone has not reached the studied section yet (in the foreground).

b) The cone is above the studied section: the filler material flows to form the weld bead.

c) The source cone has passed: the weld bead is no longer fed.

Figure 3. Formation of a weld bead on a plate.

Application and Validation: TIG Welding of Steel Plates

Among several application tests, we present here the one issued from the experimental work of Depradeux [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : validation sur cas tests de complexité croissante[END_REF]. It consists of a TIG welding processed on a thick plate of 316L stainless steel, without filler material. Figure 4 shows the geometry and the tetrahedral mesh which is refined (1 mm) along the welding path, and rather coarse in the rest of the sheet (7 mm).

Figure 4. TIG welding on a thick plate: geometry and finite element mesh.

The aim of the test is to assess the ability of the numerical simulation to predict the transient temperature field. As specified in [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : validation sur cas tests de complexité croissante[END_REF], the effective heat input is 13 MW/m 2 , and it is assumed uniformly distributed in a 5 mm source radius. The welding speed is 1 mm/s. The temperature dependent thermophysical properties of the material and the radiative and convective thermal exchanges with the environment are taken as specified in [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : validation sur cas tests de complexité croissante[END_REF]. However, the thermal conductivity is doubled for temperature exceeding the melting temperature, in order to simulate the temperature homogenization due to fluid flow in the weld pool [START_REF] Brickstad | A parametric study of residual stresses in multi-pass butt-welded stainless steel pipes[END_REF]. Figure 5 and Figure 6 show the temperature evolution on the upper face of the plate at a lateral distance of 10 mm from The comparison between the calculated weld pool and a picture of the melted zone, from the reference [START_REF] Depradeux | Simulation numérique du soudage -acier 316L : validation sur cas tests de complexité croissante[END_REF] shows again a very good agreement (Figure 7, left). It validates the heat source modeling for the test case with a weld speed of 1 mm/s. However, when taking a lower weld speed (0.67 mm/s), it can be seen in Figure 7 right that the weld pool depth is clearly underestimated in the numerical simulation, probably because the fluid convection in the weld pool is neglected in the present study. It can be thought that at lower weld speed, convection effects such as the Marangoni effect can fully develop and modify the transport of energy, and consequently the weld pool shape. In this case, a simple modification of the heat conductivity is not adequate: the solution consists then in identifying volume heat sources that could mimic such effects, or -better -in modeling concurrently the fluid flow in the pool, using a more refined finite element mesh which has to be dynamically adapted to follow the electrode trajectory. The 316L steel being entirely austenitic from room temperature to its temperature of fusion, no metallurgical transformation occurs in the heat affected zone (HAZ). Therefore, in order to test and illustrate the metallurgical-thermal coupling, the same welding test has been simulated, but using the data of A508 steel. Austenization and solid state phase change have been considered according to the formulation described above. The bainite fraction reaches a maximum value of 85% and is shown in Figure 8a. The bainitic phase is surrounded in the HAZ by a thin layer of ferrite phase as shown in Figure 8d. The parent material of a welded plate outside the HAZ consists of the initial pearlite (60%) and ferrite (40%) as shown in Figure 8c andd.

  ) in which appear different contributions: e ε & for the elastic deformation, p ε & for the deviatoric plastic deformation and th ε & for the volumetric thermal expansion, plus the contributions associated with solid state phase changes: tp ε & for the deviatoric plastic transformation and tr ε & for the volumetric change.

Figure 2 (

 2 right) illustrates the main steps of the weld joint formation.

Figure 2 .

 2 Figure 2. On the left, the meshes of the plate (in red) and of the box for weld bead formation (in blue). On the right, different steps of the weld joint formation.

Figure 5 .Figure 6 .

 56 Figure 5. On the left: temperature evolution on the upper face of the plate at two locations:x = 95 mm, y = 10 mm (calculated and measured[START_REF] Depradeux | Simulation numérique du soudage -acier 316L : validation sur cas tests de complexité croissante[END_REF]) and x = 95 mm, y = 0 mm (calculated). On the right: temperature distribution on the upper face and vertical section.

Figure 7 .

 7 Figure 7. Comparison between observed [ 11] and calculated melted zones for two weld speeds.

Figure 8 .

 8 Figure 8. Calculated volume fractions of the different phases at the end of welding.

Figure 9 (

 9 Figure9(left) shows the minor influence of phase changes on heat transfer by the release or transformation heat. In a first approach, a unique value of 76500 J/kg has been used for the transformations of austenite into ferrite, pearlite, bainite or martensite. The temperature evolution is slightly perturbed during cooling due to the bainitic transformation (110 s < t < 130 s). The transient microstructural changes at point (x = 95 mm, y = 0, z = 0) are given in Figure9b, showing the rapid austenite formation (80 s < t < 81 s) and its decomposition, yielding in this case a bainitic-ferritic structure.

Figure 9 .

 9 Figure 9. On the left: calculated temperature evolution on the upper face of the plate at location (x = 95 mm, y = 0). On the right: transient microstructural changes at the same location.

& & & & & &

Acknowledgements

This study was supported by the companies Aubert et Duval and Industeel Creusot.

Temperature (°C) Time (s)


A Pouya

M Chalhoub

Ellipsoidal anisotropy in elasticity for rocks and rock masses Elasticité à anisotropie ellipsoïdale pour les roches et les massifs rocheux

RESUME:

. Ces modèles permettent aussi une analyse simple des données expérimentales. Dans ce travail, nous montrons que ces modèles peuvent ajuster avec une bonne approximation le comportement de certaines variétés de roches et de massifs rocheux anisotropes.

) that many closed-form solutions for basic problems involving linear isotropic materials could be extended by linear transformation to cover a variety of "ellipsoidal" materials. This paper will describe two main varieties of ellipsoidal elastic models and show how well they fit the in situ data for sedimentary rocks; numerical homogenization results for several varieties of fractured rock masses will also be provided.

INTRODUCTION

In some anisotropic elasticity problems, information is available on the values of an elastic parameter in the various directions, to be used in identifying the elastic tensor. A typical case in Rock Mechanics involves deducing Young's modulus from simple compression tests in different directions or measuring the acoustic velocity in different directions on a sample [START_REF] Homand | Characterization of the Moduli of Elasticity of an Anisotropic Rock Using Dynamic and Static Methods[END_REF][START_REF] Francois | Determination of the symmetries of an experimentally determined stiffness tensor : Application to acoustic measurements[END_REF]. Based on the notion that material isotropy corresponds geometrically to the image of a sphere, an expression for anisotropy can naturally be sought through an ellipsoidal variation of a number of parameters in different directions. The uncertainty then lies in how to deduce the anisotropic tensor from this assumption. Saint Venant (1863) studied this specific question intensively by introducing the approximation of ellipsoidal indicator surfaces. The indicator surface of an elastic parameter c is the polar diagram c(n), where n is a unit vector and c(n) the value of parameter c in the material direction n. In recent years, the concept of ellipsoidal anisotropy has been adopted as a guideline for the phenomenological modelling of geomaterials such as soils, rocks and concrete (Peres Rodrigues, 1970;[START_REF] Daley | Reflection and transmission coefficients for seismic waves in ellipsoidally anisotropic media[END_REF]). Yet, the concept of anisotropic elasticity has at times been employed erroneously. For instance, Peres Rodrigues (1970) attempted using ellipsoids to fit, for several types of rocks, Young's modulus values measured along different directions. It has been shown (Pouya, 2007a) however that the Young's modulus indicator surface, i.e. the polar diagram of E(n), can never be an ellipsoid (different from a sphere), hence the parameters fitted by this author do not define any possible elasticity tensor. The correct approach calls for fitting the diagram of 4 ( ) E n by an ellipsoid; this was performed by Saint Venant in 1863.

ELLIPSOIDAL MODELS

Let's define the elastic coefficient and Young's modulus in the n direction, where n is a unit vector defining direction in the material, respectively by:

c (n) = (n⊗n): :(n⊗n) E (n)= [(n⊗n): :(n⊗n)] -1 (1)
with being the fourth-order elasticity tensor and = -1 . The main family of ellipsoidal materials considered herein has been defined by the condition that the indicator surface of 4 -c( ) n , i.e. the surface given by polar equation

r(n)= 4 -c( )
n , is ellipsoidal. This condition specifies a family of materials that depends on 12 independent parameters (Pouya, 2007a). The intersection of this family and the family of orthotropic materials in turn defines another family, denoted here by Φ 4 , that depends on 6 intrinsic parameters (c 11 , c 22 , c 33 , c 12 , c 13 and c 23 ) in the orthotropic coordinate axes, as described by the following relations between elastic coefficients: (2)

The second family to be taken into consideration in this study is the subfamily of Φ 4 , for which the indicator surface of 4 ( ) E n is also ellipsoidal. This family, denoted Ψ, is dependent upon four independent parameters c 11, c 22, c 33 and η; it is described by the following conditions:

1 η η = = - (3)
with other coefficients obtained by means of index permutation {1,2,3} and, correspondingly, {4,5,6}. This family can be equivalently defined by the four parameters E 1 , E 2 , E 3 and ν, as follows, with Ψ:

1 1 2 1 3 2 1 2 2 3 3 1 3 2 3 2 3 1 3 1 2 1 1 2(1 )
2( 1)

2( 1)

E E E E E E E E E E E E E E E E E E E 1 E E ν ν ν ν ν ν ν ν ν - - ⎡ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎢ ⎢ ⎢ - - ⎢ ⎢ ⎢ - - ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ ⎢ ⎥ ⎣ ⎦
(4) [START_REF] Pouya | Analytical solutions for foundations on anisotropic elastic soils[END_REF] remarked that some of Bohler's (1975) data on the Young's modulus of various soils could be fitted by ( 4); moreover, they demonstrated that this assumption allows simplifying foundation modelling. The theoretical features of Ψ-type materials have been thoroughly explained in [START_REF] Pouya | A transformation of elastic boundary value problems with application to anisotropic behavior[END_REF] and Pouya (2007a). For Φ 4 -type materials, Pouya (2007b) demonstrated that a closed-form expression for Green's function (displacement solution for a point force within an infinite medium) can be derived, which enables generating an explicit solution for many classical elasticity problems, as well as developing numerical Boundary Element methods for Φ-type materials.

APPLICATION TO SEDIMENTARY ROCKS

For the study of seismic wave propagation in geological layers, [START_REF] Daley | Reflection and transmission coefficients for seismic waves in ellipsoidally anisotropic media[END_REF] developed the concept of an "elliptically anisotropic" medium, as distinct from the ellipsoidal anisotropy considered in the present paper. This concept has been widely used in geophysical studies and [START_REF] Thomsen | Weak elastic anisotropy[END_REF] undertook an examination within the context of "weak anisotropy" for a large variety of sedimentary rocks. [START_REF] Thomsen | Weak elastic anisotropy[END_REF] defined four dimensionless parameters ε, δ, δ * and γ, in order to characterise transverse isotropic materials, and provided their values for a wide array of sedimentary rocks. Based on these parameters, the dimensionless elastic coefficients c * ij , defined as c * ij = c ij /c 33 , can be deduced. These coefficients are given in Table 1 for some of the samples studied by Thomsen.

Rock type and sample depth are presented in the first and second columns of the table, respectively; these two data elements allow identifying the precise corresponding material in the table established by Thomsen. The effort will now be made to compare these materials with the ellipsoidal anisotropy model (2). Since the context here is one of transverse symmetry (Fig. 1), the third condition expressed in (2) is automatically satisfied and the first two conditions become equivalent. The discrepancy between the anisotropy model of these materials and the ellipsoidal model (2) can therefore be measured by the difference between the two sides of the first equality in (2). (5)

The value of d calculated for the materials listed in Table 1 is presented in the last column of the table. The parameter c * 11 serves to describe the anisotropy. As an example, it may be noticed that for the clay shale sample at a depth of 5,858.6 m, with distinguishable anisotropy, i.e. c * 11 = 1.38, the assumption of an ellipsoidal model induces an error of just 0.3% (d = 0.003). The other rows of the table indicate that despite the pronounced anisotropy, distance to the ellipsoidal model remains relatively small. The mean value of d calculated for all sandstone, limestone, mud shale, clay shale and shale samples (about 25 in all) presented in the Thomsen table (1986) equals approximately 0.03. The ellipsoidal model therefore seems to provide a good fit for the anisotropic parameters of a variety of sedimentary rocks.

APPLICATION TO FRACTURED ROCK MASSES

Numerical homogenization represents a current method for determining fractured rock mass properties [START_REF] Pouya | Determination of rock mass strength properties by homogenization[END_REF][START_REF] Min | Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method[END_REF][START_REF] Chalhoub | Contributions of numerical homogenization methods on the rock mass classifications[END_REF]. According to this method, it proves easier to prescribe simple loads, such as simple compression or shear, along different directions and then calculate the corresponding modulus value. Fitting numerical results to an ellipsoidal model simplifies both data analysis and interpretation, in addition to reducing the number of parameters to be determined and providing effective approximate models for certain varieties of rock masses. Moreover, it yields an estimation of elastic parameter values in those directions not accessible through numerical simulation methods. [START_REF] Chalhoub | Contributions of numerical homogenization methods on the rock mass classifications[END_REF] offers a study example of a limestone rock mass slope at the border of a main road in Lebanon (see Fig. 2).

Figure 2: Studied sedimentary rock mass [START_REF] Chalhoub | Contributions of numerical homogenization methods on the rock mass classifications[END_REF] The rock mass geometry illustrated in this figure reveals the presence of two main sets of fractures, whose geometrical and mechanical properties are listed in Table 2. The geometry indicates that the rock mass exhibits orthotropic behaviour in directions 1, 2 and 3. The parameters obtained for the rock mass by homogenisation are listed in Table 3.

Using the calculated 2D elastic parameters in association with matrix transformation principles, the value of Young's modulus can be calculated in different directions. Figure 3 shows the indicator diagram of both E and 4 and allows stating that E 4 E can be well described by an ellipse. This result is promising and suggests application of the ellipsoidal model to fit numerical data. The 3D model Ψ with four parameters was chosen for this purpose.

Parameters E 1 , E 2 and E 3 , which correspond to the three orthotropic directions, can be deduced directly from Table 3 results. The Poisson's ratio ν of the model in (4) can then be deduced as an approximation, by means of the following formulae, taken from the values in Table 3:

We obtained as a result: ν =0.072. 2D numerical calculations provide all the elastic compliance parameters, with the exception of elastic compliances s 44 and s 55 [START_REF] Min | Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method[END_REF]. By using the 3D ellipsoidal model in (4), these parameters can be expressed by a combination of three Young's moduli and one Poisson's ratio ν. This result may be considered a helpful approximation of these parameters, which are inaccessible by means of 2D numerical homogenization methods.

CONCLUSION

The concept of ellipsoidal anisotropy seems to offer an attractive guideline for the phenomenological modelling of anisotropic elasticity of geomaterials, soils, rocks and rock masses; it simplifies data analysis and serves to create models with a reduced number of parameters and interesting theoretical properties. As discussed in Pouya (2007a), this concept also corresponds well to models of anisotropic elastic tensors obtained for micro-cracked materials using the micro-macro approach. The advantages offered in terms of handling geotechnical problems have already been discussed from a theoretical perspective in [START_REF] Pouya | Analytical solutions for foundations on anisotropic elastic soils[END_REF] and in [START_REF] Pouya | A transformation of elastic boundary value problems with application to anisotropic behavior[END_REF]. This paper has shown that such models could also provide a good fit for the experimental or numerical homogenization data on certain varieties of sedimentary rocks. The potential application of these models to wave propagation in rocks, for use in seismic analysis, is one of the promising avenues opened by the present work.

Figure 1 :

 1 Figure 1 : Ellipsoidal anisotropy with transverse symmetry

Figure 3 :

 3 Figure 3: Indicator diagrams of E and 4 ( )E n obtained for a rock mass studied using the numerical homogenization technique[START_REF] Chalhoub | Contributions of numerical homogenization methods on the rock mass classifications[END_REF]

Table 1 :

 1 Dimensionless parameters for several transverse isotropic sedimentary rocks, as deduced from the Thomsen (1986) data (indirect measurements), along with the distance d from the ellipsoidal model. The 4 c( )

	-	n

Table 2 :

 2 Geometrical and mechanical properties of the rock mass

	Joint filling material (limestone)
	l 1 : infinite	L 2 =45
	d 1 =19	d 2 =39
		K n =2870, K t =768
	Rock (clay)
	E=20000	ν=0.25
	l(cm): average fracture length; d(cm): average fracture
	spacing; K n (MPa/m): normal fracture stiffness;
	K t (MPa/m): fracture shear stiffness; E(MPa): Young's modulus of the rock; ν: Poisson's ratio of the rock

Table 3 :

 3 Calculated elastic compliance parameters

	Direction	E 1 (MPa)	E 2 (MPa)	E 3 (MPa)	G 12 (MPa)	ν 31 , ν 32	ν 13 , ν 12


Hervé Jégou

Cordelia Schmid

Hedi Harzallah

Jakob Verbeek

Accurate image search using the contextual dissimilarity measure

This paper introduces the contextual dissimilarity measure which significantly improves the accuracy of bag-of-features based image search. Our measure takes into account the local distribution of the vectors and iteratively estimates distance update terms in the spirit of Sinkhorn's scaling algorithm, thereby modifying the neighborhood structure. Experimental results show that our approach gives significantly better results than a standard distance and outperforms the state-of-the-art in terms of accuracy on the Nistér-Stewénius and Lola datasets.

This paper also evaluates the impact of a large number of parameters, including the number of descriptors, the clustering method, the visual vocabulary size and the distance measure. The optimal parameter choice is shown to be quite context-dependent. In particular using a large number of descriptors is interesting only when using our dissimilarity measure. We have also evaluated two novel variants, multiple assignment and rank aggregation. They are shown to further improve accuracy, at the cost of higher memory usage and lower efficiency.

Introduction

In this paper we address the problem of finding images of the same object or scene viewed under different imaging conditions. Initial approaches used simple voting based techniques [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF]. More recently they were extended based on a bag-of-features image representation [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF][START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF]. Our paper builds upon these approaches and presents methods to improve the accuracy.

The main contribution of this paper is the contextual dissimilarity measure (CDM) which takes into account the neighborhood of a vector. This measure is obtained by iteratively regularizing the average distance of each vector to its neighborhood. This regularization is motivated by the observation that a "good ranking" is usually not symmetrical in an image search system. To be more precise, if an image i is well-ranked for a query j, then j is not necessarily well-ranked for query i. Intuitively, this phenomenon yields suboptimal accuracy, as will be confirmed in this paper.

The dissimilarity measure described in this paper improves the symmetry of the k-neighborhood relationship by updating the distance, such that the average distance of a vector to its neighborhood is almost constant. This regularization is performed in the spirit of the Sinkhorn's scaling algorithm [START_REF] Sinkhorn | A relationship between arbitrary positive matrices and double stochastic matrices[END_REF]. It is also somewhat similar to a local Mahalanobis distance. Indeed, assuming all directions to be equivalent, the average distance computed over the neighborhood can be seen as a local variance.

Our CDM is learned in a unsupervised manner, in contrast with a large number of works which learn the distance measure from a set of training images [START_REF] Chopra | Learning a similarity metric discriminatively, with application to face verification[END_REF][START_REF] Frome | Image retrieval and classification using local distance functions[END_REF][START_REF] Goldberger | Neighbourhood components analysis[END_REF][START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF]. In contrast to category classification where class members are clearly defined and represented by a sufficiently large set, this does in general not hold for an image search system. Our approach is somewhat similar to the weighting schemes from text retrieval, e.g. the term frequency/inverse document frequency weighting [START_REF] Salton | Term-weighting approaches in automatic text retrieval[END_REF], which can be seen as a simple way to improve the distance [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF][START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF]. Experimental results show that the gain due to our CDM is significantly higher than the one obtained by a weighting scheme. Furthermore, these approaches can be combined. This paper also analyzes the impact of different parameters. We show that using a large number of descriptors is not always desirable, except when using the CDM. As previously shown in [START_REF] Stewénius | Object recognition benchmark[END_REF], the dataset used to generate the visual vocabulary strongly impacts the accuracy of the search. Accuracy is much higher when the vocabulary is learnt on the dataset to search, especially when using large visual vocabularies.

The first proposed variant consists in assigning each local descriptor to several visual words instead of only one. It provides a moderate but consistent improvement, especially for large visual vocabularies. As it significantly impacts the efficiency of the query, it should only be considered when very high accuracy is required. The second proposed variant is rank aggregation [START_REF] Fagin | Efficient similarity search and classification via rank aggregation[END_REF], which has not been used in the context of bag-of-features based image search before. The idea is to combine the results of several concurrent image search systems which differ in the visual vocabularies learnt on distinct subsamples of descriptors. The performance improves significantly, and increases with the number of image search systems used in parallel. However, using several image search systems has a high cost in terms of memory and computation time. We have found experimentally that the choice of three systems is a good compromise, as it provides most of the accuracy improvement.

This paper is organized as follows. Section 2 reviews the bag-of-words image retrieval approach of [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] and describes some variants. The contextual dissimilarity measure and its relationship with Sinkhorn's algorithm are de-scribed in Section 3. Section 4 presents the approach for rank aggregation. The relevance of our approach and its parameters are then analyzed in Section 5.

Overview of the image search scheme

In the following, we present the different steps of our image search framework, similar to [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF], and the tested variations.

Descriptors:

The n database images are described with local descriptors. We combine the SIFT descriptor [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] with the affine Hessian region extractor [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF]. As a variant, the 128-dimensional SIFT descriptors are reduced to 36-dimensional vectors using principal component analysis.

Visual words:

The visual words quantize the space of descriptors. Here, we use the k-means algorithm to obtain the visual vocabulary. Note that, although the generation of the visual vocabulary is performed off-line, it is time consuming and becomes intractable as the number of visual word increases (> 100000). As a variant, we use the fast hierarchical clustering approach described in [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF].

Assigning the descriptors to visual words: Each SIFT descriptor of a given image i is assigned to the closest visual word. The histogram of visual word occurrences is subsequently normalized with the L 1 norm, generating a frequency vector f i = (f i,1 , . . . , f i,V ). As a variant, instead of choosing the nearest neighbor, a given SIFT descriptor is assigned to the k-nearest visual words. This variant will be referred to as multiple assignment (MA) in the experiments.

Weighting frequency vectors:

The components of the frequency vector are weighted with a strategy similar to the one in [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF]. Denoting by n the number of images in the database and by n j the number of images containing the j th visual word, the weighted component w i,j associated with image i is given by

w i,j = f i,j log n n j . (1)
The resulting visual word frequency vector w i = (w i,1 , . . . , w i,j , . . . , w i,V ), or simply visual word vector, is a compact representation of the image.

Distance: Given the visual word vector w q of a query, similar images in the database are represented by vector(s) w i minimizing d(w q , w i ), where the relation d(•, •) is a distance on the visual word vector space. Note that the weighting scheme previously described can be seen as part of the distance definition.

Our contextual dissimilarity measure described in Section 3 operates at this stage. It updates a given distance d(•, •), e.g., the Manhattan distance, by applying two weighting factors δ i and δ j that depends on the vectors w i and w j between which the distance is computed: CDM(w i , w j ) = d(w i , w j ) δ i δ j .

(

) 2
The distance update term is computed off-line for each visual word vector of the database. The extra-storage required to store this scalar is negligible. We will show in Section 3 that computing this term is not required for the query vector.

Efficient search: The distance computation is optimized with an inverted file system exploiting the sparsity of the visual word vectors [START_REF] Zobel | Inverted files versus signature files for text indexing[END_REF]. Such an inverted file can be used for any Minkowski norm [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] when the vectors are of unit norm. For huge vocabulary sizes, the strategies proposed in [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] and [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] greatly reduce the cost of assigning the descriptors to visual words.

Rank aggregation:

The idea is to use multiple visual vocabularies, i.e., to combine the results of several image search systems where each one uses a different vocabulary. Each vocabulary is learnt on a different subset of the descriptors. The approach is described in Section 4.

Contextual dissimilarity measure

In this section, we first motivate our new measure. We then introduce the update procedure of the dissimilarity measure. This first step of this procedure, by itself, produces a new dissimilarity measure (non-iterative approach). The proposed CDM is then obtained by iterating this update step until a stopping criterion is satisfied. Finally, we underline the relationship between our approach and the projection to doubly-stochastic matrices. We also show how to efficiently compute the CDM for large datasets.

Neighborhood non-reversibility and its impact

The toy example of Fig. 1 illustrates the non-reversibility of the neighborhood for a k-nearest neighbor search. Vector 3 is a 3-nearest neighbor of vector 5, but the converse is not true. In contrast, it is (trivially) the case for an ε-search, where the neighborhood of a vector q is defined as the set of vectors x such that d(q, x) < ε.

To measure the reversibility of a neighborhood we introduce the neighborhood reversibility rate. Let us consider the neighborhood N (i) of a given visual word vector w i and #N (i) the cardinality of this set. Obviously, #N (i) = k is constant within the k-nearest neighbors framework. The neighborhood reversibility rate is then defined as follows: Figure 1: Toy example: the 3-nearest neighbors of vector 5 without (solid) and with CDM (dashed). The circles display the average distances of vector 3 (blue) and 5 (red) to their neighborhood.

1 n wi 1 #N (i) wj ∈N (i) revers(w i , w j ), (3)
where revers(w i , w j ) = 1 if w i is a neighbor of w j and w j is a neighbor of w i , 0 otherwise. The neighborhood properties of a bag-of-feature image search system are illustrated in Fig. 2. The first line shows the returned images for the query image on the left. We can observe that the three relevant images are not ranked in the first three positions. However, if we submit each of the 10 highest ranked images to the system, we can observe that the initial query is returned in the 10-neighborhood of the relevant images only, see1 columns of Fig. 2. In other words: the neighborhood reversibility is satisfied for the relevant images only. This suggests that accuracy may be improved by

• verifying for each returned image of the short-list, if the reversibility property is satisfied. However, note that this would require to perform many additional queries, i.e. one per image of the short-list2 .

• by improving the reversibility of the neighborhood.

Due to the non-reversibility of the neighborhood, we can observe that some images are returned relatively often, while others are returned only when submitting the image itself. These images are referred to as too-often-selected images and never seen images and are defined for a given neighborhood size k. Section 5.2.3 shows experimentally that the CDM significantly reduces the number of too-often-selected and never seen images. Note that the never seen images can not be retrieved even when iteratively browsing the dataset, i.e., when choosing any image in the short-list of size k as new query.

Non-iterative approach

The above mentioned problems of neighborhood non-reversibility suggest a solution which regularizes the visual word vector space. Intuitively, we would like the k-neighborhoods to have similar diameters in order to approach the reversible ε-neighborhood.

Let us consider the neighborhood N (i) of a given visual word vector w i defined by its #N (i) = n N nearest neighbors. We define the neighborhood distance r(i) as the mean distance of a given visual word vector w i to the vec-tors of its neighborhood:

r(i) = 1 n N x∈N (i) d(w i , x), (4)
where d(•, •) is a distance or dissimilarity measure, e.g. the distance derived from the L 1 -norm. The quantity r(i) is shown in Fig. 1 by the circle radii. It is computed for each visual word vector and subsequently used to define a first dissimilarity measure d * (., .) between two visual word vectors:

d * (w i , w j ) = d(w i , w j ) r r(i)r(j) , ( 5
)
where r is the geometric mean neighborhood distance obtained by

r = i r(i) 1 n . ( 6
)
This quantity is computed in the log domain. Note that the arithmetic mean can be used as well and leads to similar results. In contrast to [START_REF] Jegou | A contextual dissimilarity measure for accurate and efficient image search[END_REF], we do not use any smoothing factor denoted α in [START_REF] Jegou | A contextual dissimilarity measure for accurate and efficient image search[END_REF]. Indeed, the new update term r/ r(i)r(j) used here (instead of its square) amounts to choosing α = 0.5. For this non-iterative approach, it provides close to optimum results in terms of accuracy, see [START_REF] Jegou | A contextual dissimilarity measure for accurate and efficient image search[END_REF]. The relation d * (•, •), referred to as non-iterative contextual dissimilarity measure (NICDM), is not a distance: although the symmetry and the separation axioms are satisfied, the triangular inequality does not hold. This is not a problem in the context of finding the nearest neighbors of a given vector w i . Comparison measures that do no satisfy the properties of a distance have been used in information retrieval. For instance, the image search system of [START_REF] Örster | Image retrieval on large-scale image databases[END_REF] explores the use of the Shannon-Jenson divergence and a metric derived from a LDA model.

Note that in Eq. 5 the terms r(i) and r do not impact the nearest neighbors of a given vector. They are used to ensure that the relation is symmetric.

Let us now consider the impact of the approach on the average distance of a given vector w i to the others. This impact is formalized by the following ratio:

j d * (w i , w j ) j d(w i , w j ) = j r r(i)r(j) . (7)
Together with the observation that j r(j) = rn , we have

j d * (w i , w j ) j d(w i , w j ) = r r(i) n , (8)
which in essence means that, on average, the NICDM d * (•, •) reduces distances associated with isolated vectors (with r(i) > r) and, conversely, increases the ones of vectors lying in dense areas.

Iterative approach

The update of Eq. 5 is iterated on the new matrix of dissimilarities. The rationale of this iterative approach is to integrate the neighborhood modification from previous distance updates. Denoting with a superscript (k) the quantities obtained at iteration k, we have

d (k+1) (w i , w j ) = d (k) (w i , w j ) r(k) r (k) (i)r (k) (j) . (9)
Note that, at each iteration, the new neighborhood distances r (k) (i) are computed for each visual word vector w i .

The objective of this iterative approach is to minimize a function representing the disparity of the neighborhood distances, in other terms to optimize the homogeneousness of the dissimilarity measures in the neighborhood of a vector. This function, here defined as

S (k) = i |r (k) (i) -r(k) |, (10)
is clearly positive. Its minimum is zero and satisfied by the trivial fixed-point of Eq. 9 such that ∀i, r(i) = r.

Let us define a small quantity ε > 0. As a stopping criterion, the algorithm terminates when the inequality S (k) -S (k+1) > ε is not satisfied anymore. This ensures that the algorithm stops within a finite number of steps. In practice, for ε small enough, we observed that this criterion led r (k) (i) to converge towards the fixed-point of Eq. 11.

Let us recall that, by contrast to [START_REF] Jegou | A contextual dissimilarity measure for accurate and efficient image search[END_REF] where a smoothing factor had to be set, here Eq. 9 amounts to choosing α = 0.5. This choice does not impact the accuracy, as it is observed that for α < 0.9, the algorithm converges towards the same set of values.

At this point, we can only compute the CDM between visual word vectors of the database, due to the iterative design of this distance. In order to compute directly the CDM from the original distance, one has to maintain a cumulative distance correcting term δ

(k) i during iterations, as δ (k+1) i = δ (k) i r(k) r (k) (i) . ( 12
)
Denoting by δ i the quantity δ (k-1) i when the algorithm terminates, it is easy to show that The k-nearest neighbors of a given query q are then the minima given by NN(q) = k-argmin j d(q, w j ) δ j .

d k (w i , w j ) = d(w i , w j ) δ i δ j . ( 13
) (a) -3 -2 -1 0 1 2 3 -2 -1 0 1 2 (b)
Note that finding the nearest neighbors of a query vector q does not require the knowledge of the update term associated with q, as shown in Eq. 14. That is why we only need to compute the partial term d(q, w j ) δ j .

Hence, it is possible to find the nearest neighbors with the CDM for a vector which is not in the database. One has just to store together with a given database visual word vector w i the corresponding distance update term δ i , which in terms of storage overhead is clearly negligible. Given the original distance matrix and the parameters k (neighborhood size) and ε (convergence threshold), the pseudo-code for computing the update terms of the CDM is given by Algorithm 25. This algorithm may be advantageously implemented in the log-domain.

Relationship between CDM and projection to doubly-stochastic matrices

Here we briefly describe the projection of distance matrices to doubly-stochastic matrices, which is closely related to the CDM introduced above. As discussed the square distances between the corresponding points in a set Q, where the points in Q are obtained by a stereographic projection of the points in P . The points in Q are confined to a hypersphere of dimension d embedded in IR d+1 , where d is the dimensionality of the subspace spanned by the points in P (or, if the points in P are confined to a hypersphere, d is the dimension of that hypersphere). In Fig. 3 we illustrate the projection of a set of points P ∈ IR 2 to a sphere in IR 3 . Note that the points in Q all have the same average squared distance to other points, since all rows (and columns) of B sum to unity. Thus, if we consider the "new" distances between the i-th point and other points j, given by [B] ij = [A] ij δ i δ j , we see that they are given by a scalar correction δ j of the original distances [A] ij , and in addition we have the scalar correction δ i which is constant for all j. Clearly, CDM and projection to doubly stochastic matrices modify the distances from a point i to other points j in the same way: by multiplicative correction terms for each j.

The projection to doubly-stochastic matrices suffers from one weakness in the context of image retrieval: the projection takes into account the distances between all pairs of points. However, as high dimensional data-like our visual word frequency vectors-usually live on a (non-linear) manifold of lower dimension embedded in the vector space, only small distances are meaningful in the sense that they tend to correspond to small geodesic distances along the manifold. Large distances, however, are not indicative for the corresponding geodesic distances along the manifold: the geodesic distance may vary greatly for constant distance in the embedding space. For this reason, our CDM method regularizes pairwise distances in smaller neighborhoods instead of regularizing all pairwise square distances. Note that the flavor of the CDM method-global analysis of properties in small overlapping neighborhoodsresembles that used by many recently developed methods for non-linear dimensionality reduction inspired by ISOMAP [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] and LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. The relevance of this choice is demonstrated in 5.2.1.

CDM for very large sets

For very large datasets, the bottleneck of the CDM is the distance computation between all frequency vectors, which in theory is of quadratic complexity in the number of images. Fortunately, finding the true neighborhood of frequency vectors is not required to obtain accurate update terms. Suboptimal approximate nearest neighbor search of frequency vectors as proposed in [START_REF] Jegou | A contextual dissimilarity measure for accurate and efficient image search[END_REF] greatly improves the efficiency of the update terms' calculation. We showed that using this strategy to compute the CDM update terms moderately decreases the accuracy of the search while allowing the use of the CDM for a set of 1 million images [START_REF] Jegou | A contextual dissimilarity measure for accurate and efficient image search[END_REF].

Another possible simple method for computing approximate update terms consists in choosing a fixed set of frequency vectors, not necessarily extracted from the dataset to index, to compute the neighborhood distance. This results in a very fast computation of the update terms. Moreover, the terms do not depend on the dataset to index: they only depend on the chosen fixed set. As a consequence, this avoids having the so-called out-of-sample extension, as adding new frequency vectors to the dataset does not modify the other update terms.

Exploiting multiple vocabularies using rank aggregation

The search accuracy is improved by using t independently generated vocabularies instead of only one, i.e., by generating t different visual word codebooks. Hence, t distinct image search systems are used, each of which is implemented with an inverted file. The underlying motivation is that it is very unlikely that each system returns the same false positives, while it is very likely that true positives are returned often. Rank aggregation combines the results of t different subsystems with the method of [START_REF] Fagin | Efficient similarity search and classification via rank aggregation[END_REF] which was proposed to perform approximate nearest neighbor search of vectors in the spirit of locality sensitive hashing [START_REF] Shakhnarovich | Nearest-Neighbor Methods in Learning and Vision: Theory and Practice[END_REF]. For each retrieved image we compute its median rank over all ranked lists returned by the t sub-systems. Ties are resolved arbitrarily, but not randomly.

Example: Let us consider t = 3 different subsystems and a query for which an image of the dataset is ranked 1 st in list 1, 4 th in list 2 and 3 rd in list 3. The set of ranks obtained for this image is (1, 3, 4), hence its median rank is 3.

Note that this approach can also be applied to other quantiles, i.e., instead of the median one can choose the first quartile. However, we have observed that the median rank consistently provides good results, in contrast to other quantiles.

This approach improves the accuracy, as shown in Section 5. Its main drawback is that the storage requirements are t times higher. Moreover, except for region extraction and descriptor computation, all the other steps of the image search system, i.e., descriptor quantization and inverted file querying are computed t times instead of only one. For very large datasets, where querying the inverted file becomes the bottleneck of the algorithm, the query becomes roughly t times longer.

Experiments

Datasets and evaluation criteria

The evaluation is performed on two datasets: Figure 4: Query examples: ranked lists obtained with and without CDM. We can observe that the CDM significantly improves the ranking.

• the Nistér and Stewénius (N-S) dataset [START_REF] Stewénius | Object recognition benchmark[END_REF],

• a set of frames [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] extracted from the movie "Run, Lola, Run!".

The N-S dataset is composed of 2550 objects or scenes, each of which is imaged from 4 different viewpoints. Hence the dataset contains 10200 images. The Lola dataset is composed of 164 video frames extracted at 19 different locations in the movie.

For all the experiments, we used the Hessian-Affine detector [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF]. Except when explicitly specified, the threshold is set to 100, resulting in an average number of 2269 descriptors per image, and the L 1 distance is used.

Three datasets have been used to perform the k-means clustering: the Corel set which is uncorrelated with the evaluation sets, as well as the N-S and Lola datasets used for the evaluation. For this purpose we have extracted from these datasets subsamples of about 1 million SIFT descriptors, except for the Lola dataset where the whole set of descriptors was used.

Two different measures are used to evaluate the impact of the various parameters and variants: the average normalized rank (ANR) and, for the sake of comparison, the measure used by Stewénius and Nistér [START_REF] Stewénius | Object recognition benchmark[END_REF]. For a given query image, the ANR [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] is given by where n q is the number of queries, n is the number of dataset images and n rel (i) is the number of images which should be retrieved for image i. This measure indicates the average normalized position (between 0 and 1), in which a relevant image appears. For instance, ANR≈0.01 means that the average rank of a relevant image is approximately equal to 1000 for a dataset of 100000 images.

ANR = 1 n q nq i=1 1 n n rel(i)   n rel(i) j=1 rank(j) - n rel(i) (n rel(i) + 1) 2   , (16)
Clearly, a lower ANR signifies better accuracy. The measure proposed in [START_REF] Stewénius | Object recognition benchmark[END_REF] counts the average number of correct images among the four first images returned. This measure is meaningful because there are 4 relevant images per object in the N-S dataset.

CDM

CDM vs Sinkhorn algorithm

All the experiments in Table 1, Table 2 and Table 3 show a significant improvement when using a distance regularization method (CDM or Sinkhorn). Note that the parameters are summarized in the caption. Table 1 shows that the Sinkhorn algorithm improves the results. However, the gain due to the CDM is significantly higher, and this for all the tested parameters. Thus, regularization with local distances only is very important in our context. The relevance of the CDM is also confirmed by experiments on the preprocessed data of [START_REF] Stewénius | Object recognition benchmark[END_REF], as shown in Fig. 7.

Fig. 4 illustrates some typical queries for which the CDM significantly improves the results. For the N-S dataset (first two lines), the query with no CDM returns flowers, which are often irrelevantly returned. The capability of the CDM to reduce the impact of the too-often-selected images is clear in this context. The query on the Lola database (two last lines) is even more impressive. The first three images are correct with and without CDM. Although four images seem wrong for both queries, they are in fact correct for the CDM, as the images correspond to the same location (Deutsche Transfer Bank) observed from significantly different viewpoints.

Neighborhood size of the CDM

The only parameter of the CDM is the neighborhood size n N . Fig. 5 shows the impact of this parameter on the performance of the iterative approach. We can observe that the sensitivity to this parameter is moderate: the accuracy increases significantly in the case of very small neighborhoods and decreases gracefully when using large neighborhoods. A small neighborhood also results in lower computational cost. In the rest of this paper, the size n N is fixed to 10, although better results may be obtained by optimizing this parameter.

too-often-selected and never seen images

The impact of the CDM on the neighborhood reversibility is very important. This has been verified on the N-S dataset with a vocabulary size of 10000. For a neighborhood size of 10 the neighborhood symmetry rate (Eq. 3) increases from 0.37 to 0.62. The percentage of never seen images, see Section 3.1 for the definition, decreases from 9.7% to 0.2%. Similarly, for the 10200 queries of the N-S dataset the most frequent image is returned 54 times in the first 10 positions with the CDM, against 1062 times using the standard L 1 distance.

Impact of the parameters and variants

Clustering

We have implemented and evaluated the hierarchical clustering approach [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF].

Comparing Table 3 Exp. #1 and Table 1 Exp. #2, we can see that hierarchical clustering reduces the accuracy. However, its significantly reduces the computational cost for assigning SIFT descriptors to visual words, especially for large vocabularies. Note that the concurrent approach of [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] offers similar efficiency as [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF], but provides better accuracy (very similar to k-means).

The dataset used for the clustering may have an impact on the accuracy, as shown in Tables 1 and2. For these two datasets we compare in column "clustering dataset" k-means clustering on an uncorrelated dataset (Corel) with kmeans clustering on the evaluation dataset itself (either N-S or Lola). In both cases the results are improved by generating the visual vocabulary with a subset of the dataset on which the experiments are performed.

This confirms the observation made by Nistér and Stéwenius [START_REF] Stewénius | Object recognition benchmark[END_REF], i.e., that using the evaluation set for clustering significantly improves the results. This is particularly true when using a large vocabulary, as shown in Fig. 7 which shows some results obtained with the preprocessed dataset of [START_REF] Stewénius | Object recognition benchmark[END_REF].

When comparing the experiments #1-#4 with #5-#7 in Table hold for the Lola dataset (see Table 2). A possible explanation is that for this set the clustering was performed on the entire set of descriptors and not only on a subsample, hence emphasizing the adaptation of the visual vocabulary to the evaluation dataset.

Number of descriptors

Fig. 6 shows that the number of descriptors extracted for each image has a strong impact on the accuracy. We can observe that the accuracy increases up to a certain point only, i.e., using a too high number of descriptors decreases performance. A possible explanation is that the strongest interest points, i.e., with high cornerness values, are diluted among those with low cornerness, and that this results in noise in the frequency vectors. However, the best number of descriptors to be used depends on the other parameters. In particular, the CDM benefits from a high number of descriptors. Note that for this experiment, the interest points have been generated using a low threshold and are then filtered based on their cornernesses value. This is slightly different from the standard setup, where the suppression of nonmaxima is performed after the thresholding.

Vocabulary size

Table 1, Table 2 and Fig. 7 show that bigger vocabularies provide better retrieval accuracy. However, the gain is rather moderate, except when the visual vocabulary is learned on the evaluation set itself (see Fig. 7). In this 2.9 case the bag-of-features image search system becomes very similar to an approach which matches individual descriptors, i.e., for very large vocabularies the number of visual words is equal to the number of descriptors in the clustering dataset and all the descriptors are used as centroids. Note that search results obtained when matching individual descriptors outperform those of bag-of-features based search, but significantly increase the search complexity.

Norm

It was observed in [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] that the Manhattan distance provides better results than the Euclidean one. This observation is confirmed in our experiments for the two datasets and also holds when the CDM is used, see 7: N-S dataset: impact of the vocabulary size and of the clustering dataset (independent dataset "CD" or the N-S set). These curves have been generated using the N-S preprocessed data [START_REF] Stewénius | Object recognition benchmark[END_REF]. Fixed parameter: n N = 10.

Dimensionality reduction of SIFT with PCA

We reduce the dimensionality of the SIFT descriptor with PCA from 128 to 36 dimensions. Table 3 Exp. #3 shows that PCA marginally reduces the accuracy, while decreasing the computational cost of the visual word assignment. However, the hierarchical SIFT assignment of [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] decreases the assignment cost more significantly and at the same time obtains comparable results, see Table 3 Exp. #1 and Exp. #3. Using PCA in this context is, therefore, of limited interest.

Multiple assignment of SIFT descriptors

The MA of SIFT descriptors to visual words slightly improves the accuracy of the search (see Table 3, Exp. #4 to #7) at the cost of an increased search time, due to the impact of the method on the visual word vector sparsity. For instance, for V = 30000 visual words the number of multiplications performed is 7 times higher for MA×3 than for the simple assignment. It should be used for applications requiring high accuracy. Note that the number of assignments must be small, e.g. 2 or 3, as we have observed that the accuracy decreases for larger values.

Rank aggregation

Table 4 presents the results obtained with the rank aggregation method described in Section 4. The visual vocabularies have been generated using distinct SIFT subsamples of the N-S dataset, obtained by modifying the seed of the random number generator. The number of votes required for an image to be added to the ranking list is equal to 2 for 3 visual vocabularies, 3 for 5, 5 for 9 and 10 for 19.

The results show that rank aggregation improves accuracy. The scores are consistently improved for all sets of parameters. The best score obtained with rank aggregation is 3.68 against 3.57 for a single vocabulary. The trade-off between accuracy and efficiency can also be adjusted by choosing a smaller number of distinct visual vocabularies. Hence using only 3 distinct ranking sets is sufficient to obtain a fair improvement.

Comparison with the state-of-the-art

For the N-S dataset, the CDM approach obtains a N-S score of 3.55 (maximum 4) for a CDM computed with n N = 10 neighbors and 30000 visual words learnt on the Corel dataset. Combining the CDM with the MA improves this results to 3.61 for 50000 visual words. Our best score of 3.68 has been obtained using rank aggregation (see Table 4). The best previous score [START_REF] Stewénius | Object recognition benchmark[END_REF] most time consuming approach and a visual vocabulary learned on the N-S dataset itself. Our best ANR score for the Lola movie is 0.0046, significantly outperforming the previous best score 0.0132 [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF]. Note that, by contrast to their work, we use only one type of descriptor (in this case their best score is 0.0196) and no temporal filtering. Our approach is still better (0.0118) if the visual words are learned on uncorrelated data. In [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] the visual vocabulary was learnt on the Lola dataset.

Conclusion

This paper introduces the contextual dissimilarity measure to compare frequency vectors of a bag-of-features image representation. This new measure is based on a distance regularization algorithm in the spirit of the Sinkhorn's algorithm, which projects distance matrices on doubly-stochastic matrices. In contrast to this algorithm, our regularization uses local distances only, similar to recently proposed methods for non-linear dimensionality reduction.

The performance of our approach has been demonstrated for a bag-of-features based image search system. A large set of experiments shows that the accuracy is significantly and consistently improved by the CDM for two different datasets. We also analyze several variants and the impact of the main parameters of our image search system. Our final system significantly outperforms the state-of-the-art on both datasets.

Figure 2 :

 2 Figure 2: Illustration of the neighborhood non-reversibility of a bag-of-features image search system (L 1 distance). The top-left image is the query image. The first line displays the ordered set of images returned for this query. The three relevant images are ranked 1 st , 5 th and 9 th . Each column represents a new image search performed by querying with the top image of this column. The boxes indicates where the initial query (top-left image) is returned for the new queries.

Figure 3 :

 3 Figure 3: (a) A set of points in IR 2 with zones of different density. (b) Projection onto the sphere in IR 3 using Sinkhorn's algorithm to yield doubly stochastic distance matrix. Note the much more uniform density of the points on the sphere.

Figure 5 :

 5 Figure 5: CDM: impact of neighborhood size n N for the N-S dataset. Parameters: k-means clustering performed on corel, L 1 distance.

Figure 6 :

 6 Figure 6: N-S dataset. Impact of the number of descriptors on the retrieval accuracy for a vocabulary of size 30000 learned on the N-S dataset.

Table 1 :

 1 N-S dataset. Impact of distance regularization (CDM with n N = 10 and Sinkhorn algorithm[START_REF] Sinkhorn | A relationship between arbitrary positive matrices and double stochastic matrices[END_REF]), of the dataset used for k-means clustering (uncorrelated Corel dataset or the N-S dataset itself) and of the vocabulary size.

		clustering clustering vocab.		N-S score
		method	dataset	size	no CDM Sinkhorn CDM
	#1	k-means	corel	1000	2.90	3.05	3.45
	#2	k-means	corel	10000	3.10	3.21	3.53
	#3	k-means	corel	20000	3.12	3.25	3.54
	#4	k-means	corel	30000	3.14	3.28	3.55
	#5	k-means	N-S	1000	2.91	3.14	3.49
	#6	k-means	N-S	10000	3.16	3.25	3.57
	#7	k-means	N-S	30000	3.26	3.31	3.57

Table 2 :

 2 Lola dataset. Impact of k-means clustering dataset, vocabulary size, norm (Manhattan L 1 or Euclidean L 2 ) and number of neighbors n N used in the CDM calculation.

	the next

Table 3 :

 3 1, we can observe that the CDM is less influenced by the learning set. This remark does not N-S dataset. Impact of the variants and of the following parameters: clustering method (k-means or hierarchical[START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF]) performed on the corel dataset, vocabulary size, norm (L 1 if not specified or L 2 ), use of the PCA (36 dimensions), multiple assignment (MA) of descriptors to visual words. Fixed parameter: n N = 10.

		clustering	training vocab.	N-S score
		method	set	size	no CDM CDM variant
	#1 hierarchical	corel	10000	2.90	3.41	[3]
	#2	k-means	corel	30000	2.68	3.37	L 2
	#3	k-means	corel	30000	3.01	3.49	PCA
	#4	k-means	corel	30000	3.12	3.59	MA×2
	#5	k-means	corel	30000	3.07	3.60	MA×3
	#6	k-means	corel	50000	3.21	3.60	MA×2
	#7	k-means	corel	50000	3.18	3.61	MA×3

Table 3

 3

	Exp. #2 and

Table 4 :

 4 is 3.29 for their number of distinct vocab. N-S dataset. Rank aggregation: impact of the number of distinct visual vocabularies (1 to 19), here learned on the N-S dataset itself. Note that the three first rows (1 distinct ranking only) correspond to no rank aggregation. Fixed parameter: n N = 10.

			N-S score	N-S score
	visual vocab.	size	without CDM with CDM
		1000	2.91	3.49
	1	10000	3.16	3.57
		30000	3.26	3.57
		1000	2.93	3.54
	3	10000	3.20	3.63
		30000	3.29	3.63
		1000	2.94	3.55
	5	10000	3.21	3.64
		30000	3.31	3.65
		1000	2.95	3.56
	9	10000	3.22	3.65
		30000	3.32	3.67
		1000	2.96	3.57
	19	10000	3.22	3.66
		30000	3.33	3.68

Note that the reversibility rate for this query and a neighborhood of size 10 is equal to 0.3.

We have tested this variant and observed that its performance is inferior to the distance regularization proposed in this paper.

Acknowledgements

We would like to acknowledge J. Sivic, A. Zisserman, D. Nistér and H. Stewénius for kindly providing their datasets.

% compute their geometric mean r := ( i r(i))

end for % update the pairwise vector dissimilarity measures % diag(δ) is the diagonal matrix with diagonal δ

in the previous section the CDM rescales distances d(w q , w i ) by a scalar factor δ i . The correction factors δ i are set in such as manner that for all points w i the average distance from w i to its nearest neighbors, r(i), become similar:

The CDM is a modification of Sinkhorn's scaling algorithm [START_REF] Sinkhorn | A relationship between arbitrary positive matrices and double stochastic matrices[END_REF]. Sinkhorn's algorithm takes a positive matrix A and iteratively normalizes the rows and columns to have unit L 1 norm. The algorithm is guaranteed to converge to a unique fixed point, which is a doubly stochastic matrix: all rows and columns sum to unity. Interestingly, when applied to a squared distance matrix A there is a geometric interpretation of Sinkhorn's algorithm [START_REF] Johnson | On the diagonal scaling of Euclidean distance matrices to doubly stochastic matrices[END_REF]. In this case Sinkhorn's algorithm yields a matrix B = ∆A∆ which is doubly stochastic, and ∆ is a diagonal matrix with diagonal elements δ i . If P is the set of points that generated the square distance matrix A, then elements of B correspond to


R Bernal

Ch Tassius

F Melo

J-Ch Géminard

Elastic response and wrinkling onset of curved elastic membranes subjected to indentation test

Keywords: . 46.32.+x Static buckling and instability -81.70.Bt Mechanical testing, impact tests, static and dynamic loads -68.60.Bs Mechanical and acoustical properties

Starting from a polymeric-fluid droplet, by vulcanization of the fluid free surface, curved elastic membranes, several nanometers thick and a few millimeters in diameter, which enclose a constant fluid volume, are produced. In an indentation-type test, carried out by pushing the membrane along its normal by means of a microneedle, under some conditions, wrinkles are likely to appear around the contact region. Interestingly, we observe that the instability does not significantly alter the force-displacement relation: the relation between the force and the displacement remains linear and the associated stiffness is simply proportional to the tension of the membrane. In addition, we determine that the wrinkles develop when the stretching modulus of the membrane compares with its tension, which provides a useful method to estimate the elastic constant.

Introduction

The mechanical properties of thin elastic membranes are of wide interest to control processes such as cell motility [1,2], wound healing [3], electronic surface protection [4] and fabrication and characterization of ordered structures at the nano and micro scale [5,6]. For instance, in biological applications, it has been recently shown that hollow structures, like vesicles or capsules, have high potential for holding sensitive contents until release into the external environment is appropriate [7,8]. The properties and structure of these capsules are assessed with indentationtype tests by deforming them axi-symmetrically with the help of microcantilevers or flat surfaces. Relations between the indentation force, the membrane inner pressure and the structure deflection have been obtained [7,9]. However, due to large membrane deflection involved in such tests, wrinkling is often observed when the (compressive) stress exceeds a certain threshold [10,11,12]. In turn wrinkles structures might reveal important information on membrane properties and thickness as well [12,13,14].

Here we develop further the potential of indentation tests for the assessment of mechanical properties of thin curved membranes. Our configuration has the advantage of providing well-defined boundary-conditions, which allows us to develop simple numeric description accounting for large deflections of the membrane. We observe that the applied force depends linearly on the imposed displacement and, from the theoretical and numerical analysis of the problem, we show that the stiffness relates to the membrane tension. In addition, we observe that radial wrinkles, which are likely to form around the indenter, do not significantly affect the force-displacement relation.

We determine that the onset of the instability depends on the membrane tension, wrinkles forming when the tension compares with the stretching modulus of the membrane. Thus, measurements of the tension, thanks to the indentation test, and the determination of the instability onset provides, in a single experiment, a useful method to estimate the membrane elastic modulus.

Experimental procedure

The sample consists of a curved thin membrane attached to the upper circular edge of a cylindrical rigid cavity (Fig. 1, Top). First a homogeneous substrate is obtained by curing a layer (thickness h = 500 µm) of poly dimethyl siloxane (PDMS, Sylgard 184), previously deposited onto a piece of glass plate. The cylindrical cavity is then created by cutting a hole of diameter 2r 0 = 3 mm in the substrate.

To obtain the elastic membrane, a given volume V of PDMS polymer fluid is injected inside the cavity obtaining a drop whose free surface has a given curvature. Drop curvature is easily increased by increasing the volume V . Membranes are prepared by vulcanization, as described in Refs. [1,2], of the surface of the polymeric fluid. Fluid vulcanization is simply produced by the action of a homogeneous flame from a Bunzer-type burner [1,2,15]. During this process, and in order to obtain a more homogeneous film, the sample is quickly turned upside-down and displaced at constant horizontal velocity (a few mm/s) just above the blue part of the flame. The flame width being about 5 mm in our case, the contact-time between the flame and the fluid surface is about 1 s which insures homogeneous membrane thickness at large scale. Thus, in our experimental conditions, we produce a membrane of reproducible thickness h m = (60 ± 2) nm, at the scale of tip size, as determined by an AFM-based technique [15]. The elastic constant associated to the stretching is B = (240 ± 10) mN/m [15] with B = Eh m /(1 -σ 2 ), where σ and E stand for the Poisson ratio and Young Modulus respectively. The resulting height H and radius of curvature R of the cross-section in a vertical plane (Fig. 1) are obtained by means of a camera (Nikon DXM1200) which images the sample from the side. Thus, our procedure results in an almost spherical (radius of curvature R) membrane, stretched over a circular frame, which encloses a chosen constant-volume of liquid.

A force F is applied at the membrane center by means of a flexible micro-needle (Fig. 1, Top), whose tip is rounded to avoid damaging of the membrane. The radius of the bead tip is ρ ≃ 11 µm (Fig. 1, Bottom). In addition, in order to insure that only the bead touches the membrane, the needle tip is bent making almost a right angle with the membrane plane around the contact point. The micro-needle calibration is performed by recording the tip deflection due to the force applied by a AFM cantilever of known compliance. For the whole set of data reported below, the force F = k n δ n , where the stiffness k n = (65 ± 3) 10 -3 N/m (δ n stands for the deflection of the needle along the normal to the surface). To obtain δ n and the corresponding indentation, δ m , the membrane is imaged from the side. Good contrast in images is obtained by guiding light from a LASER inside the needle (Fig. 1, Bottom). The distance δ m is obtained directly from the images whereas the deflection δ n is obtained as the vertical deflection of the needle tip with respect to a reference line attached to the needle. Note that the total vertical displacement imposed to the needle writes δ m + δ n . Notice that a right choice of the microneedle stiffness ensures a membrane deflection that is of the same order than the tip deflection which, in turn optimizes the sensitivity of the force measurements. The sample is also imaged from the bottom by an optical microscope (Nikon Eclipse TE-300) and a high resolution SLR digital camera (Nikon D1X) using the Nomarsky contrast mode which makes possible both to determine the membrane features and to detect the wrinkles, even if the amplitude of the associated vertical displacement is small (Fig. 1, right panels).

Force-displacement relation

Qualitatively, when the membrane center is pushed downwards, a region of positive radial slope develops around the needle tip. The radius r c of this region is defined by the location of the circular apex, clearly visible as a brighter ring in the images from the bottom (Fig. 1, right panel a, dashed circle). The radius r c increases when the imposed displacement δ m is increased and, depending on the preparation of the membrane, radial wrinkles are likely to appear above a given onset (Fig. 1, right panel b). Experimentally, F is proportional to δ m in a large range (δ m ∈ [0, 2ρ], Fig. 2). The conclusion interestingly holds true even when wrinkles form.

In the following, we estimate the relation between the stiffness k and the initial tension of the membrane T 0 by calculating numerically the variation of the membrane surface energy which results from the indentation at constant inner volume V . In principle, the surface energy of the membrane has two contributions : a surface energy γ (corresponding to the liquid-membrane interface) and the energy associated to the elastic modulus B (the membrane being very thin, we neglect the energy associated with the bending). However, the surface tension γ is expected to be small (at least much smaller than the surface tension associated with the air-PDMS liquid interface, γ ≃ 23 mN/m) as the liquid and the membrane are made of the same material. Thus, in a first step, we can neglect the contribution of the surface tension and assume that the initial tension of the membrane T 0 only originates by the elastic contribution (the initial stretching of the membrane is due to the preparation process). We shall see that the experimental value of T 0 is usually large enough for the latter assumption to be well satisfied. The main conceptual steps to achieve our numerical calculations are given below, a detailed description will be presented elsewhere. Our analysis is based on the minimization of the elastic free energy per unit surface, Γ , which only depends on the strain field u and the elastic constant B. The axissymmetric membrane is described by the radial δr and the vertical δz displacements which are calculated with respect to a membrane without tension. This means that the tension T 0 itself is accounted for by an initial isotropic stretching of the membrane. Thus, in any of the considered cases, either the flat or the spherical membrane, the shape of the indented membrane results from an initial deformation produced by a homogenous tension T 0 , plus the additional deformation produced by the pushing tip. The calculation of the latter strain field, satisfying the conditions at both the pushing tip and the circular outer boundary, is the goal of the numerics. In a first approximation, the boundary condition at a spherical tip being difficult to account for, we consider instead a cylindrical tip, of radius r t , and we mimic the indentation by a vertical displacement of the membrane in a circle of radius r t . In practice, the Lagrangian is, L = S0 Γ dSλ S dV , where the first term is estimated by integrating over the initial shape of the membrane whereas the second term is estimated by integrating over its deformed shape. The differential equations that govern δr and δz are obtained from the variational transformation whose approximated solution are obtained by a shooting method. Starting from the center, where the action of the indenter is a vertical displacement of the membrane, i.e., δr = 0 and δz = -δ m in a circle of radius r t , we look the correct matching at the outer boundary, where we consider that the indentation does not produce any additional deformation; δr = 0 and δz = 0 at r = r 0 . The integration is performed with an adaptive step size control in the Runge Kutta integration method.

With the approximation above, the stiffness is found to satisfy,

k T 0 = 0.39 + 2π(r t /r 0 ) β , (1)
with β = 0.37 ± 0.04. Equation 1 holds true, to within 5%, in a large range of radius of curvature R and elastic modulus B (∈ [0, 10 T 0 ]). Note that k is thus found, in the elastic case, to increase slowly with r t [like (r t /r 0 ) β and not like log(r t /r 0 ) as found when the surface tension alone is taken into account] which explains why the forcedisplacement relation still remains linear in the case of a spherical indenter. Indeed, in this latter case, r t corresponds to the radius of the contact region between the tip and the membrane which rapidly reaches a value of the order of ρ when δ m is increased. Finally, we point out that a rather large imposed displacement, δ m , at the center does not necessarily induce large in-plane strains.

Experimentally, we observe that k decreases when the initial volume V (or equivalently 1/R) increases. However, according to Eq. 1, this dependance of k on V is not due to the change in the curvature 1/R of the membrane. The dependance of k on V is due to the thermal dilation during the preparation process which leads to an initial tension T 0 that depends on the geometry of the droplet. Indeed, during the preparation, the fluid is heated in order to produce the polymerization of the free surface. After polymerization, the temperature and, thus, the surface area of the membrane without tension and the inner volume change. The thermal effects thus tend to change the surface to volume ratio and, as a consequence, to affect the membrane tension. In order to prove that the amplitude of the thermal effects is enough for accounting for the experimental changes in the membrane tension, we considered numerically the case of only a thermal contraction of the inner volume and estimated the resulting change in the tension T 0 . Assuming that the polymerization leads to a reproducible initial tension T * 0 and considering only a relative decrease θ of the entire fluid volume, we obtained that the tension T 0 of the membrane at the center should decrease according to T 0 = T * 0 -1.66Bθ H r0 in our experimental conditions. Experimentally, the tension is expected to vanish for H/r 0 ≃ 0.54 (Fig. 2, inset) which permits to estimate the dilation θ that would account for the experimental variation of the tension. Taking into account the ratio B T * 0 ≃ 4 [15], we obtain θ ≃ 30%, which corresponds, the thermal expansion coefficient of the fluid being 3.5 10 -4 K -1 [16], to a temperature-change of about 200 K, which is a reasonable value considering that the fluid layer is thin and the temperature increase required for surface vulcanization is the same order. Thus, our analysis indicate that the membrane highly contract immediately after preparation. Note also that the thermal effects are likely to explain the large value of the tension T 0 (about 140 mN/m) which is expected in the limit H → 0.

Thus, indentation provides a reliable method for measuring the local tension T 0 of the membrane. The force is proportional to the displacement. The stiffness k is proportional to the tension T 0 , does not depend on the local curvature and only slowly increases with the size of the indenter. The result holds true even when the displacement is as large as the indenter diameter and when the membrane undergoes a wrinkling instability. However, we expect the curvature of the membrane to play a role when it compares with the typical size of the tip. Indeed, from the tip point of view, the membrane is flat as long as the radius of curvature of the membrane is much larger than the tip typical size. In additional experiments, to be published elsewhere, using latex membranes of 20 cm in diameter, inflated at a controlled pressure and indented by tips of large size, we have observed that the effect of membrane curvature remains weak.

Membrane shape and instability

We observe that the radius of the concave part around the tip, r c , increases almost exponentially with δ m (Fig. 3a),

r c = r max c 1 -exp - δ m δ typ m , (2)
and the appearance of wrinkles does not alter significantly the dependance of r c on δ m . The radius r c only slightly depends on the initial height H and, thus, on the tension T 0 in the whole experimental range.

From the exponential interpolation of the data, we determine the characteristic deflection δ typ m ≃ (2.8 ± 0.1) × 10 -3 r 0 and the asymptotic value r max c ≃ (38 ± 1) × 10 -3 r 0 . The theoretical analysis as well as the numerical results suggest that, in the limit of small δ m , r c ∝ √ Rδ m . In the inset of the figure 3a, we report r 2 c /R as a function of δ m and obtain experimentally r c ≃ Rδ m /8. Note however that the numerical results, that involve equations that are only valid in the limit of small strains, fail in accounting for the exponential growth of r c when the displacement δ m is further increased (typically, when δ m ≃ 5 × 10 -3 r 0 which corresponds to δ m ≃ 1.5 ρ, the vertical displacement being thus larger than the radius of the pushing tip).

If T 0 is small enough and δ m is increased above the critical value, δ c m , the membrane is subjected to an ax-ial wrinkling instability (Fig. 1, right panel b). We note that wrinkles extend to the border of the concave region, i.e. the length of the wrinkle, L, almost equals the radius r c (Fig. 3b). The most interesting result lies in the fact that the critical displacement δ c m drastically depends on the tension T 0 (Fig. 3b, inset). For small T 0 , the transition is nearly supercritical as δ c m almost vanishes. On the contrary, for slightly larger T 0 , the instability is subcritical as δ c m takes a finite value: the membrane is subjected to a wrinkling instability provided that the imposed displacement is larger than a given, finite, value and, at the transition, the wrinkle length jumps to a finite value [L = r(δ c m ), Fig. 3b, circles and squares]. If T 0 is initially too large, wrinkles are no longer observed.

Results above can be understood considering that the membrane instability requires that the orthoradial tension becomes negative (compression) at regions where wrinkles are observed. The variation of the tension in the orthoradial direction at the distance r from the center, σ ⊥ , is about B δr r where δr denotes the displacement of the membrane toward the center due the action of the tip so that σ ⊥ ∼ T 0 + B δr r . Due to the indentation, the membrane is locally pulled toward the center and we can estimate, assuming that the membrane is wrapping around the pushing tip, δr ∼ r t (sin α-α) at the distance r t α from the center. The angle α stands for the angle that the membrane makes with the horizontal at the contact with the indenter. Even expected to mainly depend on the ratio δ m /ρ, the displacement δr can hardly be estimated without solving the entire problem. However, the local-strain is estimated to never exceed 1-2 π ≃ 0.36 when the displacement δ m is increased. Accordingly, the membrane should be stable as long as T 0 > 0.36 B. Taking into account the value of the elastic modulus B, we thus predict, in acceptable agreement with the experimental results, that no wrinkle can form for a tension T 0 typically larger than 86 mN/m. Finally, even if such a behavior cannot be accounted by the description of the system in the small displacement limit, we would like to mention that the interpolation of the experimental data leads to δ c m = ρ(aT 0 /B) 6 , with a = (5.0 ± 0.2) for B = (240 ± 10) mN/m (Fig. 3b, inset). Although this latter relation would deserve further theoretical effort, the constant a is a universal constant for the problem as long as the radius of the frame keeps much larger than the other lengths in the system. Thus, our result can be useful for the assessment of the modulus B in similar systems.

Conclusion

We demonstrated experimentally that the measurement of the stiffness while pushing, with the help of a microneedle, a curved membrane enclosing a constant volume provides information about the tension T 0 of the membrane (or equivalently on the inner overpressure ∆P , provided the radius of curvature is known) and not about the elastic constant B. Such stiffness measurements hold even in the presence of wrinkles. Wrinkling instability occurs when the membrane tension T 0 compares with the elastic modulus B, which provides an alternative method for measuring B in small sized systems. For instance for capsules, made of a polymer network and inflated by osmotic pressure [7], one could decrease the osmotic pressure down to values at whose wrinkles develop. By reporting the behavior of the critical displacement δ c m as a function of ∆P , which is simply related to T 0 , the elastic modulus B would be obtained.

In summary, our study shows that, the force -displacement relation in the linear regime, for an indentation type measurement, only provides information about initial tension of the membrane. The assessment of elastic modulus B, can be achieved either from a nonlinear regime as in the case of reference [15] or by examination of the onset of the wrinkling pattern as in the present manuscript. This research was supported by Conicyt-Chile under Anillo project No ACT-95.

Fig. 1 .

 1 Fig. 1. Sketch of the sample. Top: Schematic of force measurement indicating the deflection, δn and the indentation δm. Bottom: The needle tip shape designed to avoid off center contact. Right panel a) A concave central part marked by a brighter ring; H/ro = 0.23. Right panel b) Radial wrinkles for H/ro = 0.40.

Fig. 2 .

 2 Fig. 2. Force F vs. displacement δm. The force F increases with δm. The result holds true without (Open symbols) or with (Full symbols) wrinkles. Inset, the stiffness k and, thus, the tension T0 decrease with the initial height H (or volume V ).

Fig. 3 .

 3 Fig. 3. Radius rc vs. displacement δm. a) rc increases almost exponentially (grey line) with δm. The result holds true without (Open symbols) or with (Full symbols) wrinkle pattern (The typical error bar is displayed on a single point). Inset: r 2 c /R vs δm shows that rc ∝ √ Rδm for small δm, with an initial slope ≈ 1/8. b) L increases almost exponentially with δm. The grey line is the exponential fit displayed in panel (a). Arrows mark the transitions to wrinkling and L is taken equal to zero below these thresholds. Inset: δ c m , drastically increases with T0. Wrinkles (W) and No Wrinkles (NW) regions.