Guinot
email: guinot@msem.univ-monpt2.fr

Caroline Sart
email: sart.caroline@laposte.net

Jean-Pierre Baume
email: jean-pierre-baume@cemagref.fr

Vincent Guinot

Adaptation of Preissmann's scheme for transcritical open channel flows Adaptation du schéma de Preissmann pour les écoulements transcritiques à surface libre

Keywords: Fluvial hydraulics, hydraulic jump, numerical model, numerical scheme, open channel flow, transcritical flow, transient flow

Adaptation of Preissmann's scheme for transcritical open channel flows

Introduction

An essential practical aim in river hydraulics is to model hydraulic networks, whether for irrigation enhancement, flood planning or management of urban-drainage. This implies simulations of lengthy events and inclusion of hydraulic structures [START_REF] Baume | SIC: A 1D hydrodynamic model for river and irrigation canal modeling and regulation[END_REF].

Even if the flow is generally subcritical, supercritical regions may appear. The different regimes as well as the transitions between these -bores and critical points -must be numerically represented. Nonetheless, there is no need for a precise localisation of a sharp front. In other words, the expense of sophisticate techniques created for the solution of dam-break problems would be neither justified nor appropriate for our above applications.

For subcritical flows, Preissmann's scheme [START_REF] Cunge | Practical aspects of computational river hydraulics[END_REF]) is one of the most widely used in industrial models. This implicit finite-difference method is unconditionally stable and extremely robust. Moreover, each discrete equation implying only two discrete spatial positions makes easy the inclusion of hydraulic structures that can be represented between two sections, without perturbing the modelling of neighbouring fluid zones.

Unfortunately, while also well suited for fully supercritical flow, it can be demonstrated that in its original usage, the scheme is invalid for transcritical flows, i.e. for flows in which the two regimes coexist [START_REF] Meselhe | Invalidity of Preissmann scheme for transcritical flow[END_REF].

A classical procedure allows the treatment of transcritical flows with a numerical method created for subcritical flows [START_REF] Kutija | On the numerical modelling of supercritical flows[END_REF]). The solution is to apply the method to a degenerated system of equations, reducing the influence of the convective momentum term as the flow becomes supercritical. But this trick does not conserve volumes.

Noticeably, circumventing Meselhe and Holly's invalidity conditions, [START_REF] Johnson | A box scheme for transcritical flow[END_REF] have proposed a method for transcritical flow based on Preissmann's scheme involving the complete Saint-Venant equations. Using the idea of splitting methods [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF]), Preissmann's equations from two adjacent cells are combined. The matrix of the resulting scheme loses the bi-diagonal structure associated with Preissmann's scheme. This structure is highly desirable as it permits solving algebraically the implicit system.

This research focuses on the direct adaptation of Preissmann's scheme for transcritical flows. First, the observations of [START_REF] Meselhe | Invalidity of Preissmann scheme for transcritical flow[END_REF] concerning the invalidity of direct use of Preissmann's method for the solution of transcritical flows using the full Saint-Venant's equations are discussed. Then the interesting features proposed by [START_REF] Johnson | A box scheme for transcritical flow[END_REF] to describe transcritical transitions with a through method (i.e. a method whereby all possible wave configurations are handled within one single mathematical formula, without the need of specific or ad hoc treatments) compatible with Preissmann's scheme are highlighted. The resulting method uses Preissmann's scheme in non-transcritical domains, which are linked by internal boundary conditions characterizing the physics of transcritical transitions. Furthermore, an automatic procedure based on discrete Froude numbers is used to solve the transcritical implicit system. It corresponds to an adaptation of the simple-and double-sweep methods which only require localising the transcritical transitions, thus keeping one of the major advantages of Preissmann's scheme. The method is presented on a transient problem involving appearance and disappearance of a supercritical region between a critical point and a bore. Finally, conclusions on this research are drawn.

2

Generalities on numerical methods for free surface hydraulics

Saint-Venant's equations

Governing equations in conservation form

To model one dimensional (1D) open channel flows [START_REF] Cunge | Practical aspects of computational river hydraulics[END_REF], the depthintegrated equations for conservation of mass and momentum can be written for the variables wetted cross-sectional area A and discharge Q as

) ( ) / ( 0 1 2 f o x t x t S S gA gI A Q Q Q A (1)
Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440.

The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

where g = gravitational acceleration, S o = bottom slope, (gI 1 ) = integral of hydrostatic pressure over the cross-section, i.e. I 1 =A 2 /(2B) for a rectangular channel of width B, friction slope S f =Q|Q| /De 2 , where De =KAR h 2/3 with K = Strickler coefficient and R h = hydraulic radius. Equation (1) corresponds to the conservative form of the equations that is desirable when dealing with transcritical flows. The vector-form for a conserved variable u=(A Q) T of flux f with source terms s is

) ( ) ( u s u f u x t (2)

Jacobean matrix and characteristic form of Saint-Venant's equations

Written in conservative form, the system of Saint-Venant's equations is of hyperbolic nature. This means that it can be diagonalised. Using the Jacobean matrix of the system

1 2 2 2 1 0 d d VDV v v c J u f (3) c v a c v a a a D ) 2 ( ) 1 ( ) 2 ( ) 1 ( ; : ; diag (4)
where 6) can be re-written using total derivatives in the (x,t) plane as

a (i) (i=1,2) = Eigenvalues of J, V = matrix of
) ( ) ( ' d d i i s t w along 2 , 1 d d ) ( i a t x i (7)
introducing the so-called characteristic curves along which perturbations propagate.

The nature of the flow regime is associated with these waves propagating in the domain, and similarly with the value of the Froude number F=v/c representing the relative importance of inertial and pressure governing forces. In a subcritical flow (F<1), pressure forces dominate and the wave velocities (a (1) and a (2) ) are of opposite signs, while a supercritical flow is velocity driven and waves can only travel in the downstream direction.

If there is a transcritical transition, the sign of the celerity of the slowest wave changes (for a positive discharge). A convergence of the associated characteristic lines would imply an over-determination (couples of characteristics issued from different zones giving different values at one point). The solution is then divided into two domains separated by a discontinuity. The Rankine-Hugoniot shock relations (Godlewski and Ravaiart 1991) connect the values on both sides of the discontinuity, defining the shock speed. On the

1 2 2 2 1 0 d d VDV v v c J u f
Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648 contrary, if shifting from a subcritical to a supercritical flow, the transition is continuous;

the direction of propagation of the characteristic associated with the critical point changes continuously from negative to positive, the fan of characteristics is centred on the critical point. There, v=c and the characteristic is steady, namely

0 d d along ' d d ) 1 ( ) 1 ( t x s t w (8)
and for the critical point

c Q A w 2 d d d ) 1 ( (9)
2.2 Classical Preissmann's method for Saint-Venant's equations Discrete system for Saint-Venant's equations Consider a temporal evolution between discrete times (t n =ndt) over a spatial domain represented by N+1 discrete locations (nodes) x 0 <x 1 < …<x N . Preissmann's method applies over a cell [x j ;x j+1 ]x[t n ;t n+1 ]. This finite difference scheme can be interpreted as an approximation of the integral form of the conservation law (2) over the cell

1 1 1 1 d d d )] , ( )) , ( [ d )] , ( ) , ( [ 1 1 j j n n n n j j x x t t j t t j n x x n x t s(x,t) t t x f(u t x f(u x t x u t x u (10)
where the spatial and temporal integrals are given by a linear interpolation using the vertexes of the spatio-temporal domain, with a weighting coefficient for values at t n+1 and =1/2 for values at x j+1 . This relation with the integral form says that the method is conservative, which is essential to obtain the correct shock speeds. Rankine-Hugoniot relations are obtained as a limit of the integral conservation law.

The equations given in Eq. ( 1) discretized with Preissmann's method are formally written as

0 2 / 1 j R (11)
where R=(R A R Q ) T is called "cell residual". These vector relations correspond to the discretized form of the conservation equations for mass ("A") and momentum ("Q")

n n j j j n n j j j j Q j n n j j j j A j j j j j j s s s s t f f f f Q Q R Q Q Q Q A A R 1 1 1 1 2 / 1 1 1 2 / 1 1 1 2 1 2 1 (12)
Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440.

The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

where u=u n+1 u n for the scheme at the known time step n, = t/ x, the flux and source terms of the momentum equation are f = Q 2 /A+gI 1 and s = gA(S o S f ).

Stability and accuracy

While the scheme is only based on two spatial positions, it is second-order accurate in both time and space if = 1/2, but is then marginally stable. Practically > 1/2 is thus used, leading to an accuracy of first order in time and second order in space. A von Neumann stability analysis also indicates that the linearized Preissmann's scheme, for the Saint-Venant's equations in the frictionless horizontal bed case, is marginally stable if F = 1 [START_REF] Meselhe | Invalidity of Preissmann scheme for transcritical flow[END_REF]. This particular configuration occurs only locally for transcritical flows.

Boundary conditions

The set of cell residual's equations (R j+1/2 = 0, j=0,…, N 1 ) has to be completed with exactly two Boundary Conditions (BC) to close the system. If the flow is supercritical both are imposed upstream, whereas for subcritical flow one is imposed upstream and one downstream. This corresponds to giving exactly the BC physically required, namely one at each end where a characteristic line is entering the domain. Other choices for the location of the given BC would not respect the physics of the travelling waves and then the numerical method would not be stable [START_REF] Macdonald | Analysis and computation of steady open channel flow[END_REF].

Note that the scheme uses an indiscriminate space-difference approximation in subcritical and supercritical flow regions, so nothing changes in the scheme whether information is coming from upstream or downstream. But each unknown is implied in two systems of equations i.e. in two boxes. Hence it is appropriate both for fully subcritical and fully supercritical flows. Giving the adapted BCs corresponds to forcing the directions of propagation of the information in the numerical model.

Solution algorithm

Owing to the two-spatial-point formulation of Preissmann's discretisation, the matrix of the implicit system for Saint-Venant's equations has a double-bidiagonal structure allowing for a simple solution algorithmic procedure, depending on the imposed boundary conditions.

Figure 1 shows a schematic representation, from where it can be seen that each node, for the new time step, needs two equations or two "arrows" to be solved.

Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

In the supercritical case, a simple sweep is sufficient. The two boundary conditions give the solution at the upstream boundary. One discrete position after the other, Preissmann's scheme gives the solution down to the downstream boundary node.

In the subcritical case, a double-sweep is necessary. The upstream boundary condition permits simplifying Preissmann's scheme for the first cell. One equation is still linking the two nodes (three discrete values) and the other giving a relation for the values of the second node. This last step is used as the subcritical upstream boundary condition for the second cell. An upper-diagonal matrix is obtained at the end of this first sweep. The downstream boundary condition then allows obtaining the solution of the complete system at the last downstream cell. A second sweep calculates all other discrete unknown values, one after the other, from downstream to upstream using, for each cell, the equation linking the nodes.

Each method reflects the properties of the corresponding flows: either only upstream information is responsible for the evolution, or both upstream and downstream information have an effect.

The in-cell arrows represent the direction of propagation of information associated with the characteristics as well as the two-cell residual's equations for the Saint-Venant system and the "way of using them" once boundary conditions are given down-and/or upstream.

Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648 Even if many methods and algorithms are available for computing numerical solutions of the shallow water equations, including transcritical flows, few have been applied to real life problems [START_REF] Alcrudo | A state-of-the-art review on mathematical modelling of flood propagation[END_REF]. No exhaustive list of all the numerical methods solving transcritical flows is presented below, yet the terms of references are described, namely ability to include hydraulic structures, to model networks, to use relatively large time steps, conservativity and robustness.

A practical method to allow for the modelling of transcritical flows uses a model built for subcritical flows but to work with an approximate system for the physical equations [START_REF] Kutija | On the numerical modelling of supercritical flows[END_REF]. The standard (implicit) method is used but gradually reducing the inertia term in the momentum equation as F 1. Then, the characteristic lines of the associated reduced system keep opposite directions. This allows keeping the algorithmic structure of the initial model, imposing boundary conditions inherent to subcritical flow. This advantageous method is used in several commercial packages combined, for example, with

Preissmann's or Abbott-Ionescu's scheme. Unfortunately, while maintaining robust solution procedures, these approximations may imply accuracy problems and numerical oscillations. A detailed study of these drawbacks in conjunction with the application to hydraulic networks was conducted by [START_REF] Djordjevic | Simulation of transcritical flow in pipe/channel networks[END_REF].

To solve the complete hyperbolic system for discontinuous solutions, the first idea was to use shock fitting methods, which explicitly calculates the evolution of each discontinuity following its position. This is known to be a difficult problem if real cases are considered [START_REF] Liu | A conservative front tracking method[END_REF]. Shock capturing methods have then been advanced. as if it was associated with a discontinuity. Either this involves an (approximate) solution of a Riemann problem as in Godunov-type methods [START_REF] Guinot | Godunov-type schemes: An introduction for engineers[END_REF], van Leer 1994[START_REF] León | Application of Godunovtype schemes to transient mixed flows[END_REF], or this calculus is avoided by taking advantage of the finite domain of influence of the discontinuity as in centred methods [START_REF] Jiang | Non-oscillatory central schemes for multi-dimensional hyperbolic conservation laws[END_REF]Tadmor 1998, Toro and[START_REF] Toro | Centered TVD-schemes for hyperbolic conservation laws[END_REF].

These are appropriate for the precise computation of extreme flows in dam break problems [START_REF] Paquier | Validity of 1D model for simulating dam-break wave[END_REF][START_REF] Zoppou | Explicit schemes for dam-break simulations[END_REF][START_REF] Burguete | Implicit schemes with large time step for nonlinear equations: Application to river flow hydraulics[END_REF]. One of their drawbacks is that they require small time steps and are highly time-consuming. Toro and Garcia-Navarro ( 2007) reviewed Godunov-type methods for free-surface flows that are now routinely used by the engineering community to compute dam break problems.

Most classical shock capturing methods are explicit. Some implicit methods have also been proposed. But, according to the Total Variation Diminishing (TVD) criterion [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF]) to avoid numerical oscillations, even if they are implicit, they are constrained by a limited time step (except, formally, for a fully implicit method). Furthermore, the implicit systems need to be linearized eventually leading to accuracy problems if large time steps are used [START_REF] Burguete | Implicit schemes with large time step for nonlinear equations: Application to river flow hydraulics[END_REF]). To our knowledge, these methods have not yet been used in practice.

The attractive features of these methods are mostly related to the ability to minimise diffusion in presence of strong gradients. Amongst shortcomings are the computation burden and the difficulty to model inline structures and complex looped network [START_REF] Mignot | Experimental and numerical modeling of symmetrical four-branch supercritical cross-junction flow[END_REF]. None of these methods is really adapted to our requirements.

For an irrigation canal, the simulation can typically last from one week to several months. The main objective is not to model accurately sharp bores, but rather to model correctly the low and medium frequencies dynamics, and the water distribution along the canals using the lowest computational cost. The scheme must also preserve water conservation. The Preissmann scheme meets the present requirements in non-transcritical cases, leading to an adaptation for all flow conditions.

2.4 Meselhe and Holly's restrictions on Preissmann's scheme [START_REF] Meselhe | Invalidity of Preissmann scheme for transcritical flow[END_REF] have demonstrated that Preissmann's method cannot be used directly to solve the complete Saint-Venant's equations for transcritical flows. The three main points of their instability arguments are discussed hereafter. The first problem if applying Preissmann's scheme for transcritical flow directly comes from its formulation that imposes exactly two BC for the system. If the flow presents one transition from super-to subcritical (shock/bore) the physics requires two BC upstream and one BC downstream, which numerically leads to one excess equation. If the transition is from subcritical to supercritical (critical point) only one upstream BC is required, which leads to one missing equation. Thus, if one transcritical transition does exist, then imposing exactly the external physical BCs results in an ill-posed numerical system.

Problems related to usage for transcritical flows

Another problem is related to the fact that Preissmann's scheme applied to Saint-Venant's equations is marginally stable if F = 1. Hence, this method cannot solve the case of an even number of transcritical transitions. Even for an apparently well-posed discrete system, the method is likely to diverge.

Flexibility according to boundary conditions

Saying that Preissmann's scheme is not sufficiently flexible to admit any kind of BC, [START_REF] Meselhe | Invalidity of Preissmann scheme for transcritical flow[END_REF] compared it with Euler's implicit scheme. They argued that the discrete system is closed for the interior nodes and is flexible according to BC as it will accept imposing any number from one to four BC. But this supposes that 'the remaining four unknowns at the boundary nodes can be evaluated by either imposition or extrapolation of BC'. Thus the advantage of flexibility is converted into the problem of defining these extrapolations.

The treatment of boundary conditions is a key point in the design of numerical methods -either external BC or internal BC around hydraulic structures. For a conservative scheme, agreement with physics depends on what enters and leaves the domain, depending on the BCs. Moreover, the error growth factor in the domain also depends on the way BCs are treated. The best ones to be imposed do correspond to the physical conditions.

Marginal stability

Concerning the problem of marginal stability, it must be noted that almost all schemes dealing with transcritical flow add an amount of diffusion in transcritical zones (adding an explicit diffusion term or through limitations of the variations). This necessity can be linked to the invalidity of the initial Saint-Venant approximation in these domains, as they are deduced under the assumptions of hydrostatic pressure and uniform distribution of the velocity along the vertical axis. Accounting for turbulence or non-hydrostatic pressure distribution in models is computationally expensive. Adding numerical diffusion permits to damp the infinite gradients that appear within the classical hypothesis framework. Neither procedure is exact for real cases or susceptible to correct numerical treatment.

Residual distribution concept

The method of residual distribution was proposed by [START_REF] Morton | Upwind iteration methods for cell vertex scheme in one dimension[END_REF] to introduce some up-winding using cell-vertex methods. The problem is approached considering the lack of agreement between the number of discrete unknowns and equations. The idea is to combine, for each section, the cell residuals of the right and left neighbouring cells (Fig.

2 a), thus obtaining one "nodal residual" associated with each section

j j j / j- j j B R R R 2 / 1 2 / 1 2 1 2 / 1 (13)
where ± . = transfer matrices and the vector B accounts for external and internal boundary conditions. For + -1/2 =0 and -N+1/2 =0, the complete system is given by

N j R j ... 0 : 0 (14)
To preserve the conservative property, transfer matrices must satisfy 2 1 2 1 / j / j (15) where identity matrix.

Splitting with characteristic decomposition

To obtain the distribution according to wave propagating in each cell, a local Jacobian matrix J j+1/2 is used as

0 2 0 2 0 0 2 ~1 0 v v c J ( 16
)
where the particular mean-values v ~ and c ~ are given by [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] approximation

Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

1 1 1 2 / 1 2 1 2 2 / 1 2 ) ( ~ j j j j j j j j j j c c c v c v v c c c (17)
This choice ensures that bores propagate with the correct velocity. Hence, the transfer matrices used by Johnson to define the nodal residuals (13) correspond to the application of up-winding on the characteristic variables defined in the Eigenspace of the system of Eqs.

(3) to (5) as

1 2 / 1 2 / 1 2 / 1 2 / 1 j j j j V V (18)
where V j+1/2 = matrix of Eigenvectors for J j+1/2 and the matrices ±

j j / j- j j B R R R 2 / 1 2 1 2 / 1 0 j supercritical, j+1 subcritical: j j j / j- j B R R R 2 / 1 2 / 1 2 1 1
Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440.

The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648 The use of nodal residuals solves the over-determination at shocks. Nonetheless, this does not ensure that the new problem is well-posed since the equations still have to be independent. In particular, if cell j 1/2 is subcritical and j+1/2 supercritical, the transfer

matrices defining R j satisfy 2 ) ( rank 2 / 1 2 / 1 j j (21)
which means that the resulting system is ill-posed. To provide an extra linearly independent equation, Johnson then adds an internal boundary condition

j j j j u B ) ( 2 / 1 2 / 1 (22)
characterizing a critical point located exactly at x j .

The non-linear implicit system is solved by an iterative method. The approximate Jacobian matrix must be given at the new time step, which means calculating it at each iteration. This is necessary to allow the critical interface to cross more than one computational cell within one time step.

Limitation

The general drawback of Johnson's method is that it requires, for each cell, calculating mean-value propagation velocities using Eq. ( 17) and testing the dimension of transfer matrices as given in Eq. ( 21). This is necessary to obtain an automatic procedure, adapted to any regime. But then, the method leading to Eq. ( 13) uses three points, even in the subcritical regions where the initial Preissmann's scheme would be appropriate. The following focuses the necessity of this three points approach.

Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440.

The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648 4 Preissmann's scheme with internal boundary conditions

Direction of solution and internal boundary conditions

To achieve stability, a solution along directions of propagation must be respected. With

Preissmann's scheme, this means that the direction of solution has to change if the regime changes. Hence an equivalent/internal boundary condition has to be prescribed at transcritical transitions to permit a change in the solution direction. These relations should describe the particular physical properties associated with transcritical transitions.

Actually, the ill-posedness of the classical implementation of Preissmann's method for both possible transitions can be related to the physical properties of the transitions. In both cases the problem is completely similar to solving the Saint-Venant system with the method of characteristics:

-For a bore, while convergence of characteristics leads to an over-determination;

Rankine-Hugoniot jump conditions prove that influences from each side must be taken into account to determine the flow evolution (to obtain right and left states and shock speed). The "over-determination" obtained with Preissmann's classical system corresponds to this fact saying this complete piece of information actually does have to be used to describe the flow at this point.

-At a critical point, the steady characteristic has a proper definition. The "missing equation" observed with Preissmann's classical system corresponds to this fact asking for supplementary information in the description of this specific point.

The idea is thus to translate the physical properties of transcritical transitions to link the different regime sub-domains, within which Preissmann's scheme is adapted. Johnson's method shows that it is possible to treat the transcritical transitions with a through method.

The mathematical description of Johnson's combination shows the two systems are exactly equivalent for entirely subcritical or supercritical domains.

Treatment of transcritical flows

Preissmann's scheme is used in each fully subcritical or fully supercritical domain. This regime, at time steps n+1 (with Newton solution, the last approximation is used, hence for the first solution, the regime at time step n is used).

Formally, the system of Preissmann's equations over a non-transcritical sub-domain can be reduced to a set of equations at the boundaries (in the idea of simple-and double-sweep methods). The attention is thus restricted to cells where the regime changes. The numerical problem is approached as reduced to giving boundary conditions for each subcritical or supercritical region, consistent with Preissmann's scheme (Fig. 3). The existence of the critical point does not compromise the conservation property over the cell. Preissmann's system is then retained, but the equation defining the critical point ( 8) is added to characterize the particular evolution due to the critical point.

If F j < 1 F j+1 a critical point exists in this cell, at point x c =x j + dx (0< 1 to be defined), such that v c =c c . The system can be simply written then

0 0 2 1 = B = R AQ / j (23)
where B AQ = 0 is the discrete equation for the steady characteristic written in terms of the conserved variables A and Q, that must be defined at point x c .

As the solution is continuous across the critical point, the position can be interpolated linearly between the two neighbouring mesh points. This is consistent with Preissmann's discretisation over the cell (with the space weighting =1/2). The point is characterised by

a (1) =0, resulting in ) 1 ( ) 1 ( 1 ) 1 ( j j j a a a ( 24
)
Using Eq. ( 9) gives

j j j j j j AQ c Q A c Q A B 2 1 2 1 1 1 (25)
where the values characterising the wave are given at the approximate new state.

If this modification is considered in splitting methods, the supplementary equation corresponds to the correction required as the simple Roe-type splitting is applied. This does not always satisfy the entropy condition without introducing appropriate modifications. This is generally achieved either by using another approximate solver (Delis 2003, Ying and[START_REF] Ying | Improved implementation of the HLL approximate Riemann solver for one-dimensional open channel flows[END_REF] or by adding artificial viscosity [START_REF] Burguete | Implicit schemes with large time step for nonlinear equations: Application to river flow hydraulics[END_REF].

As compared to classical corrections including Johnson's method, where a supplementary term is added to an equation of the initial system, the above formulation expresses the "physical need" for diffusion terms in classical schemes based on Saint-Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440.

The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

Venant's equations. It is similar to the supplementary "compatibility condition" proposed by [START_REF] Djordjevic | Simulation of transcritical flow in pipe/channel networks[END_REF].

Shock

If a supercritical-subcritical transition is encountered, the direction of shock propagation must be tested to define which sections are influenced by the shock. Roe's definition is appropriate there, since the Eigenvalue of his approximate Jacobian is the shock speed if the values on both sides are linked by a jump.

F j 1 > F j+1 , determine ) 1 ( 2 / 1 ~ j a -if 0 ~) 1 ( 2 / 1 j a
, the shock is travelling downstream (Fig. 5a), the solution at point j can be given by Preissmann's scheme in cell j 1/2, while the solution at point j+1 accounts for influences from both upstream and downstream,

-if 0 ~) 1 ( 2 / 1 j a
, the shock is travelling upstream (Fig. 5b), the solution at point j must account for influences from both upstream and downstream. To account for upstream and downstream influences at point k (=j, or j+1) if the solution at point k 1 is known, the solution there must remain indeterminate as long as the subcritical region is not solved. At the same time, information must be transmitted downstream to the subcritical region. It is also essential that conservation be respected. To

Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648 combine these requirements, and following the idea of Johnson's method, Preissmann's solutions for the two cells k 1/2 and k+1/2 are added and the equation for the downstream travelling characteristic issued from cell k+1/2 gives an additional equation in the same unknowns (linked through the definition of dw ( 2) ) as

0 0 2 / 1 ) 2 ( ) 2 ( ) 2 ( ) 2 ( 2 / 1 2 1 k x t k / k s w a w R R (26)
The conservation property inherent to Preissmann's method is globally preserved and not compromised by the addition of the second equation giving additional information to be transmitted to the subcritical domain. Furthermore, this summation cancels the flux at the intermediary boundary, which is that of the near-critical node, potentially responsible for instability.

According to the definition of dw (2) , the supplementary equation can be written with the physical increments dA and dQ. Using the matrix of Eigenvectors, a combination of the equations of conservation of mass and momentum is found as

0 2 / 1 ) 1 ( k Q A R R a (27)
which gives the third equation for the two cells implicated in the shock.

Resulting scheme

Preissmann's scheme is directly used in each fully subcritical or fully supercritical domain and specific equations are added at transcritical transitions. The external boundary conditions are then physical.

The procedure begins with computing the Froude number for each node. The system of equations is then developed by determining the equations one cell after the other, according to the regimes of the nodes, tested at n+1 as: 1) Impose the boundary conditions corresponding to the inflow regime.

2) For j = 0…N 1 -If (F j 1)(F j+1 1)>0, use Preissmann scheme for cell j+1/2

-If F j 1 > F j+1 , compute ) 1 ( 2 / 1 ~ j a if 0 ~) 1 ( 2 / 1 j a
, then use Eq. ( 28) with k=j+1

-verify F j+2 < 1, otherwise add Eq. ( 29) for l=j+1 (in cell j+3/2)

-j:=j+1 if 0 ~) 1 ( 2 / 1 j a
, previous cell (j 1/2) is not solved with Preissmann's scheme boundary conditions for the supercritical zone which can be solved by a classical simple sweep-method.

Numerical results

The method is applied on a transient problem involving appearance and disappearance of a supercritical region between a critical point and a bore. Starting from a steady subcritical state, in a canal of mild slope (0.01%) followed by a steeper (1%) slope, the outflow boundary condition is varied as dh/dt=5 mm/s. First, the depth is reduced toward a steady transcritical state, and then raised back to the original depth. The inflow discharge is constant at 50 m 3 /s. The channel is rectangular of width 5 m and has a Strickler roughness coefficient of K=50 m 1/3 /s. The solution involves space and time steps of dt=10 s and dx=10 m, and Preissmann scheme is used with =0.55.

For this problem, along with a classical implementation using Preissmann's scheme, divergence appears as soon as the transcritical zone appears [START_REF] Djordjevic | Simulation of transcritical flow in pipe/channel networks[END_REF]). The result of the proposed method is generally the same as with Johnson's method, except for the upstream moving bore, where a divergence can appear with Johnson's method.

Results are presented for different specific states of the solution. The transcritical zone appears at t=480 s in Fig. 6a. The entire subcritical domain is then separated at the location of the change of bed slope by a critical point immediately followed by a shock. As the downstream water level is continuously reduced, a shock moves downstream (Fig. 6b).

Then, the boundary conditions are kept constant and a steady transcritical state is established (Fig. 6c). As the downstream flow depth is raised back to its initial value, the shock moves upstream. Figure 6d shows the stronger variation appearing at the shock location. As the shock attains the critical point, the supercritical zone disappears and the initial state is found back.

Some details are given in Fig. 7 to show both transcritical transitions at the intermediate steady state and the evolution of the solution when the transcritical domain appears and disappears.

At the intermediate steady state the flow characteristics may be described with an inflow velocity v=2.5 ms -1 , c 6.26 ms -1 ; in the steep slope region, the normal depth is h 1.90 m, corresponding to v 5.25 ms -1 , c 4.32 ms -1 . The Courant number is thus around 9.

The solution is also shown for a high Courant number in Fig. 6c. The same problem is treated but taking dt=5dx=50. A small local disturbance in discharge is observed at the location of the steady bore (Fig. 8). This remains a problem essentially in the steady case. In the transient case, the disturbance is embedded in the strong modifications of the flow around this point. Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

The method obtained for the solution of transcritical flows directly uses Preissmann's scheme. The classical use of this two-point scheme is preserved in the major flow part.

Even though there is an explicit determination of the cells experiencing the transcritical changes, it is a shock capturing method. It can be directly related to the three-point method of Johnson whose formulation is more systematic. The regime test is not explicitly done in the definition of the algorithm but implicitly calculated in the scheme. Johnson's method is elegant but finally requires more calculation.

The problem of transcritical transitions is solved locally, either by associating the cells involved in a bore and adding an equation to characterize the information transferred in the subcritical domain, or by the addition of an internal boundary condition to characterize the expansion fan at the critical point.

The internal boundary condition is on the one hand an alternative to the introduction of entropy corrections to solve the classical problem of rarefaction-shocks linked to the use of so-called Roe-type splitting. On the other hand, the solution of shocks is primarily based on Preissmann's scheme written for non-transcritical sections, and considering the solution in this domain, the supplementary discretisation of the characteristic equation adds some diffusion to the Preissmann's system.

The required test on the discrete Froude number is not expensive as compared to other particular treatments that can be proposed for transcritical flows. This is a simple way to conserve the two-point discretisation with its robustness and adapted form for the inclusion of hydraulic structures. Besides, the double-sweep algorithm can be used. Moreover, this "regime-detector" is also a useful tool to organize external boundary conditions. It is intended to make the scheme able to support any flow evolution or any change imposed by hydraulic structures. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

  Eigenvectors, v=Q/A = mean flow velocity and c = celerity of pressure waves with c 2 =g dI 1 /dA =gA/B.

Figure 1

 1 Figure 1 Schematic representation of simple-and double-sweep solution, propagation of information in (a) supercritical state, (b) subcritical state, (c) on how necessary information is given by equations

  Bases of Johnson's scheme[START_REF] Johnson | A box scheme for transcritical flow[END_REF] proposed a method based on Preissmann's scheme, adapted to the treatment of transcritical flows and using physical boundary conditions. It combines, for each section, the two formulations obtained with Preissmann's scheme for the right and left neighbouring cells. The combination is obtained through a local decomposition (residual distribution concept) along the wave propagation directions (splitting with characteristic decomposition), in the idea of "flux vector splitting" methods[START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF].

  cells do not have the same regime (Fig. 2 b), the nodal residual is j subcritical, j+1 supercritical:

Figure 2

 2 Figure 2 Nodal residual of (a) general from, Johnson's nodal residual for (b1) supercriticalsubcritical transition, (b2) subcritical-supercritical transition

  means determining the regime at each point (what is generally implicitly done in classical methods for transcritical flows). The flow regime is tested via the local Froude number F = |v|/c at each point. The characterisation of "fully" subcritical or supercritical means, for the solution at n+1 in cell {j,j+1}, that values at both points j and j+1 must have the same Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440.The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

Figure 3 Figure 4

 34 Figure 3 Reduction of solution to neighbour of transcritical transitions for (a) supercriticalsubcritical transition, (b) subcritical-supercritical transitionCritical pointThe rarefaction fan issued from the critical point shows that, from this point, information is transferred right and left. Hence, the equation defining the critical point is able to furnish the internal boundary condition giving, at the same time, the BC "closing" the upstream subcritical flow and the second BC necessary to define the evolution of the supercritical flow downstream.

Figure 5

 5 Figure 5 Methods for shock travelling (a) downstream, (b) upstream. The sign of Roe's first Eigenvalue gives direction of shock propagation

Figure 6 Figure 7

 67 Figure 6 Solution evolution with (a) Appearance of the transcritical zone, (b) Downstream moving shock, (c) Intermediate steady state with dt=dx and dt=5dx, (d) Upstream moving shock

Figure 8

 8 Figure 8 Profile of Q(x) at intermediate steady state

Notationsa

  (i) = Eigenvalues of Jacobian matrix (i = 1 or 2) A = Wetted cross-sectional area Author-produced version of the article published in Journal of Hydraulic Research, 2010, 48(4), 428-440. The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648 B = Channel width B j = Boundary condition at node j B AQ = Scalar internal boundary condition c = (gA/B) 1/2 = Celerity of infinitesimal wave c ~= Roe's approximation for infinitesimal wave celerity dw = V -1 du = variation of Eigenvariable D = diag{a (1) ;a (2) } = Diagonal matrix formed with Eigenvalues a (1) and a (2) f = Flux of momentum equation f = General vector flux function F = Froude number g = Gravitational acceleration I 1 = Hydrostatic-pressure term J = Jacobian matrix of Saint-Venant systemK = Strickler coefficient Q = Discharge R j+1/2 = (R A j+1/2 R Q j+1/2 ) T =Preissmann's cell residual R j = Nodal residual s = Source term of momentum equation s = General vector source term s = V -1 s = Vector source term in Eigen-system S o = Bottom slope S f = Friction slope t = Time u = (A Q) T = General vector conserved variable v = Q/A = Mean celerity of flow v ~= Roe's approximation for mean celerity of flow V = Matrix of eigenvectors of J w = Vector-conserved variables in Eigensystem x = Space coordinateAuthor-produced version of the article published in Journal of HydraulicResearch, 2010, 48(4), 428-440.

  Author-produced version of the article published in Journal of HydraulicResearch, 2010, 48(4), 428-440.The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

  Author-produced version of the article published in Journal of HydraulicResearch, 2010, 48(4), 428-440.The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

  Author-produced version of the article published in Journal of HydraulicResearch, 2010, 48(4), 428-440.The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

  Author-produced version of the article published in Journal of HydraulicResearch, 2010, 48(4), 428-440.The original publication is available at http://www.tandfonline.com/ doi : 10.1080/00221686.2010.491648

	16			16
	12			12
	y (m)			y (m)
	8			8
	Bottom				Bottom
	Water level				Water level
	Initial state				Critical depth
	4			4
	0	500	x (m)	1000	0	500	x (m)	1000

Acknowledgements

The authors wish to thank the JHR editorial Board for their valuable comments improving the quality of this paper.

-use Eq. ( 28) with k=j

-If F j < 1 F j+1 compute from Eq. ( 24), add Eq. ( 29) for l=j

use Preissmann system for cell j+1/2

3) Impose the boundary conditions corresponding to the outflow regime.

Test on external boundary conditions

In the described procedure, it is supposed that the numerical regime corresponds to the one imposed by the user-prescribed boundary conditions. Tests must be added at the boundaries to verify this consistency. The user has to specify if the proposed boundary conditions are to prevail or not.

Solution of the system

This non-linear system is solved by Newton's iterative method; transcritical positions and their characteristics are computed for each iteration. The solution of the linearized system, for (half-) implicit schemes, requires inverting a matrix. One of the advantages of Preissmann's scheme is that simple algorithmic solutions can be used. Using the present transcritical system, these methods can be adopted by just localising the transcritical transitions.

Bore: Combination of sections j 1, j, j+1

The sweep in the supercritical regime gives the solution from the inflow to cell j 1.

The system (28) can then be simplified giving the equivalent of an inflow BC at j and a system of two equations for the cell j+1/2. The subcritical region is then classically solved with a double sweep method.

Critical point: Addition of one equation for sections j, j+1

The first sweep in the subcritical domain, from the inflow boundary condition to cell j 1/2, leads to one equation in (Q j , A j ). This statement simplifies the equations of cell j+1/2, giving three equations with three unknowns. It correspond to obtaining the solution at j+1, being able to perform the second subcritical sweep, and giving the required upstream


CEPN Damien Besancenot

Piecework versus merit pay: a mean field game approach to academic behavior

Keywords: Mean Field Game, Academic production, incentives, Lotka's law. JEL Classification : C72, I23, O31

 to the analysis of the researchers' academic productivity. It provides a theoretical motivation for the stability of the universaly observed Lotka's law. It shows that a remuneration scheme taking into account the researchers rank with respect to the academic resume can induce a larger number of researchers to overtake a minimal production standard. It thus appears as superior to piecework remuneration.

Introduction

Over the last twenty years, academic publications tend to become the generalized gauge of the academic performance. Would assessment be directly grounded on the number of publications, on the number of citations or on the impact factor of the journals, all bibliometric criteria give a dramatic role to the characteristics of scientific productions in the evaluation of individual researchers [START_REF] Groot | Research quality and efficiency: an analysis of assessments and management issues in Dutch economics and business research programs[END_REF]Garcia-Valderrama, 2006, Bence andOppenheim 2004). In return these publications play an increasing role in setting individual salaries and promotions [START_REF] Hamermesh | Scholarship, citation and salaries : economic rewards in economics[END_REF]. 1A vast strand of research deals with the influence of the academic resume on the researchers' careers and rewards. For instance, [START_REF] Stephan | The Economic of Science[END_REF] emphasizes that any publication contributes to the researcher's social recognition which is a main part of the academic reward. [START_REF] Diamond | What Is a Citation Worth?[END_REF] shows that the marginal worth of a citation may vary between 50 and 1300 dollars. In a same way, Swidler and Goldreyer (1998) state that the present value of the first top finance journal article is between 19 493 and 33 754 dollars. In France, since 2008, the jury of the "Concours d'Agregation pour le recrutement des professeurs d'universite" in charge of the recruitment of professors of economics is explicitly using a bibliometric evaluation rule based on the number of publications and the quality of their support (as revealed by the CNRS ranking) to assess applicants (Levy [START_REF] Garboua | Rapport sur le premier concours national d'agrégation de l'enseignement supérieur pour le recrutement des professeurs des universités[END_REF], Ellison 2010).

Despite these major changes in the governance of the research system, the distribution of the scientific production among the academic population appears as an unchanging rule of the academic production pattern. Almost one century ago, [START_REF] Lotka | The frequency distribution of scienti.c productivity[END_REF], already highlighted the inverse relationship between the number of scientific publications in an academic resume and the number of authors that achieve such a CV. Dealing with the productivity of researchers in chemistry and physics, the initial work of Lotka revealed a statistical regularity whereby the majority of authors (60%) publish a single contribution during their career while a small number contributes an essential part of the total number of scientific contributions 2 .

Lotka's law is regarded as one of the three pillars of contemporary scientometrics 3 . Its validity has been the subject of numerous empirical verifications. An illustration of this law is given by Le [START_REF] Coadic | Mathématique et statistique en science de l'information et en science de la communication: Infométrie mathématique et infométrie statistique des revues scientifiques[END_REF] to characterize the production of French researchers in chemistry. It is also verified in Economics [START_REF] Cox | Patterns of Research Out-put and Author Concentration in the Economics Literature[END_REF], finance [START_REF] Chung | Patterns of Productivity in the Finance Literature: A Study of the Bibliometric Distributions[END_REF], medicine [START_REF] Kawamura | Lotka.s law and the pattern of scienti.c productivity in the dental science literature[END_REF]... The robustness of this law is also attested by similar studies using other indexes: Combes and Linnemer (2001) focus on the number of articles weighted by the quality of journals and [START_REF] Courtault | Research in Economics and Management in France: A bibliometric study using the h-index[END_REF][START_REF] Courtault | De la réputation scientifique et de sa mesure[END_REF] considered the distribution of the h index [START_REF] Hirsch | An index to quantify an individual's scientific output[END_REF]. Both found that these indexes are distributed according to Lotka's law. 4 In these settings, the aim of the present paper is to show that Lotka's distribution may be a persistent equilibrium distribution of the scientific production in a game where the publication choices of researchers are made under alternative regimes of incentives. In this purpose, we develop a highly stylized mean field game in the line of Guéant, Lasry and Lions (2010) to 2 Lotka's law states that the number of authors making n publications is about 1/n a of those making one publication, where a is approximately equal to two. For example, the number of researchers who have ten contributions represents one percent of the number of researchers who have one publication only.

3 the two other are Bradford's law (on a specific topic, a limited number of journals publish a large number of contributions and a large number of journals publish a limited number of items) and Zipf's law (if in a given text, words are ordered according to their frequency, the product rank by frequency is approximately constant). It is easy to show that all these laws are particular cases of Pareto's law, i.e.: roughly 80% of the effects comes from 20% of the causes (Cf. [START_REF] Egghe | Introduction to Informetrics : quantitative methods in library, documentation and information science[END_REF]Rousseau 1990, Newman 2005). 4 This robustness of similar distribution has been shown in a larg number of fields in economics and finance (income and wealth, the size of cities and firms, stock market returns, trading volume, international trade etc...), Cf. [START_REF] Gabaix | Power laws in economics and finance[END_REF] consider the optimal decision of a rational academic researcher interacting with a continuum of other researchers.

Introduced by [START_REF] Lasry | Mean field games[END_REF], the Mean Field Games Theory allows the study of strategic behavior when many agents are in a situation of interaction. When the number of players is important (mainly when there is a continuum of agents), one can hardly assume that a given player is able to take into account all the other players interactions and to compute the Nash equilibrium of a game. In this case, the Mean Field Game Theory adopts the methodological approach of statistical physics in the modelization of the interaction of a great number of particles. Faced to this insurmountable computational problem, physicists consider each particle as being influenced by a mean field exerted by all other particles while simultaneously taking into account the influence of each particle on the mean field. In order to formalize the behavior of a great number of rational agents, the Mean Field Game Theory assumes each agent to be influenced by the mean field made of the distribution of other players' behavior and considers the consequences of each individual decision on this mean field. In a mean field game, an agent cannot directly influence other agents but each agent contributes to the formation of a mean field which influences in turn the behavior of each individual.

Developed in continuous time, this approach involves methods of stochastic calculus and dynamic optimization. In a Mean Field Game, the dynamics of the system is governed by two equations: a backward Hamilton-Jacobi-Bellman equation describing the optimal behavior of agents given the distribution of the other players and a forward Kolmogorov equation which takes into account the influence of each player on the mean field. The fixed point of those two equations gives the Nash equilibrium of the game : Lasry and Lions proved that the coupled system, i.e. the backward HJB equation and the forward Kolmogorov equation can be obtained as the asymptotic limit of Nash equilibrium of N players as N goes to infinity.

It is important to note that the approach of mean field games is different from what is done in standard games theory : the specific partial differential equations system (forward/backward structure) studied in MFG theory is new and does not appear in the literature before MFG, see for instance [START_REF] Khan | Non-cooperative games with many players[END_REF] or [START_REF] Aumann | Markets with a Continuum of Traders[END_REF]. Morover, one have to mention that the forward/backward structure is a real difference between MFG theory and econophysics where the anticipations process of agents is not sufficiently taken into account.

In the paper at hand, we consider infinite living researchers who decide at each point of time the level of their academic production taking into account the tradeoff between the cost of this production and the reward that it provides. The equilibrium of such a game is computed for alternatives compensation mechanisms. In particular, the game allows to show that the Paretian distribution is a stable solution under merit pay. Finally, some comparative statics allows to show that the fraction of academic researchers having a publication score above a given threshold is an increasing function of monetary incentive and a decreasing function of time preference.

In the line of Pareto (1896) and [START_REF] Allais | Classes Sociales et Civilisations[END_REF] who believed that economic policy could not influence the distribution of income and wealth (supposed to be independent even of the economic system), we found that economic policy cannot change the distribution of the scientific production. However, policy may influence both the speed at which this distribution moves and its dispersion and proves to have an impact on the overall research output.

Section 2 presents the mean field game approach of academic behavior both with and without externalities. Section 2.2, takes explicitly into account the influence of ranking on the researchers' remuneration. The analysis of the effects of this externality in the reward system is the main innovation of the paper. Section 3 interprets the economic meaningfulness of the results. A last section concludes the paper.

A mean field game approach to academic behavior)

Let a t denote the instantaneous scientific production of a researcher at date t and x t be the cumulative measure of this production at the same date. 5The relation between a t and x t is given by :

dx t = a t dt
Let there be a continuum of researchers. Define m(t, x) as the density of the distribution of x t at date t. By assumption, the support of the distribution at date 0 is [1, ∞), where we consider that the worst researcher at date 0 is a newly mint PhD student with only his PhD Thesis. At date t, the optimization problem of the researcher is given by :

U (t, x) = max (as) +∞ t f (x s , a s , m(s, x s ))e -r(s-t) ds/x t = x (1)
where f (x s , a s , m(s, x s )) is the researcher's instantaneous utility function and r takes into account the researcher's time preference. Below, we will consider several specific expressions for the function f (., ., .).

If we suppose, for instance, that m is known the solution of this optimization problem is given by the Hamilton-Jacobi-Belmann (HJB) equation 6 :

∂U ∂t (t, x) + max a f (x, a, m(t, x)) + a ∂U ∂x -rU = 0 (2)
The control is given by a = argmax(f (x, a) + a ∂U ∂x ). Heuristically one can obtain this equation using the dynamic programming principle of Bellman :

U (t, x) = max (as) f (x, a t , m(t, x))dt + e -rdt U (t + dt, x + dx t )/x t = x (3)
At time t=0, each researcher is endowed with an initial academic production. Let m 0 be the initial density of x 0 , the scientific production of the researchers at time 0. m(t, •) is the solution of Kolmogorov equation:

∂m ∂t (t, x) + ∂ ∂x (a(t, x)m(t, x)) = 0 (4)
For the reader convenience, we give the formal derivation of this equation which is also useful for the stochastic case. Let Φ be a smooth function with compact support,

m(t + dt, x)Φ(x)dx = m(t, x)Φ(x)dx + ∂m ∂t (t, x)Φ(x)dtdx
But m(t, x) is the density of the distribution of x t at date t, so

m(t + dt, x)Φ(x)dx = E(Φ(x t+dt )) = E(Φ(x t ) + Φ (x t )dx t ) = E(Φ(x t ) + Φ (x t )a t dt) = m(t, x)(Φ(x) + Φ (x)a t dt)dx = (m(t, x)Φ(x) - ∂ ∂x (a(t, x)m(t, x))Φ(x)dt)dx Hence ∂m ∂t (t, x)Φ(x)dtdx = - ∂ ∂x (a(t, x)m(t, x))Φ(x)dt)dx
and this identity implies the Kolmogorov equation (Φ is a test function).

In order to get explicit solutions we have now to consider explicit forms for f (., .).

A model without externalities

In their research activities, researchers have to consider the trade-off between the cost of carrying out their research and the specific reward of total and instantaneous scientific production (Swidler and Goldrayer 1998). This trade off between costs and remuneration may be formalized as :

f (x, a) = αx + βa -a 2 , ( 5
)
where β is the instantaneous reward that a researcher gets from an additional publication. This reward incorporates both a monetary premium (as in some scientific institutions such as business schools which pay researchers for each new publication) and a psychological reward given by the mere fact of being published by a selective journal [START_REF] Stephan | The Economic of Science[END_REF].

α represents the reward that a researcher gets given his academic resume. Once again this incorporates both monetary and social rewards. The cumulative measure of scientific production is essential in determining promotions and raises for scientific researchers [START_REF] Diamond | What Is a Citation Worth?[END_REF][START_REF] Stephan | The Economic of Science[END_REF]. Beside, the more an academic researcher publishes, the greater his readership and the higher will be the influence of his work on the academic field.

The last term in the previous equation measures the research cost associated with an instantaneous production a. This cost is assumed to be increasing at an increasing rate. For the sake of simplicity, we consider a quadratic cost function.

In this section, we ignore the researchers interactions. The externalities generated by the ranking of the academic researchers within in a field of research will be dealt with Mean Field Game in subsection 2.

Althought we do not need to use the MFG equations for this model without externalities, for the reader convenience, we revisit the solution of the model via the MFG "tools". For a traditional exposition using calculus of variation, see the appendix.

The control is given by a

= argmax(f (x, a) + a ∂U ∂x ) = 1 2 ( ∂U ∂x + β). Repl- cing in (2), we find ∂U ∂t (t, x) + f (x, 1 2 ( ∂U ∂ x + β)) + 1 2 ( ∂U ∂ x + β) ∂U ∂x -rU = 0 ∂U ∂t (t, x) + αx + β 2 4 -rU + β 2 ∂U ∂x + ( ∂U ∂x ) 2 4 = 0 (6)
It is usual, in regards with the terms of this equation, to look for a solution with the following form U (t, x) = b(t)x 2 + c(t)x + g(t).

Replacing in [START_REF] Brock | Optimal Control and Economic Dynamics[END_REF] and by identifying the coefficients of the polynomial gives :

b (t) = -b 2 (t) + rb(t) (7)
c (t) = -b(t)c(t) -βb(t) + rc(t) -α (8)
g (t) = -c 2 (t)/4 -β/2c(t) + rg(t) -β 2 /4 (9)
For the solution to be admissible (see the integrand in the maximization criterion), one must have :

x(t) = x 0 + α + βr 2r t (10)
and a(t) = α+βr 2r . Moreover, we know that m(t, x) is solution of :

∂ t m(t, x) + a(t)∂ x m(t, x) = 0 (11) m(0, x) = m 0 (x) (12)
Even this is the simplest transport equation we give the following details for the reader convenience. First, one can see that m(t, x(t

)) = cte = m(0, x(0)), indeed dm(t, x(t)) dt = ∂m ∂t + x (t) ∂m ∂x = ∂m ∂t + a(t) ∂m ∂x = 0
then (expressing x(0) as a function of x(=x(t))) we find

m(t, x) = m 0 x - α/r + β 2 t (13)
In this first model, the instantaneous production of each researchers is constant (da/dt = 0). Each researcher exhibits the same instantaneous scientific production (∂a/∂x = 0) and any initial distribution (and in particular the Pareto distribution) will be preserved over time while "traveling" (like a travelling-wave) with the following velocity : α/r+β 2 . At time t, m(t, x) is thus the translation of m 0 by α/r+β 2 t. To illustrate graphically this phenomena let us fix the parameters α = r = 0, 05 and β = 0 and take for m(0, x) a bell-type curve distribution, we obtain the following graphics (Fig. 1) for the densities at t = 0 and t = 20, i.e. m(0, x) and m(20, x) : Fig. 1 The speed of the wave describes the intensity of the publication activity of a generation since it compares the stock of publications of the same pool of researchers at two different dates. It is an increasing function of α and β and a decreasing function of the researchers' time preference. Indeed, the higher this time preference is, the lower is the discounted value of net benefit of a publication and so the lower is the quantity of publications.

Recall that β is the measure of the direct utility of an extra publication, in this model α r is the indirect utility of this extra publication through its impact on the academic resume. As an additional publication increases the stock of publications at each future date, its discounted value is 1 r .

A model with ranking

This section explicitely introduces the role of ranking in the academic reward system. The reward of the academic work is rarely strictly tight to the absolute volume of publication. Instead, most of the remuneration system rest on the comparative number of works. The objective of any researcher is to be in the highest possible place in the overall ranking of his peers. Following [START_REF] Merton | The matthew effect in science[END_REF][START_REF] Merton | Behavior pattern of scientists[END_REF] the goal of a scientist is to be the first to obtain a particular result, the proof of the priority being given by publication. The reward to priority is the recognition awarded by the scientific community to being first. This phenomenon of giving rewards to priority is particularly illustrated with Nobel prizes which are given to persons having obtained an outstanding result in the first place, however, this is not limited to Nobel Prizes, but is the rule of the promotion game for every researchers. Moreover, the competition between researchers is not limited only to the publication race. Obviously, the process of promotion is the result of a tournament in which each academic resume is assessed with respect to others.

Beside, peer recognition within a field is related to one's ranking in the academia. For instance, REPEC provides every month a ranking of the registered authors within one country or a wider geographic country. Newspapers also give periodically ranking of academic institutions or even of scientific researchers. Scientific journal publish similar results ( To implement the influence of the relative situation of a researcher in a scientific community, we assume now that the utility function depends on the ranking of a researcher. This brings an externality into the model which can be fruitfully studied with Mean Field Games.

One way to choose this dependence is to take the remuneration inversely proportional to m. 7 Indeed, the number of people having achieved a particular ranking is smaller the higher the ranking. Instead of using the cumulative distribution function as an index of ranking, it is therefore equivalent and it will prove easier to use the inverse of the distribution function.

Hence we propose the following model :

max (at) +∞ 0 α 1 m(t, x(t)) + β ẋ(t) -ẋ2 (t) e -rt dt/x(0) = x 0 (14)
where dx t = a t dt. We introduce the value function :

U (t, x) = max (as) +∞ t α 1 m(s, x(s)) + β ẋ(s) -ẋ2 (s) e -rs ds/x(t) = x (15)
The Partial Differential Equations system is given by :

∂U ∂t (t, x) + max a α m(t, x) + βa(t) -a 2 (t) + a ∂U ∂x -rU = 0 (16) ∂m ∂t (t, x) + ∂ ∂x (am) = 0 (17) m(0, x) = m 0 (x) (18)
The control is given by a = argmax α m(t,x(t)) + βa(t) -a 2 (t) + a∂ x U = ∂xU +β 2

. Replace this expression in the equations ( 16) -( 17) -( 18) we get :

∂ t U (t, x) + α m(t, x) + 1 4 (β + ∂ x U ) 2 -rU = 0 ( 19
)
∂ t m(t, x) + ∂ x ∂ x U + β 2 m = 0 (20) m(0, x) = m 0 (x) (21)
We recall that we are interested by the "persistence" of a Pareto density in this system, so as explained in the introduction we take m(0, x) = m 0 (x) =

1 x 2 1 x>1 .
Theorem 2.1. For an initial distribution of researcher given by m(0, x) = m 0 (x) = 1 x 2 1 x>1 , the system (19) -(20) admits as solution :

m(t, x) = α r t + 1 (x -β 2 t) 2 1 x>β/2t+ α r t+1 (22)
U (t, x) = α r (x -β 2 t) 2 α r t + 1 + β 2 4r ( 23
)
the control is given by :

a(t, x) = α r (x -β 2 t) α r t + 1 + β 2 (24)
and for the researcher who starts with x 0 , we find :

x(t) = ( α r t + 1)x 0 + β 2 t. (25)
Proof. Suppose that we are looking for an x(t) in the form x(t) = h(t)x 0 + v(t), then a(t, x) = h (t)x 0 + v (t) and the Kolmogorov equation :

∂ t m(t, x) + h (t) h(t) m + a(t, x)∂ x m(t, x) = 0 (26) h(t)∂ t m(t, x) + h (t)m + a(t, x)h(t)∂ x m(t, x) = 0 then for m(0, x) = 1 x 2 1 x>1 , we find m(t, x) = 1 h(t) m(0, (x -v(t))/h(t) = h(t) (x-v(t)) 2 1 x>v(t)+h(t)
. Replacing this expression in [START_REF] Groot | Research quality and efficiency: an analysis of assessments and management issues in Dutch economics and business research programs[END_REF] gives :

∂ t U (t, x) + α(x -v(t)) 2 h(t) 1 x>v(t)+h(t) + 1 4 (β + ∂ x U ) 2 -rU = 0 (27)
One can now check that

U (t, x) = α r (x-β 2 t) 2 α r t+1
+ β 2 4r is solution of ( 27) for v(t) = β 2 t and h(t) = α r t + 1. Hence

x(t) = ( α r t + 1)x 0 + β 2 t
and the control

a(t, x) = α r (x -β 2 t) α r t + 1 + β 2 ,
or in other terms :

a(t) = a(t, x(t)) = α r (x(t)-β 2 t) α r t+1 + β 2 = α r x 0 + β 2 .
Note that the solution given in Theorem 2.1 is a static Pareto's Optimum. Indeed, any two researchers with initial positions given by x 0 and x0 have identical marginal rates of substitution :

∂ x U (t 1 , x(t 1 ) ∂ x U (t 2 , x(t 2 ) = α r (x(t 1 )-β 2 t 1 ) α r t 1 +1 α r (x(t 2 )-β 2 t 2 ) α r t 2 +1 = 2 α r x 0 2 α r x 0 = 2 α r x0 2 α r x0 = 1 (28)
3 Economic interpretation : Impact of the remuneration scheme on the density profile

This section studies the impact of the remuneration scheme on the distribution of the academic production over the continuum of researchers. Fig. 2 presents the densities m(0, x) and m(5, x) for the parameters value α = r = 0, 05 and β = 0 : Fig. 2 Recall that at date 0, the worst researcher has an academic resume of 1, at date t, his scientific production will be: x(t) = ( α r t + 1) + β 2 t. Fig. 2 illustrates two important features of the distribution over time : first, for all dates, the distribution of academic resumes is always Pareto although the law deforms over time. Second, the instantaneous scientific production of any researcher is constant over time according to a(t) = α r x 0 + β 2 . this speed is an increasing function of the initial CV of the researcher. In order to maintain his initial position, a researcher is induced to produce accordingly. Over time, the gap between scientific CV of researchers is growing. Moreover, the score of publication x(t) and the instantaneous scientific production a(t) are both increasing functions of α and β, and a decreasing functions of time preference. However, the consequences of a change in these parameters are contrasted. An increase in β has the same positive influence for every researcher, be they bad or good, whereas an increase in α will foster relatively more the academic production of the best scientists [START_REF] Lasry | Mean field games[END_REF]. Similarly, the influence of r is mediated through the initial position.

In order to measure the relative influence on the dispersion of the scientific production, we introduce the following measure:

F (t, x) = +∞ x+ β 2 t+ α r t+1 m(t, y)dy (29) = +∞ x+ β 2 t+ α r t+1 α r t + 1 (y -β 2 t) 2 1 y> β 2 t+ α r t+1 dy = α r t + 1 x + α r t + 1 (30)
which represents the number of researchers who exhibit an academic resume over and above the minimum production β 2 t + α r t + 1 plus an arbitrary level x of scientific production. 8 We can interpret this measure as the fraction of outstanding researchers.

One interesting aspect of F (t, x) is that it shows that only the parameter α influences positively the number of people having an academic resume of at least x. Instantaneous remuneration β has no influence on this measure.

Although it has an influence on the speed of the wave and hence on the total productivity of the whole community, the fraction of outstanding researchers can be increased by a rise in α.

we might guess that the speed at which the equilibrium distribution moves would slow down as at each period. The older generation of researchers who have in averaged the most important academic resume is replaced by a generation of junior researchers having weaker academic achievement. The ripple effect exerted on the speed of the distribution is more likely to be reduced. Second our approach oversimplifies the remuneration scheme. For instance we do not take into account the ratchet effect implied by some academic procedure such as tenure. Once tenure is acquired, a drop of the academic production has little impact over remuneration. Besides, we did not consider either the influence of ranking with respect to the instantaneous academic production which would allow us to take into account remuneration schemes of the type of the french "Prime d'Excellence Scientifique". Finally, scientific production is not the only output of academic researchers, the scientific method requires, in fact, that many researchers can devote a significant part of their research time not to make fundamental discoveries, but to validate the innovative research and to reveal, through their citations, the best ideas produced.

Despite these limitations our paper is interesting for at least two reasons. It provides first a theoretical foundation to the empirical stability of Lotka's law. Above this, in a period where governments intend to foster applied and academic research and aim at improving the efficiency of their research institutions, this paper brings some insight on the effect of various incentive schemes on scientific production. The Mean Field Game approach developped here allows to analyze the consequences of competition between researchers and to underlines the positive effects of the race for the best ranks both on the volume of scientific production and on the visibility of the best researchers. -2ẍ(t) + r2 ẋ(t) = α + rβ

The homogeneous equation admits a general solution Ae rt + B. Given the right hand side of the equation and the fact that we are searching admissible solution, one must have A = 0, which gives x(t) = x 0 + α+βr 2r t.

  Dubois et al 2010, Courtault et al 2011, Combe et Linemer 2001, Bosquet et al 2010).

5

  Appendix : Resolution of the model without externalities using variational calculusIt should be noted though that this model can be easily solved using calculus of variations. Replace a(t) by ẋ(t) in the optimization problem (1ẋ(t) -ẋ2 (t) e -rt dt/x(0) = x 0 which can be written as max (xt) +∞ 0 J(t, x(t), ẋ(t))dt/x(0) = x 0 where J(t, x, v) = αx + βv -v 2 e -rt . The Euler-Lagrange equation for this prob-β -2 ẋ(t))e -rt ) = αe -rt , -2ẍ(t) -r(β -2 ẋ(t)) = α, i.e.

[START_REF] Boyes | Publish or Perish : fact or fiction ?[END_REF] emphasizes that research output is more important for promotion than in-house activities like teaching or administrative tasks

xt and at could represent the number of papers weighted by the quality of the scientific reviews or any other bibliometric index.

See for exampleFleming anf Soner (2005) or[START_REF] Brock | Optimal Control and Economic Dynamics[END_REF]

1/m can be used as a ranking since m is decreasing as we will show that the equilibrium distribution is a Pareto distribution.

conclusionIn this paper, we use the Mean Field Game approach to study the impact of different types of remunerations over scientific production. One of the main result of the paper is that of the motivation of Lotka's law. Indeed we showed that Pareto's distribution is persistent. Once we start with a Pareto distribution for the academic resume, this distribution stays Pareto. However, the caracteristics of this distribution evolves over time. This result is achieved through competition among the researchers who aim at maintaining their rank within the distribution. One interesting result of the paper is that remuneration of ranking with respect to the academic resume influences both the speed at which the overall distribution moves over time, in the same way as piecework remuneration do, and more importantly boosts the productivity of outstanding researchers. In order to achieve these results, the model had to make some simplifying assumptions. First, we consider that researchers have an infinite horizon, a more appropriate model should consider overlapping generations. It is difficult to predict what will be the exact influence of this assumption on the equilibrium distribution. However,[START_REF] Combes | La production d'articles de recherche en économie en France[END_REF] Taking the number of researchers exhibiting an academic resume above x only would not make sense as for any given level x there is a time t for which the proportion of people exhibiting an academic resume above x is equal to one


Sébastien Alvo
email: sebastien.alvo@isir.upmc.fr

Jérôme Dejeu

Patrick Rougeot

Stéphane Régnier
email: stephane.regnier@isir.upmc.fr

Michaël Gauthier
email: michael.gauthier@femto-st.fr.

Jér Ôme Dejeu

Overview of coupling effects on interaction forces in micro-nano-world

Keywords: Microhandling, modeling, interaction forces

come   L'archive ouverte pluridisciplinaire

Overview of coupling effects on interaction forces

in micro-nano-world

I. INTRODUCTION
Nowadays, micro-grippers are mostly not forcecontrolled so forces around one thousand times greater than the adhesion force are commonly applied to guarantee the detachment of the two micro-objects. A better prediction of these forces could improve tele-operating manipulation thanks to haptic feedback [START_REF] Bolopion | Haptic Teleoperation for 3D Microassembly of Spherical Objects[END_REF]. Moreover, modeling these forces on the microscale is necessary for automatic tasks and for creating simulators for testing new manipulation strategies. A better knowledge of physical phenomenon enables to exploit their effects rather than undergo them. For example, (i) coating or functionalizing the handling area could prevent from adhesion forces or increase it [START_REF] Dejeu | Reducing the Adhesion between Surfaces Using Surface Structuring with PS Latex Particle[END_REF]; (ii) optimization of roughness can be achieved to minimize adhesion [START_REF] Dejeu | Reduction of micro-objects adhesion using chemical functionalisation[END_REF], [START_REF] Dejeu | Adhesion Forces Controlled by Chemical Self-Assembly and pH: Application to Robotic Microhandling[END_REF] . . . Micromanipulation is mostly based on top-down strategies so microgrippers are still catching micro-objects in the same way that macrorobots e.g. by pressing objects between two compliant fingers. However their design has to be adapted to the differences between macro-and microphysics. Moreover adhesion forces on the microscale are due to phenomena which have nanoscale origins. So micromanipulations use topdown tools but are modeled with bottom-up methods. Knowledge of the two scales can be used on the microscale but requires special attention on validity domains.

Phenomena that influence interactions between micro-objects have to be studied in order to highlight their relative influence according to the scale. Models from nanophysics presented in section II will be used to bring out each phenomenon influence on the microscale. Four specific phenomena influenced by scale effects between micro-and nanoscales are presented in this paper. It is shown in the section III that electrostatic forces prevail on van der Waals forces on the microscale. Section IV deals with the impact of chemical functionalisation on interaction force. Roughness and surface structuration are studied in section V. The section VI focuses on the coupling between the deformation and the van der Waals forces.

II. MODELING THE MICROWORLD

Differences between macro-and microphysics is now well established but what about the micro-and nanophysics ? As for macro-and micromanipulation, miniaturization of micromanipulation processes is not the right way to do nanomanipulation. Miniaturization is not just a technical issue; it also must take into account differences in physics.

A. From nano phenomena . . . The main issue of the micromanipulation is to release micro-objects properly. The corresponding force is the pull-off force which is the force to apply to break contact surface between two objects. This force is usually considered as normal to the contact. In dry environment they are induced by several phenomena as van der Waals, electrostatic and capillary forces [START_REF] Fearing | Survey of Sticking Effects for Micro Parts Handling[END_REF].

Each of them depend on several parameters (materials, roughness, humidity, temperature) and are due to nanomechanisms.

At the atomic level, van der Waals force is the attraction force between two single atoms. It is commonly assumed to be an additive force so for two objects S 1 and S 2 distant from r, the total van der Waals force is the double sum of forces applied on atom of S 2 by atoms of S 1 [START_REF] Hamaker | The London -van der Waals attraction between spherical particules[END_REF], [START_REF] Israelachvili | The Nature of van der Waals Forces[END_REF], [START_REF] Parsegian | Van der Waals forces[END_REF]. Adhesion is currently modeled by considering a minimal distance z = z 0 between both objects. z 0 is usually set to 0.3 nm which is the equilibrium distance between two atoms in vacuum. All along this paper, we choose to illustrate phenomena with adhesion between a glass sphere of radius r 2 and a plane. The van der Waals force induced by this geometry is:

F vdw = - Ar 2 6z 2 , (1)
where z = z 0 for contacting objects and A is the Hamaker coefficient (A = 6.5 10 -20 J for glass). Electrostatic forces are due to surface charges added by friction or chemical treatment (cleaning, oxidization . . . ). They are modeled in the same way as the van der Waals forces replacing volume atoms distribution by a surface charge distribution. The electrostatic force applied by the plane (surface charges σ 1 ) on the sphere (surface charges σ 2 ) is:

- → F pull-in = q 2 . - → E 1 = 2πr 2 2 σ 1 σ 2 3 0 -→ n 1 ( 2
)
where 0 is the vacuum permittivity. Capillary forces F capil are due to water meniscus created thanks to ambient humidity or adsorbed layers. This paper does not deal with capillary effects.

B. . . . to the microscale modeling

The models presented above can be used to model pull-off forces on the microscale. So the most natural way to model pull-off force is to sum them:

F po = F vdw + F es + F capil . (3)
This bottom-up model of pull-off forces let us understand origins of the adhesion. Top-down models are also currently used to model adhesion. Energetic modeling derived from continuum mechanics can be used for estimating pull-off forces between objects with simple geometries. The most famous models are the JKR [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF] and the DMT [START_REF] Derjaguin | Effect of Contact Deformations on the Adhesion of Particles[END_REF] models. They are wildly used but they do not consider the influence of electrostatic and capillary forces. This kind of model is not flexible: adding a complementary force or changing geometry is complex because of the model computation method.

C. Force measurements

Pull-off forces are measured thanks to an atomic force microscope (AFM). In the present case, beads are fixed on an AFM cantilever (Fig. 1) whose relative position to the studied plane is controlled by a micropositioning device. The cantilever has a low stiffness (between 0.2 N/m and 0.3 N/m in our case) in order to measure weak non-contact forces.

10µm Bead

AFM cantilever

Fig. 1. A microbead is fixed on an AFM cantilever in order to measure pull-off forces. Force measurements have been performed on gold coated beads (see photo) and silica beads.

In this article, we chose to take advantage of both the bottom-up and the top-down approaches to model pull-off forces on the microscale. The flexibility of eq. ( 3) will let us studying the effect of each part of the pull-off force according to the nanoscale point of view. Moreover, the deformations will be taken into account too in order to stay close to macroscopic theories.

III. SCALE EFFECT ON ELECTROSTATIC FORCES

There are usually fluctuations in charges distributions on micro-nano objects so associated forces are not repeatable and then are currently considered as a perturbation. As an example, the maximum amplitude of electrostatic force for silica beads is first estimated in this section. The resulting distribution of charges is then used to illustrate the scale effect on electrostatic forces.

Two kinds of attractive non-contact forces can be observed between micro-objects: the van der Waals forces and the electrostatic forces. They can be measured when pull-in occurs. Pull-in is a jump to the contact that occurs when the bead is approached to the plane. The corresponding force is defined as the strongest attractive force observed just before the contact. One way to determine the predominant force consists in looking at the interaction range. Indeed electrostatic forces have a longer interaction range than van der Waals forces (typically several micrometers to tens of nanometers for van der Waals forces). In case of an interaction between a silicon plane and a borosilicate sphere, typical surface charges about σ = 0.15 mCm -2 have been observed experimentally [START_REF] Alvo | Analysis and Modeling of Scale effects at the Frontier between Micro and Nanohandling[END_REF]. This level of charge induced by tribolectrification is typically lower than those obtained by chemical functionalisation in liquid (see section IV).

A. Scale effect consequences on pull-off measurements

In order to evaluate the effect of electrostatic forces on adhesion phenomena, we chose to compare them to the other cause of adhesion: the van der Waals forces. In order to study the influence of scale effect on forces, it can be noticed that the van der Waals forces are proportional to the radius of the bead (1) whereas, according to (2), electrostatic forces are proportional to its square. This difference induces different scale effects on both phenomena. In the case of borosilicate glass, it can be seen on Fig. 2 that the influence of the electrostatic forces have to be considered for objects whose size is greater than 1 µm. On nanoscale, the proportion of the pull-off force due to the electrostatic forces becomes negligible compared to these due to van der Waals forces.

Fes << Fvdw Fig. 2. Electrostatic (solid line) and van der Waals (dashed line) forces between a sphere and a plane according to the sphere radius. The electrostatic force is considered as negligible as it is lower than the tenth of the van der Waals forces (dash-dot line). Considered materials are borosilicate glass, with A = 6.5 10 -20 J, r 0 = 0.3 nm and σ = 0.15 10

-3 C/m 2 .
It is possible to control the charge density on the surface by chemical treatments [START_REF] Dejeu | Reduction of micro-objects adhesion using chemical functionalisation[END_REF] (see section IV below) or using conductive materials. However, in the case of non-conductive and non-treated materials, the charge density could not be controlled and fluctuate. On microscale, the adhesion is linked to the charge density and its variations; it is therefore very difficult to make repeatable pull-off measurements with nonconductive micro-objects.

At the nanoscale, it is possible to measure the van der Waals forces between non-conductive materials without variations observed on the microscale because the electrostatic forces can be considered as negligible compared to the van der Waals forces. This explains why the van der Waals theory has been validated experimentally on the nanoscale for all types of materials, even non-conductive.

IV. CHEMICAL FUNCTIONALISATION AND ELECTROSTATIC FORCES

Based on the fact that electrostatic forces can be a predominant effect on microscale. Electrochemical effects can be exploited in order to control electrostatic surface density via chemical equilibrium between the surface and the liquid medium. Surface functionalisation of both object and gripper can be obtained by different methods (physisorption, grafting...). Two chemical functionalisations have been tested (see in figure 3):

• the silane, 3 (ethoxydimethylsilyl) propyl amine (APTES); • the silane, (3 aminopropyl) triethoxysilane (APDMES). Both chemical compounds (APTES, APDMES) used to surface functionalisation are amine functions NH 2 which can protonated or ionised to NH + 3 according to pH. In acidic pH, the anime is totally ionised, then the ionisation decreases and is null in basic pH (between pH 9 and 12). Our principle is also based on the protonation of silica, which enables to switch from SiO 2 to SiO -according to the pH. The combination of both effect enables to obtain a surface whose electrostatic charges switch from a positive value for low pH to a negative value for high pH. Pull-off force measurements were done on a functionalized plane using a cantilever with a nonfuntionalised sphere. The results of the pull-in and pulloff forces is presented in figure 4.

In this figure, we noted that the pH influences significantly the forces between the cantilever and the surface. At natural pH, a attractive pull-in is measured (near -60 nN, Figure 4(a)) with an important pull-off (-350 nN, Figure 4(b)). When the pH increases the pull-in force is inverted and becomes repulsive respectively 280 nN and 770 nN at pH 9 and 12. Moreover, the adhesion forces disappear. In order to be able to extrapolate this result to other geometries, the electrical surface density induced by the functionalisation have been studied. The repulsive electrostatic force applied by the plane (surface charges σ 1 ) on the object (surface charges σ 2 ) is defined in [START_REF] Dejeu | Reducing the Adhesion between Surfaces Using Surface Structuring with PS Latex Particle[END_REF], where surface charges have been established in [START_REF] Dejeu | Reduction of micro-objects adhesion using chemical functionalisation[END_REF], [START_REF] Dejeu | Adhesion Forces Controlled by Chemical Self-Assembly and pH: Application to Robotic Microhandling[END_REF] (see in table I). An other important parameter which modifies significantly the adhesion is the roughness. Using nanostructuration of surfaces, the roughness of the surface can be controlled. It enables to study properly the impact of roughness on interaction forces and also to control the roughness in an application case.

Let us consider the nanostructure described in figure 5 which represents the position of nanoscaled PS spheres on a surface. These nanostructurations have been built by the EMPA institute, Thunn, Switzerland [START_REF] Dejeu | Reducing the Adhesion between Surfaces Using Surface Structuring with PS Latex Particle[END_REF] using self-assembly methods. In an application case and also during force measurements the location of the sphere on the probe up the structured surface cannot be controlled precisely. When the probe r 1 is approaching, it touches the nanospheres r 1 on a non-controlled position. We have shown that the force is included between a minimum and a maximum which verifies [START_REF] Dejeu | Reducing the Adhesion between Surfaces Using Surface Structuring with PS Latex Particle[END_REF]:

F T vdwmin = Z 2 i,j A 12 r 12 6z 2 ij min . r 2 + z 0 + r 1 r 2 + z ij min + r 1 (4)
F T vdwmax = Z 2 i,j A 12 r 12 6z 2 ij max . (r 2 + z 0 + r 1 ) 2 -(4/3.r 2 ) r 2 + z ij max + r 1 (5)
where:

z ij min = ((r 2 + z 0 + r 1 ) 2 + 4r 2 2 (j 2 -ij + i 2 )) 1 2 -r 1 -r 2 (6) z ij max = ((r 2 + z 0 + r 1 ) 2 + 4r 2 2 (j 2 -ij -i -j + i 2 )) 1 2 -r 1 -r 2 (7) r 12 = r 1 .r 2 /(r 1 + r 2 ) (8)
The comparison between value predicted by the model and the measurement, plotted in 6, shows a good concordance. So 90 % of the experimental points validate the model. The other 10% of the experimental points are very near to the predicted value, just few nN below the model. The second results deals with the determination of a minimum of the interaction force which represents an optimum of adhesion reduction in the applicative field of micromanipulation. In our experimental case, the optimum radius r 2 in order to minimize the adhesion is between 45 and 100 nm. If the radius r 2 is lower than this optimum, more and more spheres should be considered in the sum thus increasing the force. This value depends of the diameter of borosilicate sphere . glue to the cantilever. The model can be extended to different diameter of sphere r 1 using ( 4) and ( 5). The proposed model can be used to determine the diameter of the optimal PS particles spheres to be placed on a gripper to minimize adhesion force with a grasped sphere S b . Experimentally, some PS spheres has been deposed on silicon grippers [START_REF] Agnus | Silicon endeffectors for microgripping tasks[END_REF] (see in figure 7).

PS spheres nanostructure

Surface in contact with the grasped object Van der Waals forces (1) are usually computed on non deformed objects but local deformations are able to significantly increase their value on microscale. Moreover, the greater part of van der Waals forces is due to the contact area, where both objects are the closest. Deformations should therefore be taken into account for the force calculation. This section shows that the coupling between van der Waals forces and deformations is also a characteristic of the microscale.

A. Coupling principle

Because of its short interaction range, the van der Waals forces are close to be a surface force. In this model, we chose to replace this local force by an equivalent external load which induces a deformation. This deformation increases the contact surface and then the global van der Waals forces too. This coupled problem can be seen as an algorithm that uses sequentially two models (Fig. 8). The first one computes van der Waals forces according to the object shape. The other one computes deformation shape according to an external load. An iterative calculation is able to converge to the physical equilibrium. Fig. 8. Algorithm proposed for calculating the adhesion force F dvdw between two objects using the coupling between deformation and van der Waals forces. The algorithm starts with non-deformed objects so the initial contact radius a 0 is set to zero. The van der Waals force Fn can be computed as the sum of a deformable dependent part F defo and a non deformable one F vdw [START_REF] Alvo | A van der Waals Force-Based Adhesion Model for Micromanipulation[END_REF]. The deformation model used comes from the Hertz theory.

This section deals with the contact between a sphere and a plane as an example. In order to highlight the influence of the scale effects on deformations, the contact geometry will be simplified. We chose to replace the deformed sphere by a truncated sphere. The high δ of the removed spherical hat is calculated thanks to the Hertz deformation model. It is then possible to calculate the van der Waals force between the truncated sphere and the plane according to δ (details of the calculation can be found in [START_REF] Alvo | A van der Waals Force-Based Adhesion Model for Micromanipulation[END_REF]) and then close the model (Fig. 8).

B. Model of the coupling

The adhesion force F dvdw calculated using this modeling principle is presented on Fig. 9 in the case of silicon objects. A critical radius can be extract from the model equations [START_REF] Alvo | A van der Waals Force-Based Adhesion Model for Micromanipulation[END_REF]:

R c = 2E 2 z 7 0 (1 -ν 2 ) 2 A 2 , (9)
where E and ν are the mechanical characteristics of the objects. In the case of glass objects, R c = 0.5 µm. If the sphere radius is smaller than this critical radius R c , the computed pull-off force tends to be the classical van der Waals force presented in eq. ( 1). So, on the nanoscale, the influence of deformation on van der Waals forces becomes negligible.

R c F ≈ F v dw F ≈ F dv dw
Sphere radius (µm) Force (µN) influence of deformations Fig. 9.

Comparison of forces computed with classical and deformable van der Waals theories on nano and microscales. The numerical solution of the deformable van der Waals (dvdw) model (solid line) matches classical van der Waals (eq. ( 1), dashed line) on nanoscale and analytical formula (eq. ( 10), dash-dotted line) on microscale.

The influence of deformations on the pull-off force has to be only considered on the microscale. We have shown in [START_REF] Alvo | A van der Waals Force-Based Adhesion Model for Micromanipulation[END_REF] that on the microscale, for r 2 R c the adhesion force tends to be:

F dvdw = - A 3 r 2 2 48z 9 0 E * 2 . ( 10
)
The difference between both force models is increasing with the radius of the sphere. This asymptotic model is based on analytic formulation. A full numerical model that considers exact shapes thanks to finite element models should rather have a R k (k > 0) scale dependence. This second scale effect highlights another characteristic of the microscale physics where deformations significantly increase van der Waals force.

VII. CONCLUSION

The micrometer seems to be a critical size where physics is changing and where couplings are appearing between several physical effects. Four features of physics on the microscale were highlighted in this paper: (i) electrostatic forces disturb micromanipulation whereas it seems to be negligible on the nanoscale, (ii) surface functionalisation enables to modify significantly interaction forces, (iii) roughness is a critical parameter of the pull-off force and (iv) the material deformations strongly increase the pull-off forces only on the microscale. Because of those scale effects, pulloff force measurements on the microscale present large variations whereas forces could be predicted easier on the nanoscale. It is well established that a top-down strategy is not well adapted to micromanipulation. We have shown that a pure bottom-up strategy based on the extrapolation of nanosciences is not a better way. Both macroscale and nanoscale effects have to be considered for microphysics modeling.

Fig. 3 .

 3 Fig. 3. Molecules used for the silica functionalisation.

Fig. 4 .

 4 Fig. 4. Force-distance curve for the APTES functionalised substrate in wet medium at different pH obtained with a tip whose spring constant is 0.3 N/m.

Fig. 5 .

 5 Fig. 5. Arrangement of the PS spheres on the substrate.

Fig. 6 .

 6 Fig. 6. Comparison between the model (minimal force, blue solid line, and maximal force, red dash line) and experimental measurements (error bar) on the structuring surface for borosilicate 20 µm in function of the nanostructuration r 1 .

Fig. 7 .

 7 Fig. 7. Structuring gripper by PS particles of 1 µm : Joint work between EMPA institute, Thunn, Switzerland, and FEMTO-ST institute, Besanc ¸on, France [2]

ACKNOWLEDGMENT

This work was partially supported by the EU under HYDROMEL (contract NMP2-CT-2006-026622): Hybrid ultra precision manufacturing process based on positional-and self-assembly for complex microproducts, and FAB2ASM (contract FoF-NMP-2010-260079): Efficient and Precise 3D Integration of Heterogeneous Microsystems from Fabrication to Assembly, and by the French National Agency (ANR) under NANOROL (contract ANR-07-ROBO-0003): Nanoanalyse for micromanipulate. We would also like to acknowledge the support of EMPA institute in nanostructuration.