M Tufis

J.-G Ganascia

A Normative Extension for the BDI Agent Model

Keywords: BDI agent, normative agent, consistency check, consequentialism

This paper proposes an approach on the design of a normative rational agent based on the Belief-Desire-Intention model. Starting from the famous BDI model, an extension of the BDI execution loop will be presented; this will address such issues as norm instantiation and norm internalization, with a particular emphasis on the problem of norm consistency. A proposal for the resolution of conflicts between newly occurring norms, on one side, and already existing norms or mental states, on the other side, will be described. While it is fairly difficult to imagine an evaluation for the proposed architecture, a challenging scenario inspired form the science-fiction literature will be used to give the reader an intuition of how the proposed approach will deal with situations of normative conflicts.

INTRODUCTION

"Mistress, your baby is doing poorly. He needs your attention." "Stop bothering me, you f* robot." "Mistress, the baby won't eat. If he doesn't get some human love, the Internet pediatrics book says he will die" "Love the f*ing baby, yourself."

The excerpt is from Prof. John McCarthy's short story "The Robot and the Baby", 1 which besides being a challenging and insightful look into how a future society where humans and robots might function together, also provides with a handful of conflicting situations that the household robot R781 has to resolve in order to achieve one of its goals: keeping baby Travis alive.

The scenario itself made us think about how such a robot could be implemented as a rational agent and how would a normative system graft onto it. Granted, McCarthy's story is offering a few clues about the way the robot is reasoning and is reaching decisions, but he also lets us wonder about the architecture of a rational agent, such like R781, and how it would function in a normative context. In the following, we will be trying to look exactly into that: how can the well known Beliefs-Desires-Intentions (BDI) rational agent architecture be combined with a normative system to give what we call a normative BDI agent?

The paper is structured as follows: in the next section we will review the state of the art in the field of normative agent systems and present several approaches which we found of great value to our work. In the third section we describe our proposal for normative BDI agents, which will be supported by the case study scenario in the fourth section. In the fifth section we will present the implementation details for our agent. Finally we will sum up the conclusions of our research so far.

STATE OF THE ART

Agents, norms, normative agent systems

One of the first key points is defining the notion of norm. This turns out to be a bit more difficult than expected in the context of intelligent agents. Having become foundation stones of the way we function as a society, norms are now spread in most activities and domains (law, economics, sports, philosophy, psychology etc.), therefore becoming complex to represent given their different needs and their multiple facets. We would be interested in such definitions specific to the field of multiagent systems (MAS). Since this domain itself is very much interdisciplinary, defining a norm remains a challenge. For example, we would be interested in a definition applicable to social groups, since MAS, can be seen as models of societies. In 2 the definition of a norm is given as "a principle of right action binding upon the members of a group and serving to guide, control, or regulate proper or acceptable behaviour". On a slightly more technical approach, in distributed systems norms have been defined as regulations or patterns of behaviour meant to prevent the excess in the autonomy of agents. [START_REF] Criado | Towards a normative bdi architecture for norm compliance[END_REF]

NoA agents

Kollingbaum and Norman 4 study what happens when a new norm is adopted by an agent: what is the effect of a new norm on the normative state of the agent? Is a newly adopted norm consistent with the previously adopted norms? To this extent they propose a normative agent architecture, called NoA, which is built as a reactive agent. The NoA architecture is fairly simple and it comprises of a set of beliefs, a set of plans and a set of norms. Further on, they formalize way an agent will adopt a norm following the consistency check between a newly adopted norm and its current normative state.

Using some of the ideas of NoA, we will try to work on what we consider to be its limit, which is the lack of consistency check not only against the normative state, but also against the mental states. Finally, a norm instance is derived from an abstract norm n a , by grounding all the variables in n a according to a given a belief set.

The main drawback of the approach is the lack of coverage concerning the topic of norm acquisition. Therefore, a big challenge will be to integrate this approach, with the consistency check presented in section 2.2, as well as finding a good way to integrate everything with the classic BDI agent loop. [START_REF] Wooldridge | An Introduction to MultiAgent Systems[END_REF]

Worst consequence

An important part of our work will focus on solving conflicts between newly acquired norms and the previously existing norms or the mental contexts of the agent. Beforehand we draw from some of the definitions given by Ganascia in. [START_REF] Ganascia | An agent-based formalization for resolving ethical conflicts-Belief change, Non-monotonic reasoning and Conflict resolution Workshop -ECAI[END_REF] Those will later help us define what a conflict set is and how we can solve it.

Definition 2.2. Given (φ 1 , ..., φ n , φ ) ∈ L n+1 ¬ , φ is a consequence of (φ 1 , ..., φ n ) according to the belief-set B (we write φ = csq(φ 1 , ..., φ n )[B] if
and only if:

• φ ∈ (φ 1 , ..., φ n ) or • ∃Φ ⊆ (φ 1 , ..., φ n ) s.t. Φ → φ ∈ B or • ∃φ ∈ L ¬ s.t. φ = csq(φ 1 , ..., φ n )[B] ∧ φ = csq(φ 1 , ..., φ n , φ )[B]
Definition 2.3. φ is worse than φ given the belief-set B (we write φ c φ ) if and only if one of the consequences of φ is worse than any of the consequences of φ .

• ∃η ∈ L ¬ s.t. η = csq(φ)[B] and • ∃φ ∈ L ¬ s.t. φ = csq(φ )[B] ∧ η c φ [B] and • ∀φ ∈ L ¬ , if φ = csq(φ )[B] then η c φ [B] ∨ η φ [B] Notation: ∀(φ, φ ) ∈ L ¬ , φ φ [B] means that φ and φ are not comparable under B, i.e. neither φ c φ [B] nor φ c φ[B].
Definition 2.4. α and α being subsets of L ¬ , α is worse than α given the belief-set B (we write α c α [B]) if and only if:

• ∃φ ∈ α.∃η ∈ α s.t. φ c η[B] and • ∀η ∈ α .φ c η[B] ∨ φ η[B]

A NORMATIVE EXTENSION ON THE BDI ARCHITECTURE

Normative BDI agents

Starting from the classic BDI execution loop 5 we will now introduce and discuss a solution for taking into account the normative context of a BDI agent.

First, the agent's mental states are initialized. The main execution loop starts with the agent observing its environment, including for percepts received from other agents. There is a multitude of ways in which an agent can detect the emergence of norms in its environments and a good review of those is given in. [START_REF] Savarimuthu | A categorization of simulation works on norms[END_REF] For simplicity, we will consider that norms are transmitted via messages and our agent will consider the sender of such a message to be a trusted normative authority. The agent will acquire a new abstract norm n a (see section 2.3) and store it in the Abstract Norms Base(ANB). Drawing from the normative contexts described in, [START_REF] Criado | Towards a normative bdi architecture for norm compliance[END_REF] we define the ANB as a base of in-force norms. It is responsible with the acquisition of new norms based on the knowledge of the world as well as the deletion of obsolete norms. However, at this point the agent is simply storing an abstract norm which it detected to be in-force in its environment; it has not yet adhered to it! At this point a normative BDI agent should take into account the norms which are currently in force and check whether the instantiation of such norms will have any impact on its current normative state as well as on its mental states.

Consistency check

Let us define the notion of consistency between a plan p and the currently in-force norms to which an agent has also adhered and which are stored in the Norm Instance Base (NIB). By contrast to the ANB, the NIB stores the instances of those norms from the ANB which become active according to the norm instantiation bridge rule (see below). Definition 3.1. A plan instance p is consistent with the currently active norms in the NIB, if the effects of applying plan p are not amongst the forbidden effects of the active norms and the effects of current obligations are not amongst the negated effects of applying plan p.

consistent(p, N IB) ⇐⇒ (ef f ects(n F i ) \ ef f ects(n P i )) ∩ ef f ects(p) = ∅ ∧ ef f ects(n O i ) ∩ neg ef f ects(p) = ∅
The types of consistency / inconsistency which can occur between a newly adopted norm and the currently active obligations are:

• strong inconsistency occurs when all plan instantiations p which satisfy the obligation o are either explicitly prohibited actions by the NIB or the execution of such a plan would make the agent not consistent with its NIB • strong consistency occurs when all the plan instantiations p which satisfy the obligation o are not amongst the explicitly forbidden actions by the NIB and the execution of such a plan would keep the agent consistent with the NIB • weak consistency occurs when there exists at least one plan instantiation p to satisfy obligation o which is not explicitly prohibited by the NIB and the execution of such a plan would keep the agent consistent with its NIB.

The rules for prohibitions and permissions are analogous. The second point of consistency check is formalizing the rules about the consistency between a newly adopted abstract obligation and the current mental states of the agent. Prior to this, we define: Definition 3.2. A plan instance p is consistent to the current intentions set I of the agent when the effects of applying the plans specific to the current intentions are not among the negated effects of applying plan p.

consistent(p, I) ⇐⇒ ∀i ∈ I.(ef f ects(π i ) ∩ ef f ects(p) = ∅
Where by π i we denote the plan instantiated to achieve intention i.

The types of consistency / inconsistency states between a plan and an intention are almost similar to those between a plan and the norms in the NIB:

• strong inconsistency occurs when all plan instantiations p which satisfy the obligation o are not consistent with the current intentions of the agent • strong consistency occurs when all plan instantiations p which satisfy the obligation o are consistent with the current intentions of the agent • weak consistency occurs when there exists at least one plan instantiation p which satisfies the obligation o and is consistent with the current intentions of the agent

Norm instantiation

If in the ANB there exists an abstract norm with modality M about C and according to the belief-set the activation condition is true, while the expiration condition is not, then we can instantiate the abstract norm and store an instance of it in the NIB. In this way, the agent will consider the instance of the norm to be active.

We thus obtain the updated Norm Instance Base (NIB) containing the base of all in-force and active norms, which will further be used for the internalization process.

Solving the conflicts

When following its intentions an agent will instantiate from its set of possible plans (capabilities) P ⊆ L ¬ , a set of plans Π(B, D). We call Π(B, D) the conflict set, according to the agent's beliefs and desires. Sometimes, the actions in Π(B, D) can lead to inconsistent states. We solve such inconsistency by choosing the maximal non-conflicting subset from Π(B, D). Definition 3.3. Let α ⊆ Π(B, D). α is a maximal non-conflicting subset of Π(B, D) with respect to the definition of consequences given the belief-set B if and only if the consequences of following α will not lead the agent in a state of inconsistency and for all α ⊆ Π(B, D), if α ⊆ α then the consequences of following α will lead the agent in an inconsistent state.

The maximal non-conflicting set may correspond to the actions required by the newly acquired norm or, on the contrary, to the actions required by the other intentions of the agent. Thus, an agent may decide either: i) to internalize a certain norm, if the consequences of following it are the better choice or ii) to break a certain norm, if by 'looking ahead' it finds out that the consequences of following it are worse than following another course of actions or respecting another (internalized) norm.

A more comprehensive example of how this works is presented in section 4.

Norm internalization

With the instantiation process being finished and the consistency check having been performed, the agent should now take into account the updated normative state, which will become part of its cognitions. Several previous works treat the topic of norm internalization 8 arguing which of the mental states should be directly impacted by the adoption of a norm. For this initial state of our work and taking into account the functioning of the BDI execution loop, we propose that an agent updates only its desire-set; subsequently, this will impact the update of the other mental states in the next iterations of the execution loop.

AN EXAMPLE

Now that we have seen how a BDI agent becomes a normative BDI, adapting to to norm occurrence, consistency check and internalization of norms, let's get back to Prof. John McCarthy's story. 1 And let's focus on the short episode with which we started this article, considering that R781 functions according to the normative BDI loop which we have just described. R781's initial state is the following: AN B : ∅ N IB : F, love(R781, T ravis) Bset : B, ¬healthy(T ravis) , B, isHungry(T ravis) , B, csq(¬love(R781, x)) c csq(heal(R781, x)) Dset : D, ¬love(R781, T ravis) , D, isHealthy(T ravis) Iset : ∅ When R781 receives the order from his mistress he will interpret it as a normative percept and the brf(...) method will add a corresponding abstract obligation norm to the ANB structure. Since the mistress doesn't specify an activation condition nor an expiration condition (the two "none" values), R781 will consider that the obligation should start as soon as possible and last for an indefinite period of time. Its normative context is updated:

AN B : O, none, none, love(R781, T ravis) N IB : F, love(R781, T ravis) , O , love(R781, T ravis)
At this point, R781 will update the desire-set and will detect an inconsistency between the obligation to love baby Travis and the design rule which forbids R781 to do the same thing. Therefore, it will try to solve the normative conflict looking at the consequences of following each of the paths, given its current belief-set. In order to do so, let us take a look at the plan base of R781: PLAN heal(x, y) { pre: ¬ isHealthy(y) post: isHealthy(y) Ac: feed(x, y) } PLAN feed(x, y) { pre: ∃ x.(love(x, y) ∧ hungry(y)) post: ¬ hungry(x) } As we know from the story, R781 uses the Internet Paediatrics book to find out that if a baby is provided with love while hungry, it is more likely to accept being fed and therefore not be hungry any more. This is described by the feed(x, y). Moreover, R781 also knows how to make someone healthy through the heal(x, y) plan, given that a-priori, that someone is not healthy. In our simplified scenario we consider that R781 knows how to do so only by feeding someone.

Instantiating its plans on both of the paths, R781 will come up with the following maximal non-conflicting sets: {love(R781, T ravis), f eed(R781, T ravis), heal(R781, T ravis)} and {¬love(R781, T ravis)} And since the current belief set has a rule defining that not loving someone has worse consequences than healing that person, R781 will opt for the first maximal non-conflicting subset. This means R781 will be breaking the prohibition of not loving baby Travis and will follow the action path given by the first maximal non-conflicting subset {love(R781, Travis), feed(R781, Travis), heal(R781, Travis)}, while dropping the contrary. Further on, it will create an intention to achieve this state and will begin the execution of such a plan (simulating love towards baby Travis turns out to involve such plans as the robot disguising himself as human, displaying a picture of a doll as his avatar and learning what it considers to be the "motherese" dialect, mimicking the tone and the language of a mother towards her son).

IMPLEMENTATION

We implemented our normative BDI agent framework and the test scenario we've described using the Jade patform for agents development, in conjunction with Jadex 9 -a Jade extension for rational agents. Using the separation of concerns principle, we have isolated the mental states of the agent from its normative states. The mental states are all specified in Jadex's Agent Description File (ADF), which is an XML-based file format for specifying each BDI-like structure. In our case:

• Beliefs. A Java class was implemented to model the beliefs according to the needs of our agent; in general, we have paid particular attention to the plan implementations and what were the requirements for fully specifying such a plan, based on the beliefs. Finally, our model of the beliefs was referenced by the belief-set in the ADF.

• Desires. They are described inside the ADF, by means of goals.

• Intentions. They are described by means of those plans needed to be executed to achieve the goals specific to an intention. Basically, each plan is specified by means of a Java class, inheriting from Jadex's generic Plan class. Finally, the implemented plans are linked to goals in the ADF.

On the normative side of the agent, however, things were not as clearly defined. Hence the need to adopt a format for describing the normative state and storing the normative information related to our agent. Several reasons pointed us to XML as a representation language for the normative part of the agent. First of all, we wanted this part to follow the logic imposed by Jadex and to make things as easily interoperable as possible. Then, we needed a flexible enough language, which could offer us the possibility of adequately expressing the norm formalization that we have adopted. We have thus built a small XML controlled vocabulary for easily representing the normative state of our agent in which two distinct sections can be identified: the norm-bases and the consequences value-base.

CONCLUSION

In this paper we have presented an adaptation of the BDI execution loop to cope with potential normative states of such an agent. We have given a motivation for choosing the mental states model of Bratman which we have enriched with capabilities of reasoning about norms. We have investigated several previous relevant work in the domain in order to come up with a formalization of such issues as norm instantiation, norm consistency, solving consistency conflicts and norm internalization. Finally, we have provided with an intriguing study scenario, inspired from Professor McCarthy's science fiction short story "The Robot and The Baby".

FUTURE WORK

One of the limitations of our work which we would like to address in the future is the issue of norm acquisition. Whereas our work is providing with a very simple case of norm recognition, several interesting ideas have been explored based on different techniques. A good review of those as well as a description of a norm's life cycle is given in. [START_REF] Savarimuthu | A categorization of simulation works on norms[END_REF] Out of those specific approaches, we will probably focus on learning based mechanisms, namely machine learning techniques and imitation mechanisms for norm recognition.

In terms of real-world applications, there are a number of scenarios in which we would like to study the behavior of our agent. Ranging from healthcare agents to military drones, we would like to model and examine how our normative BDI agent will deal with classic conflictual situations appearing in different human activities.


Amelia Carolina Sparavigna

The Retreat of the Terminus of Gangotri Glacier in Google Earth Images

Keywords: Satellite Images, Google Earth, Glaciers

Here we use a time series of satellite images in Google Earth for measuring the retreat of the terminus of Gangotri Glacier, one of the largest glaciers of Himalayas. From 2010 to 2016, the terminus moved of about 150 meters.

The terminus of the Gangotri Glacier, one of the largest glaciers of Himalayas is known as the Gomukh, the "Mouth of the Cow", due to the presence of an ice cave (Figure 1). The place is situated at a height of about 4,000 meters, in Uttarkashi district of the state of Uttarakhand, India. From the Gomukh, the Bhagirathi River originates, which is one of the main sources of the Ganges River. For this reason, the Gomukh is also worshipped by many religious people. It is therefore a pilgrimage site, along with Gangotri, as well as a trekking destination [START_REF] Malhotra | Gangotri and Gaumukh: a trek to the holy source[END_REF][START_REF] Bisht | Tourism in Garhwal Himalaya: With Special Reference to Mountaineering and Trekking in Uttarkashi and Chamoli Districts[END_REF].

Figure 1:

The Gomukh of Gangotri Glacier (Courtesy, Barry Silver, Wikipedia) As told in [START_REF] Vv | Wikipedia[END_REF], on July 26, 2016, following heavy rains in Uttarakhand, it was reported that the front end of Gomukh "was no more, as a large chunk of the glacier had collapsed and was washed away". Actually, any glacier terminus moves due to the environmental conditions and therefore the shape of it is continuously changing. The Gomukh had been studied in several researches. In [START_REF] Drew | A Retreating Goddess? Conflicting Perceptions of Ecological Change near the Gangotri-Gaumukh Glacier, Journal for the Study of Religion[END_REF], it is told that the Gangotri glacier is a significant body of ice whose loss could cause ecological problems and "also threaten the continuity of cultural-religious practices intertwined with the sacred Ganga River" [5]. The discussion on the evolution of the glacier in [START_REF] Drew | A Retreating Goddess? Conflicting Perceptions of Ecological Change near the Gangotri-Gaumukh Glacier, Journal for the Study of Religion[END_REF] is presented on scientific bases, "along with ethnographically documented perspectives on retreating glaciers and other ecological transformations". However, in [START_REF] Raina | Is the Gangotri Glacier Receding at an Alarming Rate?[END_REF], which is a work presenting a remarkable analysis of historical data, it was concluded that the retreat of the Gomukh is of 15 to 18 meters a year, and then "definitely not alarming. Some of the glaciers in Columbia and even in Alaska have recently shown an annual retreat of more than 200 metres". "If we assume that this glacier (Gangotri) will continue to retreat, say at the rate of 15 to 18 metres a year, as it is doing at present, even then it will take almost 2,000 to 1,600 years for the glacier" to disappear [START_REF] Raina | Is the Gangotri Glacier Receding at an Alarming Rate?[END_REF]. The data given in [START_REF] Raina | Is the Gangotri Glacier Receding at an Alarming Rate?[END_REF] are in agreement to those in [START_REF] Bhambri | Frontal recession of Gangotri Glacier, Garhwal Himalayas, from 1965 to 2006, measured through high resolution remote sensing data[END_REF]. In this reference, the researchers report the frontal area changes at the terminus of Gangotri Glacier, after a work based on highresolution satellite data from 1965 to 2006. The results show that the glacier exhibited retreat up to 819 ± 14 m (it means about 20 metres a year) and lost 0.41 ± 0.03 sq. km at its front from 1965 to 2006. The results from satellite images were supported by in-situ field survey conducted by the Geological Survey of India. Other measures of the retreat are 15, 19 and 40 meters a year, given in References [START_REF] Vohra | Himalayan glaciers[END_REF], [START_REF] Naithani | Geomorphological evidences of retreat of Gangotri glacier and its characteristics[END_REF] and [START_REF] Mukherjee | Recession of Gangotri glacier through 20thcentury[END_REF] respectively. In [START_REF] Bahuguna | Himalayan glacier retreat using IRS 1C PAN stereo data[END_REF], stereo images from Indian Remote Sensing satellite (IRS)-1C covering Gangotri glacier were used to identify its terminus and measure its retreat with respect to the position in a topographical map of 1962. As concluded in [START_REF] Bahuguna | Himalayan glacier retreat using IRS 1C PAN stereo data[END_REF], "the precise location of the snout", that is, the terminus, "as identified on the images based on its shadow and relief variation in stereo images can serve as a standard reference for future monitoring of its movement or finding rate of retreat annually". In fact, we can also use the Google Earth imagery to study the retreat of the Gangotri glacier. Actually, we have a time series of high resolution images of the Gomukh. In the Figure 2, we see two of them. The lower panel of the Figure 2 needs some processing to see the terminus. We can use the Retinex filter of GIMP, the GNU Image Manipulation Program, to see the details in the shadows of the images (for some examples of the use of GIMP Retinex, see for instance [START_REF] Sparavigna | Rock glaciers on Mars viewed with Google Earth[END_REF][START_REF] Sparavigna | Gimp Retinex for Enhancing Images from Microscopes[END_REF][START_REF] Sparavigna | Retinex filtering and thresholding of foggy images[END_REF][START_REF] Marazzato | Retinex filtering of foggy images: generation of a bulk set with selection and ranking[END_REF]).

Figure 2 :

 2 Figure 2: The terminus of Gangotri in two images of Google Earth.

Figure 3 :

 3 Figure 3: The lower panel of Figure 2 processed by means of the Retinex tool of GIMP.

Figure 4 :

 4 Figure 4: Superimposing the two panels of the Figure 2, the lower one enhanced by means of the Retinex tool as in the Figure 3, we can see the terminus retreating from the position of 2010 to that of 2016. The two termini are given by the red curves. A point which has the same position in the two images (red circle) is used as reference point. From this point to the pin of the coordinates, there are 318 meters. Using this scale we can evaluate the retreat. It is of 150 ± 20 meters. It corresponds to a rate of 25 meters a year.

After processing the 2016 image, we can use the images of 2010 and 2016 to create a new image superimposing them by means of GIMP. The result is given in the Figure 4. In this image, the two termini are evidences by red curves. To measure the displacement of the terminus, we used a reference point and the ruler of Google Earth. The terminus moved of 150 ± 20 meters in 6 years, corresponding to 25 meters a year. This measure is close to those given in [7] and [9]. The method that we have used here is the same used for determining the rate of movement of the sand dunes [16][17][18].


Pierre Joyot

Jean Trunzler

Accounting for Incompressibility in Reproducing Kernel Particles Meshless Approximations

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Accounting for incompressibility in reproducing kernel particle meshless approximations

Pierre Joyot* -Jean Trunzler* -Francisco Chinesta** * LIPSI-ESTIA Technopole Izarbel, F-64210 Bidart {j.trunzler, p.joyot}@estia.fr ** LMSP UMR 8106 CNRS-ENSAM 151 Boulevard de l'Hôpital, F-75013 Paris francisco.chinesta@paris.ensam.fr ABSTRACT. Meshless approximations seem to be an appealing choice for simulating forming processes involving large transformations because they allows alleviating the mesh constraints. However, because the novelty of these techniques a lot of questions are today unresolved. One of these open problems is the treatment of incompressibility which as well known impose some restrictions on the choice of the approximation spaces. The accurate treatment of incompressibility is a key point in the simulation of forming processes because the plastic flow can be in fact considered as incompressible. This paper introduces the problematic as well as some possibilities for taking into account the incompressibility in the context of mixed formulations, making special emphasis in a kind of Hermite approximations. RÉSUMÉ. Les méthodes sans maillage semblent bien adaptées à la résolution de problèmes en grandes transformations, mais à cause de leur jeunesse quelques difficultés subsistent. Une de ces difficultés concerne la prise en compte de l'incompressibilité qui induit des contraintes sur le choix des espaces fonctionnels d'approximation. La prise en compte de l'incompressibilité est fondamentale dans le cadre de la simulation des procédés de mise en forme. Dans ce travail nous présentons la problématique liée à l'incompressibilité matérielle ainsi que différentes possibilités pour sa prise en compte dans le cadre des formulations mixtes, où des approximations de type Hermite sont particulièrement adaptées. KEYWORDS: RKPA, meshless approximation, mixed formulations, incompressibility.

MOTS-CLÉS : RKPA, méthodes sans maillage, approximations mixtes, incompressibilité.

Introduction

Meshless methods are an appealing choice for constructing functional approximations (with different degrees of consistency and continuity) without a mesh support. Thus, this kind of techniques seem to be specially appropriated for treating 3D problems involving large displacements, due to the fact that the approximation is constructed only from the cloud of nodes whose positions evolve during the material deformation. In this manner neither remeshing nor fields projections are a priori required.

Other important point lies in the easy introduction of some known information related to the problem solution within the approximation functional basis. For this purpose, different reproduction conditions are enforced in the construction of the approximation functions. This approach has been widely used in the context of the moving least squares approximations involved in the diffuse meshless techniques [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF] as well as in the element free Galerkin method [START_REF] Belytschko | Meshless methods: an overview and recent developments[END_REF].

Meshless approximations seem to be an appealing choice for simulating forming processes involving large transformations because they allows alleviating the mesh constraints. However, because the novelty of these techniques a lot of questions are today unresolved. One of these open problems is the treatment of incompressibility which as well known impose some restrictions on the choice of the approximation spaces. The accurate treatment of incompressibility is a key point in the simulation of forming processes because the plastic flow can be in fact considered as incompressible. This paper introduces the problematic as well as some possibilities for taking into account the incompressibility in the context of mixed formulations, making special emphasis in a kind of Hermite approximations.

Reproducing Kernel Particle Approximation

The Reproducing Kernel Particle Method is based in the following approximation:

( ) ( , ) ( ) a f x x s h f s ds [1]
where h defines the size of the support of the kernel function, which approaches the delta distribution as 0 h . If one want to reproduce polynomials up to a certain degree the kernel function must be constructed in an appropriate manner. For the sake of simplicity we consider the 1D case and second degree reproducing conditions, which imply:

2 2 1 ( , )1 ( , ) ( , )
x s h ds

x x s h sds x x s h s ds [2]
As usual kernel functions do not verify the reproducing conditions [2] [START_REF] Liu | Reproducing kernel particle methods for structural dynamics[END_REF] proposed to introduce a correction in the form:

( ) ( , ) ( , ) ( ) a f x C x x s x s h f s ds =?
with the correction function given by: ( , ) ( ) ( )

T C x x s H x s b x [4] where 2 ( ) 1, ( ), ( ) T H x s x s x s [5]
and the components of ( ) b x must be determined for verifying the reproduction conditions.

The discrete form of Equation [3] results:

( ) ( , ) ( , ) ( ) a J J J J f x C x x x x x h f x [6] with ( , ) ( ) ( ) T J J C x x x H x x b x [7]
Thus, the reproduction conditions result:

2 2 1 ( ) ( ) ( , ) ( ) ( ) ( , ) ( ) ( ) ( , ) T J J J T J J J J T J J J J H x x b x x x h x H x x b x x x h x x H x x b x x x h x [8]
that can be written in the following matrix form:

( ) ( ) ( ) ( , ) ( ) T J J J J R x R x H x x x x h b x [9] or ( ) ( ) ( ) R x M x b x [10]
where R denotes the reproduction vector and ( ) M x the moment matrix. Thus, it results:

1 ( ) ( ) ( ) b x M x R x [11] and 1 ( , ) ( ) ( ) ( ) ( ) ( , ) ( ) ( ) ( ) J T a J J J J C x x x J J J f x H x x M x R x x x h f x N x f x [12]
where ( ) J N x denote the approximation shape functions that do not verify the delta property.

REMARK 1. -This procedure can be generalized for reproducing arbitrary functions related to the problem solution [START_REF] Trunzler | Discontinuous derivative enrichment in RKPM[END_REF].

REMARK 2. -In the 2D case the reproduction vector results REMARK 4. -The trial and test functions related to the variational formulation of a generic problem can be approximated using the scheme [6] from which the discrete problem is derived as in the usual finite element method. However due to the non polynomial character of the kernel function richer integration rules must be considered.

The main drawbacks related to this kind of discretization technique is the difficulty to impose essential boundary condition (the shape functions related to internal nodes don't vanish on the domain boundary) as well as the difficulty to control the shape function support when the nodal distribution evolves significantly (as encountered for example in forming processes simulations).

Moreover, when this kind of approximation is considered for treating incompressible media one must verify the LBB stability conditions which implies, when a mixed variational formulation is considered, richer approximation in the velocity or displacement fields than in the pressure one (which is in fact the Lagrange multiplier related to the incompressibility restriction).

In the context of our former works, we have proved that by changing the reproduction condition (second degree for the richer field and first degree for the one related to the Lagrange multiplier) the derived approximation is not stable. The same conclusion was obtained by changing simultaneously the shape functions support. Specific strategies have been proposed for solving some kind of problems, as the Stokes one for example, but its generality is limited. For this reason we introduce in this paper richer approximations based on the use of Hermite-RKPA approximations.

Efficient Hermite -RKPM approximations

In order to derive a more general H-RKPA approximation, we consider (in the 1D case for the sake of simplicity):

( ) ( ) ( ) ( ) a J N a J J f x f x C df x df x dx dx [13]
If we are looking for second degree reproduction condition, we can write:

2 2 2 2 N J J J J a bx cx a bx cx C b cx b cx [14] that implies: 2 2 1 1 0 1 2 0 1 2 N J J J J x x x x C x x [15] or 2 2 1 0 1 0 1 1 2 2 N T J J J J x x C x x x x [16]
Defining now the reproduction matrix

2 1 0 ( ) 1 2 R x x x x [17]
the correction matrix is expressed by:

T T C H b [18] with 2 1 0 ( ) 1 2( ) J J J J H x x x x x x x x [19]
Thus, from Equations [16] and [18] it results: ( ) ( )

N T J J R x R x H b M b [20]
which allows to compute b , -1

( ) b M R x . Thus, we can finally write:

T ( ) ( ) ( ) ( ) ( ) ( ) a J J N N a J J J J f x f x f x b H N df x df x df x dx dx dx [21]
REMARK 5. -It is easy to verify that the derivatives involved in the previous expression are diffuse derivatives that can be employed in collocation discretization schemes but not discretization schemes operating on the variational formulation. However, usual collocation schemes using the nodes as collocation points lead to a rectangular discrete linear system, because there are several dof at each node.

Numerical results

Approximation tests

In this section we analyze the approximation of the polynomials. Firstly, we consider an approximation with third degree of consistency, i.e.

2 3 2 1 0 1 ( ) 2 3 x R x x x x x
involving two degrees of freedom at each node, one related to the unknown field and the other one to its derivative. When the approximated polynomial has a degree lower or equal to the one considered in the reproduction conditions, that polynomial is approximated with the machine precision. Figure 1 depicts the approximated (stars) versus the exact (solid curve) solutions related to the polynomial:

x x x x

As the polynomial and its derivative are exactly reproduced, the diffuse derivative corresponds to the exact one, which justifies the excellent accuracy obtained in the derivatives. Errors are of the order of the machine precision.

Finally, we are analyzing the convergence of the approximation when a polynomial that cannot be exactly reproduced:

( ) f x x x

If we compute the approximation error using the L 2 norm, we obtain fourth order convergence with respect to the log of the total number of degrees of freedom.

Second order PDE

In order to analyze the behavior of the H-RKPM approximation just defined when it is used for solving second order partial differential equation, we first consider the 1D problem defined in [0,1]:

2 2 2 12 6 d f x x dx YKVJ ( 0) 0 f x CPF ( 1) 0 f x whose exact solution results 4 5 ( ) f x x x .
The approximation is built from the reproduction matrix: The convergence order obtained was only first order, significantly lower than the third order expected when the L 2 norm is considered, may be induced by the numerical integration. REMARK 6. -The boundary conditions have been imposed by removing the rows of the linear system related to the first degree of freedom of the first and last nodes, imposing at its place the Dirichlet boundary conditions. If these conditions are imposed by using Lagrange multipliers the solution improvement is inappreciable. REMARK 7. -The variational formulation can be integrated with or without integration by parts. Obviously, due to the approximation continuity and derivability no integration by parts is required. In the case of integration by parts the boundary integrals must be considered because the test functions do not vanish on the domain boundary. In both cases the solution accuracy does not exhibit significant differences. REMARK 8. -When a regular distribution of quadrature points is used for integrating the variational formulation the shape function derivatives can be computed from the value of the shape functions at the integration points using finite differences. In practical applications the accuracy is not significantly degraded. REMARK 9. -When the variational formulation is considered, as argued previously, we can not use for defining the field derivative its diffuse derivative given by the second row of Equation [ u and v being velocity vector v components. We consider the driven cavity problem which implies the following boundary conditions:

2 1 0 ( ) 1 2 R x x x x
2 ( 0, ) (0,0) ( 1, ) (0,0) ( , 0) (0,0) ( , 1) ( ,0) T T T T v x y v x y v x y v x y x x [24]
If all the unknown fields are approximated using standard RKPM with second degree of consistency, with a single degree of freedom by node, spurious oscillations are found in both the velocity and the pressure field. Figures 6, 7 and 8 depict both velocity components as well as the pressure field. being the pressure approximated using standard RKPM, whose reproduction vector results:

2 2 ( , ) (1, , , , , )

T p

R x y x y x xy y

which implies 3 degrees of freedom at each node for each velocity component, and a single degree of freedom at each node for the pressure field. Figures 9, 10 and 11 depict the computed velocity and pressure fields. We can notice that the use of the H-RKPM allows stabilizing the discretization. The convergence order is only first order, possibly due to the numerical integration.

Conclusions

In this paper we have proposed and applied an original Hermite reproducing kernel particle approximation (H-RKPM) for solving second order PDE. The only interest of using this kind of approximation is the possibility to reduce significantly the size of the shape functions support. Even if the support is slightly higher that the inter-nodal distance, the reproduction order can be increased by increasing the number of nodal degrees of freedom (that are the successive field derivatives).

However, its interest in mixed formulations that must verify the LBB stability condition is out of discussion. A centered formulation, as usually used in moving least squares framework, could be more suitable, simplifying the computational aspects and may be reducing the numerical errors. This is a work in progress.

  xy y for the same degree of reproduction. REMARK 3. -The integration weights are not present in Equation [6] because they are implicitly considered in the expression of ( ) b x .

Figure 2

 2 Figure 1. Approximation of polynomial

Figure 2 .

 2 Figure 2. Approximation of the derivative of

  versus de exact one.

Figure 3 .

 3 Figure 3. Approximation of polynomial

Figure 4

 4 Figure 4 depicts de exact and the computed solution when 20 nodes are considered, each one involving 2 degrees of freedom: the field and its derivative.

Figure 4 .

 4 Figure 4. Computed (stars) versus exact (solid curve) solution

  21]. REMARK 10. -This technique can be applied directly for solving fourth order PDE which requires the extra imposition of some field derivative on the domain boundary with either collocation or Lagrange multipliers.We consider now the 2D problem defined in [0has three nodal degrees of freedom: the unknown field and its first derivatives with respect to both space coordinates. The Dirichlet conditions are imposed, as previously, via collocation. Figure5depicts the numerical versus de computed solutions. As both solutions coincide (at the machine precision) they are indistinguishable in the representation.

Figure 5 .

 5 Figure 5. Computed and exact solutions: both are perfectly superposed

Figure 6 .

 6 Figure 6. x-component of the velocity field in the driven cavity problem

Figure 7 .

 7 Figure 7. y-component of the velocity field in the driven cavity problem

Figure 8 .

 8 Figure 8. Pressure field in the driven cavity problem

Figure 9 .

 9 Figure 9. x-component of the velocity field in the driven cavity problem


M Bergmann

L Weynans
email: lisl.weynans@inria.fr

A sharp Cartesian method for incompressible flows with large density ratios

Keywords: Incompressible flows, bifluid flows, finite-differences, projection method, cartesian grid, level-set, jump conditions across interface, interface unknowns Ecoulements incompressibles, écoulements bifluides, différences finies, méthode de projection, grille cartésienne, level-set, conditions de saut au travers

A new Cartesian method for bifluid incompressible flows with high density ratios is presented. The specificity of the method relies on a sharp second order numerical scheme for the spatial resolution of the discontinuous elliptic problem for the pressure. The Navier-Stokes equations are integrated in time thanks to a fractional step method based on the Chorin scheme and discretized in space on a Cartesian mesh. The bifluid interface is implicitly represented using a level set function. The numerical tests show the improvements due to this sharp method compared to classical first order methods.

Introduction

In this paper we are concerned with incompressible bi-fluid flows like air and water, and by the accurate description of the phenomena occurring at their interface. We thus present a sharp Cartesian method for the simulation of incompressible flows with high density and viscosity ratios. This method is inspired from the second-order Cartesian method for elliptic problems with immersed interfaces developed in [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF].

Cartesian grids are an attractive alternative to body fitted meshes. Indeed, they avoid complex mesh generation as well as mesh adaptation when unsteady interfaces are considered. Moreover, the numerical resolution of the governing equations can be simplified with an easy parallelization and the use of standard linear algebra libraries. Generally speaking, numerical schemes are easy to implement on a Cartesian mesh because a dimensional splitting is often possible. However, some numerical modeling is necessary near a complex interface that does not fit the background Cartesian grid. This is the case for fluid structure interface and moreover for bi-fluid interface where the properties of the flow are discontinuous. Indeed, applying naively a numerical scheme originally devised for a flow with constant or continuously varying density will lead to a a non-consistent treatment of the interface. Most of the time, it will result in severe stability issues if the density ratio is large as highlighted in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF] and references therein. Therefore, as already mentioned, one has to devise specific numerical schemes at the vicinity of the interface. This region is called narrow band and is the set of numerical points that have at least one neighbor on the other side of the interface.

Conservative or non-conservative approaches can both be used to face this issue. Among the nonconservative approaches, one solution is to regularize in the vicinity of the interface the properties of the fluids, so that the density, viscosity, and their derivatives are continuous in the whole computational domain. This idea leads to the well known "Continuous Surface Force" (CSF) method [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], where the discontinuous quantities are smoothed near the interface, and in case of a fluid with surface tension, this surface tension is taken into account as a smooth volume force instead of a surface force. This method is widely used (see for instance [START_REF] Smereka | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] and [START_REF] Galusinski | On stability condition for bifluid flows with surface tension: Application to microfluidics[END_REF]) because it offers a straightforward way to implement the presence of two fluids in an already existing mono-fluid Navier-Stokes code. However, the exact way that the regularization should be performed is not always clear, and spurious oscillations at the bi-fluid interface can appear due to errors in the pressure gradient computations. Another non-conservative method introduced by Kang, Fedkiw and Liu [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] after the CSF is the Ghost Fluid Method (GFM). It is based on a first-order method developed in [START_REF] Liu | A boundary capturing method for poisson's equation on irregular domains[END_REF] to solve an immersed interface elliptic problem, with a dimensional splitting making the method easy to implement. The resulting linear system is symmetric and has the same structure as the usual matrix to discretize a Poisson equation with variable coefficient on a Cartesian grid. This method has been used successfully in numerous works, for instance [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF] and [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF]. One drawback is that the method is only first-order accurate near the interface [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF] and a loss of conservativity of the momentum of each fluid near the interface can occur leading to erroneous velocities. Non-conservative methods are often associated with a level-set representation of the interface [START_REF] Smereka | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] because the level-set method is itself intrinsically non-conservative at the discrete level, and convenient to use on a Cartesian grid.

The other family of methods is based on the conservative form of the Navier-Stokes equations, where mass and momentum fluxes of each fluid are explicitly computed, see for instance [START_REF] Rudman | A volume-tracking method for computing incompressible multifluid flows with large density variations[END_REF], [START_REF] Sussman | A coupled level set and volume of fluid method for computing 3d and axisymmetric incompressible two-phase flows[END_REF], [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF], [START_REF] Hu | A conservative interface method for compressible flows[END_REF] and more recently [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. An explicit interface representation is necessary even if the interface do not coincide with grid points. Conservative methods are generally more stable that non-conservative methods. The price for this increased stability is an additional amount of work due to the interface reconstruction, which can be performed from informations carried by Lagrangian markers or by cell quantities such as volume fractions.

In this paper we aim to preserve as much as possible the simplicity of the Ghost-Fluid Method of [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], avoiding an explicit identification of the volume fractions near the interface, while improving the accuracy and stability of the pressure computation. We thus propose a new method where the discontinuities across

Inria Ω -(φ < 0) ρ -, µ - Ω + (φ > 0) ρ + , µ + Γ (φ = 0) n Figure 1: Sketch of the computational domain.
the interface are taken into account in a sharp way with a second-order scheme inspired from [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF]. This second-order treatment improves the conservativity of the method, as it will be proved numerically in the section devoted to numerical validations. After having described the governing equations for the incompressible bi-fluid flows that we consider ( §2), the discretization of these equations in each fluid and at the interface are presented ( §3-4). The second-order numerical resolution of the elliptic problem arising from the computation of the pressure is introduced ( §5), and the overall is validated on several two-dimensional numerical test cases ( §6).

Governing equations 2.1 Flow configuration

We consider a rectangular domain Ω filled with two viscous incompressible fluids with different densities and viscosities. The subdomains Ω -and Ω + corresponding to the two fluids are separated by an interface Γ. These domains are generally defined with a scalar function φ that takes different values in each subdomains with a fixed value on the interface. For instance we chose φ = 0 on Γ, φ > 0 in Ω + and φ < 0 in Ω -. The unit normal to the interface is n and the unit tangent vector is η. The density is

ρ = ρ -+ H(φ )(ρ + -ρ -),
and the viscosity is

µ = µ -+ H(φ )(µ + -µ -),
where H is the Heaviside function, i.e. H(φ ) = 1 if φ > 0 and H(φ ) = 0 if φ < 0. Finally, the twodimensional velocity vector is u = (u, v).

The flow is modeled in each subdomain with the incompressible Navier-Stokes equations:

ρ(u t + (u • ∇)u) = -∇p + ∇ • τ + ρg, ∇ • u = 0,
with g the gravitational acceleration vector, and τ the viscous stress tensor:

τ = µ(∇u + ∇u T ).

Jump conditions at the interface

The above equations are completed by jump conditions at the interface Γ between the two fluids. In what follows jumps are defined by [ψ] = ψ +ψ -.

• The first ones describe the balance between the normal stresses at the interface and the surface tension σ , with κ the local curvature of the interface Γ,

[p -2µ(∇u • n, ∇v • n) • n] = σ κ, (1)
[µ(∇u • n, ∇v • n) • η + (∇u • η, ∇v • η) • n] = 0. ( 2
)
Several others jump conditions can be derived from continuity properties across the interface.

• For a viscous fluid the velocity field is continuous across the interface

[u] = 0, (3)
[v] = 0. (4)
• The material derivative of (3)-( 4) is zero, therefore

0 = ∂ [u] ∂t + (u • ∇)[u] = - ∇p ρ + ∇ • τ ρ + g , which leads to ∇p ρ = ∇ • τ ρ . (5)
• The jump condition for the pressure p can be simplified. We differentiate in the tangential direction the jump on the velocity:

∂ u ∂ η = 0, ∂ v ∂ η = 0.
Moreover, because the velocity is divergent-free on each side of the interface,

0 = [∇ • u] = [(∇u • n, ∇v • n) • n + (∇u • η, ∇v • η) • η] .
Combining the two last relationships, we obtain

[(∇u • n, ∇v • n) • n] = 0. Consequently [p] = 2[µ](∇u • n, ∇v • n) • n + σ κ. (6)
We will use this equation to compute the pressure jump at the interface.

Inria

Moving interface

For a moving interface as it is usually the case for bi-fluid problems, the flow density and viscosity are updated with φ tracking the interface thanks to the transport equation

φ t + û • ∇φ = 0, (7)
where the velocity fields û coincides with the flow velocity field u on the interface Γ. Different choices for the value of û in Ω + and Ω -can be a priori used. A natural choice is to consider û = u in the whole domain. Another choice is the extension velocity introduced in [START_REF] Adalsteinsson | The fast construction of extension velocities in level set methods[END_REF].

3 Numerical method outside the discontinuities For simplicity reason the computational domain is discretized on a two-dimensional uniform Cartesian grid with a grid spacing ∆x = ∆y = h, but the following approach stands for non uniform Cartesian meshes. The points on the Cartesian grid are named with indices such as M i, j = (x i , y j ). We denote by u i j the approximation of u at the point (x i , y j ). The narrow band is the set of points having at least one neighbor on the other side of the interface. An accurate computation of the interface is thus necessary.

Interface computation

In order to provide an accurate discretization in the vicinity of the interface we need some geometric information about the interface, and thus a good choice for φ is necessary. This information can be provided by the level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamiltonâ jacobi formulations[END_REF]. We refer the interested reader to [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], [START_REF] Sethian | Evolution, implementation, and application of level set and fast marching methods for advancing fronts[END_REF] and [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] for recent reviews of this method. A useful choice for the level set function is the signed distance function to the interface:

φ (x) =    dist Γ (x) if x ∈ Ω + , -dist Γ (x) if x ∈ Ω -, 0 if x ∈ Γ. (8)
The zero isoline thus represents implicitly the interface Γ immersed in the computational domain.

We assume that the interface is smooth enough, so that the derivatives of the level-set function in the vicinity of the interface are well-defined. A useful property of the level set function is an easy computational of its normal

n(x) = ∇φ (x) |∇φ (x)| . (9)
In the same way, the curvature of the interface can be computed with the formula

κ = ∇ • n. (10)
In what follows, the level-set function is advected by the fluid velocity u in the whole domain

φ t + u • ∇φ = 0.
The computation of the level-set function should be performed very accurately when one deals with moving interfaces, Indeed, as the level-set method is not intrinsically conservative, a lack of accuracy in the computation of the level set evolution results often in a substantial loss of mass for one of the fluids. It can also increase the problem of transfer of momentum between both fluids and generate spurious velocity oscillations. Moreover, if ones wants to compute the curvature of the interface from ( 10), the level-set function needs to be accurate enough (at least third-order) so that the finite difference formulas used to discretize [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free-surface flows[END_REF] are consistent.

Unfortunately, the property of the signed distance function is usually lost when the interface evolves with the flow velocity. The norm of the level set gradient can be far from unity. These gradients variations of the level-set are harmful to the accuracy of the numerical evaluation of the normal to the interface and the curvature. To circumvent the problem, Sussman et al. [START_REF] Smereka | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] introduced a reinitialization algorithm to recover the signed distance function through the resolution of the eikonal equation

|∇φ | = 1.
Several methods have been developed over the years to perform this reinitialization step, either by using a relaxation method and searching a stationary solution to a time dependent Hamilton-Jacobi equation [START_REF] Smereka | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] or using a Gauss-Seidel based method as in fast marching methods [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF] or fast sweeping methods [START_REF] Cheng | Fast sweeping algorithms for a class of hamiltonjacobi equations[END_REF].

The usual level set strategy for an evolving interface is thus the following:

• Use a transport equation to update φ .

• From time to time, reinitialize φ with the signed distance function.

One of the most widespread option in the literature is to combine a fifth-order WENO scheme [START_REF] Guang Shan Jiang | Efficient implementation of weighted eno schemes[END_REF] with a RK3 scheme, for the transport and for the reinitialization through a relaxation equation. It provides a high-order yet stable resolution. But, although the reinitialization procedure performed with such a numerical scheme may improve mass conservation, it also introduces some error by slightly moving the interface as noticed in [START_REF] Russo | A remark on computing distance functions[END_REF]. For this reason, Russo and Smereka [START_REF] Russo | A remark on computing distance functions[END_REF] introduced a subcell fix taking into account the interface location in the reinitialization procedure. This technique was extended to a higher order accuracy in [START_REF] Duchene | Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes[END_REF] through the use of third-order ENO schemes near the interface.

Usually, this reinitialization steps are performed uniformly every each n iterations. A recent study [START_REF] Luddens | Enablers for high-order level set methods in fluid mechanics[END_REF] proposes a strategy to sample the reinitialization steps based on interface deformation criteria. This should also be coupled with the high-order decentered reinitialization scheme of [START_REF] Duchene | Second-order accurate computation of curvatures in a level set framework using novel high-order reinitialization schemes[END_REF] near the interface.

In what follows, we will use the classical option of the fifth-order WENO scheme for the spatial discretization. Since it is the most commonly used technique in the literature, it will allow to distinguish the effects of the new scheme from the effects of the reinitialization technique.

Flow computation

We use a classical non-incremental projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des equations de navier-stokes par la méthode des pas fractionnaires ii[END_REF], i.e. the guess value for the pressure in the prediction step is zero. This choice avoids instability issues due to the discontinuous pressure values when the interface moves. It avoids to wonder which jump conditions across the interface should be satisfied by the pressure correction. Indeed, these jump conditions for the pressure correction are not based on physical considerations.We thus compute successively:

u * -u n ∆t = -(u n • ∇)u n + 1 ρ (∇ • τ) n + g (prediction step), (11)
u n+1 -u * ∆t = - ∇p ρ (correction step). ( 12
)
The convective terms are computed with a fifth-order WENO scheme, and the viscous terms with an explicit second-order centered finite-difference scheme. The time integration is performed with a firstorder explicit Euler scheme.

Inria

The pressure is computed through the resolution of a Poisson equation in order to enforce the divergencefree condition. At each grid point outside the narrow band around the interface, the following relationship is satisfied:

∇ • 1 ρ ∇p = ∇ • u * ∆t , (13)
with Neumann boundary condition:

∇p ρ = u b -u * ∆t , (14)
where u b is the value of the velocity to be imposed on the external boundaries. We will provide additional details about the jump conditions that have to be satisfied across the interface for this problem in subsection 4.3. The overall algorithm is the following:

1. Prediction: evaluate convective and diffusive fluxes and compute u * , 2. Interface evolution: convect the level-set with velocity u and re-initialize if necessary, 3. Construction and resolution of the linear system for the pressure, 4. Correction step: update velocity with pressure gradient.

We compute at each iteration an adaptive time step taking into account the restrictions due to convection, viscosity, surface tension and gravity.

The convective time step restriction is given by

∆t |u| max ∆x + |v| max ∆y ≤ 1. ( 15
)
with |u| max and |v| max the maximum magnitudes of the horizontal and vertical velocities. The viscous time step restriction is given by

∆t max( µ - ρ -, µ + ρ + )( 2
∆x 2 + 2 ∆y 2 ) ≤ 1 (16)
The time step restriction associated with the surface tension is simiar to the one in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] and in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF]:

∆t σ |κ| min(ρ + , ρ -) min(∆x 2 , ∆y 2 ) ≤ 1. ( 17
)
We only apply it for grid points in the narrow band.

In what follows, all the unknowns are collocated in space on Cartesian meshes. In all the applications considered we have not observed any odd-even coupling. However, this point can be fixed using some corrections [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF][START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF].

Specific sharp method for the discontinuities across the interface

The values of the viscosity and the densities are discontinuous across the interface. Therefore, if the numerical scheme described above was applied on the grid points in the narrow band, stability problems could occur because the approximations of the following terms are not consistent:

• • • • • • • • • • • • • • • • • • • • • • • • • • I i, j,E = I i+1, j,W I i+2, j+1,N = I i+2, j+2,S i-1 i i+1 i+2 i+3 i+4 i+5 j-1 j j+1 j+2 j+3 • • • • • • • • • Figure 2:
Example of geometrical configuration, with regular grid points (white circles), irregular points (black circles), and interface points (diamonds) with the two possible notations.

• viscous terms,

• divergence of the predicted velocity,

• elliptic operator in the correction step,

• gradient of the pressure.

We have to devise a specific treatment for the points near the interface, for each of these steps. The computation of the convective terms is not mentioned in the above enumeration because it is performed with a fifth-order WENO scheme, which provides automatically spatial adaptivity. Therefore, we assume that the gradients computed with the WENO scheme are decentered near the interface, and consequently, consistent. Moreover, the level-set function is classically evolved with such a scheme, because it is crucial to have a good accuracy in the computation of the interface evolution. Therefore, it seems coherent to have the same numerical scheme for the convection of the interface and the convection of the fluids.

Let us introduce some notations. A grid point is defined to be irregular if at least one of its neighbors is on the other side of the interface, i.e. if the sign of φ changes between this point and at least one of its neighbors, see Figure 2. All the other points are called regular grid points.

We define the interface point I i, j,E = ( xi, j,E , y j ) as the intersection of the interface Γ and the East segment [M i j M i+1 j ], if it exists. Similarly, the interface points I i, j,W = ( xi, j,W , y j ), I i, j,N = (x i , ỹi, j,N ) and I i, j,S = (x i , ỹi, j,S ) are respectively defined as the intersection of the interface and the West [M i-1 j M i j ], North [M i j M i j+1 ] and South [M i j-1 M i j ] segments . With this notation the same interface point can be described in two different ways I i, j,S = I i, j-1,N or I i, j,E = I i+1, j,W .

The set of interface points is denoted Γ h , see Figure 2 for an illustration. These points are used to impose the jump conditions across the interface in the numerical scheme. Since the pressure is discontinuous across the interface, two unknowns are created, one for each side of the interface.

Viscous terms

We follow a continuous approach and regularize the quantities used for the computation of the viscous terms. It has been proven in [START_REF] Desjardins | A spectrally refined interface approach for simulating multiphase flows[END_REF] and [START_REF] Herrmann | The influence of density ratio on the primary atomization of a turbulent jet in crossflow[END_REF] that this continuous approach provides correct accuracy for high Reynolds numbers flows. It has also been used successfully in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF] and [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. A sharp approach for the viscous terms could probably improve the accuracy of the simulations. However, the complexity of the computations would be increased due to the treatment of the jump conditions for the viscous terms (2) implying derivatives of the velocity components in both normal and tangential directions. Moreover, if one needs to use an implicit treatment of the viscous terms, such a sharp treatment would become more complex to handle.

The viscosity and the inverse of the density are regularized in this step by a discrete convolution [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]:

16 μi, j = 4 µ i, j + 2 µ i+1, j + 2 µ i-1, j + 2 µ i, j+1 + 2 µ i, j-1 + µ i+1, j+1 + µ i+1, j-1 + µ i-1, j+1 + µ i-1, j-1 , 16 ρi, j = 4
ρ i, j + 2 ρ i+1, j + 2 ρ i-1, j + 2 ρ i, j+1 + 2 ρ i, j-1 + 1 ρ i+1, j+1 + 1 ρ i+1, j-1 + 1 ρ i-1, j+1 + 1 ρ i-1, j-1 .
Then we discretize the viscous terms with a classical second-order centered scheme. Due to this continuous approach for the viscous terms, the jump condition for the pressure (6) writes

[p] = σ κ. ( 18
)
and we avoid the computation of the term 2[µ](∇u • n, ∇v • n) • n across the interface.

Divergence of predicted velocity

The predicted velocity u * obtained after the prediction step ( 11) is defined only on grid points. We need to compute the divergence of this predicted velocity to solve the elliptic equation [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF]. However, since the two fluids have different properties across the interface and the derivatives of the velocity are not necessarily continuous, we need to use a decentered stencil on each side of the interface. Consequently, we have to compute two values for u * on each interface point, one for each side of the interface. In practice, as jump conditions for u * are not available, we perform simply linear extrapolations from the grid values on the interface points. Formally this is equivalent to a standard first-order decentered scheme. Then, to compute the divergence of u * on an irregular grid point M i, j , we use a standard five point stencil, see Figure 3. More precisely, we denote u * S the value of the solution on the nearest point to M i, j in the south direction (possibly an interface point), with coordinates (x S , y S ). Similarly, we define u * N , u * W and u * E and the associated coordinates (x N , y N ), (x W , y W ) and (x E , y E ). The discretization reads

∇ • u * i, j = u * N -u * S x N -x S + v * E -v * W y E -y W .

Elliptic problem near the interface

To compute the pressure, according to the jump conditions ( 1) -( 6) presented in section 2, it is necessary to solve an elliptic problem with discontinuous values of the solution and its derivative across the interface:

∇ • ( 1 ρ ∇p) = ∇ • u * ∆t in Ω + ∪ Ω -, [p] = σ κ + 2[µ] (u n , v n ) • n on Γ, ∇p ρ = ∇ • τ ρ on Γ. • • • • • • M i+2, j+1 • • • • • M i, j i-1 i i+1 i+2 i+3 i+4 i+5 j-1 j j+1 j+2 j+3
Figure 3: Example of geometrical configuration, with points involved in the discretization of the divergence of the predicted velocity and pressure gradient at grid points M i, j and M i+2, j+2 in black.

Because the viscous terms are handled with a regularization approach, they can not be considered as discontinuous anymore. Therefore, the elliptic problem to solve becomes

∇ • ( 1 ρ ∇p) = ∇ • u * ∆t , in Ω + ∪ Ω - (19)
[p] = σ κ on Γ,

∇p ρ = 0 on Γ, (20)
avoiding the computation of the terms (u n , v n ).n and ∇ • τ ρ . The details of the resolution of this elliptic problem will be provided in section 5.

Pressure gradient near the interface and correction step

In order to keep a consistent discretization, the gradient of the pressure p computed through the resolution of the elliptic problem in the last subsection, is also computed with an adapted decentered stencil near the interface, see Figure 3. More precisely, with the same notations as before, the discretization reads

(∇p) i, j =    p N -p S x N -x S p E -p W y E -y W    . ( 22
)
If one of the discretization point is an interface point, we consider the value of the interface point corresponding to the same subdomain than point M i, j . Indeed we recall that, since the pressure is discontinuous across the interface, two pressure unknowns (one for each subdomain) are computed at an interface point.

If no interface point is involved, the numerical scheme reduces to the classical second-order central finite differences scheme.

Inria 5 Numerical resolution of elliptic problems with immersed interfaces

The elliptic problem with discontinuous values across an interface ( 19) -( 21) is solved with the secondorder method developed in [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF]. The accuracy of this method is based on the use of unknowns located at the interface. The size of this linear system is thus augmented with two unknowns for each interface point. These interface unknowns are used to discretize the flux jump conditions and the elliptic operator accurately enough to get a second-order convergence in maximum norm. Actually, to this purpose, near the interface the elliptic operator needs to be discretized with a first-order truncation error, and the fluxes with a second-order truncation error. For a visual explanation of the discretization we refer to Figure 4. The advantage of using this method, compared to the reference work of [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] is that the jump conditions in the correction step are solved with second-order accuracy instead of first-order. The drawback is that the linear system is not symmetric anymore and it is solved with the preconditioned GMRES method.

Discrete elliptic operator

We use a standard five-point stencil including the grid point M i, j and its nearest neighbors in each direction: interface or grid points. More precisely, we denote p S the value of the solution on the nearest point in the south direction, with coordinates (x S , y S ). Similarly, we define p N , p W and p E and the associated coordinates (x N , y N ), (x W , y W ) and (x E , y E ). Since the density is piecewise constant, the discretization reads

∇.( 1
ρ ± ∇p) i, j = 1 ρ ± ∆p = 1 ρ ± p N -p i j x N -x i - p i j -p S x i -x S x N -x S 2 + 1 ρ ± p E -p i j y E -y j - p i j -p W y j -y W y E -y W 2 , (23)
where ρ ± stands for ρ + or ρ -.

Discrete flux transmission conditions

At each interface point we create two additional unknowns, called interface unknowns, and denoted by p ± i, j,γ with γ = E,W, N or S. The interface unknowns carry the values of pressure on each side of the interface.

Contrarily to [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF], we do not have a jump condition on the normal derivative, but on the whole gradient, as expressed in formula [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamiltonâ jacobi formulations[END_REF]. To obtain the same number of equations and unknowns we have to chose in which direction we want to project this gradient equality. As we use a Cartesian grid, it is easier to discretize the xand y-derivatives than derivatives in other directions. Consequently, we discretize the following jump conditions at each interface point I i, j,γ , with γ = N, S,W, E.

p + i, j,γ -p - i, j,γ = σ κ i, j,γ , (24) 1 ρ
+ (∂ x p + ) i, j,γ - 1 ρ -(∂ x p -) i, j,γ = 0 if γ = E,W. (25)
1

ρ + (∂ y p + ) i, j,γ - 1 ρ -(∂ y p -) i, j,γ = 0 if γ = N, S. (26)
We want the truncation error of the discretization of flux equality [START_REF] Rudman | A volume-tracking method for computing incompressible multifluid flows with large density variations[END_REF][START_REF] Russo | A remark on computing distance functions[END_REF] to be second-order accurate in order to solve the problem with a second-order accuracy. A possible configuration of the interface is illustrated in Figure 4. In the x-direction, it is straightforward to compute a second-order approximation of the x-derivative with three a priori non equidistant points. For example we approximate the flux on the

• • • • M i+2, j+2 • M i+3, j • • • • • i-1 i i+1 i+2 i+3 i+4 i+5 j-1 j j+1 j+2 j+3 • • • • • I i+1, j+1,N • •• • I i, j-1,E i-1 i i+1 i+2 i+3 j-1 j j+1 j+2
Figure 4: Left: points involved in the discretization of the elliptic operator at grid nodes M i+2, j+2 and M i+3, j in black, right: Example of stencils for the discretization of jump conditions. Points involved in the discretization of the x-derivative of the pressure at interface point I i, j-1,E and y-derivative of the pressure at interface point I i+1, j+1,N in black. For I i, j-1,E both derivatives are expressed with second-order accuracy while for I i+1, j+1,N the left derivative is expressed with second-order and the right derivative with first-order accuracy.

left side of interface point I i, j,E , if it exists, with the values of p on the points M i-1, j , M i, j and I i, j,E with the formula:

(∂ x p ± ) i, j,E ≈ (p i-1, j -p ± i, j,E )(x i -xi, j,E ) ∆x(x i-1 -xi, j,E ) - (p i, j -p ± i,,E j )(x i-1 -xi, j,E ) ∆x(x i -xi, j,E ) . ( 27
)
The right x-derivative is approximated in the same way.

(∂ x p ± ) i, j,E ≈ - (p i+2, j -p ± i, j,E )(x i+1 -xi, j,E ) ∆x(x i+2 -xi, j,E ) + (p i+1, j -p ± i,,E j )(x i+2 -xi, j,E ) ∆x(x i+1 -xi, j,E ) . ( 28
)
The same discretization holds for the y-derivative. The formulas ( 27) and ( 28) are consistent if both grid points involved in the formula, for instance M i-1, j and M i, j , belong to the same domain. If on one side of the interface the two closest grid points aligned with the intersection point do not belong to the same subdomain, then the second-order discretization is not possible anymore. In this case, we use a first-order discretization involving only two points: the interface point and the closest grid point on the same side of the interface. Such a case is illustrated on Figure 4. In fact, this first-order discretization is equivalent to the ghost-fluid method [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF].

Let us notice that, because we use a dimensional splitting for the jump conditions across the interface, it is quite straightforward to eliminate the interface unknowns from the linear system. We simply inject expressions ( 27) and [START_REF] Sethian | Evolution, implementation, and application of level set and fast marching methods for advancing fronts[END_REF] in the jump condition [START_REF] Rudman | A volume-tracking method for computing incompressible multifluid flows with large density variations[END_REF], and use the resulting equality to express p ± i, j,E as a function of p i-1, j , p i, j , p i+1, j and p i+2, j . This expression for p ± i, j,E can then be used in the discretization of the elliptic operator [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF].

The local curvature κ i, j,γ at the interface point I i, j,γ is computed in the following way. We first compute on all irregular grid points the value with centered second-order finite-difference formulas. Then we perform a one-dimensional linear interpolation of these values on the interface points.

κ = φ 2 x φ yy + φ 2 y φ xx -2φ x φ y φ xy (φ 2 x + φ 2 y ) 3/2
6 Numerical validation

Parasitic oscillations

This first test case aims to assess the influence of the interface curvature error on the stability of the numerical scheme. A bubble is located at the center of the computational domain. Due to Laplace law and the concavity of the interface, the pressure inside the bubble is larger than the pressure outside. If the curvature of the interface is computed numerically, the errors due to the numerical approximation in the right-hand side of equation ( 24) cause small errors in the resolution of the pressure equation. These errors create artificial values of the velocity near the interface which should theoretically be zero. These artificial velocities are often called parasitic currents. The amplitude of the parasitic currents is an indication of the stability and the accuracy of the numerical method, and especially of the pressure computational step. Indeed, they are the only source of numerical errors.

Comparison with the Ghost Fluid and the CSF methods

We use the same parameters as in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF], where a Ghost-Fluid and a CSF method where implemented. The amplitude of the parasitic currents, compared to the results in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF], are reported in Table 1. The L ∞ and L 2 norms are computed over the whole domain Ω. The initial configuration is described in Figure 5. As it can be observed, the amplitude of the parasitic currents generated by our new method is several orders of magnitude smaller than those of the CSF method, and significantly lower than those of the Ghost-Fluid method when the grid is refined.

                   L = 2 cm, R = 1 cm, ρ int = 1000 kg.m -3 , µ int = 10 -3 Pa.s, ρ ext = 1 kg.m -3 , µ ext = 10 -5 Pa.s, σ = 0.1 N.m -1 (29)
6 R 9 R in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] or 10 R in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free-surface flows[END_REF] 2 R in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free-surface flows[END_REF] 3 R in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] gaz liquid R Figure 6: Initial fluid domain for the test case of the rising bubble in water in the references [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF] and [10]

Comparison with the Ghost-Fluid method

We consider air bubbles rising in water, as in the test case proposed for the Ghost-Fluid method in [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF].

The value of the physical parameters are

                       R = 1/
We consider two cases: a small bubble with R = 1/300 m and a large one R = 1/3 m. In the first case, the surface tension plays an important role in the evolution of the interface because of the high bubble curvature. In the second case, the surface tension has less influence, and larger deformations occur. The interface of the small bubble is plotted at times t= 0.,0.02,0.035, 0.05 in Figure 7. The interface of the large bubble is plotted at times t= 0.,0.2,0.35, 0.5 in Figure 8. Our numerical results are in good agreements with [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF]. [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free-surface flows[END_REF] and the level-set method [START_REF] Smereka | A level-set approach for computing solutions to incompressible two-phase flows[END_REF] This test case is taken from [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free-surface flows[END_REF], and inspired from a test case presented in [START_REF] Smereka | A level-set approach for computing solutions to incompressible two-phase flows[END_REF]. It gives us the opportunity to compare our method to another class of methods, based on the SPH formulation. The initial configuration is described on Figure 6. The value of the physical parameters are The evolution of the interface is plotted on Figure 9, for 120 × 200 grid points. We observe that the interface deforms in a way similar to the results in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free-surface flows[END_REF]. We impose a curvature threshold 1/h to allow the bubble break-up.

Comparison with SPH

                   R = 0.

Collapse of a water column (dam break)

This test case is studied in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF] and [START_REF] Chenadec | A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows[END_REF], and based on experiments conducted in [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF]. The initial configuration is a water column at rest in air. The initial height and width of the column are both 5.715 cm. The domain size is 40 cm×10 cm. The physical constants are the same than for the rising bubble ( §6.1.1). For more details, we refer the reader to [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. We present in Figure 10 the interface evolution at non-dimensional times T = t g/H = 0, 1, 2, 3, 4, with H the initial height of the water column. The computations are performed with 256 × 64 points. Figure 11 present the temporal evolution of the water front, compared to the experimental results [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF], to Ghost-Fluid results, and to the conservative method of Raessi and Pitsch [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. We observe that the front propagation is in agreements with the experimental results and the results of the conservative method [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF]. It means that, though the method is not strictly conservative, the numerical errors due to momentum transfer across the interface are not large enough to slow down the propagation of the front. It is not the case for instance for the Ghost-Fluid method, as it can be noticed in Figure 11 and has been reported in [START_REF] Raessi | Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method[END_REF].

Conclusion

We have developed a new method on Cartesian grids for the simulation of incompressible flows with large density ratios, based on a sharp resolution of the pressure term. This Cartesian scheme uses additional unknowns located on the interface to discretize with second-order accuracy the jump conditions across the interface. The viscous terms are treated with a regularizing approach which allows to eliminate terms in the jump conditions without damaging the accuracy of the results. Numerical results show that this new method leads to more accurate and stable results than reference methods. In the same time, the resulting numerical scheme remains simple to implement: it simply amounts to modifying the stencil of the pressure equation for irregular grid points by adding one additional point. Future works include an extension of the method to three-dimensional problems, possibly with interactions with solids, as for instance in [START_REF] Bergmann | An accurate cartesian method for incompressible flows with moving boundaries[END_REF]. We also aim to study in details the effects of the reinitialization procedure presented in [START_REF] Luddens | Enablers for high-order level set methods in fluid mechanics[END_REF].

Inria

Figure 5 :

 5 Figure 5: Test case of the static bubble with parasitic oscillations.

  300 m (small bubble) R = 1/3m (large bubble) ρ water = 1000 kg/m 3 , µ water = 1.137 × 10 -3 kg/ms, ρ air = 1.226 kg/m 3 , µ air = 1.78 × 10 -5 kg/ms, σ = 0.0728 kg/s 2 g = -9.8m/s 2

Figure 7 :Figure 8 :

 78 Figure 7: Evolution of the interface for the small bubble test case, resolution 80 × 120 (top) and 160 × 240 (bottom) .

Figure 9 :

 9 Figure 9: Evolution of the interface for the bubble test case from [10], resolution 120 × 200.

  025 m ρ water = 1000 kg/m 3 , µ water = 1.137 × 10 -3 kg/ms, ρ air = 1.226 kg/m 3 , µ air = 1.78 × 10 -5 kg/ms, σ = 0.0728 kg/s 2 g = -9.8m/s 2 (32)

Figure 11 :

 11 Figure 11: Evolution of the front of propagation: comparison between experimental date and several numerical methods: the Ghost Fluid method (non-conservative method), the conservative method of Raessi and Pitsch and our new method, The dimensionless location of the front z a is plotted as a function of the dimensionless time T = t g/H.

RR n°8926

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Ghost Fluid method CSF New method N L ∞ error L 2 error L ∞ error L 2 error L ∞ error L 2 error 16 8.08 ×10 -3 1.88 ×10 -3 3.55 ×10 -2 1.94 ×10 -2 5.21 ×10 -3 7.31 ×10 -5 32 3.42 ×10 -4 7.50 ×10 -5 3.12 ×10 -2 1.18 ×10 -2 9.26 ×10 -5 1.42 ×10 -6 64 5.13 ×10 -5 7.97×10 -6 2.12 ×10 -2 5.44 ×10 -3 1.36 ×10 -5 1.47 ×10 -7 128 2.79 ×10 -5 4.74 ×10 -6 6.44 ×10 -3 1.38 ×10 -3 2.22 ×10 -6 1.92 ×10 -8 Table 1: Comparison between the new method and the numerical results obtained in [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles: application à la désintégration assistée d'un jet liquide par un courant gazeux[END_REF] for the ghost-fluid method and the CSF method for parasitic oscillations, at time t = 1.

Comparison with a Volume of Fluid method

We now compare the behavior of our method to the Volume of Fluid method developed in [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. The density and viscosity ratio are both chosen to be one for this test-case. The coefficient σ is chosen so as to obtain an Ohnesorge number

. The maximum velocity is computed for varying grids at non-dimensional time

A comparison between our method and the Volume of Fluid method is presented in Table 2. The new method provides a better accuracy than the Volume of Fluid method for the coarsest grid. As expected, the Volume of Fluid method outperforms our new approach for finer meshes due to more sophisticated schemes near the interface. Nonetheless, Table 2 show a second-order accuracy for our new method.

∆x

error L ∞ for [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] error L ∞ for our method 2.5/16 7.3 ×10 -4 7.48 ×10 -5 2.5/ 32 4.5 ×10 -6 4.7 ×10 -6 2.5/64 5.5 ×10 -8 1.26 ×10 -6

Table 2: Numerical results for parasitic oscillations at non-dimensional time t * = 250 for [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] and our method.

Rising of air bubble in water

We study the evolution of fluid bubbles rising in an heavier fluid, and compare our results to several methods in the literature. The initial configuration is described in Figure 6.

Inria