Changyi Xu

Eric Niel
email: eric.niel@insa-lyon.fr

Nicolae Brinzei
email: nicolae.brinzei@univ-lorraine.frchangyi.xu

Discrete Event System Formal Approaches Contribution onto Global Reliability Markov Chain Generation Contribution des approches formelles des systèmes à événements discrets pour la génération des modèles markoviens de fiabilité

published or not. The documents may come

Discrete event system formal approaches contribution onto global reliability Markov chain generation

Background

Model Based Safety Assessment (MBSA) is an inherent task in most safety-critical engineering projects, it is a guarantee for the safe-design. Appreciating the effectiveness of a system through its faulty behavior model, by considering both the operating properties and system architecture, is a common basis for MBSA [START_REF] Niel | Maitrise des risques et surete de fonctionnement des systemes de production[END_REF]. Conventionally the establishment of MBSA requires professional skills for engineers and contains four major steps: system modeling, simulation for validation, fault logic expression integration, performances assessment [START_REF] Joshi | Model-based safety analysis[END_REF]. Linked to Model Based System Engineering (MBSE) for which developments are objects of interconnected models MBSA considers various criteria conform to operating capacity and causality [START_REF] Wymore | Model-based systems engineering[END_REF]. Fault propagation relies to functional and operational structure are then totally taken into consideration. The location arrangement of the local components establishes logic connections (basic structure as parallel and series operating structure). Each faulty behavior of one local component directly influence the global behavior of the whole system through that connection. So the global system behavior results of diversity influenced both by the local component performance and system architecture. Moreover, the whole system structure can be appreciated into series, parallel or more complex operating relationship. Indeed, related to granularity each local component operates in collaboration following an admissible principle. Also related to those logic connections, when several local components break leading to the total failure of the whole system. This principle defines the fault logic expression of this faulty system. So analyzing the system architecture to fulfil the system logic expression and then to establish a model strictly satisfying all the logic expression belongs to MBSA approach [START_REF] Ford Jr | Network flow theory[END_REF]. The proposal aims to produce a formal as well as automatic oversimplified approach coherent to MBSA with application to Markov chain reliability model generation.

Research Introduction

Thanks to Reliability Block Diagram (RBD) [START_REF] Čepin | Reliability block diagram[END_REF] and Fault tree Analysis (FtA) [START_REF] Ericson | Fault tree analysis[END_REF], the capacity and complexity of the faulty system can be analyzed and the Boolean Value fault logic expression formula can be elaborated, in which several Boolean Value semantics are connected by conjunction and disjunction logic symbols. Markov chain model is used for representing the behavior of local components (operating state, failure state, repair transform, failure transform, etc). Each local Markov chain can be translated into event driven local automata for the purpose of composition and satisfying deterministic requirement for dependability analysis [START_REF] Cassandras | Introduction to discrete event systems[END_REF]. A set of Boolean Value state properties is established consistent with the local component behaviors and a formal associate function is established for the automata states and these properties [START_REF] Baier | Principles of model checking[END_REF].

GFA is achieved as the automata composition combining faulty behavior. Via the automata composition, the local Boolean Value state properties are generated to be the composed state properties Boolean value in GFA, which are the representatives of the combining faulty behavior [START_REF] Cassandras | Introduction to discrete event systems[END_REF]. A judgment calculating operation is used to check each Boolean Value property of each state in GFA whether to be a sufficient condition of the fault logic expression formula (the judgment operation can be performed by tools, for example Reduced Ordered Binary Decision Diagram) [START_REF] Weaver | Extending existential quantification in conjunctions of BDDs[END_REF], if the result is true, it implies that this state belong to TBS. The TBS identification contribution is achieved by this judgment operation applied to all the states in GFA. System operation requirements expression are issued by the synchronization of GFA and event trajectory specifications (ETS). The function of ETS is to specify the GFA following the operation requirements, and ETS state evolution trajectory is formally planned by several specification automata work function establishment based on the requirements [START_REF] Cassandras | Introduction to discrete event systems[END_REF]. TBS treatment lies in the blocking operation applied to the identified TBS, automatically prosecuted by the tool e Congrès de Maîtrise des Risques et Sûreté de Fonctionnement λµ21 Reims 16-18 octobre 2018 Supremica [START_REF] Fabian | Supremica--An integrated environment for verification, synthesis and simulation of discrete event systems[END_REF]. After synchronizing GFA and ETS, GRA respecting the system operation requirement is established by the TBS identification and TBS treatment. The global reliability Markov chain model is by a further translation issued from GRA and mean time indicators could be then prosecuted [START_REF] Niel | Maitrise des risques et surete de fonctionnement des systemes de production[END_REF]. Figure 1 is the work flow of this research and the main contributions will follow 5 major steps: 1. Transform Markov chains to event driven automata. The Markov chains model of the local faulty components includes physical performance as states (operating state, failure state, etc) and behaviors as transitions (failure event transition, repair event transition, etc). This Markov chain can be translated into local dependability event driven automata. This translation satisfies the deterministic dependability analysis requirement. The state Boolean Value properties is established consistent with physical performance. 2. Cartesian product [START_REF] Eilenberg | Automata, languages, and machines[END_REF] of automata to produce GFA. In order to establish GFA, the Cartesian product is used for the composition operator of all the local dependability models and state Boolean Value properties is generated via composition. 3. TBS identification. By analyzing the system architecture using RBD and FtA, the system failure logic Boolean Value expression formula is established. Judgement operation of state Boolean Value properties comparing to System Failure Logic Expression Formula is applied for TBS identification [START_REF] Baier | Principles of model checking[END_REF][10] 4. Operation Requirement expression. Establishing formal work function of ETS according to the operation requirement, the ETS is able to specify event evolution trajectory in GFA to be requirement respecting [START_REF] Cassandras | Introduction to discrete event systems[END_REF]. 5. GRA issuing. After synchronizing ETS and GFA, the TBS treatment is applied to the model and GRA is established with TBS blocking operation.

Markov chains and Event Driven Automata Transformation

Shown in Figure 2 left, a Continuous Time Markov Chain representing a fail-repair component with a failure rateλf1 , a repair rate μ r1 , initial state S0 and failure state S1.

According to Arden's lemma, Continuous Time Markov Chain is able to be transformed as a sequence of unique events. In equation {1}, A and B are the universe of event, Li is one of the sub languages representing a system, A and B are the universe of system unique events. If there is no empty sequence in A, the solution of {1} is unique Li=BA*. If an empty sequence belongs to A, the equation solution is Li=(B+C)A*, there C⊆E* and E is the event set.

Li= LiA + B {1}

Normally, the Markov chain transition is labelled by an evolution rate and as the Markov chain can also be expressed as a sequence of unique events, so the transition labelled by happening rate and unique event is shown in figure 2 middle: λ f1 is presented by failure behavior event f1 with a failure rate λ; μr1 is presented by repair behavior event r1 with a repair rate μ [START_REF] Ionescu | Evaluation quantitative de séquences d'événements en sûreté de fonctionnement à l'aide de la théorie des langages probabilistes[END_REF]. Moreover, regardless the happening rate, event is established as component behavior, the focus on the event function results the event driven automaton like figure 2 right side.

� L1 = L2f1 L2 = L1r1
So, L(1) = L1f1r1 and the solution is:

� L1 = (f1r1) * L2 = (f1r1) * f1
The alphabet Σmc expalined by the work function is expressed:

� δmc (S0, f1) = S1,

GFA Establishment 1 Local faulty component model establishment

As it is introduced, principle for modeling a component is: performs to be automata states, behaviors to be automata transition events. Figure 3 shows two local faulty components by the help of the transform between continuous time Markov chains and Event driven automata modeling for component G1 and component G2. To consummate the component perform property, Boolean value state property is adopted to improve the state significance description. Boolean value state property is in fact Boolean value symbols of state characteristic, for example when there is a need to denote one state on the failure characteristic, a Boolean value symbol < F > can be created to associate with this state, and if there is a need to

L3(A,C)=L1(A)∧L2(C)= <¬ F1∧¬ F2> L3(A,D)=L1(A)∧L2(D)= <¬ F1∧ F2> L3(B,C)=L1(B)∧L2(C)= <F1∧¬ F2> L3(B,D)=L1(B)∧L2(D)= <F1∧ F2>
The GFA is achieved by the model generation of local faulty automata. Consummating the expression of system performs, Boolean value state property is associated with each local state and via the model generation of local automata, the global states in GFA are also associated with Boolean value state property, this state property not only has a full expression of the GFA system perform but also is the key factor for further TBS identification and TBS treatment.

TBS identification 1

System failure logic expression formula obtaining Thanks to basic system dependability analysis technologies, Reliability Block Diagram (RBD) [START_REF] Čepin | Reliability block diagram[END_REF] and Fault tree Analysis (FtA) [START_REF] Ericson | Fault tree analysis[END_REF], system architecture is able to be analyzed. Based on functional dependency, the simplest structure are parallel or series in the application of RBD and FtA. System Failure Logic Expression Formula provides the system Boolean Value failure property in logic connection form will help for the TBS identification aim. There exists direct translation RBD safe structure between FtA by the relation that parallel structure in RBD is equal to 'and' gate in FtA and series structure in RBD is equal to 'or' gate in FtA. By connecting the local system Boolean Value failure property of the system with combination of "and" gate and "or" gate, from the FtA, we can get the System Failure Logic Expression Formula of this global system with the help of computer calculation software, such as GRIF.

Figure 5. System Failure Logic Expression Formula

System Failure Logic Expression Formula is expressed as the form that a set of Boolean Value failure property connected with the symbol '∨'( logic connection "or" ) and with the symbol '∧' ( logic connection "and" ), showing how the local components breakdown contribute to the globally total breakdown of this architecture. Suppose system architecture of our two components G1 and G2 are working in the hot redundancy form: only if G1 and G2 are both breakdown, the whole system breaks down. In figure 5, as the parallel safe structure in RBD, in FtA there is an 'and' gate connected F1 and F2, here as introduced before F1 and F2 are denoted as Boolean Value failure property of two local components. Based on the calculation principle of FtA: 'and' gate result in a conjunction '∧' of the Boolean Value properties; 'or' gate results in a disjunction '∨' of the Boolean Value properties, the System Failure Logic Expression Formula is expressed: Ψ=F1∧F2 Ψ is prepared for the TBS identification reference which is compared to the state Boolean Value property in this approach. And also for the case of operation requirement expression designing specification automata ETS, as the event trajectory ETS will be possible formalized by a sequence of failure events exciting Ψ happening.

TBS identification judgment method

System Failure Logic Expression Formula is established through RBD and FtA or a given formula by designer, identification judgment method detail is issued by the compare of the Boolean value state property of each state in GFA with System Failure Logic Expression Formula, the states whose significance satisfys the total breakdown condition will result that its Boolean value state property is a sufficient condition of System Failure Logic Expression Formula and this state is operated to be extracted from the global state set into a TBS state set. Judgment method mathematical formula exposition is as following:

Given a automaton G = (Q; E; δ; q0; qm; AP; L) The Boolean value state property of one state satisfying the System Failure Logic Expression Ψ contributes to the TBS state set QTBS.

QTBS={x| x∈Q, L(x) ⊨Ψ} {8}

Here, the satisfaction symbol '⊨' is applied to denote L(x) is a sufficient condition of Ψ and QTBS is a subset of Q: QTBS⊆Q. The Boolean value state property satisfaction judgment is applied in the verification technique of almost automata model checking tools. In our approach, a brief judgment calculation tool named BDDs (Reduced Ordered Binary Decision Diagrams) is applied for the TBS identification of GFA. In BDDs, the state property is coded as a combination of Boolean value property elements, the identification formula is coded as disjunction or conjunction of system Boolean value property elements like formula Ψ , where state property value and identification formula are both consistent as what they are defined. The sufficient condition judgment satisfaction operation is implemented by the BDDs code order operation 'imp' which means 'imply'. Submitting BDDs code, if the result only show an assignment 'T', this means the state property value is a sufficient condition of the identification formula, which indicates that this state is identified to be TBS and this state should be extracted into QTBS. If BDDs result shows not 'T', this state is not a satisfying state (The result will shows a diagram telling how Boolean Value develops will be 'T' or how it will be 'F', 'F' meaning false, this result sort is not our concerning knowledge) [START_REF] Weaver | Extending existential quantification in conjunctions of BDDs[END_REF].

Operation Requirement Expression

Normally the system operation is required to respect several principles, thus an operation requirement context is established to instruct the system operators. And there are various sorts of operation requirements, for example: when there happens a series of breakdown accidents onto the local components, the repair action should be prosecuted, then a problem of repair action before and after order is available, for one certain case "One important component has the priority to be first repaired" and other rule that "FIFO: the component firstly breaks down and then following another component breaks down secondly, requirement context rules that respect the time order first breaks first repaired". Operation Requirement Expression means GFA behaviors should respect the operation requirement context, which is issued by several formal specification automata (ETS) synchronizing with GFA, specifying GFA transition evolutions.

For the problem that requirement context various as different cases various, it is not possible to establish all the ETS in the world, but a formal work function of ETS respecting a sequence events happing can be established. Two cases are used to explore our approach.

"One important component has the priority to be first repaired". Assume that our faulty component G1 and G2, component G1 should be first repaired whenever G1 breaks down, so this requirement context should be established in ETS by a sequence of failure event (denoted by SFE) happening: SFE1=" f1 " This means that if f1 happens, the operation requirement is excited to be effective. Given the definition ETS1= (QSP1; ESP1; δSP1; q0SP1; qmSP1), QSP1 is the state set {S0, S1}; ESP1is the event set ESP1=E1∪E2={f1,f2,r1,r2}; q0SP1 is initial state S0; qmSP1 is marked state S1, the transition sequence route from initial state to marked state is SFE1=" f1 "; Shown in Figure 7, the specification ETS1 is designed like: from the initial S0 state after SFE1 S1 is reachable and requirement context is effective at S1, where S1 has the selfloop {f1,f2} with meaning only r1 is allowed (r2 is forbidden here), the self-loop of S0 is {r2, r1, f2} is of the meaning all the events except SFE1="f1" happening makes no requirement context excited at S0. By SP1 whenever f1 happens, r1 is only allowed, which is equal to "component G1 should be first repaired whenever G1 breaks down". Moreover the simulation result proof our method of the result "whenever f1 and f2 have both happened, r1 is priority to be excited than r2" FIFO: Assume that our faulty component G1 and G2, G1 first breaks down and then following G2 breaks down, G1 is of priority to repair. So this requirement context should be established in ETS by a sequence of failure event (denoted by SFE) happening:

SFE2=" f1 f2 "

This means that f1 happens and then f2 happens, the operation requirement is excited to be effective. Shown in Figure 8, the specification ETS2 is designed like: from the initial S0 state after SFE2 S2 is reachable and requirement context is effective. By the self-loop designing of S2 the only repair event allowed is r1, which is equal to "G1 first breaks down and then following G2 breaks down, G1 is of priority to repair" and the simulation result shows the correctness in figure 8 lower part.

TBS treatment for GRA establishment

The establishment principle of a reliability model is to convert break down state as an absorbing state thus the TBS treatment lies in the GRA establishment by converting TBS output transitions deleted, operated on GFA. Thanks to a new application of automata tool Supremica, after TBS is identified, the TBS treatment is implemented automatically and formally.

According to equation {8}, with the hot redundancy working form system architecture,the TBS set is Q3TBS={ (B,D) } Assume that the FIFO operation requirement: "G1 first breaks down and then following G2 breaks down, G1 is of priority to repair" is respected, so ETS2 shown in figure 8 is used for the final GRA establishment.

Based on the work flow in figure 1, we synchronize GFA with ETS2: G3//ETS2=G1//G2//ETS2 shown in figure 8. Up to this step, our GFA model is already applied by Operation Requirement Expression.

Conclusion

The related contributions belongs to MBSA by combining requirement and faulty structures. The major steps are local behavior model, requirement expression, TBS identification and TBS absorbing. They all satisfy the incompatibility properties for dependability analysis. By the help of system assessment tools, such as GRIF, the mean time indicators assessment is possible. Furthermore, all models can be simulated and are able to validate the design steps.

Figure 1 .

 1 Figure 1. Work Flow

Figure 2 .

 2 Figure 2. Markov chain and Event Driven automata transform According to equation {1}, shown in figure 2 <Transition rate with unique event association>, Markov chain model is defined by CTMC=( Qmc; Σmc; δmc; q0), Qmc is the state set{S0,S1}, Σmc is the alphabet which is labelled by two parameters rate{λ,μ} and the event name{ f1, r1}, qo is the initial state{S0}; δmc is the transition function. The proof of formalizing CTMC is based on equation {1} solution :

  by the rate λ δmc (S1, r1) = S0, by the rate µ Event driven automaton transformed from the Continuous Time Markov Chain is defined: G=( Q; E; δ; q0; qm ): Q is the state set: Q= Qmc={ S0, S1} ; E is the event set: E= { f1,r1 }; Initial state q0= S0; Marked state qm= S1; And the work function compared to δmc is expressed:� δ (S0, f1) = S1 δ (S1, r1) = S0By this transform method, the faulty component model previously delivered by continuous time Markov chains is able to be delivered by associated event driven automata. The purpose of this method is to present transition rate to component behavior marked by visual system event with the ability prepared for model composition, by the reason composition operation of event driven automata offers a formal solution for model generation, moreover after the global faulty automaton model has been established the associated Markov chain is able to be achieved for global system time indictor assessment.

2

  Model generation by Cartesian productModel generation is issued by the local automata composition operation theory to establish a global automaton. Composition operation is in fact the Cartesian product, assume two automata for our components: G1 = (Q1; E1; δ1; q01; qm1) and G2 = (Q2; E2; δ2; q02; qm2), we use symbol '//' to denote the automata composition and symbol '×' to denote the Cartesian product, the composition operation of two automata is defined[START_REF] Cassandras | Introduction to discrete event systems[END_REF][11]: G3= G1//G2 = Ac(Q1×Q2;E1∪E2;δ3; q01×q02;qm1×qm2) {3 } Ac stands for taking the accessible part.

Figure 4 .

 4 Figure 4. Automata Composition And the work function δ3 issued by Cartesian product is:

Figure 6 .

 6 Figure 6. state (B,D) from G3 BDDs method

  Then based on SFE1 the work function δSP1 is defined: Evolution function: δSP1 (So,f1)=S1 δSP1 (S1,r1)=S0 Self-loop function: δSP1 (S0, ESP1\f1)=S0 δSP1 (S1, ESP1\r2)=S1 here, '\' denotes for complement set operation.

Figure 7 .

 7 Figure 7. specification ETS1

  Given the definition ETS2= (QSP2; ESP2; δSP2; q0SP2; qmSP2) QSP2={S0,S1,S2}; ESP1=E1∪E2={f1,f2,r1,r2}; q0SP2=S0; qmSP2=S2; Then based on SFE2 the work function δSP2 is defined: Evolution function: δSP2 (So,f1)=S1 δSP2 (S1,r1)=S0 δSP2 (S1,f2)=S2 δSP2 (S2,r1)=S0 Self-loop function: δSP2 (S0, ESP2\f1)=S0 δSP2 (S1, ESP1\r1\f2)=S1 δSP2 (S2, ESP1\r1\r2)=S2 here, '\' denotes for complement set operation.

Figure 8 .

 8 Figure 8. specification ETS2

Figure 9 .

 9 Figure 9. GRA establishment Continuing operating on GFA, normally the TBS identification is issued by the judgment of state Boolean value property with System Failure Logic Expression Formula: Ψ=F1∧F2 of the GFA shown in figure 8 lower part like what it is introduced in the TBS identification part. But considering ETS2 is only established by event trajectory specifying without any new Boolean value property added into it (in ETS2, S0 and S1 have no association with state Boolean value property), there is no local state Boolean value property generation via the composition of G3 and ETS2. So to be brief, the automaton state Boolean value property reminded to be the same before or after specification automata composition, state (B,D,S1) and state (B,D,S2) share the same state Boolean value property: L(B,D,S1)= L(B,D,S2)= L3(B,D)=<F1∧F2>.


Mohamed El-Demerdash

Pedro Lealdino Filho

Christian Mercat

KINESTHETIC PROMOTION OF FUNCTION GRAPH RECOGNITION AT UNIVERSITY LEVEL

Keywords: Demerdash, Lealdino, Mercat

published or not. The documents may come

INTRODUCTION

The emergence of motion-controlled technologies within the increasing usage of embodied cognition and augmented reality environments opens new doors for students to use their different senses into the learning process [START_REF] Mercat | La diffusion: un lieu pour une mathématique plus humaine? Paper presented at EMF[END_REF]. The goal besides making them more motivated to pay attention to the information presented and better retain the information [START_REF] Tsai | The Development and Evaluation of a Kinect Sensor Assisted Learning System on the Spatial Visualization Skills[END_REF] is to let students experience an alternative way to build concepts.

The concept of function is unifying idea in mathematics. It plays an important role throughout mathematics curriculum [START_REF] Nctm | Standards for school Mathematics[END_REF]. Like many concepts in mathematics, developing a thorough understanding of mathematical functions is challenging for many students [START_REF] Clement | What do students really know about functions[END_REF]. Students have learning difficulties with mathematical functions include operations with functions, recognizing functions in different context. Tall and Viner (1981) identified some difficulties that students faced while using graphical methods on functions. [START_REF] Aspinwall | Uncontrollable mental imagery: Graphical connections between a function and its derivative[END_REF] asserted that in many cases the graphical representations can cause cognitive difficulties, because the perceptual analysis and synthesis of mathematical information presented implicitly in a diagram often makes greater demands on a student that any other aspect of a problem. For students to develop a holistic understanding of the concept of mathematical functions, they have to be able to identify the connection elements of a functional dependency and to combine these [START_REF] Nitsch | Students'competencies in working with functions in secondary mathematics educationempirical examination of a competence structure model[END_REF]. The ability to translate between various functional representations is essential for this purpose [START_REF] Ainsworth | Examining the effects of different multiple representational systems in learning primary mathematics[END_REF][START_REF] Seufert | Supporting coherence formation in learning from multiple representations[END_REF]. [START_REF] Smith | Learning angles through movement: Critical actions for developing understanding in an embodied activity[END_REF] found that students who focused on static representations of angles experienced less learning gains than those who participated in the movement-based lessons.

In order to provide different ways to perceive mathematical functions, we developed a game which is the subject of this research, using embodied motion as a fundamental strategy.

El-Demerdash, Lealdino, Mercat

1 -2
The contents of the kinesthetic digital materials include a wide range of mathematical functions aimed at freshmen students: Constant function, linear function, quadratic function, cubic function, absolute value function, fractional function, trigonometric functions, exponential function, logarithm function.

An achievement test is developed to assess the impact of this experiment in enhancing students' recognition and understanding of mathematical function graphs.

METHODOLOGY The Development of Kinesthetic Digital Resources

In order to get kinesthetic feedback from the students, we used the Microsoft Kinect sensor and some programming, using Unity 3D Development Kit, Microsoft Kinect (SDK) -Software Development Kit -and CindyScript the programming language associated with Cinderella -Dynamic Geometry Software, linked through a UDP (User Datagram Protocol) connection. The software recognizes students' body gestures as input representations for mathematical function graphs.

As shown in the figure below, Cinderella software plots two function graphs, the red one corresponds to body gestures, and the green one corresponds to the goal graph which the student should perform to obtain the positive score during the game.

The Development of Achievement Test

Since the purpose of the present study is to promote students' recognition and understanding of mathematical functions, we designed an achievement test (It is available as a download from: https://dl.dropboxusercontent.com/u/5761672/AchievementTest.pdf), that can be used in assessing students' performances before and after administering the prepared kinesthetics digital resources.

The following table shows the test specification in terms of mathematical function types that have been covered in the digital resources with respect to recognition and understanding levels of mathematical achievement. The test consists of 40 items that are equally distributed among the two levels of achievement: recognition and understanding. 20 items are designed at each level in isomorphic pairs that will be administered randomly as whether pretest or posttest.

The items type for the test are multi-choice questions of 4 choices to eliminate the guessing issues from the students' side. The test items were written in verbal and nonverbal, symbolic and graphical ways.

The total score of this test is set to be (40) marks. One-point mark is given to each correct item response.

Procedure of Experimental Study

We are going to set up a challenge where volunteer teams of freshmen students will have to compete by achieving highest scores and be rewarded by some goodies (T-shirts, caps…). The learning gain is going to be computed and correlated with the scores in performing the kinesthetic game in order to look for evidence of a relation between their engagement and achievements.

Figure 1 :

 1 Figure 1: Cinderella plots the function graph in accordance to body gesture detected by the Kinect sensor. On the left side of the figure above we show Cinderella plotting the functions in accordance with what the player (right side) is performing in front of Kinect sensor.

Table 1 :

 1 Specification of the Achievement Test

	Achievement Level	Recognition	Understanding
	Function Type	No. of Items	Items Numbers	No. of Items	Items Numbers	Sum	%
	Constant function	2	11, 12	2	23, 30	4	10%
	Linear function	2	9, 10	2	21, 24	4	10%
	Quadratic function	2	1, 15	2	27, 37	4	10%
	Cubic function	2	2, 16	2	28, 38	4	10%
	Square root function	2	3, 6	2	32, 36	4	10%
	Absolute value function	2	17, 18	2	29, 35	4	10%
	Fractional function	2	7, 8	2	22, 31	4	10%
	Trigonometric function	2	19, 20	2	39, 40	4	10%
	Exponential function	2	4, 13	2	25, 34	4	10%
	Logarithm function	2	5, 14	2	26, 33	4	10%
	Sum	20		20		40
	%	50%		50%			100%

ACKNOWLEDGEMENTS

This study was carried out in the framework of the MC Squared project (http://mc2-project.eu/) funded by the European Commission, under the 7th Framework Programme for Research, Technological Development and Demonstration (ICT-2013.8.1 "Technologies and scientific foundations in the field of creativity", Project 610467).

1 -4 CAPES -Coordenação de Aperfeiçoamento de


Sergey Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

Henri Gouin
email: henri.gouin@univ-amu.fr

Dynamic boundary conditions for membranes whose surface energy depends on the mean and Gaussian curvatures

Keywords: 68.03.Cd, 68.35.Gy, 02.30.Xx. Mathematics Subject Classification (2000) 74K15, 76Z99, 92C37

, the dependence on the mean and Gaussian curvatures was considered in statics and under a restrictive assumption of the membrane inextensibility. The authors derived the shape equation as well as two scalar boundary conditions on the contact line.

In this paper -thanks to the principle of virtual working -the equations of motion and boundary conditions governing the fluid membranes subject to general dynamical bending are derived without the membrane inextensibility assumption. We obtain the dynamic 'shape equation' (equation for the membrane surface) and the dynamic conditions on the contact line generalizing the classical Young-Dupré condition.

Introduction

The study of equilibrium, for small wetting droplets placed on a curved rigid surface, is an old problem of continuum mechanics. When the droplets' size is of micron range the droplet volume energy can be neglected. The surface energy of the surface S can be expressed in the form :

E = S σ ds,
where σ denotes the energy per unit surface. Two types of surfaces are present in physical problems:

rigid surfaces (only the kinematic boundary condition is imposed) free surfaces (both the kinematic and dynamic boundary conditions are imposed)

We will see the difference between the energy variation in the case of rigid and free surfaces.

The simplest case corresponds to a constant surface energy σ, but in general, σ also depends on physical parameters (temperature, surfactant concentrations, etc. [START_REF] Gouin | Interfaces endowed with non-constant surface energies revisited with the d'Alembert-Lagrange principle[END_REF][START_REF] Rocard | Thermodynamique[END_REF][START_REF] Steigmann | Energy minimizing states of capillary systems with bulk, surface and line phases[END_REF]) and geometrical parameters (invariants of curvature tensor). The last case is important in biology and, in particular, in the dynamics of vesicles [START_REF] Alberts | Molecular biology of the cell, 4th edn[END_REF][START_REF]Structure and dynamics of membranes[END_REF][START_REF] Seifert | Configurations of fluid membranes and vesicles[END_REF]. Vesicles are small liquid droplets with a diameter of a few tens of micrometers, bounded by an impermeable lipid membrane of a few nanometers thick. The membranes are homogeneous down to molecular dimensions. Consequently, it is possible to model the boundary of vesicle as a two-dimensional smooth surface whose energy per unit surface σ is a function both of the sum (denoted by H) and product (denoted by K) of principal curvatures of the curvature tensor :

σ = σ(H, K).
In mathematical description of biological membranes, one often uses the Helfrich energy [START_REF] Helfrich | Elastic properties of lipid bilayers: theory and possible experiments[END_REF][START_REF] Tu | Geometry of membranes[END_REF] :

σ(H, K) = σ 0 + κ 2 (H -H 0 ) 2 + κK, (1)
where σ 0 , H 0 , κ and κ are dimensional constants. Another purely mathematical example is the Wilmore energy [START_REF] Willmore | Riemannian Geometry[END_REF] :

σ(H, K) = H 2 -4K.
This energy measures the "roundness" of the free surface. For a given volume, this energy is minimal in case of spheres. One can also propose another surface energy in the form :

σ = σ 0 + h 0 (H 2 -H 2 0 ) 2 + k 0 (K -K 0 ) 2 ,
where σ 0 , h 0 , H 0 , k 0 and K 0 are dimensional constants. This kind of energy is invariant under the change of sign of principal curvatures, (i.e. the change of sign yields H → -H, K → K). It can thus describe the 'mirror buckling' phenomenon : a portion of the membrane inverts to form a cap with equal but opposite principal curvatures. It is also a homogeneous function of degree four with respect to principal curvatures. The equilibrium for membranes (called "shape equation" by Helfrich) is formulated in numerous papers and references herein [5-7, 14, 16, 17]. The "edge conditions" (boundary conditions at the contact line) are formulated in few papers and only in statics. In particular, in [START_REF] Rosso | Adhesive borders of lipid membranes[END_REF] the shape equation and two boundary conditions are formulated for the general dependence σ(H, K) under the assumption of the membrane inextensibility. However, the boundary conditions obtained do not contain the classical Young-Dupré condition for the constant surface energy. In the case when the energy depends only on H the generalization of Young-Dupré condition was obtained in [START_REF] Gouin | Vesicle Model with Bending Energy Revisited[END_REF].

The aim of our paper is to develop the theory of moving membranes which are in contact with a solid surface. The surface energy of the membrane will be a function both of H and K. We obtain a set of boundary conditions on the moving interfaces (membranes) as well as on the moving edges.

The motion of a continuous medium is represented by a diffeomorphism φ of a three-dimensional reference configuration D 0 into the physical space. In order to analytically describe the transformation, variables X = (X 1 , X 2 , X 3 ) T single out individual particles corresponding to material or Lagrangian coordinates, subscript " T " means the transposition. The transformation representing the motion of a continuous medium occupying the material volume D t is :

x = φ (t, X) or x i = φ i (t, X 1 , X 2 , X 3 ) , i = {1, 2, 3},
where t denotes the time and x = (x 1 , x 2 , x 3 ) T denote the Eulerian coordinates. At t fixed, the transformation possesses an inverse and has continuous derivatives up to the second order (the dependence of the surface energy on the curvature tensor will regularize the solutions, so the cusps and shocks do not appear).

At equilibrium, the unit normal vector to a static surface ϕ 0 (x) = 0 is the gradient of the so-called signed distance function defined as follows. Let

d(x) =    min|x -ξ|, if ϕ 0 > 0, 0, if ϕ 0 = 0, -min|x -ξ|, if ϕ 0 < 0, (2)
where the minimum is taken over points ξ at the surface, and | | denotes the Euclidien norm. The unit normal vector is :

n = ∇d(x).
In dynamical problems, the main difficulty in formulating boundary conditions comes from the fact that one cannot assume that for all time t the unit normal vector to the surface is the gradient of the signed distance function.

Indeed, if the material surface is moving, i.e. the surface position depends on time t, the surface points of the continuum medium are also moving and they will depend implicitly on x. Let ϕ(t, x) = 0 be the position of the material surface at time t. Its evolution is determined by the equation :

ϕ t + u T ∇ϕ = 0, ( 3
)
where u is the velocity of particles at the surface. Equation ( 3) is the classical kinematic condition for material moving interfaces. Let us derive the equation for the norm of ∇ϕ. Taking the gradient of Eq. ( 3) and multiplying by ∇ϕ, one obtains :

(|∇ϕ|) t + n T ∇ u T ∇ϕ = 0, (4)
where n = ∇ϕ |∇ϕ| is the unit normal vector to surface ϕ(t, x) = 0. It follows from Eq. ( 4) that, even if initially |∇ϕ| = 1 (i.e. unit normal n is defined at t = 0 as the gradient of the signed distance function), this property is not conserved in time.

The following definitions and notations are used in the paper. For any vectors a, b, we write a T b for their scalar product (the line vector is multiplied by the column vector), and a b T for their tensor product (the column vector is multiplied by the line vector). The last product is usually denoted as a ⊗ b. The product of a second order tensor A by a vector a is denoted by A a. Notation b T A means the covector c T defined by the rule c T = (A T b) T . The identity tensor is denoted by I. The divergence of A is covector divA such that, for any constant vector h, one has (divA) h = div (A h),

i.e. the divergence of A is a row vector, in which each component is the divergence of the corresponding column of A. It implies div (Av) = (divA) v + tr A ∂v ∂x , for any vector field v. Here tr is the trace operator. If f is a real scalar function of x, ∂f ∂x is the linear form (line vector) associated with the gradient of f

(column vector) : ∂f ∂x = (∇f ) T .
If n is the unit normal vector to a surface, P = Inn T is the projector on the surface with the classical properties : P 2 = P, P T = P, Pn = 0, n T P = 0.

For any scalar field f , the vector field v and second order tensor field A, the tangential surface gradient, tangential surface divergence, Beltrami-Laplace operator, and tangent tensors are defined as :

v tg = Pv, A tg = PA, ∇ tg f = P∇f, div tg v tg = tr P ∂v tg ∂x , ∆ tg f = div tg (∇ tg f ) ,
and for any constant vector h,

div tg (A tg h) = div tg (A tg ) h.
The following relations between surface operators and classical operators applied to tangential tensors in the sense of previous definitions are valid :

div tg v tg = divv tg + n T ∂n ∂x T v tg , (5)
div tg v tg = n T rot (n × v tg ) , (6)
div tg A tg = divA tg + n T ∂n ∂x T A tg , (7)
div tg (f v tg ) = f div tg v tg + (∇ tg f ) T v tg , (8)
div tg (f A tg ) = f div tg A tg + (∇ tg f ) T A tg , (9)
where rot denotes the curl operator. The proof is straightforward. Indeed, since

∂(n T v tg ) ∂x = n T ∂v tg ∂x + v T tg ∂n ∂x = 0, one has div tg v tg = tr P ∂v tg ∂x = divv tg -n T ∂v tg ∂x n = divv tg + n T ∂n ∂x T v tg ,
which proves relation [START_REF] Biscari | Impermeability effects in three-dimensional vesicles[END_REF]. To prove relation [START_REF] Capovilla | Stresses in lipid membranes[END_REF], one uses the following identity valid for any vector fields a and b :

rot (a × b) = a divb -b diva + ∂a ∂x b - ∂b ∂x a.
We apply this identity to the vectors a = n and b = v tg . Multiplying on left by n T , one obtains relation [START_REF] Capovilla | Stresses in lipid membranes[END_REF]. Relations [START_REF] Fournier | On the stress and torque tensors in fluid membranes[END_REF], [START_REF] Gavrilyuk | A new form of governing equations of fluids arising from Hamilton's principle[END_REF], [START_REF] Gavrilyuk | Multiphase Flow Modeling via Hamilton's principle[END_REF] are direct consequences of relation [START_REF] Biscari | Impermeability effects in three-dimensional vesicles[END_REF].

Curvature tensor

The unit normal vector being prolonged in the surface vicinity, we can directly obtain the expression of its derivative :

∂n ∂x = P ϕ |∇ϕ| ,
where ϕ is the Hessian matrix of ϕ with respect to x. One obviously has

n T ∂n ∂x = 0.
However, since in dynamics n is not the gradient of the signed distance function, we cannot have the property :

∂n ∂x n = 0. (10)
The curvature tensor is defined as :

R = -P ϕ |∇ϕ| P = - ∂n ∂x P.
Hence, in dynamics

R = - ∂n ∂x .
Let us note that the derivation of the shape equation and boundary conditions in statics always uses property [START_REF] Germain | The method of the virtual power in continuum mechanics. -Part 2: microstructure[END_REF] and the curvature tensor coming from the definition of the signed distance function. In dynamics, we cannot use these properties and new tools should be developed.

Tensor R is symmetric and has zero as an eigenvalue :

R = R T , Rn = 0.
In the eigenbasis, tensor R is diagonal :

R =   c 1 0 0 0 c 2 0 0 0 0   ,
where c 1 , c 2 are the principal curvatures. The two invariants of curvature tensor R are :

H = c 1 + c 2 , K = c 1 c 2 .
Invariant H is the double mean curvature, and invariant K is the Gaussian curvature. They can also be expressed in the form :

H = tr R = -tr ∂n ∂x , 2K = (tr R) 2 -tr R 2 = tr ∂n ∂x 2 -tr ∂n ∂x 2 .
Lemma 1 The following identities are valid :

div tg P = H n T , div tg R = ∇ tg T H + H 2 -2 K n T , R 2 = H R -K P.
Proof : First, let us remark that P = P tg , R = R tg . One can apply Eq. ( 7) to obtain :

div tg P = -div n n T + n T ∂n ∂x T P = -(div n) n T -n T ∂n ∂x T + n T ∂n ∂x T I -n n T = -(div n) n T ,
which proves the first relation. The proof of the second relation is as follows :

div R = -div ∂n ∂x + div ∂n ∂x n n T = - ∂(div n) ∂x + div ∂n ∂x n n T + n T ∂n ∂x 2 T = - ∂(div n) ∂x + div ∂n ∂x n n T + tr ∂n ∂x 2 n T + n T ∂n ∂x 2 T = ∂H ∂x P + tr ∂n ∂x 2 n T -n T ∂n ∂x T R.
Consequently,

div tg R = ∂H ∂x P + tr ∂n ∂x 2 n T .
Using tr ∂n ∂x

2 = tr R 2 = H 2 -2 K
, we obtain the second relation of the lemma. Now, the curvature tensor satisfies the Cayley-Hamilton theorem :

R 3 -H R 2 + K R = 0.
The minimal polynomial is :

R 2 -H R + K P = 0,
which proves the third relation.

Virtual motion

Let a one-parameter family of virtual motions

x = Φ (t, X, λ)
with scalar λ ∈ O, where O is an open real interval containing zero and such that Φ (t, X, 0) = φ (t, X) (the motion of the continuous medium is obtained for λ = 0). The virtual displacement of particle X is defined as [START_REF] Gavrilyuk | Multiphase Flow Modeling via Hamilton's principle[END_REF][START_REF] Serrin | Mathematical principles of classical fluid mechanics[END_REF] :

δx(t, X) = ∂Φ(t, X, λ) ∂λ | λ=0 .
In the following, symbol δ means the derivative with respect to λ at fixed Lagrangian coordinates X and t, for λ = 0. We will also denote by ζ(t, x) the virtual displacement expressed as a function of Eulerian coordinates :

ζ(t, x) = ζ (t, φ(t, X)) = δx (t, X) .

Variational tools

We assume that D t has a smooth boundary S t with edge C t . We respectively denote D 0 , S 0 and C 0 the images of D t , S t and C t in the reference space (of Lagrangian coordinates). The unit vector n and its image n 0 are the oriented normal vectors to S t and S 0 ; the vector t is the oriented unit tangent vector to C t and n = t × n is the unit binormal vector (see Fig. 1). F = ∂φ(t, X)/∂X ≡ ∂x/∂X is the deformation gradient. For the sake of simplicity, we will use the same notations for quantities as F, n, etc. both in Eulerian and Lagrangian coordinates.

Lemma 2 We have the relations :

δ det F = det F div ζ , (11)
δn = -P ∂ζ ∂x T n , (12)
δ

F -1 n = -F -1 ∂ζ ∂x n + F -1 δn , (13)
δ ∂n ∂x = ∂δn ∂x - ∂n ∂x ∂ζ ∂x . (14)
Proof of Rel. ( 11):

The Jacobi formula for determinant is :

δ(det F) = det F tr F -1 δF . Also, δF = δ ∂x ∂X = ∂δx ∂X . Then tr F -1 δF = tr ∂X ∂x ∂δx ∂X = tr ∂δx ∂X ∂X ∂x = tr ∂ζ ∂x = div ζ.
Proof of Rel. ( 12): Surface ϕ(t, x) = 0 is a material surface. It can be represented in the Lagrangian coordinates as ϕ(t, x) = ϕ 0 (X) which implies that δϕ = 0. Also,

δ ∂ϕ ∂x = δ ∂ϕ ∂X F -1 = ∂δϕ ∂x - ∂ϕ ∂x ∂ζ ∂x = - ∂ϕ ∂x ∂ζ ∂x .
Here we used the following expression for the variation of F -1 coming from the relation

F -1 F = I : δF -1 = -F -1 ∂ζ ∂x .
One also has :

δ|∇ϕ| = (∇ϕ) T δ∇ϕ |∇ϕ| .
Finally, taking the variation of n = ∇ϕ |∇ϕ| , one can obtain

δn = n T n -I ∂ζ ∂x T n = -P ∂ζ ∂x T n.
Proof of Rel. ( 13):

δ F -1 n = δ F -1 n + F -1 δn = -F -1 ∂ζ ∂x n + F -1 δn.
Proof of Rel. ( 14):

δ ∂n ∂x = δ ∂n ∂X F -1 = ∂δn ∂X F -1 + ∂n ∂X δF -1 = ∂δn ∂x - ∂n ∂x ∂ζ ∂x .
We denote by σ the energy per unit area of surface S t . The variation of σ is δσ. This variation depends on the physical problem through the dependence of σ on geometrical and thermodynamical parameters. For now, we do not need to know this variation in explicit form, the variation will be given further. The next lemma gives the variation of the surface potential energy [START_REF] Gouin | Interfaces endowed with non-constant surface energies revisited with the d'Alembert-Lagrange principle[END_REF][START_REF] Gouin | Vesicle Model with Bending Energy Revisited[END_REF].

Lemma 3 Let us consider a material surface S t of boundary edge C t . The variation of surface energy

E = St σ ds is δE = St δσ -∇ T tg σ + σH n T ζ ds + Ct σ n T ζ dl,
where ds, dl are the surface and line measures, respectively 1 .

1 It is interesting to remark that the combination δσ = δσ -∇ T tg σ ζ is the variation of σ at fixed Eulerian coordinates. Indeed, since the symbol δ means the variation at fixed Lagrangian coordinates, and δ is the variation at fixed Eulerian coordinates, this formula is a natural general relation between two types of variations (cf. [START_REF] Gavrilyuk | A new form of governing equations of fluids arising from Hamilton's principle[END_REF][START_REF] Gavrilyuk | Multiphase Flow Modeling via Hamilton's principle[END_REF]).

Proof : We suppose that the unit normal vector field is locally extended in the vicinity of S t . For any vector field w one has :

rot(n × w) = n div w -w div n + ∂n ∂x w - ∂w ∂x n.
From relation n T n = 1, we obtain n T ∂n ∂x = 0. Using the definition of H, (H = -div n), we deduce on S t :

n T rot(n × w) = div w + H n T w -n T ∂w ∂x n. (15)
The surface energy is given by :

E = St σ |d 1 x ∧ d 2 x|,
where

d i x = ∂x ∂s i ds i (i = 1, 2
) and s i are curvilinear coordinates on S t . This integral can also be written as :

E = St σ det (n, d 1 x, d 2 x) = S0 σ det(FF -1 n, Fd 10 X, Fd 20 X).
Here d i0 X = ∂X ∂s i0 ds i0 and s i0 are the corresponding curvilinear coordinates on S 0 . Finally,

E = S0 σ (detF) det(F -1 n, d 10 X, d 20 X) = S0 σ det((detF) F -1 n, d 10 X, d 20 X).
Let us remark that (detF) F -1 n is the image of n and is not the normal vector to S 0 because F is not an orthogonal transformation. One has :

δE = S0 δσ det F det (F -1 n, d 10 X, d 20 X) + S0 σ δ det F det (F -1 n, d 10 X, d 20 X) .
Using Lemma 2, one gets :

S0 σ δ det F det (F -1 n, d 10 X, d 20 X) = St σ div ζ det(n, d 1 x, d 2 x) + σ det (δn, d 1 x, d 2 x) -σ det ∂ζ ∂x n, d 1 x, d 2 x = St div(σ ζ) -(∇ T σ) ζ -σn T ∂ζ ∂x n ds. Relation (15) yields div (σ ζ) + σH n T ζ -n T ∂(σ ζ) ∂x n = n T rot (σ n × ζ). It implies S0 σ δ det F det (F -1 n, d 10 X, d 20 X) = St -σH n T + (∇ T σ) P ζ ds + St n T rot (σ n × ζ) ds. Since P ∇σ ≡ ∇ tg σ, one has St n T rot (σ n × ζ) ds = Ct det(t, σ n, ζ) dl = Ct σ n T ζ dl,
and we obtain Lemma 3.

Lemma 4 Let σ be a function of curvature tensor R, or equivalently, a function of H and K. Then,

∂σ ∂R = a I + b R with a = ∂σ ∂H + H ∂σ ∂K and b = - ∂σ ∂K , (16)
where for the sake of simplicity, we indifferently write σ(R) or σ(H, K). In particular, this implies :

n T ∂σ ∂R ∂n ∂x = 0. ( 17
)
Proof : Since H = tr R, 2K = (tr R) 2 -tr R 2 , and ∂tr R k ∂R = k R k-1 , one gets ∂σ ∂R = ∂σ ∂H + H ∂σ ∂K I - ∂σ ∂K R. Since R = - ∂n ∂x P and ∂σ ∂R = a I + bR, (18)
we obtain

n T ∂σ ∂R ∂n ∂x = a n T ∂n ∂x -b n T ∂n ∂x 2 = 0. 5 Variation of σ
This is a key part of the paper. The variation of the surface energy per unit area is obtained in the general case σ = σ(H, K). The membrane is determined by a surface S t having a closed contact line C t on a rigid surface S = S 1 ∪ S 2 (see Fig. 1). The dependence on other parameters such as concentrations of surfactants on the membranes can further be taken into account as in [START_REF] Gouin | Interfaces endowed with non-constant surface energies revisited with the d'Alembert-Lagrange principle[END_REF][START_REF] Steigmann | Energy minimizing states of capillary systems with bulk, surface and line phases[END_REF].

Lemma 5 The variation of surface energy σ(R) is given by the relation :

δσ = -div tg ∂σ ∂R R ζ + P ∂σ ∂R δn + div tg ∂σ ∂R R ζ + div tg P ∂σ ∂R δn. (19)
Proof : Using Lemma 2, we have :

δR = -δ ∂n ∂x P = - ∂δn ∂x - ∂n ∂x ∂ζ ∂x P + ∂n ∂x δ nn T .
By taking account of Eq. ( 12) and δ nn T = δn n T + n δn By using Eq. ( 6) and Lemma 3 associated with the Stokes formula, and property n T P = n T , we obtain :

δE = St div tg ∂σ ∂R R -div tg P ∂σ ∂R R + div tg P div T tg P ∂σ ∂R n T -σ H n T -∇ T tg σ ζ ds + Ct n T σI - ∂σ ∂R R -div T tg P ∂σ ∂R n T ζ - ∂σ ∂R δn dl.
From Lemma 2 we deduce :

-n T ∂σ ∂R δn = n T ∂σ ∂R P ∂ζ ∂x T n,
which proves Lemma 7.

Equations of motion and shape equation

The vesicle occupies domain D t with a free boundary S t which is the membrane surface, and S 1 which belongs to the rigid surface S = S 1 ∪ S 2 . S 1 denotes the footprint of D t on S, and C t is the closed edge (contact line) between S 1 and S 2 (see Fig. 1). We denote n 1 the external unit normal to S 1 along contact line C t . Then denoting t 1 = -t, one has :

n 1 = t 1 × n 1 = n 1 × t.
The surface energy of membrane S t is denoted σ. Solid surfaces S 1 and S 2 have constant surface energies denoted σ 1 and σ 2 . The geometrical notations are shown in Fig. 1.

A n 1 n n' n 1 ' t S 2 C t S t S 1
Fig. 1 A drop lies on solid surface S = S 1 ∪ S 2 ; St is a free surface; n 1 and n are the external unit normal vectors to S 1 and St, respectively. Contact line Ct separates S 1 and S 2 , t is the unit tangent vector to Ct on S. Vectors n 1 = n 1 × t and n = t × n are the binormals to Ct relatively to S and St at point A of Ct, respectively.

One can formulate the virtual work principle in the form [START_REF] Germain | The method of the virtual power in continuum mechanics. -Part 2: microstructure[END_REF][START_REF] Gouin | The d'Alembert-Lagrange principle for gradient theories and boundary conditions[END_REF]:

δA e + δA i -δE = 0,
where δA e is the virtual work of external forces, δA i is the virtual work of inertial forces, and δE is the variation of the total energy. The energy E is taken in the form :

E = Dt ρ ε dv + St σ ds + S1 σ 1 ds,
where specific internal energy ε is a function of density ρ. As we mentioned before, one can also include in this dependence several scalar quantities which are transported by the flow (specific entropy, mass fractions of surfactants, etc.). From Lemma 2, Eq. ( 11) and the mass conservation law :

ρ det F = ρ 0 (X),
we obtain the variation of the specific energy and density at fixed Lagrangian coordinates in the form :

δε = p ρ 2 δρ with δρ = -ρ divζ,
where p is the thermodynamical pressure. Consequently, the variation of the first term is [START_REF] Berdichevsky | Variational Principles of Continuum Mechanics: I. Fundamentals[END_REF][START_REF] Gavrilyuk | Multiphase Flow Modeling via Hamilton's principle[END_REF][START_REF] Serrin | Mathematical principles of classical fluid mechanics[END_REF]:

δ Dt ρ ε dv = δ D0 ρ 0 ε dv 0 = D0 ρ 0 δε dv 0 = Dt ρ δε dv = - Dt p div ζ dv.
The variation of the surface energy is given in Lemma 3. The third term is the surface energy of S 1 with energy σ 1 per unit surface. The virtual work of the external forces is given in the form :

δA e = Dt ρ f T ζdv + St T T ζds + Ct σ 2 n T 1 ζds,
where ρ f is the volume external force in D t , T is the external stress vector at the free surface S t , and σ 2 n 1 is the line tension vector exerted on C t . The last term on the right-hand side comes from Lemma 3 which can be also applied for rigid surfaces. Finally,

δA i = - Dt ρ a T ζdv
is the virtual work of inertial force, where a is the acceleration. The virtual work of forces δT applied to the material volume D t is defined as :

δT = Dt -ρ a T + ρ f T -∇ T p ζ dv + S1 (p + H 1 σ 1 ) n T 1 ζ ds + St -div tg ∂σ ∂R R + div tg P ∂σ ∂R R (20)
-

div tg P div T tg P ∂σ ∂R n T + (p + Hσ) n T + ∇ T tg σ + T T ζ ds - Ct (σ 1 -σ 2 ) n T 1 + σ n T -n T div T tg P ∂σ ∂R n T -n T ∂σ ∂R R ζ + n T ∂σ ∂R P ∂ζ ∂x T n dl.
As usually, H 1 and H are the sum of principle curvatures of surfaces S 1 and S t , respectively. Terms on D t , S 1 , S t are in separable form with respect to the field ζ. Expression (20) implies the equation of motion in D t and boundary conditions on surfaces S 1 , S t [START_REF] Schwartz | Théorie des distributions[END_REF]. Virtual displacement ζ must be compatible with conditions of the problem; for example, S 1 is an external surface to domain D t and consequently ζ must be tangent to S 1 . This notion is developed in [START_REF] Berdichevsky | Variational Principles of Continuum Mechanics: I. Fundamentals[END_REF]. They are presented below.

Equation of motion

We consider virtual displacements ζ which vanish on the boundary of D t . The fundamental lemma of virtual displacements yields :

ρ a + ∇p = ρ f , (21)
which is the classical Newton law in continuum mechanics.

Condition on surface S 1

Due to the fact that the surface S 1 is -a priori -given, the virtual displacements must be compatible with the geometry of S 1 . This means that the non-penetration condition (slip condition) is verified :

n T 1 ζ = 0. ( 22
)
Constraint ( 22) is equivalent to the introduction of a Lagrange multiplier P 1 into (20) where ζ is now a virtual displacement without constraint. The corresponding term on S 1 will be modified into

S1 (p + H 1 σ 1 -P 1 ) n T 1 ζ ds.
Since the variation of ζ on S 1 is independent, Eq. ( 20) implies :

P 1 = p + H 1 σ 1 . (23)
This is the classical Laplace condition allowing us to obtain the normal stress component P 1 n 1 exerted by surface S 1 .

Extended shape equation

Taking account of Eqs. ( 21) and ( 23), for all displacement ζ on moving membrane S t , one has from Eq. ( 20) :

St -div tg ∂σ ∂R R + div tg P ∂σ ∂R R -div tg P div T tg P ∂σ ∂R n T + (p + Hσ) n T + ∇ T tg σ + T T ζ ds = 0.
It implies :

p + Hσ -div tg P div T tg P ∂σ ∂R n +∇ tg σ -div T tg ∂σ ∂R R + R div T tg P ∂σ ∂R + T = 0. ( 24
)
Equation ( 24) is the most general form of the dynamical boundary condition on S t . Due to the fact that surface energy σ must be an isotropic function of curvature tensor R, i.e. a function of two invariants H and K, we obtain (for proof, see Appendix) that the following vector

∇ tg σ -div T tg ∂σ ∂R R + R div T tg P ∂σ ∂R
is normal to S t and consequently T writes in the form :

T = -P n.
Here scalar P has the dimension of pressure. One obtains from Eq. (44) (see Appendix) :

Hσ -∆ tg a -b ∆ tg H -∇ T tg b ∇ tg H -div tg (R ∇ tg b) + 2K -H 2 ∂σ ∂H -HK ∂σ ∂K = P -p. (25)
Relation ( 25) is the normal component of Eq. [START_REF] Tu | Geometry of membranes[END_REF].

It is important to underline that equation ( 24) is only expressed in the normal direction to S t . This is not the case when surface energy σ also depends on physico-chemical characteristics of S t , as temperature or surfactants. In this last case, Marangoni effects can appear producing additive tangential terms to S t . Using Lemma 1 (second equation) and expressions of scalars a and b given by Eq. ( 16), we get the extended shape equation:

H σ -K ∂σ ∂K + 2K -H 2 ∂σ ∂H -∆ tg ∂σ ∂H -H ∆ tg ∂σ ∂K -∇ T tg H ∇ tg ∂σ ∂K + div tg R ∇ tg ∂σ ∂K = P -p. (26)
Equation ( 26) was also derived in [START_REF] Rosso | Adhesive borders of lipid membranes[END_REF] under the hypothesis [START_REF] Germain | The method of the virtual power in continuum mechanics. -Part 2: microstructure[END_REF] and the assumption of inextensibility of the membrane. Our derivation does not use these hypotheses. For example, the inextensibility property is not natural even in the case of incompressible fluids (at fixed volume, the surface of a 3D body may vary).

Helfrich's shape equation

The Helfrich energy is given by Eq. ( 1). The shape equation ( 26) immediately writes in the form :

σ 0 H + κ 2 (H -H 0 ) [4K -H (H + H 0 )] -κ ∆ tg H = P -p , (27)
which is the classical form obtained by Helfrich 2 .

7 Extended Young-Dupré condition on contact line C t Let us denote by θ = n , n 1 = π + n, n 1 [mod 2π] the Young angle between S 1 and S t (see Fig. 2). Due to the fact that C t belongs to S 1 , the virtual displacement on C t is in the form :

ζ = αt + βn 1 , (28)
where α and β are two scalar fields defined on S 1 . Let us remark that condition (28) expresses the non-penetration condition ( 22) on S 1 . Moreover, since n, n 1 , n 1 belong to the normal plane to C t at A (see Fig. 2), one has :

n = n 1 sin θ -n 1 cos θ. (29)
2 Let us note that Helfrich considered the vesicle as an incompressible fluid. He also assumed that the membrane has a total constant area. Then, the virtual work can be expressed as

δT = D ρ f T ζ dv + S T T ζ ds -δ S σ ds + λ 0 δ S ds + δ D p div ζ dv,
where the scalar λ 0 is a constant Lagrange multiplier and p is a distributed Lagrange multiplier. The 'shape equation' is similar to (27).

A

n 1 n n' n 1 ' ¡ C t t
Tangent plane to S t Tangent plane to S 1

But relation ζ

T n 1 = 0 implies P ∂ζ ∂x T n 1 + P ∂n 1 ∂x T ζ = 0.
Replacing (29) into (20) one has :

δT = - Ct (σ 1 -σ 2 ) n T 1 + σ n T -n T div T tg P ∂σ ∂R n T -n T ∂σ ∂R R ζ + n T ∂σ ∂R P ∂ζ ∂x T n dl = - Ct (σ 1 -σ 2 ) n T 1 + σ n T -n T div T tg P ∂σ ∂R n T -n T ∂σ ∂R R + cos θ n T ∂σ ∂R P ∂n 1 ∂x T ζ + sin θ n T ∂σ ∂R P ∂ζ ∂x T n 1 dl = 0. ( 30
)
We choose now the virtual displacement in the form ζ = β n 1 . One has :

∂ζ ∂x = n 1 (∇β) T + β ∂n 1 ∂x , ∂ζ ∂x T = ∇β n T 1 + β ∂n 1 ∂x T .

Since ∂n 1 ∂x

T n 1 = 0, it implies :

∂ζ ∂x T n 1 = ∇β.
The integral (30) becomes :

Ct (σ 1 -σ 2 ) n T 1 + σ n T -n T div T tg P ∂σ ∂R n T -n T ∂σ ∂R R + cos θ n T ∂σ ∂R P ∂n 1 ∂x T n 1 β + sin θ n T ∂σ ∂R P ∇β dl = 0. ( 31
)
Since β and the components of ∇β can be choosen as independent, relation (31) implies two boundary conditions. The first condition on line C t is :

sin θ n T ∂σ ∂R P = 0. ( 32
)
The second condition is :

(σ 1 -σ 2 ) n T 1 + σ n T -n T div T tg P ∂σ ∂R n T -n T ∂σ ∂R R + cos θ n T ∂σ ∂R P ∂n 1 ∂x T n 1 = 0.
(33) The case sin θ = 0 all along C t is degenerate. If θ = 0, this corresponds to a hydrophobic surface (the contact line is absent). If θ = π, this corresponds to a complete wetting. In the last case n 1 = -n , n 1 = n, and the condition (33) becomes trivial : σ 1 -σ 2 -σ = 0.

The general case corresponds to the partial wetting (sin θ = 0). Due to Eq. ( 18),

n T ∂σ ∂R P ≡ n T (a I + b R)P ≡ a n T + b n T R ≡ n T ∂σ ∂R .
Hence, Eq. (32) yields

n T ∂σ ∂R = 0. ( 34
)
Equation ( 34) implies (see Lemma 4) :

n T ∂σ ∂H + H ∂σ ∂K I - ∂σ ∂K R = 0.
Consequently, n is an eigenvector of R. We denote c n the associated eigenvalue c 2 . Then

∂σ ∂H + H ∂σ ∂K = c n ∂σ ∂K . (35)
Due to the fact that t is also eigenvector of R with eigenvalue c t = c 1 (t and n form the eigenbasis of R along C t ), we get H = c t + c n and the equivalent to the boundary condition (35) in the form :

∂σ ∂H + c t ∂σ ∂K = 0. ( 36
)
From Lemma 4, Eq. ( 16), we immediately deduce :

div tg P ∂σ ∂R = ∇ T tg a + a H + b H 2 -2 b K n T + ∇ T tg b R + b ∇ T tg H. ( 37
)
Due to the fact that n T n = 0, we obtain :

n T div T tg P ∂σ ∂R = n T [∇ tg a + R ∇ tg b + b ∇ tg H] = n T [∇a + R ∇b + b ∇H]
Consequently, one obtains the second condition on C t in the form :

σ 1 -σ 2 + σ cos θ-sin θ n T (∇a + b ∇H + R ∇b) = 0. ( 38
)
This is the extended Young-Dupré condition along contact line C t between membrane S t and solid surface S (3 ).

∂ζ ∂x = t(∇α) T + α ∂t ∂x , ∂ζ ∂x T n 1 = α ∂t ∂x T n 1 . Since ∂t ∂x = c N t T ,
where N is the principal unit normal and c is the curvature along Ct, one obtains :

∂ζ ∂x T n 1 = α c t N T n 1 and sin θ n T ∂σ ∂R P ∂ζ ∂x T n 1 = α c sin θ n T ∂σ ∂R t N T n 1 ,
which is equal to zero thanks to Eq. (34). Moreover, thanks to Eq. (34), we immediately obtain that term

(σ 1 -σ 2 ) n T 1 + σ n T -n T div T tg P ∂σ ∂R n T -n T ∂σ ∂R R + cos θ n T ∂σ ∂R P ∂n 1 ∂x T t α
is vanishing. Hence, new boundary conditions do not appear on Ct.

In the case of Helfrich's energy given by relation [START_REF] Alberts | Molecular biology of the cell, 4th edn[END_REF], we obtain the extended Young-Dupré condition (38) in the form :

σ 1 -σ 2 + σ cos θ-κ sin θ n T ∇H = 0. (39)
This last condition was previously obtained in [START_REF] Gouin | Vesicle Model with Bending Energy Revisited[END_REF].

8 Surfaces of revolution

Shape equation for the surfaces of revolution

Along a revolution surface, the invariants of the curvature tensor depend only on s which is the curvilinear abscissa of meridian curve denoted by Γ [START_REF] Aleksandrov | Intrinsic Geometry of Surfaces[END_REF] :

H = H(s), K = K(s).
One of the eigenvectors, denoted e 1 , of the curvature tensor R is tangent to meridian curve Γ (see Fig. 3). Let us remark that for any function f (s), one has :

∇ tg f = df ds e 1 , ∆ tg f = d 2 f ds 2 .
Indeed, the first equation is the definition of the tangential gradient. The second equality is obtained as follows : The Frénet formula was used here :

de 1 ds = c 1 n. Also, div tg (R∇ tg f ) = div tg df ds Re 1 = div tg df ds c 1 e 1 = d ds c 1 df ds .
For surfaces of revolution the shape equation ( 26) becomes :

H σ -K ∂σ ∂K + 2K -H 2 ∂σ ∂H - d 2 ds 2 ∂σ ∂H -H d 2 ds 2 ∂σ ∂K - dH ds d ds ∂σ ∂K + d ds c 1 d ds ∂σ ∂K = P -p. A C t n 1 n'=e 1 t=e 2 n n 1 '
Tangent plane to S 1

Meridian plane

Fig. 3 The case of a revolution domain. The line Ct (contact edge between St and S 1 ) is a circle with an axis which is the revolution axis collinear to n 1 . The meridian curve is denoted Γ ; normal vector n and binormal vector n are in the meridian plane of revolution surface St. We have n = e 1 and t = e 2 , corresponding to the eigenvectors of the curvature tensor R at A.

Extended Young-Dupré condition for surfaces of revolution

One has along C t , t = e 2 , n = e 1 . It implies n T R t = 0. Also, one has :

n T (∇a + b ∇H + R ∇b) = da ds + b dH ds + c 1 db ds .
The Young -Dupré condition (38) becomes :

σ 1 -σ 2 cos θ-sin θ da ds + b dH ds + c t db ds = 0. Since a = ∂σ ∂H + H ∂σ ∂K , b = - ∂σ ∂K ,
one finally obtains :

σ 1 -σ 2 cos θ-sin θ d ds ∂σ ∂H + c n d ds ∂σ ∂K = 0.
For the Helfrich energy (1) this expression yields :

σ 1 -σ 2 cos θ-κ dH ds sin θ = 0.

Conclusion

Membranes can be considered as material surfaces endowed with a surface energy density depending on the invariants of the curvature tensor : σ = σ(H, K). By using the principle of virtual working, we derived the boundary conditions on the moving membranes ("shape equation") as well as two boundary conditions on the contact line. In limit cases, we recover classical boundary conditions. The "shape equation" and the boundary conditions are summarized below in the non-degenerate case (see ( 26), ( 36), ( 38)) as the equation for the moving surface S t :

• H σ -K ∂σ ∂K + 2K -H 2 ∂σ ∂H -∆ tg ∂σ ∂H -H ∆ tg ∂σ ∂K -∇ T tg H ∇ tg ∂σ ∂K + div tg R ∇ tg ∂σ ∂K = P -p.
the clamping condition on the moving line C t :

• ∂σ ∂H + c t ∂σ ∂K = 0,
Also, (t, n, n ) -which is the Darboux frame -are the eigenvectors of curvature tensor R.

dynamic generalization of the Young-Dupré condition on C t :

• σ 1 -σ 2 + σ cos θ-sin θ n T ∇ tg ∂σ ∂H + (HP -R) ∇ tg ∂σ ∂K = 0.

In the case of Helfrich's energy the generalization of Young-Dupré condition is reduced to equation (39):

σ 1 -σ 2 + σ cos θ-κ sin θ n T ∇ tg H = 0.
The last term, corresponding to the variation of the mean curvature of S t in the binormal direction at the contact line, can dominate the other terms. It could be interpreted as a line tension term usually added in the models with constant surface energy (cf. [START_REF] Babak | Stability of the lenslike liquid thickening (the drop) on a solid substrate[END_REF]). It should also be noted that the droplet volume has no effect in the classical Young-Dupré condition. This is not the case for the generalized Young-Dupré condition since the curvatures can become very large for very small droplets (they are inversely proportional to the droplet size). The clamping condition for the Helfrich energy fixes the value of H on the contact line :

H = H 0 -c t κ κ .
The new shape equation and boundary conditions can be used for solving dynamic problems. This could be, for example, the study of the "fingering" phenomenon appearing as a result of the non-linear instability of a moving contact line. This complicated problem will be studied in the future.

Appendix

Since σ = σ(H, K), we get :

∇ tg σ = ∂σ ∂H ∇ tg H + ∂σ ∂K ∇ tg K. ( 40
)
From Eq. ( 16), we obtain :

div tg P ∂σ ∂R = ∇ T tg a + a H + b H 2 -2 b K n T + ∇ T tg b R + b ∇ T tg H. (41)
Also, one has : From relations (40), ( 41), (42), we deduce :

∇ tg σ -div T tg ∂σ ∂R R + R div T tg P ∂σ ∂R = 2 a K + 3 b H K -a H 2 -b H 3 n + ∂σ ∂H ∇ tg H + ∂σ ∂K ∇ tg K -R ∇ tg (a + b H) -(a + b H) ∇ tg H + ∇ tg (b K) +R ∇ tg a + a H + b H 2 -2 b K R n + R 2 ∇ tg b + b R ∇ tg H + T = 0.
Using relations R n = 0, Eq. (1) 3 and expressions of a and b given by Eq. ( 16), we obtain :

∂σ ∂H ∇ tg H + ∂σ ∂K ∇ tg K -R ∇ tg (a + b H) -(a + b H) ∇ tg H + ∇ tg (b K) +R ∇ tg a + a H + b H 2 -2 b K R n + R 2 ∇ tg b + b R ∇ tg H = 0.
where all tangential terms disappear in the boundary condition on S t .

Fig. 2

 2 Fig. 2 Tangent planes to membrane St and solid surface S 1 : n 1 and n are the unit normal vectors to S and St, external to the domain of the vesicle; contact line Ct is shared between S and St and t is the unit tangent vector to Ct relatively to n; n 1 = n 1 × t and n = t × n are binormals to Ct relatively to S and St at point A, respectively. Angle θ = n , n 1 . The normal plane to Ct at A contains vectors n, n , n 1 , n 1 .

e 1 ⊗ e 1 = tr d 2 f ds 2

 12 Pe 1 ⊗ e 1 + c 1 (s) df ds n ⊗ e 1 = d 2 f ds 2 .

2 .

 2 div tg ∂σ ∂R R = div tg (a R) + div tg b RDue to[START_REF] Gavrilyuk | Multiphase Flow Modeling via Hamilton's principle[END_REF], one has :div tg (a R) = (∇ T tg a) R + a ∇ T tg H + a H 2 -2 K n T , div tg b R 2 = div tg [b (HR -KP)] = ∇ T tg (bH) R + b H ∇ T tg H + H 2 -2 K n T -∇ T tg (bK) -b KHn T . Consequently, div tg ∂σ ∂R R = ∇ T tg (a + b H) R (42) + (a + b H) ∇ T tg H -∇ T tg (bK) + (aH 2 + bH 3 -2aK -3bHK)n T .From relations (40), (41), (42), we deduce :∇ tg σ -div T tg ∂σ ∂R R +R div T tg P ∂σ ∂R = 2 aK + 3 bHK -a H 2 -b H 3 n.Using (41), one obtains :P div T tg P ∂σ ∂R = ∇ tg a + R ∇ tg b + b ∇ tg H.One deduces : div tg P div T tg P ∂σ ∂R = ∆ tg a +div tg (R ∇ tg b)+ b ∆ tg H +∇ T tg b ∇ tg H. (43)

= 2

 2 aK + 3 bHK -a H 2 -b H 3 n.Finally, using (43), one obtains :p + Hσ -∆ tg a -b ∆ tg H -∇ T tg b ∇ tg H -div tg (R ∇ tg b) + 2 a K + 3 b H K -a H 2 -b H 3 n + T = 0,

The virtual displacement taken in the most general form (28) does not produce new boundary conditions. Due to the linearity of the virtual work, to prove this property it is sufficient to take ζ = α t. We obtain

Acknowledgments. The authors thank the anonymous referees for their noteworthy and helpful remarks that greatly contributed to improve the final version of the paper.


Lucas Bourreau

Laurent Gaillet
email: laurent.gaillet@ifsttar.fr

Véronique Bouteiller

Franck Schoefs

Benoit Thauvin

Julien Schneider

Better understanding of tide's influence on half-cell potential and electrical resistivity measurements for reinforced concrete in marine environment

Keywords: reinforced concrete, bridge, marine environment, corrosion diagnosis, concrete durability, half-cell potential

Main damage occurring on reinforced concrete (RC) structures can be attributed to the corrosion of the rebar. Chloride ions penetration and/or carbonation, as major causes of RC corrosion, are well understood for onshore structures. When dealing with structures in marine environment, more parameters must be considered, as the different exposure zones (tidal, splash and atmospheric) will have a significant influence on the corrosion behavior. This article deals with the corrosion diagnosis of the piers of the ile de Ré bridge. The objective is to provide an effective method to identify the different exposure zones based on the non-destructive measurements (half-cell potential) that are modelled taking into account the influence of the tide.

Introduction

Corrosion is the main cause of degradation of reinforced concrete structures [START_REF] Broomfield | Corrosion of steel in concrete -Understanding, investigation and repair[END_REF][START_REF] Cigna | COST 521 -Corrosion of steel in reinforced concrete structures -Final Report: European Communities[END_REF][START_REF] Poupard | Corrosion diagnosis of reinforced concrete beams after 40 years exposure in marine environment by non destructive tools[END_REF]. This later is relatively well known for atmospheric structures but many other parameters are involved in marine conditions, including the different exposure zones [START_REF] Costa | Chloride penetration into concrete in marine environment -Part I: Main parameters affecting chloride penetration[END_REF][START_REF] Costa | Chloride penetration into concrete in marine environment -Part II: Prediction of long term chloride penetration[END_REF][START_REF] Montemor | Chloride-induced corrosion on reinforcing steel: from the fundamentals to the monitoring techniques[END_REF]. Non-destructive techniques for corrosion evaluation of rebars are more difficult to apply in this case in order to establish a corrosion diagnosis [START_REF] Polder | Test methods for on site measurement of resistivity of concrete -a RILEM TC-154 technical recommendation[END_REF][START_REF] Andrade | RILEM TC 154-EMC:Electrochemical Techniques for Measuring Metallic Corrosion -Recommendations -Test methods for on-site corrosion rate measurement of steel reinforcement in concrete by means of the polarization resistance method[END_REF][START_REF] Breysse | Accounting for variability and uncertainties in NDT condition assessment of corroded RC-structures[END_REF].

In this article are presented results from half-cell potential measurements obtained on the piers of the Ile de Ré Bridge. An innovative data processing method is used to identify the different exposure zone heights taking into account the influence of the tide on the results.

Material and methods

The corrosion measurements were performed on a set of piers of the ile de Ré bridge on the Atlantic coast. Two sides of the piers (Side C and side G) were considered corresponding to different environmental conditions. Assessment zone included 3 different exposure zones (tidal, splash and atmospheric) with 1m height each and 1.85 m width. A complete corrosion diagnosis methodology was established consisting in different stages [START_REF] Bourreau | Diagnostic de corrosion sur ouvrage -Fiabilité et aide à la décision[END_REF]. Corrosion diagnosis was performed using a vessel. Half-cell potential measurements were performed using a Canin+ device (: copper/saturated copper sulphate electrode)

Results

Corrosion diagnosis

Performing a corrosion diagnosis is not obvious on reinforced structures in a marine environment and the interpretation of the results can be tricky. Fig 1 represents the half-cell potential mapping of the investigated area analysed with the ASTM C876-91 standard [11, all the values done on the pier's zones with heights lower than a given value would refer to a high probability of corrosion (Fig. 1) as potential values are less than-350 mV/CSE and therefore can be attributed to a probability of corrosion superior to 90% (Fig. 1). For the same experimental data, the interpretation based on COST-509 or RILEM TC 154 documents also lead to erroneous conclusions that were not supported by the visual inspection observations [7] [12]. This half-cell potential map (fig 1) can be explained by the fact that high negative values correspond to rebar situated in the tidal zone where the concrete is saturated of sea water and there is a lack of oxygen [START_REF] Astm C | Standard test method for halfcell[END_REF][START_REF]Standard test method for halfcell potentials of uncoated reinforcing steel in concrete[END_REF][START_REF] Gonzalez | Considerations on reproducibility of potential and corrosion rate measurements in reinforced concrete[END_REF][START_REF] Hussain | Underwater half-cell corrosion potential bench mark measurement of corroding steel in concrete influenced by a variety of material science and environmental engineering variables[END_REF]. Above the tidal zone, concrete is less saturated and more oxygenated and therefore potential values are in a medium range.

These very negative potentials can be attributed to an active corrosion but also to the marine environment for structures in submerged or in tidal zones. So, a method is proposed to withdraw the tide's effect on the measured potentials.

Enhanced corrosion diagnosis

To study and quantify the tide's effect on the potential measurements, mean trends of a dataset of potential values have been computed. This 40th IABSE Symposium, 19-21 September 2018, Nantes, France. Tomorrow's Megastructures 3 selected dataset corresponds to pier sides where no corrosion has been detected.

Tide's effect on potential values

The mean value (μEcorr) and the standard deviation (σEcorr) have been calculated for each horizontal measurement lines. In order to delete the outlier values from the main trend, the absolute values which were not included in the range defined by μ +/-3*σ were deleted of the dataset. Figure 2 presents the mean potential trend for the selected pier sides versus water height. Different evolution's behaviour (dashed coloured lines) can be observed on Figure 2 that can be related to the exposure zones.

Figure 2: Mean potential trend versus water height (based on four selected pier sides of the bridge).

A first statement is that a change of behaviour in the mean trend represents the transition from an exposure zone to another, which corresponds to a boundary between two exposure zones. A data classification method was used to identify the transitions based on an analysis that focused on the slope's values computed from two adjacent potential mean values. The dataset consisted in the two sides of the piers as exposure zones locations are assumed not to be dependent on surface concrete conditions.

After dataset has been normalized, K-means classification algorithm has been applied to determine a modification of the value's evolution into mean trends. The optimal number of cluster was determined using the Caliński-Harabasz criterion [16]. The number of cluster (k) is optimal when the Caliński-Harabasz criterion is maximal. This analysis has been reiterated for all seasons. An overall optimal number of clusters of 3 is determined and will be considered in the following.

Boundaries between two exposure zones are located between two slopes values. Using the classification results, an identification and location of these boundaries is achievable. Boundaries between the high and median zones are located at water height 6.12 -5.64 m CD and 5.02 -4.62 m CD for potential values.

These calculated values for exposure zones and their boundaries have been compared to "physical" values for steel structures in marine conditions according to literature (Tab. As a result of this comparison; the low, median and high zones are attributed to the tidal, splash and atmospheric exposure zones respectively.

Tide's effect quantification

As an assumption, concrete's exposition to the wetting/drying cycles is considered as the origin for tide's effect on half-cell potential data. From this point of view, the two sides of the pier have different behaviours and must be studied separately. A side C (most exposed to wetting/drying cycles) and G (less exposed to wetting/drying cycles) were labelled.

The same methodology as previously described was used to determine the mean trend of the potential values. Results are presented in Fig. 3 where the boundaries between the exposure zones are also indicated.

The modelling method was based on a fitting of every exposure zones by a linear regression in order to determine potential's evolution profiles as illustrated in Fig.

4 From Figure 4, The potential values are more negative for side G than side C and this can be explained because side C is more exposed to wetting/drying cycles This methodology of correction of initial potential data to withdraw the tide's effect on the corrosion potential has been tested on potential dataset previously analysed in paragraph 2 (Fig. 2).

Reference corrosion potential calculated from the developed methodology are subtracted for every experimental potential value. The equicontour potential line plots from initial data and after using the developped methodology of correction are presented on Figure 5. Interpretation using the ASTM C876-91 standard is now different. In the earlier stage, a high probability of corrosion for lower heights and an 40th IABSE Symposium, 19-21 September 2018, Nantes, France. Tomorrow's Megastructures undefined risk of corrosion for the others zones were identified. After correction, the whole investigated area is associated to a low probability of corrosion.

Conclusions

Corrosion diagnosis of a reinforced concrete structure in the sea can result difficult because of the marine environment. In order to achieve a more reliable diagnosis based on half-cell potential measurements, a model was developed. First of all, the boundaries between tidal, splash and atmospheric zones were characterized and compared to the ones of metallic structures. Then, the influence of the tide was analysed and it was found that the evolution of the potential profiles versus water heights depended on the sides C or G of the piers and therefore of the exposure (more or less windy and humid). Finally, the develop methodology was efficient in obtaining corrected potential profiles (by substracting the tide effect) that were consistent with the standards or recommendation.

Figure 1 :

 1 Figure 1: Half-cell potential map of a passive pier's side, using the ASTM recommendation.

Figure 3 :

 3 Figure 3: Mean potential trend versus water height for the 2 different sides of the piers of the bridge Parameters of the linear regression have been determined and fits are significant with quite good R² coefficients (0.66-1) for the potential data.

Figure 4 :

 4 Figure 4: Evolution of potentialprofiles versus water height

Figure 5 :

 5 Figure 5: Half-cell potential maps of a passive pier's side, equicontour potential line plot using ASTM's recommendation with initial data (a) and corrected data (b).

Acknowledgements

The authors would like to thank the team of the DéCoF-Ré Project: B. Godart, W. Traverst, P. Boujard, N. Coulaty-Chin, A. Orcesi, Y. El Rabbih, X. Derobert, M. Sissoko (IFSTTAR); R. Queguiner, S. Pasquiet, P. Boulaire, M. Rebours, C. Naudat (Cerema Ouest); F. Landrin and V. Queyrat (Cerema Ile de France), J-F. Barthélémy (Cerema Dtech ITM); M. Roche, O. Amiri (Nantes University); M. Brouxel, S. Naar (SIXENSE-Concrete), A. Audouin-Dubreuil, F. Lavoute, M. Barbier (Département de la Charente-Maritime); F. Gazet (Dekra). They also thank the marine team for their help with the vessel.


Nicolas Girard

Guillaume Charpiat

Yuliya Tarabalka

NOISY SUPERVISION FOR CORRECTING MISALIGNED CADASTER MAPS WITHOUT PERFECT GROUND TRUTH DATA

Keywords: Noisy supervision, Multi-modal alignment, Ground truth annotation correction, Optical images

In machine learning the best performance on a certain task is achieved by fully supervised methods when perfect ground truth labels are available. However, labels are often noisy, especially in remote sensing where manually curated public datasets are rare. We study the multi-modal cadaster map alignment problem for which available annotations are misaligned polygons, resulting in noisy supervision. We subsequently set up a multiple-rounds training scheme which corrects the ground truth annotations at each round to better train the model at the next round. We show that it is possible to reduce the noise of the dataset by iteratively training a better alignment model to correct the annotation alignment.

INTRODUCTION

One of the main tasks in remote sensing is semantic segmentation. The supervised approach needs good ground truth annotations most often in the form of class-labeled polygons outlining objects of the image. However, these good annotations are hard to come by because even if they exist (for example OpenStreetMap (OSM) annotations [START_REF]Planet dump retrieved from[END_REF]), they can be misaligned due to human error, imprecise digital terrain model or simply a lack of precision of the original data (scanned cadaster maps from local authorities). Each object can be misaligned in a different way compared to surrounding objects and the misalignment can include complex deformations such as slight stretching and rotation.

The aim of this paper is to correct the alignment of noisy OSM annotations when only these annotations are available. Several related works tackle the noisy label problem. Some use special losses to explicitly model the label noise [START_REF] Mnih | Learning to label aerial images from noisy data[END_REF] which penalize erroneous outputs less if they could be due to label noise. Others perform simultaneous edge detection and alignment [START_REF] Yu | Simultaneous edge alignment and learning[END_REF] which can handle small displacements in an unsu-The authors would like to thank ANR for funding the study. pervised manner. The task of aligning OSM annotations has already been tackled in [START_REF] Vargas-Muoz | Correcting misaligned rural building annotations in open street map using convolutional neural networks evidence[END_REF], using convolutional neural networks for building segmentation and a Markov Random Field for aligning buildings onto the building segmentation image. However, the neural network has to be trained on a small dataset of building image with corresponding good groundtruth annotations.

We propose in this paper to use the self-supervised multi-task multi-resolution deep learning method for aligning cadaster maps to images of [START_REF] Girard | Aligning and updating cadaster maps with aerial images by multi-task, multi-resolution deep learning[END_REF] in the noisy supervision setting. The dataset used for training that method has misalignment noise and still the model learned to align. We will explore this interesting behavior and experiment a kind of unsupervised learning to correct noisy annotations with a model trained on these noisy annotations. See Fig. 1 for an example of results. It leverages the natural tendency of neural networks to be robust to a certain amount of noise and does not require any special loss function.

METHODOLOGY

We provide here a short description of the self-supervised multi-task multi-resolution deep learning alignment method Fig. 2: One step of the base alignment method, applied repeatedly at increasing resolutions for the final alignment.

of [START_REF] Girard | Aligning and updating cadaster maps with aerial images by multi-task, multi-resolution deep learning[END_REF] (referred to by "base alignment method" from now on), focusing on the most relevant parts. Its code is available here: https://github.com/Lydorn/mapalignment. Mathematical modeling. Given two images I and J of same size H × W , but of different modalities, e.g. with I an RGB image (picture from a satellite) and J a binary image (cadaster, indicating for each pixel whether it belongs to a building or not), the alignment problem aims at finding a deformation, i.e. a 2D vector field f defined on the discrete image domain [1, H] × [1, W ], such that the warped second image J • (Id + f ) is well registered with the first image I. To do this, in a machine learning setting, we consider triplets (I, J, f gt ) consisting of two images together with the associated ground-truth deformation f gt . Image pairs (I, J) are given as inputs, and the model's estimated deformation f is optimized to be close to the ground truth deformation f gt .

Displacement map cost function. The displacement map loss function is the mean squared error between the predicted displacement map f and the ground truth displacement map f gt . The actual loss used by the base alignment method is a little more complex but for the purpose of this paper we can consider the simplified loss:

L disp ( f ) = x∈[1,H]×[1,W ] f (x) -f gt (x) 2 2 (1)
Model. The neural network used by the base alignment method is a transformed U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] with 2 image inputs and 2 image-like outputs for the displacement map and the segmentation image, see Fig. 2 for a schema of the model. The segmentation output is only used during training, having its own cross-entropy loss function. The input image I has 3 channels, with real values normalized to [-1, 1], standing for RGB. The input misaligned polygon raster J also has 3 channels, with Boolean values in {0, 1}, corresponding to polygon interior, edge and vertices. The output displacement map has 2 channels with real values in [-4 px, 4 px], standing for the x and y components of the displacement vector. The model uses a multi-resolution approach by applying a neural network at increasing resolutions, iteratively aligning polygons from a coarse to fine scale. The scales used are 1 8 , 1 4 , 1 2 and 1. Thus displacements of up to 32 px can be handled.

EXPERIMENTAL SETUP

Dataset. The model is trained on a building cadaster dataset consisting of the two available following ones: Inria Aerial Image Labeling Dataset [START_REF] Maggiori | Can semantic labeling methods generalize to any city? the Inria aerial image labeling benchmark[END_REF], and "Aerial imagery object identification dataset for building and road detection, and building height estimation" [START_REF] Bradbury | Aerial imagery object identification dataset for building and road detection, and building height estimation[END_REF]. The resulting dataset has 386 images (with the majority being 5000 × 5000 px) of 16 cities from Europe and USA. Each image has in average a few thousand buildings. The building footprints were pulled from OSM for all images. these polygon annotations are inconsistent across images alignment-wise (see Fig. 1 and3). Some are perfect, and some are misaligned by up to 30 px. However the base alignment method [START_REF] Girard | Aligning and updating cadaster maps with aerial images by multi-task, multi-resolution deep learning[END_REF] assumes perfect annotations in its formulation.

Self-supervised training. The model needs varied ground truth labels (displacements) in order to learn, while the dataset is assumed to be made of aligned image pairs only (f = 0). The dataset is thus enhanced by adding random deformations in the form of 2D Gaussian random fields for each coordinate with a maximum absolute displacement of 32 px. The polygon annotations A are then inversely displaced by the generated displacements, to compute the misaligned polygons, which are then rasterized. We obtain training triplets of the form (I, J, f ) with J = rast(A • (Id + f ) -1 ). For the multiresolution pipeline, 4 different models are trained independently with downscaling factors 8, 4, 2 and 1 (one per resolution). For clarification, labels are the ground truth displacements at each pixel, i.e. a 2D vector, and annotations are the polygons outlining objects.

Multiple-rounds training. We train the base alignment method with the same hyper-parameters as [START_REF] Girard | Aligning and updating cadaster maps with aerial images by multi-task, multi-resolution deep learning[END_REF], which were selected to avoid overfitting. We perform multiple rounds of training on the whole dataset to achieve our goal of aligning the whole dataset. It consist of iteratively alternating between training the alignment model on the available dataset (see Alg. Robustness to noise. In an additional experiment (Noisier) we misaligned all original annotations further with random zero-mean displacements up to 16 px. We then applied our alignment method for correcting these noisier annotations to study its robustness to more noise.

RESULTS

As annotations of our dataset are noisy they cannot be used as ground truth to measure quantitative results. We can first visualize qualitative results in Fig. 3. In order to measure the effectiveness of the multiple rounds training, we manually aligned annotations for one 5000×5000 px image (771 buildings) to get a good ground-truth. We chose the blooming-ton22 image because it has severe misalignment. To measure the accuracy of an alignment, for any threshold τ we compute the fraction of vertices whose ground truth point distance is less than τ . In other words, we compute the Euclidean distance in pixels between ground truth vertices and aligned vertices, and plot the cumulative distribution of these distances in Fig. 4 (higher is better) for all experiments and rounds.

DISCUSSION AND CONCLUSION

After the first round of training, the annotations are on average better aligned than the original annotations, but in some cases the polygons are pushed into the wrong direction, resulting in poorer accuracy for some threshold levels (see the blue curve sometimes under the red curve in Fig. 4). However after the second round of training, the annotation alignment has been significantly improved upon the first round (error divided by more than 3 for any quantile, cf. green curve compared to blue curve). The 3rd round does not bring any significant improvement in this case.

Note that a perfect alignment score cannot be expected, because of the ambiguity of the "perfect" ground truth. Indeed, when manually aligning bloomington22's annotations, we observed that the majority of buildings are annotated by a coarse polygon that does not outline the building precisely. Best aligning such a coarse polygon to a real, more complex building becomes an ill-posed problem, with multiple equally-good solutions, which creates ground truth ambiguity. See Fig. 5 (left) for an illustration of this problem, especially the building on the top-right. Fig. 5 (right) shows an example of a mistake of our approach. The left building was successfully aligned (through a slight vertical and horizontal squashing), but the adjoining building on the right was not, because the model only learned smooth displacement maps. A more well-designed displacement map generation allowing discontinuities could solve such problems.

The first ablation study shows the importance of aligning the original annotations in the second step of Alg. 2 as it achieves better accuracy. Indeed the aligned annotations after round 1 can be worse than the original annotations (see some blue polygons of Fig. 3 and the blue curve of Fig. 4), and consequently more difficult to align. The second ablation study shows that the re-training step in round 2 is very important, as skipping it does not improve the alignment compared to round 1.

An explanation of how this method is able to align misaligned annotations by training on these misaligned annotations could be that the dataset contains enough perfect ground truth annotations to steer the gradient descent in the right direction, while being mildly affected by noisy labels (even if the noise is not zero-mean) if overfitting is avoided. However the last experiment (Noisier) invalidates this explanation because in that case the fraction of well-aligned ground truth is negligible and still the model was able to align noisier annota-tions virtually as well as it did original annotations (it however needs a 3rd round to do so). Our current tentative explanation is that ground truth labels have a zero-mean noise (without bias). For the alignment task, the network tries to minimize the average error it makes. As such it tends to predict the mean value of the labels when it cannot do better. This is the case if the label noise is independent of the input, and if overfitting noisy labels is avoided. The network will learn the mean alignment, which corresponds to the underlying perfect ground truth. This explanation is further supported by a recent work on image restoration without clean data [START_REF] Lehtinen | Noise2noise: Learning image restoration without clean data[END_REF], where noisy images are de-noised by training a de-noising network on noisy images only.

In conclusion, even noisy/misaligned annotations are useful. Our model can be iteratively trained on them and align these annotations through a multiple-round training scheme.

Fig. 1 :

 1 Fig. 1: Qualitative alignment results on a crop of an image of Bloomington from the Inria dataset. Red: initial OSM annotations; green: aligned annotations.

2 . 1 :

 21 1) and correcting the alignment of the training dataset to provide a better ground truth for the next training round. The multiple-rounds training is explained in Alg. Algorithm Alignment training [5] Input: Images I = {I, ...} and corresponding annotations A = {A, ...} Build dataset with random deformations: D = {(I, J rand , f rand ), ...} with J rand = rast(A • (Id + f rand ) -1 ); Train multi-resolution model M to perform this mapping: (I, J rand ) → f rand ; Output: Trained model M Algorithm 2: Multiple-rounds training Input: Original annotations A 0 , number of rounds R for r = 1 to R do 1. Get model M using Alg. 1 with input A = A r-1 ; 2. Apply M r on the original annotations A 0 : A r = M r (A 0 ); Output: Aligned annotations A R Ablation studies. To justify the design choices of the multiple-rounds training, we performed ablation studies. The first ablation study (AS1) changes the second step of Alg.2 by applying the model on the previous corrected annotations instead of the original annotations: A r = M r (A r-1 ) in order to test whether it is better to iteratively align annotations. The second ablation study (AS2) trains the model only once on the original annotations, and applies it R times to iteratively align the annotations (as in AS1). This is implemented by additionally replacing step 1 of Alg. 2 by M r = M 1 for r > 1 and leaving it as is for r = 1. This is meant to test the usefulness of re-training.

Fig. 3 :

 3 Fig. 3: Qualitative alignment results on a crop of blooming-ton22 from the Inria dataset. Red: initial dataset annotations; blue: aligned annotations round 1; green: aligned annotations round 2.

Fig. 4 :

 4 Fig. 4: Accuracy cumulative distributions measured with the manually-aligned annotations of bloomington22 from the Inria dataset.

Fig. 5 :

 5 Fig. 5: Left: ambiguity of the perfect ground truth annotations. Right: alignment failure case. Magenta: manually aligned annotations; red: original dataset annotations; green: aligned annotations round 2.