Ingrid Falk

Fabienne Martin

Aspectual Flexibility Increases with Agentivity and Concreteness A Computational Classification Experiment on Polysemous Verbs

Keywords: lexical aspect, aspectual flexibility, French lexicon

We present an experimental study making use of a machine learning approach to identify the factors that affect the aspectual value that characterises polysemous verbs under each of their readings. The study is based on various morpho-syntactic and semantic features collected from a French lexical resource and on a gold standard aspectual classification of verb readings designed by an expert. Our results support the tested hypothesis, namely that agentivity and abstractness influence lexical aspect.

Introduction

It is generally accepted that the event type of a sentence and in particular its aspectual value result from a complex interplay between lexical features of the predicate on one hand and its linguistic context on the other. In previous work by [START_REF] Siegel | Learning Methods to Combine Linguistic Indicators: Improving Aspectual Classification and Revealing Linguistic Insights[END_REF], [START_REF] Zarcone | Computational Models for Event Type Classification in Context[END_REF] and [START_REF] Friedrich | Automatic prediction of aspectual class of verbs in context[END_REF], the event type of corpus utterances is predicted on the basis of syntactic and semantic features of the main predicate and the linguistic context thereof. In contrast, in this study, we investigate lexical (inner) aspect through a French lexical resource, "Les Verbes Français" [START_REF] François | Présentation de la classification des Verbes Français de Jean Dubois et Françoise Dubois-Charlier[END_REF], a valence lexicon of French verbs providing a detailed morpho-syntactic and semantic description for each reading (use) of a verb. Since the basic lexical units are the readings of the predicate rather than the (ambiguous) verbs, this resource allows one to study lexical aspect on an intermediate level between the coarse-grained predicate level and the very fine-grained corpus utterance level. Thanks to its detailed syntactic and semantic description, it also enables us to investigate the effect of particular complex semantic dimensions, as e.g. agentivity or abstractness, on the lexical aspect of verbs across their different readings.

Outline

In this work, we focus on the way the lexical aspect of a same verb varies across its different readings, using the lexical resource "Les Verbes Français" [START_REF] François | Présentation de la classification des Verbes Français de Jean Dubois et Françoise Dubois-Charlier[END_REF], henceforth LVF). We address this question through a supervised classification task. For this, we first extracted verbal readings from the LVF for a set of 167 verbs, chosen in such a way that each of the four Vendlerian aspectual classes are roughly equally represented. We manually annotated each of these readings on the basis of a finegrained aspectual classification based on eight aspectual classes (see below for a rough description). This manual annotation provides the gold standard for our classification experiments. For each annotated reading, we then automatically collected morpho-syntactic and semantic features from the LVF. Based on these features, we trained classifiers to automatically predict the aspectual class of the (manually annotated) readings. We then evaluated the dif-ferent classifiers used by comparing their accuracy in 10fold cross validation on the basis of the gold annotation. We conducted experiments with three feature sets and their combinations. The first set of features ('syn/sem' features) provides information on the morphosyntactic make-up of the verb, its argument structure and potential adjuncts, and indicates which suffixes are used in nominalisations and adjectives derived from the verb, when those are available (note that derived nouns and adjectives are coupled to a particular reading only if they preserve the meaning of the verb on this particular reading). The second feature set indicates whether the reading is agentive or nonagentive, and the third one whether the reading is concrete or abstract. The paper is organized as follows. Section 3 presents the resource we make use of. We first briefly review its design and theoretic principles and then describe the data used in our experiments. We detail our choices with respect to the analyzed instances (verbal readings) and their aspectual annotation. Finally, we list the features recorded in the lexical resource that we consider relevant for the aspectual classification. In Section 4, we present our experiments. We firstly describe the experimental approach and its technical implementation, and then discuss the results obtained.

The Data

The Resource -LVF

The LVF -"Les Verbes Français" [START_REF] François | Présentation de la classification des Verbes Français de Jean Dubois et Françoise Dubois-Charlier[END_REF] is a detailed and extensive lexical resource providing a systematic description of the morpho-syntactic and syntacticosemantic properties of French verbs. The basic lexical units are readings of the verbs, determined by their defining syntactic environment (argument structure, adjuncts) and a semi-formal predicate decomposition (provided through lexical templates called opérateurs, that use roughly the same inventory of labels and features as in the lexical templates found in [START_REF] Pinker | Learnability and Cognition: The Acquisition of Argument Structure[END_REF] or [START_REF] Jackendoff | Semantics and Cognition[END_REF]). The LVF covers roughly 12 300 verbs with a total of 25 610 usages. Each reading is associated with an elaborate morpho-syntactic and semantic description. In the following Table 1, we use sample entries for the verb élargir 'widen' to illustrate LVF's basic layout and give an idea of the underlying lexicographic and representation choices. Syntactic description (Table 1b). Each reading of a verb is coupled with a schematic representation of its acceptable syntactic frames. In principle, a verbal reading can be coupled with a transitive use, a reflexively marked use and an intransitive use unmarked by the reflexive, or be coupled with only a subset of these syntactic frames. Note that the reflexive is not systematically disambiguated between its different readings (anticausative, passive, purely reflexive, etc). However, the reflexive use is systematically associated with a particular semantic class, that indirectly provides some information on the range of readings the reflexive can have. Besides defining the argument structure (and potential adjuncts), the syntactic description also specifies some thematic features of the main arguments (e.g. whether the subject and direct object are animate and/or inanimate, whether the indirect object refers to a location, etc). Semantic description (Table 1a). Each entry in the LVF is also characterized by a semi-formal semantic decomposition, called opérateur, which provides a rough approximation of the verbal meaning instantiated on each reading. Each entry is therefore paired with a finite set of primitive semantic features and labels which, in a subsequent step, enable the classification of verbal readings into 14 wellestablished semantic classes (eg. communication verbs, psych-verbs, etc.). The semantic features and labels used in the semantic decomposition provide other cues about the type of verbs which is instantiated across each reading. For instance, for the reading 01 of élargir in Table 1a, 'r/d +qt [p]' roughly corresponds to BECOME(more(p)) ('r/d' stands for '(make) become'; '+qt' stands for an increase along a scale). From this, one can safely infer that under reading 01, élargir is a 'degree achievement' verb (a change of state verb with a multiple point scale in Beavers' (2008) terminology). The semantic decomposition may also contain abstract prepositions -see e.g. 'VRS' ('TOWARDS') in reading 04 of élargir -that stand for relations between arguments. To summarize, each reading is associated with two semantic descriptions: the semantic decomposition (the opérateur) and the semantic class. Additionally to the syntactic and semantic description, the LVF also indicates when a verb is formed through a derivational process, and in the positive case, provides information about the category of the verbal root, thus enabling one to identify deadjectival or denominal verbs. Finally, for each entry is specified which suffix is used for the available reading-preserving deverbal nominalisations and adjectives (-ment, -age, -ion, -eur, -oir, -ure or zero-derived nominalisations, and -able, -ant, -é or zero-derived adjectives).

The Annotation

We extracted 1199 entries (verbal readings) for the selected 167 verbs mentioned earlier. On average, each verb has roughly 15 readings, while 50% have more than 13. There are only 16 unambiguous verbs, while the largest number of readings is 37 and 19 (for 1 respectively 2 lemmas). Most lemmas (29) have 3 readings. Interestingly, the average number of 15 readings per verb very closely matches the number of event categories per verb obtained in the experiment reported by [START_REF] Marvel | Event categorization beyond verb senses[END_REF], who propose a new method of automatically categorizing event descriptions. 1The extracted verbal readings were manually annotated according to a fine grained aspectual classification on a 'telicity scale' of eight values. At the bottom of the scale are readings that are unambiguously stative (i.e. for which any other aspectual value is excluded), rated with 1 (STA). At the top are found achievement readings for which any other aspectual value is excluded, rated with 8 (ACH). At the middle of the scale are found predicates of variable telicity that can equally take a telic or atelic value within the context of the evaluated reading, rated with 4 (ACT-ACC). For instance, élargir under reading 01 is rated with 4, because, among others, both forand inadverbials are equally compatible with élargir on this reading. These predicates of variable telicity, that show no preference for the telic or the atelic reading on the given reading, are distinguished from 'weak accomplishment' readings, rated with 5 (W-ACC). Weak accomplishment readings are primarily telic, but nevertheless share some properties with activity readings (see [START_REF] Piñón | Weak and strong accomplishments[END_REF]). For instance, remplir 'fill' under its reading 01 (Pierre remplit le seau d'eau 'Pierre is filling the bucket with water') is analysed as a weak accomplishment rather than a verb of variable telicity (and therefore rated with 5 rather than 4). This is justified by the fact that under reading 01, remplir is more natural with an in-adverbial than with a for-adverbial. However, it is nevertheless still acceptable with for-adverbial under its 'partitive' reinterpretation (described e.g. by [START_REF] Smollett | Quantized direct object don't delimit after all[END_REF] or [START_REF] Champollion | The scope and processing of foradverbials: A reply to Deo and Piñango[END_REF]). This makes it a 'weak' accomplishment rather than a 'strong' one. 'Strong' accomplishments -like tuer in its reading 01 (Le chasseur tue un lapin 'The hunter is killing a rabbit') --, which are incompatible with the partitive interpretation of for-adverbials, are rated with 6 (S-ACC). Those accomplishments that share some properties with achievements are rated with 7 (ACC-ACH). The annotator evaluated each entry with a definite internal argument, in order to abstract away from the role of the determiner in the identification of the aspectual value of the VP (see e.g. [START_REF] Verkuyl | A Theory of Aspectuality : the Interaction between Temporal and Atemporal Structure[END_REF]). Additionally, sentences were systematically evaluated with a present tense, so as to avoid the possibility that lexical (inner) aspect is modified through grammatical (outer) aspect by coercion (see e.g. De Swart (1998)). 2In this first study, however, we used a coarser grained aspectual scale focusing on telicity. We therefore grouped the verb readings into the following classes: ATElic (rating 1-3), with VARiable telicity for those with rating 4, and TELic (rating 5 or more).

Table 2 gives an overview of the distribution of the aspectual ratings. The first interesting finding is that verbs display a considerable aspectual variability across readings.

The aspectual value of 2/3 of the 151 verbs with more than one reading (107) varies with the instantiated reading (on the 8 value scale). The greatest variation of 7 values on the telicity scale is attested for one verb. With respect to the coarser grained scale, roughly half of the verbs (69) has readings of only one of the three overarching aspectual classes.

The Features

For the design of the features for our classification we started from the features used in previous work [START_REF] Zarcone | Computational Models for Event Type Classification in Context[END_REF][START_REF] Siegel | Learning Methods to Combine Linguistic Indicators: Improving Aspectual Classification and Revealing Linguistic Insights[END_REF][START_REF] Friedrich | Automatic prediction of aspectual class of verbs in context[END_REF]) and looked for similar information in the LVF. Many of the features relevant for the aspectual value of utterances in context are obviously not encoded in the lexicon. This is in particular true of features related to outer/grammatical (i.e. non-lexical) aspect conveyed by tenses. This is also mostly true of the information conveyed by temporal adverbials, which also importantly contribute to the final aspectual value of an utterance in context. However, the LVF connects each verbal reading with specific morphological, syntactic and semantic features. Among such features, those that influence the lexical aspect of the verb in context are known to be pervasive: verbs encoding a result in their event structure are mostly accomplishment verbs; those encoding a manner but no result are mostly activity verbs (see e.g. Rappaport Hovav and Levin (1998) and subsequent work), verbs that subcategorise an indirect object are mostly bi-eventive, i.e. encode a result (see e.g. [START_REF] Pylkkänen | Introducing Arguments[END_REF]) and are therefore accomplishments, intransitive verbs with no transitive use with a -eur nominalisation are most likely to be unergative verbs and therefore atelic, etc. Luckily, most of these features relevant for lexical aspect are coded in one way or another through the LVF. So for instance, if, in the semantic decomposition of the LVF, 'BECOME' ('r/d') is present, we know that a 'pure' activity reading is excluded; if a movement verb has a prepositional phrase that encodes a direction rather than a location, we know that it has at least an atelic reading, etc. Also, in some cases, the LVF's semantic decomposition even contains a temporal adjunct that determines the aspectual value of the verb under the given reading. For instance, 'tps long' or 'tps court' in the opérateur seems to systematically define an atelic reading. Finally, some semantic classes give very clear hints to the lexical aspect of its members. For instance, readings instantiating the class of 'verbs of physical states and behaviors' are atelic, those instantiating 'enter/exit verbs' are telic, etc.

Overall we found that to a large extent, the information conveyed through the LVF features is comparable to [START_REF] Siegel | Learning Methods to Combine Linguistic Indicators: Improving Aspectual Classification and Revealing Linguistic Insights[END_REF] and [START_REF] Zarcone | Computational Models for Event Type Classification in Context[END_REF] linguistic indicators for lexical aspect. Most of the relevant linguistic information, in particular about the presence or absence of syntactic arguments and the nature of subject and/or objects could be derived from the descriptions. Further semantic cues are conveyed by the semantic decomposition (the opérateur), in particular its main component (BECOME, DO, ITER, etc.) and features indicating the encoding of a scale or a manner component. Also, as we mentioned above, the semantic class to which the verb belongs on a particular reading sometimes suffices to aspectually disambiguate the verb in context.

The first theoretical hypothesis we aimed to test through the experiments reported here is that agentivity and abstractness influence lexical aspect. The second (related) one is that under their agentive and/or concrete (literal) readings, verbs are aspectually more flexible (i.e. show a higher range of aspectual values) than under their nonagentive and/or abstract (figurative) readings. Consequently, we were also interested in the interplay of lexical aspect with two other groups of semantic features, namely the agentivity of the external argument and the abstractness of the verb under its particular readings. The LVF does not allow to directly determine the agentivity of an argument or the abstractness of a reading. However, it does provide some valuable cues that help to rate verbs along these two dimensions. Agentivity is for example suggested by the nature of the subject (animate or not). In addition, some of the semantic primitives used in the semantic decomposition as well as the semantic classes themselves suggest preferences for more or less agentivity.

Similarly, the nature of the internal argument also contributes to determine whether the verb is abstract or concrete for some verbs at least, that inanimate objects tend to be more abstract than animate ones (compare kill Mary and kill time/ the bottle). Also, the LVF systematically indicates whether the reading is literal and figurative, and often, figurative language positively correlates with abstractness. Finally, the semantic class is sometimes also indicative of more or less abstractness (for instance, when tuer 'kill' has a reading instantiating the class of psych-verbs, it is more likely to have abstract rather than a concrete sense, etc.).

To sum up, we extracted the following sets of features: syn/sem "classic" features based on syntactic/semantic diagnostics. These include information about the syntactic frames, the argument structure and the potential adjuncts. These features also comprise the main components of the semantic decomposition and the general semantic class, as well as the information about the morphological make-up of the verb (including the information on the category of the verbal root and the suffixes used in derived nominalisations and adjectives), and finally some information about the degree of polysemy of the verb (through the relative number of entries for one lemma).

agent Features indicating (more or less) agentivity. Besides the animacy of the subject, we also use some of the main components in the semantic decomposition, as well as semantic classes positively or negatively relevant for agentivity. Additionally, we also consider the availability of -age and -eur nominalisations as indicative of agentivity.

abs Features relating to abstractness or concreteness.

Here, we make use of the LVF's systematic classification in literal vs. figurative meanings, assuming that figurative readings are mostly abstract, while literal ones are mostly concrete. In addition, some components in the semantic decomposition indicate a preference for a concrete or abstract reading and some semantic classes (e.g. psych-verbs) indicate abstractness by themselves.

Aspectual Classification

For the aspectual classification, we proceeded as follows. As described in Section 3, we extracted 1199 readings from the LVF. On one hand, these were manually assigned gold aspectual values on the telicity scale of eight values (as described in the previous section). On the other hand, we automatically collected features from the LVF for each of these readings. Based on these features and the gold aspectual values, we trained classifiers to predict the aspectual values of the readings, with regard to a coarser grained telicity scale of three values: ATElic, VARiable telicity and TELic. Finally, we evaluated the results by looking at the effects of the features on the overall classification accuracy and by analyzing their contribution to the identification of each of the classes.

Classifiers We assumed a supervised learning setting and applied the classifiers3 shown in Table 4 with the implementation provided by Weka [START_REF] Hall | The weka data mining software: An update[END_REF] 4 to predict the 3-way aspectual class of the 1199 verb readings (ATE, VAR or TEL). We measured the performance of the classifiers by assessing the accuracy in 10-fold cross-validation and compare it to the accuracy of a baseline classifier which always assigns the majority class (rules.ZeroR). We see that all feature combinations and every classifier outperform the baseline where all readings are labeled with the majority class (the telic class). The best accuracy is achieved using the syn/sem features and lazy.kstar classifier. Using agentivity features and abstractness/concreteness features alone allows one to achieve similar results of 6-10 points above the baseline, suggesting that these features do indeed influence the aspectual value. However, the results improve considerably (about 8 points) when using the syn/sem diagnostic features. Again, when combining syn/sem features with agentivity and/or concreteness/abstractness features, we obtain results close to those using syn/sem features only.

Results

With respect to the classification algorithms, the best results were achieved for the lazy.kstar classifier. On the other hand, libsvm and naivebayes classifier appear less appropriate for this task since, in contrast to j48 and jrip, they did not outperform the other classifiers for any feature set. From these first set of experiments, we draw the following conclusions. Firstly, we found that we could classify the LVF verb readings with respect to their aspectual value based on features extracted from the LVF with reasonable accuracy. Secondly, our initial hypothesis that lexical aspect varies along the dimensions of agentivity and abstractness/concreteness is confirmed. However, agentivity and/or abstractness/concreteness alone are not sufficient to determine the aspectual value reliably: unsurprisingly, the syn/sem diagnostic features are a decisive factor for more accurate predictions.

Features favouring the overall classification performance

In the next experiment, we used the Wrapper [START_REF] Kohavi | Wrappers for feature subset selection[END_REF] method to select the feature subset producing the best classification accuracy. In our setting, this methodology consisted in using the prediction performance of the naivebayes learning algorithm to assess the relative usefulness of feature subsets obtained by forward selection5 .

In the following paragraph, we briefly describe the features selected this way for each of the three basic feature sets and the combination of all of these three sets.

syn/sem 7 selected features: temporal and manner complement or adjunct, evaluative meaning component, reflexive subcategorisation frame besides a transitive frame.

agent 4 selected features: (in)animacy of subject, components in the semantic decomposition suggesting more or less agentivity.

abs 5 selected features: opérateur suggesting abstractness and concreteness respectively, literal reading, (in)animacy of direct object.

syn/sem-agent-abs 10 selected features: temporal and manner complement or adjunct, iterative meaning component, subject always inanimate, components in the semantic decomposition indicating agentivity, features indicating abstractness.

Table 5 presents the accuracy in 10-fold cross validation based on the reduced feature sets. The results show that by using the reduced feature sets, the accuracy could be improved by 2-3 points. The best performing feature set and classification algorithm are the same as before, namely the syn/sem diagnostic features and the lazy.kstar algorithm (in this setting, the accuracy is of 67.56%). This leads us to the conclusion that the selected features are the most important for an accurate classification. As described above, among these selected features, many are related to agentivity and abstractness or concreteness. This, again, supports our initial hypothesis. However, while these experiments confirm the influence of agentivity and abstractness/concreteness on aspectuality, they do not lead to any conclusion about the orientation of this influence. We return to this issue in the course of the following section, dedicated to the feature analysis by class.

Feature analysis by class

Presentation. In this section, we present our application of a feature filtering method proposed in [START_REF] Lamirel | Enhancing classification accuracy with the help of feature maximization metric[END_REF], based on a feature maximization metric called FMC 6 . FMC allows us to analyze which features are the most helpful at the identification of specific classes. More concretely, it gives an indication of which features are mostly indicative of atelicity, telicity or variable telicity. Being "class oriented", this method also allows one to better tease out in which way agentivity and abstractness/concreteness influence lexical aspect. [START_REF] Lamirel | Enhancing classification accuracy with the help of feature maximization metric[END_REF] precision and recall measures. The precision of a feature f measures how discriminating f is for a class c, when compared to other features. Feature recall shows how characteristic f is of class c. Feature precision and recall are combined into feature f-measure (the harmonic mean). Only features with an f-measure above a threshold are selected. This method is independent of the classification algorithm and allows to rank the features by f-measure (or recall or precision) with regard to each class. In order to be able to take into account feature values for non-numeric features (so that we can analyze, for example, if a feature with a particular value is more characteristic of a class), we transform non-numeric features into binary features. The following example shows how this is done. For each initial feature F having e.g. two values V 1 and V 2for example the feature subj-always-inanim, which can be TRUE or FALSE --, two new features are introduced, i.e. F = V 1 and F = V 2 (subj-always-inanim=TRUE and subj-always-inanim=FALSE in our example).

For an instance, the new feature F = V 1 gets assigned the value TRUE (or 1) if for this instance F has the value V 1 and FALSE (or 0) otherwise. Respectively, F = V 2 is TRUE if for this instance F has the value V 2 and FALSE otherwise. The feature selection is then performed on these transformed features. FMC works as follows. First, the local Recall (F R f c ) and the local Precision (F P f c ) of a feature f in a class c are defined as follows:

F R f c = |i f c | |I f | , F P f c = |i f c | |I c |
where i f c is the set of instances (readings) having feature f in c, I c the set of instances in c and I f , the set of instances with feature f . The local F-measure F F f c is the harmonic mean of local recall and local precision. A feature f is said to maximize a class c iff its feature F-measure is higher for that class than for any other class and F c designates the set of features maximizing c. Finally, the feature F-measure F F c of a class c is defined as the average of the feature F-measure of the maximizing features for c:

F F c = f ∈Fc F F f c |F c |
The set S c of features that are considered characteristic of a given class c ∈ C results in:

S c = f ∈ F c |F F f c > F F (f ) and F F f c > F F
where

F F (f ) = c ∈C F F f c |C f | and F F = f ∈F F F (f ) |F |
C f represents the set of classes for which the feature f is present for some instances and F the set of features.

The two thresholds F F (f ) and F F can be interpreted as follows. F F (f ) is a threshold pertaining to a feature f and represents an average over the classes of the "usefulness" of f for each of the classes, where the "usefulness" is given by the feature F-measure F F f c . The second threshold F F is an average of the local thresholds over the whole set of features. According to this, a feature f is considered characteristic of a class c if, on one hand, its feature F-measure is more "useful" for this class than the average over classes and if, on the other hand, its contribution is greater than the average contribution of all features. Finally, the set of selected features is the union of the features selected for each class:

S = c∈C S c
Results and discussion. Table 6 presents the accuracy of the classifications with the features selected this way. The best performing configuration in this setting are the syn/sem features used together with the agentivity features and the rules.jrip classifier, achieving an accuracy of 64.80% similar to the accuracy in our initial configuration (Table 4). We can therefore conclude that it is meaningful to investigate the respective relevance of the selected features for the aspectual classes through FMC. In the best configuration, 42 binary features were selected. Among the selected syn/sem features, the most noteworthy ones are the following.

1. Availability of a canonical (non reflexive) passive.

2. Several features related to the derivational family; in particular, features indicating that the verb reading is not coupled with derived adjectives/nouns formed with particular suffixes, as well as features indicating that the verb is itself formed by derivation.

3. Presence or absence of a direct or indirect object. Among the features related to agentivity, the most important one is an 'at-least-one-true' agentivity feature, that was satisfied as soon as at least one feature of agentivity was satisfied. Another selected feature of agentivity is the non-agentive feature (satisfied whenever any of the agentivity features indicates non-agentivity). 7 Again, these results confirm the importance and effectiveness of the syn/sem features as well as the influence of the agentivity features. However, they again do not provide further insights about the direction of this influence.

In order to investigate the features most characteristic for each of the three coarse-grained aspectual classes (ATE, VAR, TEL), we looked at those features which were selected as relevant for one of those classes only (recall that the FMC selected features are based on the union of the features considered most characteristic for each class). These features are presented in Table 7. Line ATE in Table 7 shows that many of the features most helpful at predicting the atelic class are also those which suggest less agentivity. These features are the absence of any argument other than the subject (i.e. the presence of a unary predicate) and the absence of a passive form. This result is not surprising, given that non-passivable intransitive forms are mostly unergative or stative verbs. Among the best indicators of variable telicity (line VAR in Table 7) are the inanimacy of the subject, the availability of a deverbal adjectival derivation and the fact that the verb is itself formed by derivation. This last feature may simply be an experimental confirmation of the well-known correlation between verbs displaying variable telicity and deadjectival verbs (see widen, deepen, etc., cf. [START_REF] Kennedy | Measure of change: The adjectival core of verbs of variable telicity[END_REF] and references therein). The TELic class (line TEL in Table 7) seems the hardest to predict, since only two features were found to be characteristic for this class only, namely the absence of an indicator for non-agentivity and the presence of an indirect object. We do not have an explanation for the first feature, but the second one is less surprising, given that verbs subcategorizing an indirect object are very often result (causative) verbs (see e.g. [START_REF] Pylkkänen | Introducing Arguments[END_REF]) and that result verbs are very often accomplishment verbs (see e.g. Rappaport Hovav and Levin (1998)).

In Table 8, we look at those features which are found to be characteristic for two aspectual classes and therefore indicative for a greater degree of aspectual flexibility. Among the seven features characteristic for two classes, one is re-lated to agentivitynon-agentive=FALSE (Table 8a) -and two -conc=TRUE and dobj-inanim=TRUE (Table 8b) reflect abstractness/concreteness. The agentivity feature helps to identify two opposite classes (TELic and ATElic). This confirms in some way our hypothesis that agentivity positively correlates with aspectual flexibility (or, in other words, that predicates are aspectually more flexible under their agentive than under their nonagentive reading). On the other hand, features indicative of concreteness are characteristic both of the VARiable and TELic classes. This, also, confirms our hypothesis that verbs are aspectually more flexible under their concrete readings than under their abstract ones. There were no features indicative of both the ATElic and VARiable class. These preliminary findings support the hypothesis that verbs show indeed more aspectual elasticity under their agentive and concrete readings than under their nonagentive and abstract readings. In future work, we plan to refine the features reflecting agentivity and abstractness/concreteness in order to better investigate finer grained aspectual variation across readings of the same verb.

Conclusion

In this work, we focused on the way the aspectual value of French verbs varies across their different readings, as listed and described in a syntactic-semantic lexical resource. We used a machine learning approach based on features extracted from the lexicon. We showed that the extracted information enabled a fairly accurate prediction of the aspectual class. In addition, our results based on feature selection suggest that agentivity and abstractness also contribute to the selection of the correct aspectual value. Our results can also be analyzed as confirming the additional hypothesis that aspectual flexibility is higher when verbs are used under their agentive and/or concrete (literal) readings than under their nonagentive and/or abstract (figurative) readings.

Table 1 :

 1 id example a semantic decomposition sem. primitive sem. class 01 On élargit une route/ La route s'élargit. r/d+qt large become Transformation 02 Cette veste élargit Paul aux épaules. d large become Transformation 03 On élargit ses connaissances. r/d large abs become, figurative Transformation 04 On élargit le débat à la polititque étrangère. f.ire abs VRS directed move, figurative Enter/Leave (a) The four readings illustrated by sample sentences and their semantic description a Literal translations -01: One widens a road. 02: This jacket widens Paul's shoulders. 03: One widens one's knowledge. 04: One extends the debate to foreign policy. LVF entries for élargir

	id schema encoded information
	01 A30	intransitive with adjunct, inanimate subject
	T1308	transitive, human subject, inanimate direct object, instrumental adjunct
	P3008	reflexive, inanimate subject, instrumental adjunct
	02 N1i	intransitive, animate subject, prep. phrase headed by de (of )
	A90	intransitive with adjunct, subject human or thing
	T3900	transitive, inanimate subject, object human or thing
		(b) Syntactic descriptions

Table 2 :

 2 Aspectual distribution of 1199 annotated verb readings

	1	2	3	4	5	6	7	8	1-3	4	5-8
	S-STA ACT-STA S-ACT ACT-ACC W-ACC S-ACC ACC-ACH S-ACH	ATE VAR TEL
	182	67	175	195	172	227	29	152	424	195	580
				(a) 8-value scale				(b) 3-value scale

Table 3

 3

		# sample feature
	syn/sem 28 subcategorisation of an indirect object
	agent	9 subject always inanimate
	abs	10 literal or figurative sense
		(a) Basic feature sets
			#features
		agent-abs	19
		syn/sem-agent	37
		syn/sem-abs	38
		syn/sem-agent-abs 47
		(b) Combinations

presents some more details.

Table 3 :

 3 Feature sets.

Table 4 :

 4 Table4reports the classification results for the different feature sets (Table3) and classifiers. The accuracy for the best performing feature set is marked in italics for each classifier (best accuracy in row) and in bold face for the best performing classifier for each feature set (best accuracy in column). Classification results for LVF entries, using different classification algorithms and features sets. Italics mark best performing feature set for each classifier, bold face best performing classifier for each feature set.

	Algorithm	syn/sem agent	abs agent-abs syn/sem-agent syn/sem-abs syn/sem-agent-abs
	trees.j48	62.72 56.96 57.80	58.72	61.80	63.05	63.64
	rules.jrip	63.47 54.55 54.63	56.88	63.80	63.97	62.80
	lazy.kstar	64.64 57.72 56.63	59.05	62.89	62.30	60.55
	functions.libsvm	63.30 56.13 55.21	56.21	62.72	63.05	62.30
	bayes.naivebayes	60.80 55.38 54.96	55.80	60.22	60.00	59.30
	baseline	48.37 48.37 48.37	48.37	48.37	48.37	48.37

Table 5 :

 5 method defines how characteristic each feature is for each class by introducing class based feature Classification results for LVF entries, using different classification algorithms and reduced features sets. Italics mark best performing feature set for each classifier, bold face best performing classifier for each feature set.

	Algorithm	syn/sem agent	abs agent-abs syn/sem-agent syn/sem-abs syn/sem-agent-abs
	trees.j48	63.55 56.00 58.13	56.55	63.00	63.72	62.47
	rules.jrip	59.80 55.30 55.38	56.21	61.55	60.55	61.13
	lazy.kstar	67.56 55.88 55.38	57.38	67.47	65.56	65.39
	functions.libsvm	59.30 55.80 55.47	55.96	61.13	59.88	61.22
	bayes.naivebayes	64.55 56.55 56.71	58.55	65.80	64.97	66.14
	rules.zeror	48.37 48.37 48.37	48.37	48.37	48.37	48.37

Table 6 :

 6 Classification results based on feature selection with the FMC method. Italics mark the best performing feature set for each classifier, bold face the best performing classifier for each feature set

	Algorithm	syn/sem agent	abs agent-abs syn/sem-agent syn/sem-abs syn/sem-agent-abs
	trees.j48	59.55 56.96 57.72	57.72	61.72	62.97	59.88
	rules.jrip	63.64 55.21 55.13	56.55	64.80	64.72	62.97
	lazy.kstar	60.55 57.72 56.63	59.13	60.80	58.47	57.80
	functions.libsvm	55.00 56.21 55.13	56.46	58.04	58.88	60.13
	bayes.naivebayes	51.80 55.38 54.96	56.13	52.63	52.21	54.13
	rules.zeror	48.37 48.37 48.37	48.37	48.37	48.37	48.37

Their finding that a verb is roughly associated to more than 15 event categories has been further confirmed by a subsequent rater study contrasting event categories and a native speaker judgment task validating the results of this rater study; see[START_REF] Marvel | Microsenses: exploring multi-word event categorization[END_REF].

2 However, in order to check the compatibility of forand inadverbials, the annotator switched to perfective sentences, since those adverbials are not compatible with the present tense in their relevant interpretations.

We chose roughly one classifier from each class in the Weka framework.

For libsvm (the SVM implementation) we used a linear kernel and normalization.

For sets consisting of less than 10 features, we performed an exhaustive search; for the larger feature sets, this turned out to be too time consuming.

Feature Maximization Clustering or Classification (FMC) was initially developed to assess the quality of clusterings.

However, for this feature, there is no preference with regard to the value (both the features non-agentive=TRUE and non-agentive=FALSE were selected).

class features selected ATE dobj=FALSE, passive=FALSE, subj-only=TRUE, literal=FALSE TEL indobj=TRUE, non-agentive=FALSE VAR der-ment=TRUE, is-derivation=TRUE, subj-inanim=TRUE, has-adj-deriv=TRUE, is-derivation=TRUE


Mai Ta

Franck Pigeonneau
email: franck.pigeonneau@saint-gobain.com

Pierre Saramito
email: pierre.saramito@imag.fr

An implicit high order discontinuous Galerkin level set method for two-phase flow problems

Keywords: Level set method, discontinuous Galerkin FEM, high order methods

An implicit high order time (BDF) and polynomial degree discontinuous Galerkin (DG) level set method is presented in this talk. The major advantage of this new approach is an accurate mass conservation during the convection of the level set function, thanks to the implicit method. Numerical experiments are presented for the Zalesak and the Leveque test cases. The convergence rates versus time and space are investigated for both BDF and DG high orders. The capture of the zero level set interface is then improved by using an auto-adaptive mesh procedure. The problem is approximated by using the discontinuous Galerkin method for both the level set function, the velocity and the pressure fields.

Introduction

Let us denote, at any time t 0 by Ω(t) ⊂ R d , d = 2, 3, the bounded moving domain and Γ(t) = ∂Ω(t) its boundary. A level set function φ is defined for any time in a bounded computational domain denoted by Λ ⊂ R d , and containing Ω(t) at any time, such that Γ(t) = {x ∈ Λ; φ(t, x) = 0}. Since every point belonging to the boundary Γ(t) will still continue to belong to it for any time, we have: d dt (φ(t, x)) = 0 ⇐⇒ ∂φ ∂t + u.∇φ = 0 in ]0, +∞[×Λ [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF] where u = ẋ denotes the velocity field. The normal to the boundary Γ(t) writes ν = ∇φ/|∇φ| and then deformations of Γ(t) are only due to the normal component u • ν of the velocity on Γ(t). For a given velocity field u, the problem is to find φ defined in ]0, +∞[×Λ and satisfying [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF] together with an initial condition φ(t = 0) = φ0 where φ0 is given. This is a linear hyperbolic problem. The level set method was introduced in 1988 by Osher and Sethian [START_REF] Osher | Front propaging with curvaturedependent speed: agorithms based on Hamilton-Jacobi formulations[END_REF] (see also [START_REF] Sethian | Level set methods and fast marching methods[END_REF]). Notice that several choices are possible for the function φ: the only requirement being that a fixed isocontour of φ coincides with the front at each time t. A common choice is the signed distance from the front: e.g. Ω(t) is the part where φ(t, .) is negative and φ(t, x) = -dist(Γ(t), x) for all x ∈ Ω(t).

Observe also that, assuming a divergence free velocity field u, the volume of Ω(t) remains constant at any time. This property is not strictly maintained after discretization by most numerical methods suitable for hyperbolic problems and this problem is often referred as a mass loss, as the mass conservation of the fluid is violated. Several cures to this problem have already proposed. In 1999, Sussman and Fatemi [START_REF] Sussman | An efficient, interface preserving level set re-distancing algorithm and its application to interfacial incompressible fluid flow[END_REF] proposed, in the context of a finite difference method, to add a constraint in order to improve the mass conservation. This idea was extended in 2010 by Laadhari et al [START_REF] Laadhari | Improving the mass conservation of the level set method[END_REF] to the finite element context (see also [START_REF] Laadhari | Computing the dynamics of biomembranes by combining conservative level set and adaptive finite element methods[END_REF]). In 2006, Di Pietro et al. [START_REF] Di Pietro | Mass preserving finite element implementations of the level set method[END_REF] and, independently Marchandise et al. [START_REF] Marchandise | A quadrature-free discontinuous Galerkin method for the level set equations[END_REF] proposed to discretize the level set equation (1) by using a discontinuous Galerkin finite element method. Note that the discontinuous Galerkin method applies also to a closely related fluid interface problem, the Cahn-Hilliard equations (see e.g. [START_REF] Pigeonneau | Discontinuous Galerkin finite element method applied to the coupled Stokes/Cahn-Hilliard equations[END_REF]).

In this paper, we revisit this second approach and propose some improvements that dramatically decrease the mass error.

Numerical method

Discontinuous space approximation

Let T h be a triangulation of the computational domain Λ. We introduce the following finite dimensional space:

X h = {ϕ h ∈ L 2 (Λ); ϕ h|K ∈ P k , ∀K ∈ T h }
where k 0 is the polynomial degree. The variational formulation of the semi-discretized problem writes (see e.g. [11, p. 8] or [1, p. 57]):

(F V ) h : find φ h such that φ h (t) ∈ X h , ∀t > 0 and Λ ∂φ h ∂t ϕ h dx + ∂Λ max(0, -u.n) φ h ϕ h ds + K∈T h K u.∇φ h ϕ h dx + S∈S (i) h S [[φ h ]] |u.n| 2 [[ϕ h ]] -u.n { {ϕ h } } ds = 0, ∀ϕ h ∈ X h (2)
together with the initial condition φ h (0) = π h (φ0) where π h denotes the Lagrange interpolation operator on X h . This formulation introduces some notations that are usual in the context of the discontinuous Galerkin method. The second term involves the outer unit normal n on the boundary ∂Λ of the computational domain. The last term involves a sum over S (i) h , the set of internal sides of the mesh T h . Each internal side S ∈ S (i) h has two possible orientations: one is chosen definitively. On this internal side, let n denotes the normal to the oriented side S: as S is an internal side, there exists two elements Kand K+ such that S = ∂K-∩∂K+ and n is the outward unit normal of Kon ∂K-∩ S and the inward unit normal of K+ on ∂K+ ∩ S, as shown on Fig. 1. For all φ h ∈ X h , recall that φ h is in general discontinuous across the internal side S. We define on S the inner value φ - h = φ h|K -of φ h as the restriction φ h|K -of φ h in Kalong ∂K-∩ S. Conversely, we define the outer value φ + h = φ h|K + . We also denote on S the jump

K + S K - n = n -= -n + on S n - n + φ - h φ + h
[[φ h ]] = φ - h -φ + h and the average { {φ h } } = (φ - h + φ + h
The case k = 0, i.e. a piecewise constant approximation, coincides with the popular upwinding finite volume scheme.

Time discretization by an implicit scheme

The final discrete problem is obtained from the semi-discrete one (2) by using a time discretization scheme. For instance, Di Pietro et al. [START_REF] Di Pietro | Mass preserving finite element implementations of the level set method[END_REF] and Marchandise et al. [START_REF] Marchandise | A quadrature-free discontinuous Galerkin method for the level set equations[END_REF] used an explicit Runge-Kutta (RK) method [START_REF] Gottlieb | Strong stabilitypreserving high-order time discretization methods[END_REF]. Here, we propose to use a high order implicit scheme based on the following backward differentiation formula (BDF):

∂φ ∂t = 1 ∆t p k=0 α k φ(t -k∆t) + O (∆t p )
where ∆t > 0 is the time step, p 1 is the order of the scheme and (α k ) 0 k p are the p + 1 coefficients of the formula. When p = 1, the scheme coincides with the usual backward Euler method. For 1 p 2, the scheme is inconditionally stable, for 3 p 6, the scheme is almost inconditionally stable while when p > 6, the scheme is unstable and cannot be used (see e.g. [13, p. 349]). The coefficients when p 6 are given in Table 1. Let us introduce the following bilinear form, defined for all φ, ϕ ∈ X h by

p α0 α1 α2 α3 α4 α5 α6 1 1 -1 2 3/2 -2 1/2 3 11/6 -3 3/2 -1
a(φ, ϕ) = α0 ∆t Λ φ ϕ dx + ∂Λ max(0, -u.n) φ ϕ ds + K∈T h K u.∇φ ϕ dx + S∈S (i) h S [[φ]] |u.n| 2 [[ϕ]] -u.n { {ϕ} } ds
At a time step tn = n∆t, when n p, the fully discrete problem writes (P ) h : find φ

(n) h such that a φ (n) h , ϕ h = - p k=1 α k ∆t Λ φ (n-k) h ϕ h dx
When the time step 0 < n < p, the scheme is started by using the BDF(n) formula while when n = 0 we set φ (0) h = π h (φ0). At each time step, the problem reduces to a linear system in finite dimension, which is solved by a direct method. The simulations are performed with the Rheolef C++ finite element library [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF] developed by one of the authors.

Tests and discussion

The Zalesak rotating disk test

The Zalesak slotted disk in rotation [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF] is a widely used test for comparing the performances of interface capturing methods. In this example, the slotted disk is rotated around the center of the computational domain Λ = [0, H] × [0, H] with H = 4 and a constant angular velocity 0.5, such that u(x, y) = (0.5(y -2), -0.5(x -2)). The disk of radius 0.5 is initially centered at (2, 2.75). The width of the slot is 0.12 and the maximum width of the upper bridge, that connects two parts of the disk, is 0.4. Notices that the slotted disk returns to its initial position after a period T = 4π.

Numerical computations are performed on the time interval [0, T ] with a time step ∆t = T /N where N is the number of time steps. Fig. 2 shows the effect of increasing the polynomial degree k of the discontinuous Galerkin method while using an implicit BDF(k +1) scheme. The mesh is the same for all these computations (5878 elements, 3040 vertices). The comparison between the initial disk position and the final one after one period shows a dramatic improvement when using high order polynomials and schemes: for k 3, the changes are no more perceptible and the disk after one period fits the initial shape. Fig. 3 shows the effect of an uniform mesh refinement when k = 2. Observe that the details of the slotted region, with sharp angles, becomes more precise. In order to quantify the difference between the approximate solution and the exact one, let us introduce the following L 1 error between the two shapes:

φ -φ h = max 0 n N 1 L H (φ(tn)) -H φ (n) h L 1 (Λ)
where H denotes the usual Heaviside function and L is the perimeter size of the initial interface. Fig. 4 shows, for this error measurement, the convergence of the approximate solution to the exact one when the time step tends to zero. Remark that, for each h, the error tends to a constant that is independent of the time step. We are looking for a representation of the error as the sum of two terms, one that depends only upon h and the other only upon ∆t as φ -φ h = O (h α + ∆t p ). Since the solution is regular versus t, the error is expected to depend optimally upon ∆t, asymptotically as ∆t p . Conversely, the shape is poorly regular in space, due to sharp corners, and the convergence properties Varying the mesh refinement with a fixed polynomial degree k = 2 and using an implicit BDF(3) scheme for the Zalesak test. Superposition of the initial disk (black) and after a period (red).

k = 1, ∆t = 4π/500 k = 2, ∆t = 4π/1000 k = 3, ∆t = 4π/2000
will be investigated numerically by looking for the α k + 1 power index. Fig. 5 shows the error versus h for various mesh h and polynomial degrees k. The time step ∆t has been chosen sufficiently small for the error to depend only upon h and k and not upon ∆t (see also Fig. 4). Observe that the approximate solution converges to the exact one with mesh refinement with a power index α ≈ 2 that appears to be independent upon k 1. Table 2 provides the error data for the purpose of comparison: observe that the L 1 error is of about one order of magnitude lower than those in [START_REF] Marchandise | A quadrature-free discontinuous Galerkin method for the level set equations[END_REF]. Let us choose for Ω a full disk, which is regular, instead of a slotted one. Fig. 6 shows the error versus h after one complete revolution of the full disk. It is estimated that α ≈ k/2 + 1. Thus, the power index α is strongly dependent upon k when the shape Ω is regular and increasing k decreases dramatically the error. 10 -1

10 -4 10 -3 10 -2 10 -1 10 -2 10 -1 h/H φ -φ h α = 2 k = 1 2 3 4
10 -2 10 -1 h φ -φ h α = 1 k = 1 2 3
Figure [START_REF] Marchandise | A quadrature-free discontinuous Galerkin method for the level set equations[END_REF] for the error vs mesh refinement h and polynomial degree k.

Notice that the mass of the slotted disk at time tn expresses:

m (n) h = Λ 1 -H φ (n) h dx
Remark that the mass error me = max

0 n N |m(φ(tn)) -m (n) h |
is related to the previous error in L 1 norm me L φ -φ h . The evolution of the relative mass error is reported on Fig. 7. Observe that when k = 1, the mass error on the coarsest mesh is of about 5% of the total mass during the whole period t = 4π. Conversely, on the finest mesh, it is dramatically decreased and remains bounded. As in the previous computation, the time step ∆t has been chosen sufficiently small for the error to depend only upon h and k. While increasing the polynomial degree k, these fluctuations strongly decreases in amplitude and for k = 4 the relative error is bounded on all meshes.

The vortex-in-box Leveque test

The second test, proposed in [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF], considers a disk of radius 0.15 at (0.5, 0.75) in a unit square domain and the velocity field: u(t, x, y) = -sin 2 (πx)sin(2πy)cos(πt/T ) sin 2 (πy)sin(2πx)cos(πt/T ) for 0 t T = 8. The initial disk is compressed and become very stretched: its shape tends to become very thin. The shape reaches its maximum deformation at time t = T /2, and has returned to its initial state at time T. This test is also investigated with different polynomial degrees k and with the BDF(k+1) scheme. Fig. 8 shows the L 1 error versus the mesh refinement h and polynomial degree k. Observe that the convergence rate α increases with k until k = 4 and then appears to be bounded by 3. This is due to the poor regularity of the shape, especially at the half period. From the previous paragraph, we can induce that α = min((k + 1)/2, 3). The behavior of the error on the mass conservation is similar to those of the previous test case and is not showed here.

An auto-adaptive mesh procedure is performed at each time step of the evolution of the interface. In the time interval [t n , t n+1 ], starting from φ n , the mesh adaptation loop iterates five times, generating an adaptive mesh, solving the transport problem and generating successive predictions φ * ,i , 1 i 5, to finally set φ n+1 = φ * ,5 . The anisotropic adaptive mesh is generated accordingly to the metric associated to the Hessian of the criterion χ = δ (φ n ) + δ (φ * ,i ), where δ is the Dirac-delta function. See the Rheolef library documentation [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF] for more about mesh adaptive techniques. Fig. 9 shows some steps of the revolution. The minimum edge length in the mesh was hmin = 1/200.

-0.05 0 0.05 The final solution is reported in Fig. 10 and compared with the corresponding solution obtained with a uniform mesh with h = 1/80.

0 π 2π 3π 4π tn m (n) h -m (0) h m (0) h k = 1, p = 2 h/H = 0.04 h/H = 0.02 h/H = 0.01 -0.05 0 0.05 0 π 2π 3π 4π tn m (n) h -m (0) h m (0) h k = 2, p = 3 h/H = 0.04 h/H = 0.02 h/H = 0.01 -0.05 0 0.05 0 π 2π 3π 4π tn m (n) h -m (0) h m (0) h k = 3, p = 4 h/H = 0.04 h/H = 0.02 h/H = 0.01 -0.05 0 0.05 0 π 2π 3π 4π tn m (n) h -m (0) h m (0) h k = 4, p = 5 h/H = 0.04 h/H = 0.02 h/H = 0.01

Conclusion

In this contribution, a level set transport is investigated with an implicit high order time (BDF) and polynomial degree discontinuous Galerkin finite element method. Using three well-known cases, i.e. a rotation circle, and the Zalesak and Leveque tests, we establish that our method present a nice mesh convergence. For the Zalesak case, the numerical solution converges toward the exact solution in h 2 whatever the polynomial degree mainly due to the non-regularity of the exact solution. For the Leveque test, the mesh convergence is better but stays limited at high polynomial degree due to the sharp interface observed during the advection process. For the rotating circle, up too oour tests, there is no a priori limit for the convergence rate versus h and we can conclude that the limitation of the convergence rate appears when the region Ω has poor regularity.

It is noteworthy that the mass loss becomes very small when the polynomial degree increases. This means that it is not necessary to introduce artificially the mass conservation. This is an important feature for the future work when the transport scheme will be coupled with the Navier-Stokes equations.

Figure 1 :

 1 Figure 1: Notations for the discontinuous Galerkin method.

Figure 4 :

 4 Figure 4: Zalesak test: error vs the time step ∆t for fixed mesh refinement h and polynomial degree k with scheme BDF(k+1).

Figure 2 :Figure 3 :

 23 Figure2: Varying the polynomial degree k with a fixed mesh (5878 elements, 3040 vertices) and using an implicit BDF(k+1) scheme for the Zalesak test. Superposition of the initial disk (black) and after a period (red).

Figure 5 :

 5 Figure 5: Zalesak test: error vs mesh refinement h and polynomial degree k.

Figure 8 :

 8 Figure 8: Leveque test: error vs uniform mesh refinement h and polynomial degree k.

Figure 7 :Figure 9 :Figure 10 :

 7910 Figure 7: Zalesak test: relative mass error vs time for various polynomial degree k and mesh refinement h.

Table 1 :

 1 Coefficients of the BDF(p) schemes, 1 p 6.

				/3
	4	25/12 -4	3	-4/3	1/4
	5 137/60 -5	5	-10/3 5/4 -1/5
	6 147/60 -6 15/2 -20/3 15/4 -6/5 1/6

Table 2 :

 2 6: Test of circle: error vs mesh refinement h and polynomial degree k. Zalesak test: data and comparison with

	order mesh L 1 error
	k	h/H from [7] present
	1	0.04		0.02429
	2	0.04		0.00555
	3	0.04	0.12	0.00446
	4	0.04	0.05	0.00237
	5	0.04	0.04	0.00148
	1	0.02		0.01125
	2	0.02		0.00253
	3	0.02	0.019	0.00096
	4	0.02	0.011	0.00052
	5	0.02	0.008	0.00049
	1	0.01		0.00264
	2	0.01		0.00057
	3	0.01		0.00033
	4	0.01		0.00017
	5	0.01		9.72696e-05

May 22nd -27th 2016, Firenze, Italy ICMF-2016 -9th International Conference on Multiphase Flow May 22nd -27th 2016, Firenze, Italy


Franck Pigeonneau
email: franck.pigeonneau@saint-gobain.com

Pierre Saramito
email: pierre.saramito@imag.fr

Discontinuous Galerkin finite element method applied to the coupled Navier-Stokes/Cahn-Hilliard equations

Keywords: two-phase flows, Cahn-Hilliard equation, phase-field theory, discontinuous Galerkin finite element

Two-phase flows driven by the interfacial dynamics is studied with a phase-field model to tract implicitly interfaces. The phase field obeys the Cahn-Hilliard equation. The fluid dynamics is described with the Stokes equations with an additional source term in the momentum equation taking into account the capillary forces. A discontinuous Galerkin finite element method is used to solve the coupled Stokes/Cahn-Hilliard equations. The Cahn-Hilliard equation is treated as a system of two coupled equations corresponding to the advection-diffusion equation for the phase field and a non-linear elliptic equation for the chemical potential. First, the variational formulation of the Cahn-Hilliard equation is presented. A numerical test is achieved showing the optimal-order in error bounds. Second, the variational formulation in discontinuous Galerkin finite element approach of the Stokes equations is recalled in which the same space of approximation is used for the velocity and the pressure with an adequate stabilization technique. Finally, numerical simulations describing the capillary rising in a tube is presented.

Introduction

The dynamics of two-phase flows driven by the capillary forces can be adressed by various numerical methods which can be sorted as a function of the method to track interfaces. A first class of methods tracks interfaces explicitly using boundary integral method in the limit of vanishing Reynolds number [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF][START_REF] Pigeonneau | Low-Reynolds-Number gravity-driven migration and deformation of bubbles near a free surface[END_REF] or front tracking method [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]. This approach is based on the Lagrangian description of interfaces. Despite of the high level of accuracy which can be reached by these techniques, the topological changes are difficult to take into account. A second class of methods tracks interfaces implicitly with a volume indicator like the "Volume of Fluid" technique initially developed by Hirt and Nichols [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] or by a "level-set" method [START_REF] Sethian | Level set methods and fast marching methods[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF]. They are based on the Eulerian description and can straightforward describe topological changes like coalescence or break-up of bubbles or drops. In the two previous approaches, interfaces are considered as a sharp-interfaces without volume. An alternative method considers interface with a small thickness. This technique employs the concept of "phase-field" introduced initially in statistical physics to describe phase change or spinodal decomposition [START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF]. The development of the phase-field method in fluid mechanics has been summarized by Anderson et al. [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF]. This technique is particularly well-adapted to study local effects like triple line dynamics [START_REF] Jacqmin | Contact-line dynamics of a diffuse fluid interface[END_REF] or more recently in [START_REF] Yue | Sharp-interface limit of the Cahn-Hilliard model for moving contact lines[END_REF] or coalescence of drops [START_REF] Yue | Diffuse-interface simulations of drop coalescence and retraction in viscoelastic fluids[END_REF]. Despite of many contributions, find a numerical method to solve efficiently the coupled equations of phase-field model and Navier-Stokes equations is still a research topic in numerical analysis.

In this present work, and since the phase-field method is devoted to study the fluid dynamics at small scale, the fluid inertia can be neglected. Consequently, we investigate the numerical method to solve the coupled Stokes/Cahn-Hilliard equations with a discontinuous Galerkin finite element technique. We start by the problem statement in section 2 in which the Navier-Stokes/Cahn-Hilliard equations will be recalled. The numerical method to solve only the phase-field will be presented in section 3 with a test to determine the numerical errors. Section 4 is devoted to the numerical method of the coupled Stokes/Cahn-Hilliard equations with an example of numerical simulation of a capillary rising.

Two-phase flows modelling with a phase-field formulation

We consider a domain Ω ⊂ R d (with d = 2 or 3) in which two fluids are present characterized by the density ρ1 and the dynamic viscosity η1 for the fluid 1 and by the density ρ2 and the dynamic viscosity η2 for the fluid 2. The fluid phase at any material point with a position x and at each time t is described by an "order parameter" φ(x, t). By convention, the fluid 1 is given by the order parameter φ = 1 and the fluid 2 by φ = -1. The function φ can be seen as a volume fraction.

The phase field method considers that the shift between the two phases occurs over a thin layer equal to ζ corresponding to the "diffuse interface" thickness. Under the actions of external forces, boundary conditions and interaction between the two phases, the media changes in space and time requiring balance equations. Moreover, the dynamics between the two phases has to be considered by writing an equation on the phase field φ.

The phase-field model applied for binary fluid is designated as the "model H" by Hohenberg and Halperin [START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF]. The details of the derivation of the Navier-Stokes/Cahn-Hilliard equations are not provided in this article. For more information, reader can read references [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF]. In the following, the Navier-Stokes/Cahn-Hilliard equations are written under dimensionless form with a characteristic length L, velocity U , viscosity η1 and density ρ1 which will be provided in the numerical examples presented in section 4. In the following, we assume that the fluid is incompressible. The gravity force is taken into account with a direction given by a unit vector g normalized by the magnitude of the gravity acceleration, g. The coupled system of equations is given by

div u = 0, (1)
Re ρ(φ) Du Dt = -∇P + div [2η(φ)D(u)] + Bo Ca ρ(φ)g + 3 2 √ 2 Ca Cn µ∇φ, (2)
Dφ Dt = 1 Pe ∇ 2 µ(φ), (3)
µ(φ) = φ(φ 2 -1) -Cn 2 ∇ 2 φ, (4)
in which the two first equations are the Navier-Stokes equations with D the rate-of-strain tensor and the two lasts are the Cahn-Hilliard equation. The coupling between the phase-field is achieved by the capillary force given by the last term of the right-hand side of the momentum equation, eq. ( 2) and also by the chance of the density and viscosity as function of φ. The Cahn-Hilliard equation has been written in two parts in order to have two equations at the second order useful for the numerical methods which are presented below. The last equation defines the chemical potential in which the first term is the derivative respect to φ of the double-well potential admitting two local minima in φ = ±1 and the second term takes into account the energy due to interface. At the equilibrium, the chemical potential must be equal to zero. Seven dimensionless numbers arise in the problem for which the first six numbers are defined by

Re = ρ1U L η1 , Bo = ρ1gL 2 σ , Ca = η1U σ , (5)
Cn = ζ L , ρ = ρ2 ρ1 , η = η2 η1 . (6)
The first one is the so-called Reynolds number which will be assumed much lesser than one in the following in order to neglect the fluid inertia. The Bond number Bo measures the ratio of gravity to surface tension forces characterized by the surface tension σ while the capillary number Ca compares the viscosity effect to the surface tension one. The Cahn number Cn is the ratio of the diffuse-interface thickness ζ to the length scale. The seventh dimensionless number Pe given without definition is the Péclet number comparing the diffusion time scale of the chemical potential to the convective time scale which is always greater than one. The normalized density and dynamic viscosity are given by an arithmetic average:

ρ(φ) = 1 + ρ 2 + 1 - ρ 2 φ, (7)
η(φ) = 1 + η 2 + 1 - η 2 φ. (8)
Initial and boundary conditions will be provided in the following when the numerical method will be presented.

Numerical resolution of Cahn-Hilliard equation

We start by the presentation of the numerical method to solve the Cahn-Hilliard equation with a discontinuous Galerkin finite element method. Clearly, the Cahn-Hilliard equation is a nonlinear fourth-order partial differential equation needing a high level of regularity. The application of discontinuous Galerkin method to solve the phase-field model has been initiated by Feng and Karakashian [START_REF] Feng | Fully discrete dynamic mesh discontinuous galerkin methods for the cahn-hilliard equation of phase transition[END_REF] among others. They established an optimalorder error bound by solving the fourth-order equation when the polynomial degree is greater or equal to two. Here, according the previous work of Kay et al. [START_REF] Kay | Discontinuous Galerkin finite element approximation of the Cahn-Hilliard equation with convection[END_REF], the problem is set as a system of two second-order elliptic equations. In this section, the solenoidal velocity field u is assumed known.

Discontinuous Galerkin finite element formulation

The continuous problem is Problem 1 (Cahn-Hilliard problem) find φ(x, t) and µ(x, t)

defined in Ω × [0, T ] such as ∂φ ∂t + u • ∇φ - 1 Pe ∇ 2 µ(φ) = 0, in Ω × [0, T ], (9)
-µ(φ) + φ(φ 2 -1) -Cn 2 ∇ 2 φ = 0, in Ω × [0, T ], (10) φ
(x, t = 0) = φ0(x), in Ω, ( 11
) ∂φ ∂n = f (φ), on ∂Ω, (12)
∂µ ∂n = 0, on ∂Ω.

Moreover, the normal component of the velocity u•n is equal to zero on ∂Ω for n the exterior unit normal on ∂Ω. The boundary condition on φ, eq. ( 12), has been taken in general form because in applications, the interaction with a wall for instance involves a condition on normal derivative of φ to impose a certain contact angle [START_REF] Jacqmin | Contact-line dynamics of a diffuse fluid interface[END_REF][START_REF] Yue | Sharp-interface limit of the Cahn-Hilliard model for moving contact lines[END_REF].

To solve numerically this problem, a discontinuous Galerkin method initially developed by Kay et al. [START_REF] Kay | Discontinuous Galerkin finite element approximation of the Cahn-Hilliard equation with convection[END_REF] has been chosen. The main property of the formulation is the optimal-order in error bounds as it will be verified below for a polynomial degree greater or equal to one.

Let T h be a triangulation of the domain Ω formed by finite elements K with meshsize h = diam(K) such as T h = {K} and Ω = K∈T h K. Furthermore, we consider a subset of faces F formed by interfaces between two distint mesh elements K+ and Ksuch as F = ∂K-∩ ∂K+ and boundary faces given by F = ∂K ∩ ∂Ω. For internal face, orientation has to be precised. The unit normal nF is outward of Kand inward of K+. If F i h gathers all interfaces and F b h all boundary faces, the set of faces is

F h = F i h ∪ F b h .
Finally, the local length scale at the face F noted hF is defined according to Di Pietro and Ern [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF], chap. 4 page 125.

In the below formulation, some usual notations in the context of the discontinuous Galerkin method have to be precised. First, since for all φ h ∈ X h is discontinuous, the broken gradient ∇ h is defined by [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]

(∇ h φ h ) |K = ∇ φ h|K , ∀ K ∈ T h . (14)
Moreover, at each interior face F , let define φ - h = φ h|K -the inner value and φ + h = φ h|K + the outer value. We define the jump by

[[φ h ]] = φ - h -φ + h and the average { {φ h } } = (φ - h + φ + h )/2.
The discontinuous finite element space is defined by

X h = {v h ∈ L 2 (Ω); v h|K ∈ P k , ∀ K ∈ T h }, (15)
with k the polynomial degree. This space belongs to the broken Sobolev space

H 1 (T h ) = {v ∈ L 2 (Ω); v |K ∈ H 1 (K), ∀ K ∈ T h }. ( 16
)
The variational formulation of the Cahn-Hilliard problem takes the following form

Problem 2 Find φ h , µ h ∈ X h such that Pe -1 a(w h , µ h ) + b(w h , φ h ) = 0, for w h ∈ X h , (17)
-(q h , µ h ) + (q h , φ 3 h -φ h ) + Cn 2 a(q h , φ h ) = l(q h , φ h ), for q h ∈ X h , ( 18
) in which a(w h , µ h ) = Ω ∇ h w h • ∇ h φ h dV - F ∈F i h F [{ {∇µ h • nF } }[[w h ]] + { {∇w h • nF } }[[µ h ]] -cF [[µ h ]][[w h ]]] dS, (19) b
(w h , φ h ) = ∂φ h ∂t , w h + Ω ∇ h φ h • u h w h dV + F ∈F i h F 1 2 |u h • nF |[[φ h ]][[w h ]]- u h • nF [[φ h ]]{ {w h } }) dS, (20)
l(q h , φ h ) = Cn 2 F ∈F b h F f h (φ h )q h dS, for q h ∈ X h . ( 21
)
The details of the derivation of this variational formulation can be found in [START_REF] Kay | Discontinuous Galerkin finite element approximation of the Cahn-Hilliard equation with convection[END_REF]. The penalty parameter cF is equal to the product β F with β = (k + 1)(k + d)/d (d space dimension) and F by [START_REF] Shahbazi | An explicit expression for the penalty parameter of the interior penalty method[END_REF] F = meas(∂K) meas(K)

for F ∈ K ∩ ∂Ω, max meas(∂K + meas(K + ) , meas(∂K - meas(K -) for F ∈ K-∩ K+. (22)
The convective term of the phase-field has been written according to the development of Di Pietro and Ern [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF] with an upwinding flux approximation. The main feature of the discontinuous Galerkin method is to have an accurate mass conservation applied to convective equation. This behavior has been controlled carefully in the paper [START_REF] Ta | An implicit high order discontinuous galerkin level set method for two-phase flow problems[END_REF] presented in this conference and is not addressed in the present contribution.

For the time discretization, a Backward Differential Formula (BDF) at order p ≤ 6 is used. If ∆t is the time step, the temporal derivative at the time t = n∆t of the φ n h at the order p is given by

∂φ n h ∂t = 1 ∆t p l=0 α pl φ n-l h + O(∆t p ), (23)
for which coefficients α pl can be found in the book of Süli and Mayers [START_REF] Süli | An Introduction to Numerical Analysis[END_REF] (chap. 12, page 349). This time derivative obtained by the previous formula is introduced in b h (w h , φ n h ) leading to a source term obtained with a combination of φ n-l h with 1 ≤ l ≤ p following equation (23). Moreover, the non-linearities of the Cahn-Hilliard problem are solved using a Newton algorithm.

Numerical tests

The previous problem has been implemented in the Rheolef C++ finite element library [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF]. To test the numerical solver, the numerical experiment proposed by Kay et al. [START_REF] Kay | Discontinuous Galerkin finite element approximation of the Cahn-Hilliard equation with convection[END_REF] is used for which the Cahn-Hilliard equation is solved with a velocity field given by

u = f (r)(y, -x) T , ∀(x, y) ∈ [-1; 1] 2 , (24) with f (r) = 1 + tanh[β(1 -3r)] 2
, and r = x 2 + y 2 , (25) in a domain Ω = [-1; 1] 2 . An exact solution of the phase-field given by

φe = t cos(πx) cos(πy), ∀(x, y) ∈ [-1; 1] 2 and t ∈]0; 1], (26)
is imposed by adding an adequate source term in eq. ( 9). The range of time is taken in [0; 1]. The domain has been discretized with a regular triangular elements with 16, 32 and 64 elements over each Cartesian coordinate.

We perform the numerical experiment with Cahn number equal to 10 -1 and for a Péclet number equal to 50. The temporal numerical scheme at the second order (BDF-2) has been used with a time step equal to 10 -2 . Three polynomial degrees have been tested with k = 1, 2 and 3. Errors between the numerical solution and the exact solution have been computed with L 2 and H 1 norms defined as follows

||v|| 2 L 2 ,h = Ω v 2 dV, (27)
||v|| 2 H 1 ,h = Ω ∇ h v • ∇ h vdV + F ∈F h F F [[v]] 2 dS. (28)
Figures 1 and 2 depict errors integrated over the time range [0, 1] as a function of the mesh size both for the phase field and the chemical potential. For each polynomial degree, errors behave as O(h n ). The exponent n has been provided in Figures 1 and2 for each polynomial degree. As expected, errors computed in L 2 norm behave approximately as O(h k+1 ) both for the phase field φ and for the chemical potential and as O(h k ) in H 1 norm showing that the numerical implementation is optimal-order in error bounds.

||φ h -π h (φ e )|| L 2 ,h dt k = 1 k = 2 k = 3 n = 1.9 n = 3.17
||µ h -π h (µ e )|| L 2 ,h dt k = 1 k = 2 k = 3 n = 1.93 n = 2.78 n = 4.
||µ h -π h (µ e )|| H 1 ,h dt k = 1 k = 2 k = 3 n = 1.17 n = 2.11 n = 3.05 Figure 2: Errors 1 0 ||µ h -h (µe)|| L 2 ,h dt and 1 0 h π h (µe)|| H 1 ,h
dt as a function of mesh size for the three polynomial degrees.

Numerical resolution of Cahn-Hilliard and Stokes equations

We now turn on the numerical resolution of Cahn-Hilliard and Navier-Stokes equations. Here, the Reynolds number is assumed much lesser than one in order to neglect the inertia of the fluid. In this limit, the Navier-Stokes equations are reduced to the Stokes equations. The continuous problem is similar to the set of equations [START_REF] Abramowitz | Handbook of mathematical functions[END_REF][START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF][START_REF] Cockburn | Local discontinuous galerkin methods for the stokes system[END_REF] given in § 2 in which the inertia term in the momentum equation is removed. Initial and boundary conditions for φ and µ can be written in general form as in equations [START_REF] Kay | Discontinuous Galerkin finite element approximation of the Cahn-Hilliard equation with convection[END_REF][START_REF] Pigeonneau | Low-Reynolds-Number gravity-driven migration and deformation of bubbles near a free surface[END_REF][START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF]. In order to take into account the boundary conditions for which velocity or stress are imposed, the frontier of Ω is shared between ∂ΩD for which

u = uD, (29)
and ∂ΩN where A discontinuous Galerkin finite element method is also used to solve the Stokes equations. As it will be detailed below velocity and pressure are approximated at the same polynomial degree m. So, the discontinuous finite element spaces are defined by

σ • n = tN , (30)
V h = {u h ∈ L 2 (Ω) d ; u |K ∈ Pm, ∀K ∈ T h }, (31)
Q h = {q h ∈ L 2 (Ω); q |K ∈ Pm, ∀K ∈ T h }. ( 32
)
To solve the problem as a function of time, a BDF-p scheme is used to solve the phase-field equation in which the velocity field is extrapolated at the same order that the BDF scheme using the forward difference formula [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]

u * ,n h = p l=1 p l (-1) l-1 u n-l h + O(∆t p ). (33)
At the time step n, φ n h and µ n h are determined by solving the problem 2 in which the convective equation is solved using u * ,n h . Once the phase-field problem is solved, the Stokes equations have to be solved taking into account the capillary source term given by the last term of the right-hand side of (2).

Discontinuous Galerkin finite element formulation of Stokes problem

Since the Stokes equations are in quasi-steady state, the reference of the time is removed in the following in order to simplify the notation. The variational formulation follows the method initially introduced by Cockburn et al. [START_REF] Cockburn | Local discontinuous galerkin methods for the stokes system[END_REF] which has been studied theoretically by Di Pietro and Ern [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]. The heterogeneity of the viscosity needs to generalize the formulation according to the previous developments achieved in heterogeneous diffusion [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF] (chap. 4). The discrete variational formulation writes:

Problem 3 (Stokes problem) find u h ∈ V h and P h ∈ Q h such that α(u h , v h ) + β(v h , P h ) = λ(v), ∀v h ∈ V h , (34) β
(u h , q h ) -γ(p h , q h ) = 0, ∀q h ∈ Q h , (35) with α
(u h , v h ) = Ω 2ηD h (u h ) : D h (v h )dV + F ∈F i h ∪F b,D h F [cF ηF [[u]] • [[v]]- { {2ηD h (u h ) • nF } } ω [[v h ]] - { {2ηD h (v h ) • nF } } ω [[u h ]] dS, (36)
β(u h , q h ) = - Ω q h ∇ h • u h dV + F ∈F i h F [[u]] • nF { {q h } }dS, (37)
λ(v) = Ω Bo Ca ρ h g + 3 2 √ 2 Ca Cn µ h ∇ h φ h ] • v h dV + F ∈F b,D h F u h,D • [cF v h -2ηD h (v h ) • n] dS + F ∈F b,N h F t h,N • v h dS, (38)
γ(p h , q h ) = F ∈F i h F hF [[p h ]][[q h ]]dS. (39)
The viscosity at the face F and the weighted average { {•} } ω are defined by [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]

ηF = 2η + h η - h η + h + η - h , (40)
{ {v h } } ω = η + h v - h + η - h v + h η + h + η - h . (41)
Clearly, when the viscosity is constant, the usual arithmetic average is found.

Numerical tests: capillary rising

When a liquid is introduced in a capillary tube with a diameter D, the liquid rises due to the wetting of the liquid on the tube wall. At the equilibrium, the rising height depends on the wetting angle θ, the surface tension σ, the liquid density ρ1, the gravity and the tube diameter. According to de Gennes et al. [START_REF] De Gennes | Gouttes, bulles, perles et ondes[END_REF], the height over which the liquid rises can be written as follows

hcap D = 4 cos θ Bo , (42)
in which the Bond number is defined with D as a length scale. When θ < π/2 corresponding to "wetting" condition, hcap is positif while if θ < π/2 ("non-wetting" condition), hcap is negatif meaning that the liquid goes down occurring for a liquid metal like mercury, for instance. The dynamics of rising can be studied with the Stokes/Cahn-Hilliard equations for which the wetting condition can be introduced easily. We perform with Rheolef library a numerical simulation in an axisymmetric geometry depicted in Figure 3 corresponding to the one half of the tube. The problem is normalized by a length scale equal to the tube diameter. The velocity scale is taken by writing the balance between gravity and viscous forces meaning that the capillary and Bond numbers become equivalent. In the limit of vanishing Reynolds number, five dimensionless numbers have to be considered. For the phase-field modelling, the physical length scale of interface is around few nanometers meaning that the ratio of this interface thickness and the tube diameter is too small. Fortunately, the "sharp-interface" behavior is obtained for a larger Cahn number [START_REF] Yue | Sharp-interface limit of the Cahn-Hilliard model for moving contact lines[END_REF]. Consequently, in the following, the Cahn number is taken equal to 10 -2 . The Péclet number has to be taken sufficiently large to limit the diffusion. Here, Pe is set at 10 2 . Ratios ρ and η are taken respectively at 10 On ∂ΩD, the wetting condition has been written according to the previous works of Jacqmin [START_REF] Jacqmin | Contact-line dynamics of a diffuse fluid interface[END_REF] and Yue et al. [START_REF] Yue | Sharp-interface limit of the Cahn-Hilliard model for moving contact lines[END_REF]. Moreover, the non-slip condition is applied on ∂ΩD. On the top and the bottom of the tube, we apply a pressure with stress free condition.

At the bottom, the pressure arises from a simple static equilibrium used as initial condition in order to start from the static condition. Initially, fluid 1 is located below the plane z = 0 in such of way that φ is given by the exact solution of the Cahn-Hilliard equation in one-dimension:

φ0(z, r) = -tanh z √ 2 Cn . ( 46
)
Starting with a plane interface, we impose a contact angle equal to θ = 17π/36. The numerical simulations have been achieved with BDF-2 and a time step 10 -2 . Due to the contact angle lesser than π/2, the fluid interface curves close to the wall leading to a pressure jump. Figure 4 depicts the pressure field for three times equal to 2, 20 and 60. Clearly, at the beginning, the pressure decreases close to the wall due to the curvature and Laplace pressure. Consequently, the liquid goes up. At large time, the pressure obeys at the hydrostatic solution leading to the immobilization of the liquid column. Pressure P is plotted as a function of z on the symmetric axis when t = 100 in Figure 7. The hydrostatic solution in the liquid is very well found. Pressure in the second phase is quasi-uniform since the density of fluid 2 is very small. Moreover, at the interface, the pressure jump due to the interface curvature is around 0.3 while a simple hydrostatic computation gives a value equal to 0.34 on the symmetric axis. Nevertheless, remark that the numerical solution does not reach the asymptotic static state.

Conclusion

In this contribution, we develop a numerical method to solve the coupled Stokes/Cahn-Hilliard equations with a discontinuous Galerkin finite element method. The scheme on the Cahn-Hilliard is optimal-order in error bounds and this fact has been verified numerically. A numerical method for Stokes/Cahn-Hilliard equations has been also developed for which the Stokes solver with heterogeneous viscosity has been established.

In order to study the dynamics of the triple line physics, we have solved the problem of the capillary rising in a circular tube. The dynamics is very well reproduced and the asymptotic behavior of the static pressure equilibrium is captured.

This first numerical test prove the relevance of this kind of model to describe the dynamics of triple line accurately.

n = 4

 4

Figure 1 : Errors 1 0 2 ,h dt and 1 0

 1121 Figure 1: Errors10 ||φ hh (φe)|| L 2 ,h dt and 1 0 ||φ h -π h (φe)|| H 1 ,h dt as a function of mesh size for the three polynomial degrees.

18

 18

  such as ∂ΩD ∩ ∂ΩN = ∅ and ∂Ω = ∂ΩD ∪ ∂ΩN . The stress tensor is given by -P 1 + 2η(φ)D(u) with 1 the unit tensor. Boundary faces, F b h , are now shared in two sub-sets F b,D h and F b,N h corresponding to ∂ΩD and ∂ΩN respectively.

Figure 3 :

 3 Figure 3: Geometry of a liquid rising in a tube with a radius equal to 1/2.

  -3 and 10 -2 which are the typical values for water (fluid 1) and air (fluid 2). The boundary conditions on ∂ΩD and ∂ΩN depicted in Figure 3 are the following u ∀x ∈ ∂ΩN, bottom. (45)

Figure 4 :

 4 Figure 4: Pressure field in the tube at time t = 2, 20 and 60 for a contact angle equal to 17π/36.

Figure 5

 5 Figure5presents the superposition of the phase-field and velocity field at t = 2, 20 and 60. At the beginning of the rising the velocity is relatively important and decreases when the liquid goes up. When t = 60, the magnitude of the velocity is one order of magnitude lesser than the value observed at t = 2.

Figure 5 :

 5 Figure 5: and field in the tube at t = 2, 20 and 60 for a contact angle equal to 17π/36.

Figure 6 :Figure 7 :

 67 Figure 6: Position of the triple line, zTL as a function of time for a contact angle equal to 17π/36.

May 22nd -27th 2016, Firenze, Italy ICMF-2016 -9th International Conference on Multiphase Flow May 22nd -27th 2016, Firenze, Italy