Thierry Biben

Chaouqi Misbah

Anne Leyrat

Claude Verdier

; A Leyrat

An advected -field

Keywords: 82, 70, Kj -Emulsions and suspensions, 47, 11, +j -Computational methods in fluid dynamics, 47, 55, Gf -Multiphase and particle-laden flows

come

Introduction. -Since the pioneering experimental and theoretical work of Taylor [1,2], the problem of droplet breakup, in view of its industrial relevance, has motivated a huge number of studies in various scientific communities like general physics, fluid mechanics, rheology and computational physics. The challenge is to understand and therefore give means to control the size of the fragments that are produced after breakup in out of equilibrium regimes. Using roller devices, Taylor was able to produce two-dimensional shear and elongational flows, and applied them to investigate droplet deformation in a suspending fluid. The relevant dimensionless parameters for this problem are the capillary number Ca = η out Rγ/Σ (where η out is the suspending fluid viscosity, R is the initial droplet radius, γ is the deformation rate and Σ is the interfacial tension), which is related to the droplet deformability, and the viscosity ratio p = η in /η out (η in is the viscosity of the drop), as demonstrated experimentally by Grace [3]. Critical capillary numbers for breakup are indeed determined as functions of p, in shear and elongational flows over ten decades, as well as for mixed flows [4]. Analytical models [5,6] based on the perturbation theory have been developed since, as well as numerical models, which are more appropriate for large deformations [7] [8][9][10][11]. While the latter are based on the boundary integral method [8][9][10][11], the level set method [12] has been introduced recently to cope with several difficulties encountered in the boundary integral formulation, like the front tracking problem and the treatment of breakup itself. The key point in the level set method is to regularize the macroscopically infinitely thin interface by introducing a smoothing function which is constructed from the distance field to the interface. Applications to binary fluids have been considered very recently [13,14]. Similar methods like the volumeof-fluid method have been used as well (see [15]). In the same spirit, Phase-Field models were introduced twenty years ago in the field of crystal growth, with the major difference that the auxiliary phase field has a physical meaning: it accounts for the solid-liquid phase equilibrium, and is then constructed from thermodynamical considerations. These Phase-Field methods should then naturally apply to binary fluids, with the same accuracy as level set methods, but should incorporate more physical ingredients, such as short range interactions between the interfaces (correlations due to the broad microscopic nature of the interfaces) which are absent in level set approaches. Indeed, Loewenberg and Hinch [16] have shown that in the absence of short range attractions (like the Van der Waals force) coalescence cannot be accounted for, which renders the determination of the rheological properties of concentrated emulsions difficult. Phase Field models incorporate the short range physics of the interface naturally. Coalescence can occur due to the short range interactions between the interfaces and breakup is controlled by the free energy functional of the Phase-Field itself. Van der Waals forces can be added to account for polarizability effects when necessary. Although Phase Field approaches usually give a good description of the physical phenomena (see for example [17] and references therein), they often fail to provide quantitative results since in practice the broadness of the interface is overestimated. The scope of this letter is to present a modified Phase Field approach which allows for quantitative predictions. This method should be numerically more efficient than traditional boundary integral formulations for concentrated emulsions since the resolution time mainly depends on the box volume which scales like the number N d of droplets, at a given concentration, whereas boundary integral methods scale like N d 2 in the most favorable situtation (no viscosity contrast). We will however focus here on the ability of the method to give quantitative predictions for a single droplet, where many results are available, which is certainly the worst situation from the computational point of view in comparison to other methods.

In a first part we discuss the Phase-Field formulation for binary fluids, and show that slight modifications are necessary to account for the proper macroscopic physics. As a result, the new method is intermediate between the pure Phase-Field formulation and the levelset method. We will now use the name "Advected-Field" method (AF) for this modified formulation, which will be justified later. Comparison of the predictions of the AF approach with available experimental and theoretical data for the breakup problem are presented in the second part.

The advected field approach. -A binary fluid can be described by an order parameter field which identifies the two components. Using a Landau-Ginzburg description, a purely isotropic fluid can be described by a scalar order parameter field φ(r, t) where r denotes the spatial coordinates while t is time. The simplest free energy functional that can be written to describe a binary coexistence is the usual square gradient theory:

F [φ] = V olume W (φ) + ǫ 2 2 [∇φ] 2 dr (1)
where W (φ) is a double-well potential describing the bulk coexistence of the two phases, and the square gradient contribution sets the width of the interface to √ 2ǫ. A simple prescription for W (φ) is W (φ) = (1 -φ 2 ) 2 /4 which determines the coexisting values of the field, ±1. In a Cahn-Hilliard spirit, a dynamical equation can be written for the temporal evolution of the order parameter field which is advected (i.e. transported) by the flow:

∂φ ∂t + v.∇φ = -Γ dW dφ -ǫ 2 ∆φ ( 2
)
where v is the macroscopic velocity field. While the left contribution expresses the hydrodynamic transport of the field by the flow, the right term is responsible for the chemical (diffusive) relaxation of the field to its equilibrium value. Γ is a constant which sets the relaxation time of this process. A conservative dynamics ( i.e. φdr = cste ) can be built by substituting -Γ∆ to Γ. This Cahn-Hilliard formulation (2), very similar to the formulation already given by Jasnow and Viñals [17], cannot be used however to describe the dynamics of a mesoscopic droplet (radius R ∼ 100µm) in a flow. A first practical reason is due to the value of ǫ which should be of the order of 10 -4 R to describe the actual width of a binary interface. This value is out of reach at the moment for present computers and a typical accessible value is 10 -2 R (0.04R in the present work). The second and more serious problem comes from the free energy formulation itself. The velocity field appearing in ( 2) is the macroscopic velocity field, which expresses the transport of matter. The velocity of the interface must then correspond to the local velocity field. This is not indeed the case with the above formulation, which can be illustrated quite simply: if we cancel the velocity field in (2), and start from an initial arbitrary shape for the drop (i.e. the zero isocontour of the order parameter field φ), according to the velocity field the shape should not change. Solving (2) however will lead to a spherical shape if the dynamics is conserved or the drop will simply disappear if the dynamics is not conserved, which simply expresses the minimization of the free energy functional (1).

The basic reason for this problem is due to surface tension of the order parameter field, which is not coupled to the velocity field, and thus induces a relaxation even in the absence of a flow. A simple prescription which cancels this effect is to add the counter term suggested by Folch et al. [18] in their study of the Saffman-Taylor instability, and to solve the following equation:

∂φ ∂t + v.∇φ = -Γ dW dφ -ǫ 2 [∆φ + c |∇φ|] (3)
where "c", appearing in the counter term, is the local curvature field defined everywhere in space, which corresponds to the curvature of the isocontours of φ. The precise expression of this counter term can be justified by an analysis of the asymptotic behaviour of this formulation when ǫ goes to zero (the "sharp interface limit" in the Phase-Field terminology). We cannot reproduce this detailed analysis here, but this term simply cancels to first order in ǫ the curvature dependence of the interfacial hyperbolic tangent equilibrium profile of the interface, φ(r) = tanh(r/ √ 2ǫ), where r is the algebraic distance to the φ = 0 isocontour. In other words, the laplacian appearing in (2) is replaced by ∂ 2 φ/∂r 2 and thus the Cahn Hilliard prescription only plays a role in the normal direction. The interface is thus treated as locally flat by the microscopic Cahn Hilliard theory. The global shape evolution of the interface (i.e. in the parallel direction) is entirely contained in the advection term, and is thus controlled by the Stokes equation alone. The normal and the parallel directions have thus been decoupled, keeping the physical ingredients of the Cahn Hilliard theory for a microscopically quasi-planar interface in the normal direction, and removing the unphysical coupling with the global macroscopic shape relaxation induced by a strong overestimation of the interfacial width due to numerical considerations. With this formulation, the order parameter field φ becomes a quasi-passive field (i.e. the φ = 0 isocontour is transported by the flow), the Advected Field, but the shape of the interface is fixed to its hyperbolic tangent shape thanks to the new restoring term (r.h.s. of (3)).

Obviously, surface tension is a fundamental physical ingredient for the description of droplets dynamics, and it has to be incorporated in the dynamical equation for the velocity field:

ρ ∂v ∂t + v.∇v = ∇.σ -∇P + F interf ace ( 4
)
where ρ is the density of the fluid, σ is the extra stress tensor of the fluid, P is the local pressure and F interf ace is the force exerted by the interface on the fluid. This last term can be expressed by F interf ace = d cΣδ(r)n, where d is the interface dimensionality, c is the mean curvature field introduced above, Σ is the actual interfacial tension, r is the algebraic distance to the interface (the φ ≡ 0 isosurface) and n is the normal outgoing vector at the interface. We used the convention that c = -∇.n/d, which leads to a negative value for a sphere. δ(r) is the Dirac distribution in the sharp interface limit, which can be regularized at finite values of ǫ by setting δ(r)n → ∇φ/2 , where we have used the convention that φ < 0 inside the drop. Since the AF goes from -1 to +1 while crossing the interface, we can easily check the normalization condition whatever the shape of the profile. The final expression for F interf ace is then:

F interf ace = d 2 cΣ∇φ (5)
σ, the stress tensor, is chosen to obey the local Newtonian law :

σ ij = η(φ) (∂ i v j + ∂ j v i ).
We have emphasized in this expression the spatial dependence of the viscosity η which accounts for a viscosity contrast between the internal fluid (i.e. the drop) and the external medium.

The simple prescription η = η in (1 -φ)/2 + η out (1 + φ)/2 allows the viscosity to vary smoothly across the interface to interpolate between the internal value η in and the external viscosity η out . Note that the method can also allow for a density contrast, but it is not the purpose here. Since we consider droplets of typical size about 100µm, viscosities around 1 P a.s or more, and rates in the range [10 -2 -10 2 s -1 ], the Reynolds number Re = ρa 2 γ/η out should be smaller than 10 -3 and inertial effects are negligible, the momentum equation then simply reduces to the Stokes equation where ρ has been set to zero. The pressure field can easily be eliminated in Fourier space thanks to the incompressibility condition ∇.v = 0. The solution of these dynamical equations can be obtained numerically by solving the AF equation ( 3) in the direct space using a finite difference method while the Stokes equation (which we choose to solve using a relaxation scheme, keeping ρ∂v/∂t) is best solved in the Fourier space, using an implicit scheme. For this purpose we used a regular mesh with a square or cubic unit cell, depending on dimensionality, and took special care of anisotropy problems due to the discretization of the differential operators. We used in practice operators which are isotropic up to order three in the lattice spacing. The lattice spacing is chosen to be exactly ǫ, whose value is fixed to 0.04R, where R is the radius of the drop at rest. When an external flow is imposed, the total flow can, without any restriction, be separated into two components v = v ext + u where v ext is the applied external flow, a shear or an elongational flow, and u is the flow induced by the droplet. The boundary conditions at the edges of the mesh have been chosen to be periodic for the AF as well as the component of the flow u, v ext being treated analytically.

Results. -To illustrate the ability of the method to account for the breakup process we show in Fig. 1 a breakup sequence in 3D for an initially elongated drop (spherocylindrical geometry) relaxing in a quiescent fluid. To test quantitatively the method, we compared the data provided by the AF method to the analytical study made in 2D by Buckmaster and Flaherty [6] for a drop placed in an elongational flow defined by: v x ext = γ x and v y ext = -γ y, where x is the coordinate along the elongation axis. Interestingly these authors were able to calculate the deformation of a drop in an elongational flow up to the burst transition. The two dimensionless relevant parameters are the viscosity ratio p and the capillary number Ca as defined in the introduction. Below a certain critical value Ca c the drop deforms under the elongational constraint to reach a steady shape, which can be characterized by its "elongation" 1 -B/L where B is the width and L is the length. Above Ca c no stationary shape exists and the drop elongates until it breaks. Such a dynamic behaviour corresponds to a saddle-node bifurcation, and thus an unstable branch should exist that merges with the stable branch at the critical capillary number. This is precisely the scenario found by Buckmaster and Flaherty [6]. This type of bifurcation seems to be quite generic in the physics of deformable objects under hydrodynamic constraint, it also occurs for vesicles (tank-treading to flipping transition) [START_REF] Biben | [END_REF]. We can see in Fig. 2 a comparison of the AF method with the analytical theory of [6]. This figure shows the predicted elongation of the drop at a given value of the capillary number for a viscosity ratio p = 1 in 2D. At a given value of Ca two shapes satisfy the stationary condition below the critical capillary number Ca c ∼ 0.18 (or equivalently 4π 3/2 Ca ∼ 4 from Fig. 2): convex shapes when 1-B/L < 0.72 and non-convex shapes above 0.72, and no stationary shape exists above Ca c . While the less elongated shape corresponds to the actual stationary shape, the most elongated one is unstable and thus is not physical. The AF model seems to provide quite good quantitative agreement with the analytical theory, even close to the critical capillary number since the approximation presented in Fig. 2 deviates from the full theory [6]. The value of the critical capillary number Ca c = 0.18 obtained with the AF method agrees perfectly well with the theoretical expectation.

The 3D situation can be investigated as well. It allows for a more precise comparison with experimental data. In Fig. 3 we present the results obtained for the critical capillary number as a function of the viscosity ratio in the particular case of an axisymmetric geometry, which allows for a reduction of the problem to a two-dimensional resolution. One must note that the 3D axisymmetric elongational flow is defined by: v x ext = γ x and v r ext = -γ 2 r, where r denotes the radial coordinate normal to the elongation axis x. This differs from the usual experimental situation corresponding to the four roll apparatus [4] which indeed corresponds to a two-dimensional elongational flow. The pure 2D AF resolution is indicated in Fig. 2 for comparison. The dimensionality seems to have a strong influence on the results, although we can note a good qualitative overall agreement between models and experimental findings. Interestingly, the largest deviations are observed in the low viscosity ratio regime, where both the small deformation theory and the AF method fail to reproduce the observed growth of the critical capillary number. This failure seems to be due to the particular shape of the droplets in these extreme regimes. Indeed, as can be seen in Taylor's paper [2], droplets tends to develop pointed ends, that can hardly be accounted for by a small deformation theory, but also by the AF approach. This failure of the AF approach is due to the width ǫ of the interface, which is a natural cutoff for the local curvature radius. Pointed ends cannot be observed in such situations, rather, the drop shape is seen to converge to a limiting shape when the viscosity ratio p decreases to zero. However, this problem only occurs for very low values of the viscosity ratio p < 0.01, and depends strongly on the choice of ǫ. This problem does not occur in the high viscosity ratio regime since shapes are relatively ellipsoidal in this regime, so that the curvature radius is always large compared to ǫ.

In Fig. 4 we show a comparison (log-log scale) between 2D-AF results obtained for a pure shear flow corresponding to v x ext = γ y and v y ext = 0 with experimental data. Also shown are the AF results (2D and 3D) for elongational flows. A comparison between the 2D AF results for a shear flow and experimental data can only be qualitative as results from the previous discussion, but we can observe the same trend in the large viscosity ratio regime: for p larger than 3.5 droplets cannot be broken anymore. The shear situation differs then significantly from the pure elongational case. In fact, a shear flow can be separated into two components, a purely elongational flow which axes are oriented at 45 degrees from the shear direction, and a purely rotational flow. While the elongational component tends to elongate the drop, the rotational component deflects it from the main elongation axis. Since the rotational torque felt by the drop increases with the viscosity ratio (a purely rigid ellipsoid would have a flipping motion in a shear flow [START_REF] Jeffery | Proc. R. Soc. Lond. A[END_REF]) elongation becomes very inefficient, and breakup does not occur anymore. This trend is perfectly visible in Fig. 4 for both the 2D AF model and experimental results. Full 3D simulations in the shear case should give better quantitative agreements with experimental data, as it is the case in the elongation situation, but although the numerical implementation is not sensitive to dimensionality, the computational time scales like N d where N is the number of grid points in a given direction and d the dimensionality. A systematic study then requires much more computational time and is not presented here.

Conclusion. -Apart from the very specific regime where pointed ends can occur (very low viscosity ratios) , all these results show that the AF model captures perfectly well the physical features involved in the drop deformation problem. First the burst transition which separates the existence of stationary shapes from the unsteady regime where the drop elongates is well described by the numerical model, as shown in 2D in Fig. 2. Moreover the method is able to account for topological changes an gives access to the full dynamics of the breakup process. Comparison with experimental data gives a fair agreement, in view of the scattering of experimental data and the dimensionality difference and sensitivity. We considered here a single droplet submitted to an external flow, and the same model can be used without modification to model a set of droplets (only the initial condition needs to be changed). It must be noted that this method also captures the dynamics of the fusion process, as will be shown elsewhere, since the short range attraction induced between the interfaces by the AF method is able to overcome the hydrodynamic repulsion, which is a problem in pure hydrodynamic approaches. Unlike integral formulations, which explicity take advantage of the linearity of the Stokes equation for Newtonian systems, the AF method can be used to model non-Newtonian fluids, and we believe that this method will be helpful to model problems like the dynamics of multiphase flows. We must also mention that the AF method is not restricted to binary systems, the AF method has already been applied to investigate the dynamics of more complex interfaces such as phospholipid membranes. An application to vesicle dynamics can be found in [START_REF] Biben | [END_REF].

Fig. 1

 1 Fig. 1

  Fig. 2

Fig. 1 -

 1 Fig.1-Breakup sequence of an initially elongated droplet in 3D for p = 1, the AF method allows to follow the full dynamics of the breakup. The grid size is 800x200.

Fig. 2 -

 2 Fig. 2 -Comparison between the AF method and the theoretical predictions by Buckmaster and Flaherty in 2D for drop deformation in an elongational flow. The grid size is 300x300

  Fig.3

Fig. 4

 4 Fig.4

Fig. 3 -Fig. 4 -

 34 Fig. 3 -Critical capillary number as a function of the viscosity ratio for an elongational flow. Squares and triangles correspond to experimental data while the continuous line is the AF 3D prediction, and the dashed line is the 2D prediction. The small deformation theory corresponds to the dot-dashed line. The grid size is 300x300 Fig. 4 -Same figure for elongational (bottom curves) and shear (top curves) flows. Squares and triangles are experimental data for elongational flow, while the AF predictions are the straight line in 3D and the dashed line in 2D. The dot-dashed line indicates the average location of experimental data for shear flows, while the dashed line with circles corresponds to the AF prediction in 2D.

c EDP Sciences


INTRODUCTION

The knowledge of the adhesive properties requires the determination of its stickiness using standard tests. This can be done using a double cantilever beam test [START_REF] Boucher | Enhanced adhesion between polypropylene and polyamide-6 : role of interfacial nucleation of the β-crystalline form of polypropylene[END_REF], a modified JKR test [START_REF] Deruelle | Adhesion at the solid-elastomer interface : influence of the interfacial chains[END_REF] and other kinds of experiments like tack and peel. The tack test consists in bringing a solid probe into contact with a pressure sensitive adhesive (P.S.A.) for a short time and then withdrawing it [START_REF] Zosel | Adhesion and tack of Polymers: Influence of Mechanical Properties and Surface Tensions[END_REF]. By varying the pulling velocity, the pressure and the contact time, this test allows to measure the maximum force or energy to break the bond. Peel tests, on the other hand, are carried out after forming a bond by setting an adhesive (with backing) onto a pretreated surface. Then the adhesive/backing compound is pulled away at constant angle (usually 90°) and the energy restitution rate (force per unit width) is obtained. Special peel tests can provide data over a large range of velocities or reduced velocities [START_REF] Kaelble | Theory and analysis of peel adhesion : mechanisms and mechanics[END_REF], as long as the time temperature superposition principle holds for the adhesive considered Surface treatments are also necessary before carrying out tests, to avoid pollution and to provide reproducible surfaces. Different types of treatments are available such as mechanical, chemical, flame, plasma or UV treatments, all depending on the nature of the surfaces. Indeed, chemical treatments will be more adequate for metallic surfaces while Corona discharge or flame treatments are preferred for solid polymer surfaces [START_REF] Brewis | Industrial adhesion problems[END_REF]. As a result, various techniques are used to evaluate the effect of surface changes after treatment: X-Ray Photoelectron Spectroscopy (XPS) is a method for investigating the different components or elements present at the surface, as well as Auger electron spectroscopy, ion scattering spectroscopy or even acoustic scanning microscopy [START_REF] Wu | Polymer interface and adhesion[END_REF]. Free surface energies are also very important and can be determined by conventional wetting measurements with different liquids, using dispersive and polar decomposition [START_REF] Fowkes | Additivity of intermolecular forces at interfaces. I. : Determination of the contribution to surface and interfacial tensions of dispersion forces in various liquid[END_REF][START_REF] Owens | Estimation of the surface free energy of polymers[END_REF]) or acid-base decomposition [START_REF] Van Oss | On the predominant electron-donicity of polar solid surfaces[END_REF]. Many works, theoretical as well as experimental, have focused on peel adhesion. In particular, the general shape of the curve showing the energy release rate G as a function of velocity V depends on the type of adhesive and upon the polymer degree of cross-linking.

Generally, one finds one branch for cross-linked adhesives, corresponding to interfacial peeling (low surface energy) or cohesive peeling (high surface energy) [START_REF] Verdier | Peeling of acrylic pressure sensitive adhesives : cross-linked versus uncross-linked adhesives[END_REF]. In the case of uncrossed-linked adhesives, which is more common, there are two branches, a cohesive and an interfacial one.

At low velocities or high temperatures, i.e. the cohesive regime, Benyahia et al. [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF] observed the mechanisms of peeling, and obtained the general shape of peeling master curves in a log-log plot (Fig. 1). One can recognize three regimes: at very low velocities (regime -3 -C0), G is a roughly a constant. After a typical velocity V1, a second regime (C1) shows increasing G with a slope of 0.6, corresponding to a two-dimensional flow of the adhesive.

Finally, regime C2 starts at velocity V2, with G increasing with a slope 0.8, and this threedimensional flow is characterized by the formation of regularly spaced filaments.

At very low speeds (regime C0), several authors tried to relate the plateau for G with the physico-chemical aspects. Gent and Schultz [START_REF] Gent | Effect of wetting liquids on the strength of adhesion of viscoelastic materials[END_REF] showed that the peel strength G is a product of two terms, one depending on surface energies, the other one being a rheological function φ( ) a V T of the adhesive due to the energy dissipated irreversibly in the bulk. This function should go to 1 when a T V goes to 0 :

( ) G G a V T = 0 φ( ) (1)
where a T is a function of temperature and a T V is the reduced peeling velocity. When failure is at the interface, G 0 is close to the Dupré work of adhesion ω γ γ γ is a molecular dissipation factor related to the irreversible deformation of bonds in the chains between two cross-links and M c is a homogeneous distribution of cross-links. One of the questions which remains is the generality of this formula, as well as the zero velocity limit, if measurable. Finally, the determination of φ( ) a V T is still under way but remains a major difficulty.

In regime C1, there are not many studies dealing with the prediction of cohesive flow.

The major work is by Kaelble [START_REF] Kaelble | Theory and analysis of peel adhesion : mechanisms and mechanics[END_REF], using Bikerman's works [START_REF] Bikerman | Theory of peeling through a Hookean solid[END_REF][START_REF] Bikerman | Rheology of peeling in a Newtonian liquid[END_REF], who proposed a theoretical model when the adhesive is Newtonian.

When two-dimensional peeling becomes unstable (regime C2) at velocity V2, the instability gives rise to regularly spaced ribs. This regime is comparable to hydrodynamic instabilities previously described [START_REF] Pitts | The flow of thin liquid films between rollers[END_REF][START_REF] Mill | Formation of the ribs on rotating rollers[END_REF], when investigating the flow of Newtonian fluids between counter rotating rollers. Empirical formulae predict the velocity of appearance of the ribs as well as their periodicity. Such formulae may be applied to the adhesive peeling case.

Regime C2 has been the subject of numerous works. Kaelble predicted experimentally the distribution of cleavage stresses [START_REF] Kaelble | Peel adhesion : rate dependence of micro fracture processes[END_REF] at the adhesive-substrate interface and showed that the contribution to peel force by filamentation is significant. Niesiolowski and Aubrey [START_REF] Niesiolowski | Stress distribution during peeling of adhesive tapes[END_REF] introduced the filamentation effect at the peeling front and measured the backing curvature to determine peel force. Zosel observed the filamentation mechanisms in a tack experiment [START_REF] Zosel | Adhésion et tack des polymères. Influence des propriétés mécaniques et de la structure moléculaire[END_REF] and pointed out the importance of the elongational properties of adhesives. Gent and Petrich [START_REF] Gent | Adhesion of viscoelastic materials to rigid substrates[END_REF] approximated the deformation of an adhesive (during peel) to a uniaxial elongation, and predicted peel force from extensional rheometry data. Verdier et al. [START_REF] Verdier | Influence des propriétés élongationnelles dans le pelage des polymères[END_REF] finally proposed a model to determine the peeling energy using both shear and elongational rheometrical data as well as physical observations.

In this paper, we compare experimental peeling data with an analytical model based on lubrication assumptions, allowing the different parameters to vary (adhesive thickness, viscosity, surface tension, and backing elasticity). Then we investigate the transition between C1 and C2. In the first part, we derive a theoretical relationship between the energy restitution rate G and the peeling velocity V as well as the other parameters (regime C1). In the next part, we give details of the experimental setup and results, in particular peeling curves.

Finally, we use our results and previous data to obtain a reduced peeling master curve based on our theoretical analysis. The criterion for predicting the transition from C0 to C1 is also presented.

THEORETICAL PEELING MODEL

A schematic view of the adhesive joint in the peeling zone is represented in Fig. 2. Then we will make the following hypotheses:

a) The backing is a flexible elastic strip (much longer and wide than thick) whose deformation is assumed to be one-dimensional and governed by standard beam theory.

b) The adhesive is Newtonian and undergoes two-dimensional cohesive failure at constant temperature. c) Lubrication theory for the flow of the adhesive is assumed, i.e. the velocity distribution u(y) is one-dimensional and depends on the y-coordinate only.

Under these assumptions, a rather lengthly calculation can be carried out, whichis presented in the appendix 1, and one recovers the force F, or the energy restitution rate G, as a function of the relevant parameters

        η = = - 7 4 7 4 7 6 7 9 7 3 V e h E K b F G (2)
where b is the adhesive strip width, E the elastic modulus of the backing, h the half-thickness of the backing, η the adhesive viscosity, V is the peeling velocity and K is a constant.

-5 -Further comparisons with other previous works can now be made. The first attempt to predict peeling force was made by Bikerman [START_REF] Bikerman | Theory of peeling through a Hookean solid[END_REF], who proposed a relationship linking the force required to peel a Hookean adhesive between a rigid plate and a flexible ribbon. He found:

4 1 4 3 4 1 Y E h e 0.3799 b F       = σ ( 3
)
where Y is the Young modulus of the ribbon and σ the tensile strength of the adhesive layer. This is of course rather out of the scope of our analysis, but still an interesting result. Kaelble indeed used Bikerman's equation and wrote a general relation for the peel force when the adhesive is Newtonian [START_REF] Kaelble | Theory and analysis of peel adhesion : mechanisms and mechanics[END_REF]:

F K t b e EI m =       12 1 2 η ( 4
)
where K t m is a peel rate factor and t m , a typical time, is determined by the rate at which the viscous fluid can follow the movement of ribbon. This equation may also be rewritten (using the formula for I) :

F b E h e K t m =             2 2 3 1 2 1 2 1 2 η (5)
It is now in a convenient form to be compared with (2), after we make an assumption for K t m .

From [START_REF] Kaelble | Theory and analysis of peel adhesion : mechanisms and mechanics[END_REF], it seems reasonable to assume that t m varies like V e , the typical experimental time scale. This would lead to a power law dependence of (ηV) with exponent ½ (compared to 4/7 in Eq. ( 2)) and for the thickness e, the exponent would be -1 (to be compared with -6/7 in Eq.

(

) 2
). Other parameters also give reasonable comparisons. To validate our model, we will now compare our scaling law with the experimental results.

EXPERIMENTAL MEASUREMENTS

Experimental setup

We start here with a brief description of the materials which are used throughout this study. We used two PDMS (polydimethylsiloxane), A and B, as adhesives. Measurements using another adhesive (C) will also be reported and are taken from a previous study [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF].

Regarding backings, three different ones are used, one made of Polyester, and two out of steel with different thicknesses. The characteristics of the adhesives and backings are reported respectively in Tables 1 and2. The main advantage of using PDMS is that it behaves as a Newtonian fluid in the range of velocities considered here. Moreover, precise rheometrical data has been obtained previously, as given by dynamic shear master curves covering ten decades in frequency [START_REF] Longin | Dynamic shear rheology of high molecualr weight polydimethylsiloxanes :comparison of rheometry and ultrasound[END_REF].

Adhesives are coated using a blade device. This apparatus allows us to control thickness. PDMS is applied first on the backing, placed horizontally. Then the backing is pulled slowly under a blade. This blade is held fixed on calibrated magnets, which allow maintaining a well-determined gap. After coating is achieved, we measure the weight of the assembly to check the adhesive thickness. Adhesive strips, obtained using this process, are then smoothly applied onto a Pyrex™ surface. The latter one has been treated carefully by successive degreasing, sulfochromic etching, rinsing and vacuum drying. Note that this procedure is standard in our tests but is not fundamental here since failure is cohesive, i.e. not located at the interface.

The 90° peeling experiments have been described precisely by Benyahia et al. [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF]. The system is made of two micrometric tables moving at the same controlled velocity V in the horizontal and vertical directions. The substrate is fixed on the horizontal table and the force transducer (measuring the force F) is fixed on the vertical one. The adhesive strip is set on the substrate under controlled pressure and attached to the transducer so that the peeling front is always fixed in the laboratory reference frame. The whole system is put into a chamber where temperature (T) is controlled as well as the relative humidity (RH).

For the very low speeds, we had to modify our instrumentation because of the low values of the force. We replaced the force transducer (sensitivity 0.1g) by an electronic weighting unit (sensitivity 0.001g) for increased accuracy.

Peeling curves

All experiments are carried out at a temperature of 20°C and a relative humidity around 10 %. After setting up the adhesive strip, we wait for two hours until equilibrium is reached.

Peeling starts by moving the two tables at the same moderate velocities, in order to create the adhesive free surface at the peeling front. The system is then stopped. After a short time, we start a new peeling at the required velocity V and follow the force signal. The latter one shows a transient regime first, which corresponds to deformation of the free surface, and then it reaches a constant value representing the peel force F(V). The time required to reach this steady state regime is inversely proportional to the peeling velocity and depends on the initial conditions i.e. the shapes of the peeling front and backing. At very low speeds, when long times are required, the force signals are not constant but increase slightly linearly due to the increase of adhesive/backing weight, so we need to apply a weight correction to the measured force. The energy release rate G is plotted as a function of velocity V :

G(V) b F G = = (6)
where b is the width of the adhesive strip.

As discussed above, the conditions at the beginning of a peeling experiment are very important especially at very low speeds when the energy restitution rate is of the same order as the free surface energy of the adhesive. In order to validate this principle, two types of experiments with different initial conditions were conducted (Figs 3a-b). In case (a), we take into account the creation and deformation of the free surface at the peeling front whereas in case (b), we only have a deformation of this free surface. The same behaviors for cases (a) and (b) (Fig. 4) are observed for the high velocities, as expected. At low velocities, on the other hand, there are differences, due to the fact that the creation of the free surface (a) is energy consuming. In the experiments that will follow, we will retain case (b) which corresponds to reproducible free surfaces obtained by peeling at the same moderate speed for each series of tests, before measuring the peel force at the desired velocity. Let us note that the value of the magnitude of the so-called "plateau" in the low velocity regime (Fig. 1) is hard to obtain precisely [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF] but is roughly in the range of the adhesive cohesive energy, i.e. 42 mN/m in this case. A different method is still required to be able to reach the plateau limit.

Influence of adhesive viscosity (η η η η)

We first achieve two series of measurements using a polyester backing (P23 in Table 2)

and adhesives A and B (thickness obtained from coating e = 80 ± 7 µm). Figure 5 shows the peeling curves of PDMS on a Pyrex™ substrate. Measurements have been carried out at T = 20 ± 0.1 °C with a relative humidity RH = 10 ± 5 %. For all these experiments, failure takes place in a cohesive manner.

Let us look at the behavior of adhesive A. First of all, at very slow peeling velocities (1.10 -6 m.s -1 ≤ V ≤ 5.10 -5 m.s -1 ), we observe two-dimensional flow at the peeling front (C1) as described previously [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF], this regime is associated with a slope of 0.6 in a logarithmic scale.

As we move along the peeling curve, increasing peel rate, we notice a transition region, i.e. an unstable zone where the two-dimensional flow tries to become a three-dimensional one without reaching steady state. Moreover, the measured force is not constant but decreases with the appearance of ribs and increases when ribs disappear, thus leading to an oscillatory regime of force. As we increase the speed rate (1.10 -4 m.s -1 ≤ V ≤ 8.10 -3 m.s -1 ), beyond the transition area, we clearly observe three-dimensional flow (C2). This flow is associated to a slope of 0.8 in a logarithmic scale.

Let us now observe the results for adhesive B. Regime C1 takes place at peeling velocities between 5.10 -7 m.s -1 and 8.10 -6 m.s -1 (slope of 0.6) and the three-dimensional regime (C2) at peel rates between 2.10 -5 m.s -1 and 2.10 -3 m.s -1 (slope 0.8). Between these two regimes we find a similar transition region. The differences between these two adhesives obviously come from their different viscosities.

Influence of the backing rigidity (EI)

In this part, we study the influence of the backing rigidity and deal with experiments involving adhesive B and different types of backings (see Table 2). Figure 6 shows the peeling curves on a Pyrex™ substrate. Measurements have been carried out at T = 20 ± 0.1 °C and RH =10 ± 5 %. In these peeling tests, the failure mode is still cohesive.

The behavior obtained with the S20 backing is somewhat different from the other one (P23). At low velocities, we observe the typical C1 regime, but the transition has moved.

After a short unstable regime, the adhesive does not obey a linear slope in the G(V) diagram.

The linear regime with slope 0.8 is only obtained at higher velocities. This behavior will be discussed later in §4.

Influence of adhesive thickness (e)

The influence of the adhesive thickness (e) is studied by carrying out experiments using adhesive A coated on the polyester backing P23 using different thicknesses. In addition to our previous peel tests with thickness e=80µm, we used another thickness (e=160µm). These two curves (with a Pyrex™ substrate) are shown in Figure 7. Measurements have been carried out at T = 20 ± 0.1 °C and RH = 10 ± 5 %. The behavior of the adhesive strip with thickness e=160µm is the same as the adhesive strip e=80µm but corresponds to lower levels of forces.

Usually, it is found that the peeling energy increases when decreasing thickness (see Eq. ( 5)), in the case of Newtonian adhesives. This is precisely why small thicknesses are used with common glues in order to increase adhesion.

In the next section, we will now compare theory with experiments, in order to come up with a master curve involving the different parameters above.

Analysis and discussion

Peeling master curve

Figure 8 shows the peeling master curve which is constructed by superposition of peeling curves obtained with adhesives A, B with different backings and different thicknesses. We have also included previous data, obtained by Benyahia et al. [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF]. These measurements were carried out with adhesive C (see Table 1) using the same polyester backing but with a different thickness (e=100µm) at a temperature of 50°C. We know that adhesive C has a Newtonian behavior in the low velocity range at 50°C. This is not a problem to superpose different results at different temperatures, because in this case, the only important parameters which depend on temperature are the viscosity (η) and the surface tension (γ) of the adhesive. Therefore all we need to do is to use the relevant parameters (η,γ) at the required temperatures.

In this peeling master curve, we plotted ξ 3/7 G as a function of ηV, after using a new variable ξ η =

V e E h 2 3 . ξ represents the viscous effects of adhesives versus the backing rigidity effects. This parameter can be introduced naturally when writing the equilibrium equation for the backing, considered as a beam, submitted to a vertical force, and in contact with the adhesive. The curves obtained superpose quite well in the C1 regime (ηV<10 -2 N/m), but in the second regime C2 (ηV>10 -2 N/m) they exhibit differences. Regime C1 differs from regime C2 by a different slope. Regime C1 is characterized by a slope of 1 in a logarithmic scale, according to our model. Indeed, Eq. ( 2) may be rewritten:

ξ 3/7 G = K η V (7)
Concerning regime C2, superposition of the three adhesive experiments is still good but our scaling law is no longer useful because of the three-dimensional flow and the model hypotheses. In this case, we should take into account the rib formation, where the elongation effects become dominant [START_REF] Verdier | Influence des propriétés élongationnelles dans le pelage des polymères[END_REF]. The constant K may actually be a function of the adhesive surface tension γ in this case. This is the only parameter which does not show in the model, and it will be introduced now. Moreover, to obtain a dimensionless representation, we have introduced the surface energy through the capillary number (Ca). For a careful analysis of the relevant dimensionless parameters, one can refer to Piau et al. [START_REF] Piau | Influence of rheology and surface properties in the adhesion of uncross-linked pressure sensitive adhesives[END_REF]. We can then propose:

K(Ca) ηV G ξ 7 3 = ( 8
)
where K is now a function of Ca V = η γ , the capillary number. This representation is shown -10 -in Figure 9. In the high velocity regime, the dependence on Ca is complex, involving minima, so we will limit ourselves to determine the average constant K in the low velocity regime (regime C1): K= 0.2 is the appropriate number. Finally, for the small velocities corresponding to our low-velocity model, we find:

( ) 7 4 7 3 2 3 ηV e Eh 0.2 G       = (9)
It is now needed to write this relationship in terms of the dimensionless numbers ξ and Ca.

Taking into account the previous observations showing the importance of the creation of a free surface (Fig. 3), we may assume that this creation (two surfaces in fact) requires an energy per unit surface equal to w c =2γ, so that the total energy G t becomes [START_REF] Verdier | Peeling of acrylic pressure sensitive adhesives : cross-linked versus uncross-linked adhesives[END_REF] In this final expression, we are showing the two relevant dimensionless numbers of our problem, ξ and Ca. This formula is valid in the C1 regime and holds as long as there are no cross-linking effects or physico-chemical interactions, which might create changes in the value of w c , as was suggested previously [START_REF] Carré | Polymer-aluminium adhesion II Role of adhesive and cohesive properties of the polymer[END_REF]. We believe that this relation is universal.

                      + = γ ηV ηVe Eh 0.1 1 w G 7 3 2 3 c t

Critical conditions between C1 and C2

Let us now pay attention to the transition velocity from C1 to C2. Pitts and Greiller [START_REF] Pitts | The flow of thin liquid films between rollers[END_REF] measured conditions under which the ribs start to appear, in the case of a Newtonian fluid moved by counter-rotating rollers (radius R). They changed speeds (V), viscosities (η), gap width (2e) and used different liquids to investigate the effect of surface tension (γ). They obtained the following empirical formula at criticality, relating the two relevant dimensionless numbers, i.e. the capillary number Ca and a geometrical parameter e/R: R e 62 γ ηV =

Mills and South [START_REF] Mill | Formation of the ribs on rotating rollers[END_REF] studied theoretically this problem, and found that the onset of ribbing was governed by the relation:

0.75 R e 10.3 γ ηV       = (12)
Both [START_REF] Benyahia | The mechanisms of peeling of uncross-linked pressure sensitive adhesive[END_REF] and [START_REF] Gent | Effect of wetting liquids on the strength of adhesion of viscoelastic materials[END_REF] were obtained using constant radius R. In our case, we may consider that the problem is similar, because it would correspond to the flow created by two rollers with different radii, one equal to R, and the other one being infinite. R is then the proper length.

Furthermore, R is not a constant value here and depends on the different parameters, in particuler the velocity V. R can be obtained experimentally by measuring the radius of curvature for ω=0 in Fig. 2. A circle, parabola (Eq. (A5) in appendix 1) can be used to fit the backing as well as the analytical value (Eq. (A2) in appendix 1). These methods all give similar numbers for R close to ω=0 [START_REF] Ravilly | Adhésion et pelage bidimensionnel des polymères[END_REF]. From theses values, at each different velocity, one can compute e/R and plot Ca= ηV/γ as a function of e/R, as shown in Fig. 10. This is done using a log-log scale and a resulting slope of about 0.5 is obtained :

0.5 R e 1.8 γ ηV       = (13)
This relation differs slightly from ( 11) and ( 12) but the exponent is similar and one has to realize that, in the experiments performed earlier, it has been rather difficult to have an accurate observation of the appearance of these ribs [START_REF] Coyle | Stability of symmetric film-splitting between counter-rotating cylinders[END_REF]. Therefore this may induce small changes in the critical velocity. Finally we conclude that our criterion (13) seems rather well adapted to the appearance of ribs in cohesive peeling of Newtonian adhesives.

Conclusion

In this paper, we have focused on the effects of adhesive viscosity, thickness and surface tension, as well as backing rigidity, on the peeling energy required to carry out a 90°peel test of a Newtonian adhesive. In the low velocity regime, we found a master curve, which contains two dimensionless parameters relevant to the problem. The higher velocity case is not so simple due to the appearance of regularly spaced ribs and depends on the capillary number in a more complex manner.

The transition from the two-dimensional regime to the three-dimensional one, associated with the formation of ribs, is governed by a simple criterion relating the capillary number to a geometrical parameter, and compares well with experiments carried out with counter-rotating rollers.

These relationships may be useful to predict peel adhesion between Newtonian adhesives and may be considered as a starting point in the adhesion of more complex polymeric adhesives.

The surface flow rate of the adhesive layer is:

Q V e u y dy dp dx H x V H x H x = = =       -         + ∫ ( ) ( ) ( ) ( ) 0 3 1 2 6 η (A7)
This requires the pressure gradient to vary with H(x) according to the relation:

dp dx V H x e H x = -         12 3 η ( ) ( ) (A8)
Integrating (A8) from 0 to L, after using (A5), we obtain:

( ) ( ) ∆p V RL L eR e eR L R e e Arc L eR = - + +                 3 2 2 2 2 2 2 2 η tan (A9)
Assuming that L eR >> 2 (usually true in the experiments), (A9) simply reads:

∆p k V R e e = η 2 (A10)
where k is a constant, and V is the peeling velocity. Combining the equations (A2)-(A3)-(A10), we obtain a scaling law, linking force and peeling rate:

        η = = - 7 4 7 4 7 6 7 9 7 3 V e h E K b F G (A11)
with K a constant. This equation, which links the energy release rate (G) with the peeling rate (V), takes into account the geometrical parameters, including the thickness of the adhesive layer (e) and the half thickness of the backing (h), the Young modulus of the backing (E) and the viscosity of the adhesive (η).

Figure captions

  γ s are respectively the free energies of the adhesive and substrate, and γ as the interaction energy. In cohesive failure, ω a is replaced by the cohesive energy ω γ c a = 2 . Carré and Schultz[START_REF] Carré | Polymer-aluminium adhesion II Role of adhesive and cohesive properties of the polymer[END_REF], after studying the separation of aluminium/elastomer model assemblies in air and liquid media, proposed to replace G 0 by ( ) reversible energy of cohesion measuring the energy required to break all physical and chemical bonds per unit area or entanglements). g(M c )

Figure 1 :

 1 Figure 1: General shape of the energy release rate G vs. peeling speed V in the cohesive failure mode.

Figure 2 :

 2 Figure 2: Schematic cross-section of an adhesive joint in the peeling zone.

Figure 3 :

 3 Figure 3: Schematic situations of the peeling front at the beginning of a peel test: (a) Creation and deformation of free surface. (b) Deformation of free surface only.

Figure 4 :

 4 Figure 4: Peeling curves of adhesive A (e = 80 µm) with different initial conditions (a) and (b) on Pyrex™ substrate using polyester backing P23 (T=20°C, RH=10%).

Figure 5 :

 5 Figure 5: Comparison of peeling curves of adhesives A and B (e = 80 µm) on Pyrex™ substrate using polyester backing P23 (T=20°C, RH=10%).

Figure 6 :

 6 Figure 6: Peeling curves of adhesive B on Pyrex™ substrate with different backings (T=20°C, RH=10%).

Figure 7 :

 7 Figure 7: Peeling curves of adhesive A with different thicknesses on Pyrex™ substrate using polyester backing P23 (T=20°C, %RH=10%).

Figure 8 :

 8 Figure 8: Peeling master curve (T=20°C, RH=10%) ξ 3/7 G vs. ηV

Figure 9 :Figure 10 :

 910 Figure 9: Other form of peeling master curve: ξ 3/7 G/ ηV vs. Ca (T=20°C, RH=10%)

  FIGURE 1

FIGURE 3 - 20 -

 320 FIGURE 3

  FIGURE 4

  FIGURE 5

  FIGURE 6

  FIGURE 7

  FIGURE 8

  FIGURE 9

  FIGURE 10

Table 2 :

 2 Backings characteristics

	TABLES

Appendix 1 : Peeling model for a Newtonian adhesive

We consider a strip of adhesive (b is the width and L the length of the front), constituted by the backing (thickness 2h) and the adhesive layer (thickness e). At one end we apply a traction force F under a constant angle of 2 π (Fig. 2). It follows from beam theory that in this case:

are respectively the modulus of elasticity and the moment of inertia of the backing, m is the moment arm of the peel force, d ds ω is the curvature at curvilinear abscissa s. Following [START_REF] Kaelble | Theory and analysis of peel adhesion : mechanisms and mechanics[END_REF], we obtain the radius of curvature R, at ω = 0°:

We consider a piece of flexible backing and equilibrate pressures on each side with the traction force F, thus we find : Integration of (A4) leads to:

( )


Anne Raingeard

Frédéric Lacquement

Bruno Tourlière

Hélène Julien Baptiste

Methodology for the predictive mapping of flat terrains within a catchment, by semi-automatic analysis of the Digital Elevation Model

Methodology for the predictive mapping of flat terrains within a catchment, by semi-automatic analysis of the Digital Elevation Model

Anne Raingeard 1 , Frédéric Lacquement 1 , Bruno Tourlière 1 , Julien Baptiste 1 , Hélène Tissoux 1,2 1 BRGM, Orléans, France. 2 UMR7194, Paris, France

Abstract

The knowledge of surface geological formations over the metropolitan French territory is very much heterogeneous: sometimes well represented for their contours and descriptions, sometimes absolutely absent when they can be observed on the field. The implementation of systematic and continuous acquisitions programs at appropriate scales all over the territory cannot be intended considering vastness of the country. Therefore, finding solutions to help the mapping of subsurface geological objects in order to increase the knowledge is likely to be a relevant strategy in which the BRGM has been committed for some years.

One of the research approach consists in defining methodologies for predictive mapping by image processing and analysis with validation through field feedbacks. As part of the cartographic identification of alluvial formations, we developed a methodology based on the Digital Elevation Model (25m precision from IGN (French geographical service)) and its declinations (curvature, slope, its difference with base level etc.). The implemented workflow provides a predictive map representation of morphological terrace levels as well as the colluvium from the valley bottom or lower hillside.

We developed this methodology of semi-automatic mapping by topography analysis on a selected test area in the Pyrenean foothills (sub-watershed of the Gave de Pau). A field campaign on the basins of Baïse and Ousse rivers allowed us to compare the results of the predictive cartography and the geological reality. The first results show that the proposed analysis is relevant at the scale of a multi-kilometer catchment, not affected by faults. The methodology provides valuable support for the mapping of these superficial formations, in particular by significantly saving time when drawing the map. It was noted that the estimation of the alluvial system base level was a key point in our modeling, then several interpolation methods were compared.

The development of a graphical interface allowing all the catchment analysis to be carried out in a single application will make it possible to produce predictive maps quickly and to help with the cartographic identification of alluvial stepped formations


T Guillot
email: guillot@obs-nice.fr

The composition of transiting giant extrasolar planets

Keywords: numbers: 95.10.Gi, 96.12.Ma, 96.15.Bc, 97.10.Tk, 97.82.-j, 97.82.Fs

In principle, the combined measurements of the mass and radius a giant exoplanet allow one to determine the relative fraction of hydrogen and helium and of heavy elements in the planet. However, uncertainties on the underlying physics imply that some known transiting planets appear anomalously large, and this generally prevent any firm conclusion when a planet is considered on an individual basis. On the basis of a sample of 9 transiting planets known at the time, Guillot et al. A&A 453, L21 (1996), concluded that all planets could be explained with the same set of hypotheses, either by large but plausible modifications of the equations of state, opacities, or by the addition of an energy source, probably related to the dissipation of kinetic energy by tides. On this basis, they concluded that the amount of heavy elements in close-in giant planets is correlated with the metallicity of the parent star. Furthermore they showed that planets around metal-rich stars can possess large amounts of heavy elements, up to 100 Earth masses. These results are confirmed by studying the present sample of 18 transiting planets with masses between that of Saturn and twice the mass of Jupiter.

Introduction

The discovery of exoplanets in transit in front of their parent star has led to the birth of a new field in science: exoplanetology, defined as the study of the characteristics, formation and fate of planets outside our solar system. I will not venture into trying to detail the possibility offered by these transiting systems, but they are numerous. One will be of importance for us, the possibility to measure relatively precisely both the mass and radius of these planets and thus armed with theoretical models, to constrain their compositions.

The first transiting exoplanet detected, HD209458b (Charbonneau et al. 2000, Henry et al. 2000), was shown to have a radius of 1.35 R J for a mass of 0.65 M J (R J = 71492km is Jupiter's equatorial radius at 1 bar, and M J is Jupiter's mass). This large radius implied naturally that the planet was a gas giant, mainly formed with hydrogen and helium (Guillot et al. 1996, Burrows et al. 2000), and a sign that our understanding of planet formation, badly shaken by the discovery of 51 Peg b, was not completely off-track.

However, HD209458b was then shown to be anomalously large (Bodenheimer et al. 2001, Guillot & Showman 2002), i.e. larger than predicted by the evolution of a solar composition gas with no addition of any heavy elements. In contrast, when other transiting planets were found, they appeared to be "normal" within the error bars (Baraffe et al. 2005, Laughlin et al. 2005), although some indications that they could possess relatively large amounts of heavy elements were found (Guillot 2005). This is until HD149026b was shown to be significantly smaller than expected, requiring the presence of a large amount ∼ 70M ⊕ of heavy elements in its interior (Sato et al. 2005, Fortney et al. 2006).

On the basis of the 9 transiting planets known at the time, Guillot et al. (2006) proposed to explain the structure of all planets with the same hypotheses. This required two elements: invoking a "missing physics" (e.g. energy dissipation, a modification of equations of state or opacities), with the same rules for every planet, and allowing for a variable amount of heavy elements in their interior. Their results indicated a likely correlation between the mass of heavy elements in the planets and the metallicity of the parent star. This was later confirmed by Burrows et al. (2007), on the basis of a sample that had grown to 14 members, and using as "missing physics" an increased opacity to slow the contraction of the planets.

This paper updates and extends the previous work by considering a sample of 18 planets. I purposely leave out two transiting planets with masses that are quite far from the Saturn-Jupiter mass range: GJ436 (Gillon et al. 2007), a ∼ 20M ⊕ planet, and HD147506 (aka HAT-P-2) (Bakos et al. 2007), with its 8.2 M J mass. Other planets are shown in table 1. Section 2 describes the principle of the calculations. I then show on the basis of standard evolution calculations that anomalously large planets are common, with strong implications for which mechanisms can be responsible for the slowing of the contraction. Section 4 examines several hypotheses for the "missing physics" and investigates the mass of heavy elements -stellar metallicity correlation.

Principle and hypotheses of the calculations

Once they have formed, giant planets contract monotonically and quasi-statically. The model used to calculate their evolution is described elsewhere (see Guillot 2005). It assumes that the planets are made of a central core which contributes relatively little to the global cooling, and of an extended, solar-composition envelope which accounts for most of the planetary size and internal energy.

An important problem behind models of the physical evolution of giant planets lies in the initial condition used and in their orbital evolution. Fortunately, initially extended planets cool and contract rapidly when the stellar irradiation is not dominant, so that initial conditions are generally forgotten over timescales of 10 7 to 10 8 years (see Marley et al. 2007). The situation may be more complex for Pegasids because of HAT-P-4 3.6 -6.8 0.24( 8) 1689(60) 0.680(40) 1.270( 50)

[Ko07]
The numbers in parenthesis represent the uncertainties on the corresponding last digits. M Jup = 1.8986112 × 10 30 g is the mass of Jupiter. R Jup = 71, 492 km is Jupiter's equatorial radius.

References: (Charbonneau et al. 2000, Konacki et al. 2003, Bouchy et al. 2004, Pont et al. 2004, Torres et al. 2004, Alonso et al. 2004, Sozzetti et al. 2004, Sato et al. 2005, Bouchy et al. 2005, Winn et al. 2005, O'Donovan et al. 2006, Collier Cameron et al. 2006, Knutson et al. 2007, Gillon et al. 2006, Charbonneau et al. 2006, Holman et al. 2006, Shporer et al. 2007, Winn, Holman & Fuentes 2007, Winn, Holman & Roussanova 2007, Bakos et al. 2007, Burke et al. 2007, O'Donovan et al. 2007, Mandushev et al. 2007, Torres et al. 2007, Pont et al. 2007, Gillon et al. 2007, Minniti et al. 2007, Winn, Holman & Fuentes 2007, Kovacs et al. 2007); The discovery papers are in brackets. The table is taken from F. Pont's site: http://obswww.unige.ch/∼pont/TRANSITS.htm.

their rather slow initial contraction (Guillot et al. 1996), and the possibility that they may have formed at larger distances before being brought close to their stars (Burrows et al. 2000).

To make things simple, I assume that all planets have an initial radius of 2 R J , have formed in situ, and underwent a constant irradiation from their parent star. Although this hypothesis can certainly be questioned, it has a limited impact on the results (e.g. Ikoma et al. 2006). Most importantly, applying the same set of hypotheses consistently to all the transiting planets considered will eventually enable testing them when a statistically sufficient number of such planets will have been discovered.

The assumption of solar composition in the envelope may also be surprising in regard to the fact that our giant planets appear to all be significantly enriched (by factors of ∼ 3 for Jupiter to 30 for Uranus and Neptune over the Solar value). At this point, this is a clear simplification, to be further discussed in § 4.2.

Important ingredients entering the calculations are the equations of state. The well-known EOS by Saumon et al. (1995) is used in the envelope. This EOS includes only hydrogen and helium, but a slightly larger mass fraction of helium (Y = 0.30) is used to mimic the presence of a solar composition of heavy elements. The dense core is assumed to be made of pure rocks, and the EOS is based on the experimental fit by Hubbard & Marley (1989). Unless otherwise stated, the mass of the core is the only parameter to be modified in the calculations.

The opacities are based on a table calculated by Allard et al. (2001). These opacities do not account for the possible presence of clouds. The progressive growth of an external radiative zone governs most of the evolution (in particular contraction) of irradiated giant planets (Guillot et al. 1996). I will come back on how modifications to the opacities can change the analysis of the compositions of Pegasids.

The atmospheric boundary condition used in these calculations is simplified. In the case of the first 9 transiting planets, I use specific atmospheric models calculated on the basis of Iro et al. (2005). For the remaining planets, a simple relation is used:

T 10 = 1.25T eq,0 , (1)
where T 10 is the temperature at the 10 bar level, and T eq,0 is the equilibrium temperature at the planet for a A = 0 albedo (perfect absorption). Figure 1 shows how this approximate relation compares to the more elaborate but 1D radiative transfer calculations.

The difference between the observed transit radius (at ∼mbar levels) and the modelled radius (at 10 bar) has to be taken into account (see Burrows et al. 2003). As in Guillot et al. (2006), we include this radius difference as:

∆R =

RT eq,0 µg ln

P model P transit , (2)
where P model = 10 bar, we assume P transit = 1 mbar, R is the gas constant, µ is the mean molecular weight (we use µ = 2.2) and g is the atmospheric gravity. This relation uses a mean atmospheric temperature equal to T eq,0 . Clearly eq. 2 is an approximation, both from the point of view of the transit pressure, and from that of the atmospheric temperature structure. However, what is most important is that whatever the error made here, it applies to all planets in the same way. Furthermore, I claim that it is small compared to the inherent uncertainties related to atmospheric models (including the possible presence of clouds and their structure, the inhomogeneous heating and the resulting dynamical effects).

Anomalously large planets are common!

Figure 2 shows the masses of heavy elements that are needed to explain the observed radius measurements, using a standard evolution model. The results are plotted as a

Figure 1. Temperature at the 10 bar level as a function of the equilibrium temperature, as calculated with a 1D equilibrium model (Iro et al. 2005) (cases 1 to 9) and when assuming T 10 = 1.25T eq,0 (cases 10 to 18).

function of the heavy element enrichment in the parent star (i.e. 10 [Fe/H] , where [Fe/H] is the metallicity).

As already observed by many authors (Bodenheimer et al. 2001, Guillot & Showman 2002, Chabrier et al. 2004), HD209458b is anomalously large. In figure 2, this implies that the planet appears with a negative core mass ‡, which is of course physically impossible. Of interest to us is the fact that a rather large number of planets (HD209458b, TrES-2b, WASP-1b, TrES-3b, TrES-4b, HAT-P-4b, and probably also OGLE-TR-56b, HD189733b, XO-1b, HAT-P-1b), 6 to 10 over 19, appear in the negative core mass side of the diagram.

Clearly, this rules out mecanisms that are statistically unlikely as responsible of the large radii. The small but non-zero eccentricity explanation (Bodenheimer et al. 2001) is already excluded for HD209458b (Deming et al. 2005, Laughlin et al. 2005). The spin-orbit trapping into a Cassini state (Winn & Holman 2005), has been shown to be unprobable (Levrard et al. 2007, Fabrycky et al. 2007). On the other hand, the semiconvection mechanism put forward by Chabrier & Baraffe (2007) cannot be ruled out, but one would expect it to work more efficiently at large metallicities, which seems not Figure 2. Mass of heavy elements in the planets as a function of the metal content of the parent star relative to the Sun. The mass of heavy elements required to fit the measured radii is calculated on the basis of standard evolution models. Negative core masses are required in some cases, implying that some significant physical input is missing (see text). Horizontal error bars correspond to the 1σ errors on the [Fe/H] determination. Vertical error bars are a consequence of the uncertainties on the measured planetary radii and ages.

to be compatible with the results of fig 2.

Altogether, this strengthens the case of a mechanism at work for all planets, as discussed in previous papers (Guillot 2005, Guillot et al. 2006). The rest of the article will be based on this assumption: whatever happens to HD209458b or TrES-4b must also be occuring for all other pegasids.

Explaining the observations

A hard equation of state model

A first possibility is that the hydrogen helium equation of state (EOS) used is not quite right and overestimates the density for a given pressure, temperature and composition (we need hydrogen to be "harder", i.e. less compressible). Rather large uncertainties (∼ 10 -20%) of the density predicted by different equations of state are possible at pressures and temperatures relevant to giant planet interiors (Saumon & Guillot 2004). A 30% softening of the hydrogen EOS compared to the SCVH EOS has been obtained as a result of ab-initio simulations matching high-pressure experiments (Militzer & Hubbard, personnal communication; see also [START_REF] Militzer | AGU Fall Meeting Abstracts[END_REF]). Although this last result goes in the wrong direction for us, it shows that the EOS is an important uncertainty source.

In order to test the possibility that the EOS may be the culprit in the underestimations of theoretical planetary radii, I present in fig. 3 models calculated with the SCVH EOS, but with helium removed. This is similar in essence to a ∼ 30% hardening of the EOS.

As a result, we find that this large change to the EOS is sufficient to explain the present radii of all planets except TrES-4b and HAT-P-4b. However, these two planets have been discovered recently, and still have relatively large error bars. Future measurements are needed for these objects.

Globally, the trend between M Z and [Fe/H] that was seen in the standard model is mostly lost: this effect is an artificial one that is related to the different adiabats of hydrogen and helium. As a result, the change in radius between a pure H planet and a H/He mixture is greater for planets that receive less flux from their star. It shows however that both changes in the density and specific entropies have to be carefully accounted for in new EOSs.

An increased opacity scenario

Another possibility to explain the anomalously large planets is a serious underestimation of the Rosseland opacities. The growth of the external radiative region (down to ∼kbar pressure levels) controls the cooling and hence the contraction of all pegasids (Guillot et al. 1996, Guillot & Showman 2002). This region is in a difficult, relatively cold but dense region, and generally, opacities are notoriously famous for the large uncertainties associated to them.

Figure 4 presents results obtained with Rosseland opacities that have been arbitrarily multiplied by a factor 30 compared to the fiducial values. This value corresponds approximatively to the increase that is needed to explain the radius measured for HD209458b.

It can be seen that this approach leads to physically sound values of the heavy element masses for all planets, with two exceptions: TrES-4b and HAT-P-3b. It may be premature however to rule out this explanation because these discoveries are quite recent, and experience has shown that the inferred radii which are very dependent on stellar properties have often moved a bit more than would be expected from their 1 -σ error bar. It is clear that TrES-4b and HAT-P-3b are extremely interesting planets to test these models. More discoveries will of course also help to settle the matter.

Interestingly, we again observe a clear correlation between planetary and stellar "metallicity". This result, first found by Guillot et al. (2006), was later confirmed by Burrows et al. (2007). A preliminary version of the last paper, put onto astroph, claimed that this opacity increase was related to an oversolar enrichement of the envelope, a simple, logical explanation when considering that the atmospheres of our giant planets are enriched in heavy elements. However this explanation was based on a model for which only opacities had been increased, and not the mean molecular weight. A consistent treatment, unfortunately, rules out this simple explanation as a solution to the anomalously large planets problem.

The reason for that is a balancing effect between the slower cooling and the increased mean molecular weight that implies that enriching the envelopes in heavy elements does not generally make the planets larger, and often results in the oposite effect. Figure 5 shows that the mean molecular weight effect eventually wins, and that this occurs over a timescale similar to the ages of the systems that are observed.

Figure 5. Evolution of giant planets in terms of radius vs. time, for different irradiation levels, and 2 assumed compositions: solar, and 6 times solar. (This calculation ignores second order effects as modifications of the adiabatic temperature gradient and non-linear effects in the opacity calculation, and more importantly modifications of atmospheric properties.) [From Guillot (2005)].

Could the atmosphere or the upper part of the envelope be enriched in heavy elements while the rest would not? This is regarded as extremely unlikely, because a double-diffusive instability ("salt-fingers") should set in (e.g. Stevenson 1985).

Finally, there is an interesting possibility that the envelopes of giant planets orbiting metal-rich stars are more enriched than for planets orbiting metal-poor stars. This case would lead to a stronger correlation between heavy elements in the planet and stellar metallicity.

The kinetic energy mechanism

The third mechansim that would be in action for all the planets, but at a level that depends on the irradiation that they receive from the parent star is the one proposed by Showman & Guillot (2002): a fraction (∼ 1%) of the energy recieved in form of stellar irradiation by the planet is converted into kinetic energy in the atmosphere. This energy is transported to deeper levels, and dissipated probably due to tidal interactions with the star. Although simulations and observations on Earth have shown this amount of kinetic energy to be generated in the atmosphere, several problems remain with this scenario: first the turbulent viscosity needs to be small otherwise dissipation takes place in the upper atmosphere, with little effect on the planet's contraction (Burkert et al. 2005); Second kinetic energy must be transported to the 10-100bar level at least; Third, it must be dissipated there (or at deeper levels).

Figure 6 shows results obtained with this model, assuming a 0.5% dissipation of the incoming stellar flux at the center of the planet. This level of dissipation is chosen so as to yield a small but positive value for the mass of heavy elements inside HD209458b. (Larger values lead to systematically larger masses of heavy elements derived in all transiting planets.) As previously mentioned (Guillot et al. 2006), the relation between the mass of heavy elements in the planet and the stellar metallicity is qualitatively similar than when increasing opacities.

Conclusions

With the 18 transiting exoplanets with masses between that of Saturn and twice the mass of Jupiter known thus far, we are able to confirm the conclusions derived by Guillot et al. (2006):

• All planets can be explained in the framework of a single scenario, when accounting for some missing physics which can either be due to modifications in the EOSs, opacities, or to an additional heating source.

• Giant planets can possess a large mass in heavy elements, up to 100 Earth masses (a value which is itself model-dependent).

• There is a correlation between the mass of heavy elements in the planets and the metallicity of their parent stars, so that the most massive planets generally orbit the most massive stars. This correlation is probably not a single relation, but shows, as would be expected, some scatter. However, there is presently an absence of planets with small masses of heavy elements around metal-rich stars, planets that appear to be predicted by formation models (see Alibert et al. 2005). 2, but the mass of heavy elements required to fit the measured radii are calculated on the basis of evolution models including an additional heat source slowing the cooling of the planet. This heat source is assumed equal to 0.5% of the incoming stellar heat flux (Showman & Guillot 2002).

• Planets around stars that have metallicities comparable to that of our Sun generally tend to have very small cores. There is a clear lack of giant planets at short periods (P < 10 days) around metal-poor stars (Fressin et al. 2007). This seems to imply that metal-poor stars are not able to form giant planets and have them migrate very close-in.

It should be noted however that a planet such as TrES-4 may pose a problem for these models if its large radius is confirmed. It is presently too large to be explained by the models, but the error bars are currently quite large and prevent any definitive conclusion to be made. Future detections of both large and small transiting pegasids are extremely important to test the models and better understand how giant planets formed.

Figure 3 .

 3 Figure3. Same as fig.2, but the mass of heavy elements required to fit the measured radii are calculated on the basis of evolution models of a planet with no helium, to mimic a hard hydrogen EOS (see text).

Figure 4 .

 4 Figure 4. Same as fig. 2, but the mass of heavy elements required to fit the measured radii are calculated on the basis of evolution models with opacities arbitrarily increased by a factor 30.

Figure 6 .

 6 Figure6. Same as fig.2, but the mass of heavy elements required to fit the measured radii are calculated on the basis of evolution models including an additional heat source slowing the cooling of the planet. This heat source is assumed equal to 0.5% of the incoming stellar heat flux(Showman & Guillot 2002).

Table 1 .

 1 Transiting planets included in this study

	Name	Age [Ga]	[Fe/H]	T eq,0 [K] M p [M J ]	R p [R J ]	Refs.
	HD209458	4 -7	0.02(3) 1469(120)	0.69(20)	1.320(25)	[Ch00]Wi05/Kn07
	OGLE-TR-56	2 -4	0.25(8)	2214(72)	1.29(12)	1.300(50)	[Ko03]To04/Po07
	OGLE-TR-113 0.7 -10	0.15(10)	1340(80)	1.35(22)	1.090(30)	[Bo04]Gi06
	OGLE-TR-132	0.5 -2	0.37(7) 2110(100)	1.19(13)	1.180(70)	[Bo04]Gi07
	OGLE-TR-111 1.1 -10	0.19(7)	1033(36)	0.52(13)	1.010(40) [Po04]Sa06/Wi07/Mi07
	OGLE-TR-10	1.1 -5	0.28(10)	1509(80)	0.61(13) 1.122(100)	[Ko05]Po07/Ho07
	TrES-1	2 -6	0.06(5) 1156(140)	0.75(7)	1.081(29)	[Al04]So04/Wi07
	HD149026	1.2 -2.8	0.36(5) 1740(150)	0.33(3)	0.726(64)	[Sa05]Ch06
	HD189733	0.5 -10	-0.03(4) 1199(100)	1.15(4)	1.154(32)	[Bo05]Po07
	XO-1	0.65 -8 0.015(40) 1255(100)	0.90(7)	1.184(25)	[Mc06]Ho06/M06
	HAT-P-1	2.6 -4.6	0.13(2)	1316(61)	0.53(4)	1.203(51)	[Ba07]Wi07
	TrES-2	2.8 -7.8 -0.15(10)	1459(90)	1.28(9)	1.220(45)	[OD06]So07
	WASP-1	0.3 -3	0.26(3)	1849(80) 0.867(73)	1.443(89)	[Ca06]Sh06/Ch06
	WASP-2	0.3 -10	? 1293(136)	0.88(7)	1.038(50)	[Ca06]Ch06
	XO-2	1 -3	0.45(2)	1316(21)	0.57(6)	0.973(30)	[Bu07]
	TrES-3	0.3 -10	?	1644(90)	1.92(23)	1.295(81)	[OD07]
	TrES-4	0.3 -10	?	1757(89)	0.84(20)	1.674(94)	[Ma07]
	HAT-P-3	0.3 -10	0.27(4)	1150(38) 0.599(26)	0.890(46)

[To07]

‡ Note that I did not attempt to calculate hydrostatic models using negative core masses: this is a result of the interpolation (in this case extrapolation) technique used to obtain the core mass. As a result, the precise value of M Z is of little significance when negative.

Acknowledgments

This research was supported by the Programme National de Planétologie. The author thanks an anonymous referee for helpful comments that have improved the manuscript.