Federica Adobbati

Carlo Ferigato

Stefano Gandelli

Adrián Puerto Aubel

Two Operations for Stable Structures of Elementary Regions

The set of regions of a condition/event transition system represents all the possible local states of a net system the behaviour of which is specified by the transition system. This set can be endowed with a structure, so as to form an orthomodular partial order. Given such a structure, one can then define another condition/event transition system. We study cases in which this second transition system has the same collection of regions as the first one. When it is so, the structure of regions is called stable. We propose, to this aim, a composition operation, and a refinement operation for stable orthomodular partial orders, the results of which are stable.

Introduction

This work studies the interrelations between local states, locally observable properties, and events of distributed systems. To this aim, its framework is set in the relation between elementary Petri net systems, and the labelled transition systems expressing their behaviour, their case graphs. Indeed, labelled transition systems are commonly used models for verification of properties of the specified system [START_REF] Clarke | Model checking[END_REF].

We focus on the particular case in which local states act like Boolean variables, forbidding executions by carrying true or false values. Thus, the framework is narrowed down to elementary and condition/event models, and their close relation to the theory of regions. In [START_REF] Ehrenfeucht | Partial (set) 2-structures. part I: basic notions and the representation problem[END_REF][START_REF] Ehrenfeucht | Partial (set) 2-structures. part II: state spaces of concurrent systems[END_REF], regions were shown to be the key to solving the synthesis problem: Given a labelled transition system, determine whether it is the case graph of some elementary net system. The extension of a local state of a net system in its case graph is the set of global states at which it holds the value true. In a labelled transition system, a region is a set of global states having a consistent incidence with respect to the actions, so as to be the extension of a local state. Given an elementary, or condition/event transition system, a net system can be constructed with regions as local states, such that its case graph will be precisely the specified transition system. It was further shown [START_REF] Desel | The synthesis problem of Petri nets[END_REF] that given an elementary (condition/event) system, it suffices to consider the sub-collection of regions that is minimal with 1 respect to set inclusion [START_REF] Bernardinello | Synthesis of net systems[END_REF]. The remaining regions correspond to local states which are redundant in the net system.

Regions, as subsets of global states of a system, carry an algebraic structure [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF]. They can be ordered by set inclusion, in a relation that expresses implication. Negation, disjunction, and conjunction can be defined as set operations, only partially in the latter two cases. The resulting structure was shown to be a prime and coherent orthomodular partial order [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF], omp for short.

Orthomodular partial orders typically represent, in an algebraic fashion, systems about which the acquisition of information depends on the point of view, or experiment. They were first introduced by Birkhoff and Von Neumann [START_REF] Birkhoff | The logic of quantum mechanics[END_REF] in order to formalise the logic of properties of quantum systems. Foulis and Randall [START_REF] Foulis | Operational Statistics. I. Basic Concepts[END_REF][START_REF] Randall | Operational statistics. II. Manuals of operations and their logics[END_REF] extended their approach to provide a generalized model.

In the view that omp's are a suitable specification for the observable properties of a system, we are concerned with a second synthesis problem: Given an omp determine whether it is the regional structure of some elementary transition system. This problem was first posed in [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF], and a solution was pointed at. When an omp is prime, one can make use of a representation theorem to find a set of states such that the elements of the omp are subsets of it. By considering a complete graph whose vertices are these states, one can label each edge so that all elements of the omp can be found as regions of this new transition system. When the structure of regions of this new system is isomorphic to the specified omp, then the latter is called stable.

Not every omp is stable, and the full characterization of the class of stable omp's is an open problem which represents a long term objective. It is conjectured that every regional omp, the structured collection of regions of some elementary transition system, is stable. The present work includes itself in a series of papers [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF][START_REF] Bernardinello | Synthesis of transition systems from quantum logics[END_REF][START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF][START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF] the scope of which is to study this stability of omp's, in an attempt to prove the conjecture.

The present work extends [START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF], in which some classes of stable omp's were presented. Here, composition and refinement operations are defined, which preserve stability, thus proving the stability of regional omp's for a wide range of systems.

Background

This work implicitly deals with elementary net systems (ENS for short) [START_REF] Nielsen | Elementary transition systems[END_REF]. Even though these models will not be formalised here, they underlie the main ideas in the presented results. ENS are a class of pure, and simple Petri net systems with a particular firing rule which ensures that the system remains 1-safe in case of contact. Conditions of ENS can only be in one of two states, marked or un-marked, and as such they can be seen as Boolean variables. Moreover, out of Petri's extensionality principle, two events having the same incidence on the set of conditions, and viceversa for conditions in relation to events, should be considered the same [START_REF] Petri | General net theory. computing system design[END_REF].

Condition/Event net systems (CENS for short) are structurally identical to ENS with the addition of a backward firing rule which provides symmetry to the reachability relation amongst markings; reachability becomes an equivalence relation, and reachability classes substitute the initial marking required for ENS. To the left, a CENS. To the right, its sequential case graph. The state q 1 corresponds to the marking {b 1 , b 4 }, at which e 1 and e 2 are concurrently enabled. To the right, the subset of states {q 1 , q 2 } is in evidence; it corresponds to the extension of b 1 , namely the set of markings containing b 1 .

Transition Systems and Regions

The behaviour of CENS can be expressed by means of a labelled transition system. In such a model, states correspond to the reachable markings of the net system. A marking enabling an event will correspond to a state from which an arrow, labelled with this event, points to the state which is reached when the event is fired. Such a labelled transition system is the sequential case graph of the underlying net system. Definition 1 (Transition System). A transition system is a structure A = (Q, E, T ), where Q is a set of states, E is a set of events and T ⊆ Q × E × Q is a set of transitions carrying labels in E, such that:

1. the underlying graph of the transition system is connected; 2. ∀(q 1 , e, q 2 ) ∈ T q 1 = q 2 ; 3. ∀(q, e 1 , q 1 )(q, e 2 , q 2 ) ∈ T q 1 = q 2 ⇒ e 1 = e 2 ; 4. ∀e ∈ E ∃ (q 1 , e, q 2 ) ∈ T .

All the structures we consider are finite. A fundamental notion used in this contribution is the one of region. Regions were introduced in [START_REF] Ehrenfeucht | Partial (set) 2-structures. part I: basic notions and the representation problem[END_REF][START_REF] Ehrenfeucht | Partial (set) 2-structures. part II: state spaces of concurrent systems[END_REF] as the fundamental tool for solving the so called synthesis problem for Petri Nets.

Definition 2 (Region). A region of a transition system

A = (Q, E, T ) is a subset r of Q such that: ∀e ∈ E, ∀(q 1 ,
e, q 2 ), (q 3 , e, q 4 ) ∈ T :

1. (q 1 ∈ r and q 2 / ∈ r) implies (q 3 ∈ r and q 4 / ∈ r) and 2. (q 1 / ∈ r and q 2 ∈ r) implies (q 3 / ∈ r and q 4 ∈ r).

Regions are subsets of states which have a consistent orientation with the labelling of transitions. As such, they can be assigned an incidence with respect to each event, and interpreted as Boolean conditions which hold at the states composing them.

Example 1. With reference to Figure 1, the extension of b 1 is a region. As a subset of states, all occurrences of e 1 exit it and e 4 enters it. All occurrences of the remaining events do not cross the border of b 1 .

The set of regions of a transition system A will be denoted by R(A) and, given a state q of A, R q will denote the subset of the regions in R(A) which contains the state q. Given a transition system A and its set of regions R(A), the pre-set and post-set operators can be defined: Definition 3 (Pre-and Post-sets). • ( . ), the pre-set operator, and ( . ) • , the post-set operator, are defined on the events E and on the regions R(A) of the transition system A as follows. Let r be a region in R(A)

1. • r = {e ∈ E | ∃ (q 1 ,
e, q 2 ) ∈ T such that q 1 / ∈ r and q 2 ∈ r};

2. r • = {e ∈ E | ∃ (q 1 ,
e, q 2 ) ∈ T such that q 1 ∈ r and q 2 ∈ r};

3. • e = {r ∈ R(A) | e ∈ r • }; 4. e • = {r ∈ R(A) | e ∈ • r}.
We are concerned with a specific class of transition systems, the condition/ event transition systems (cets for short). cets can be defined as those labelled transition systems which are the sequential case graph of some CENS. An alternative definition, in terms of regions and the separation axioms, follows from the results of [START_REF] Ehrenfeucht | Partial (set) 2-structures. part I: basic notions and the representation problem[END_REF][START_REF] Ehrenfeucht | Partial (set) 2-structures. part II: state spaces of concurrent systems[END_REF].

Definition 4 (CETS -Condition/Event Transition System).

A cets is a transition system in which the following conditions hold:

1. ∀ q 1 , q 2 ∈ Q R q1 = R q2 ⇒ q 1 = q 2 ; 2. ∀ q 1 ∈ Q ∀ e ∈ E • e ⊆ R q1 ⇒ ∃ q 2 ∈ Q such that (q 1 , e, q 2 ) ∈ T ; 3. ∀ q 1 ∈ Q ∀ e ∈ E e • ⊆ R q1 ⇒ ∃ q 2 ∈ Q such that (q 2 , e, q 1 ) ∈ T .
This definition requires of a transition system, that its set of regions is sufficient to fully determine the incidence of each event with respect to each state. Under these conditions, a CENS can be constructed, which has the set of regions for conditions, and the transition labels as events, and such that its case graph is isomorphic to the given transition system. This tight relation between CENS and cets was thoroughly formalised in [START_REF] Nielsen | Elementary transition systems[END_REF], in the equivalent framework of elementary systems.

The following basic properties of the regions R(A) of a cets A are proven in [START_REF] Bernardinello | Synthesis of net systems[END_REF]:

1. ∅, Q ∈ R(A); 2. r ∈ R(A) ⇒ Q \ r ∈ R(A); 3. r 1 , r 2 ∈ R(A) ⇒ (r 1 ∩ r 2 ∈ R(A) ⇔ r 1 ∪ r 2 ∈ R(A)).

Orthomodular Partially Ordered Sets

Orthomodular partially ordered sets, indicated as omp in what follows, are well known algebraic structures in the literature on Quantum Logics. The following definition can be found as the finite version of Definition 1.1.1 in [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF]. Note that ∧ and ∨ are the usual meet and join operations on a partial order, denoting, respectively, the greatest lower bound, and the least upper bound.

Definition 5 (OMP -Orthomodular Partially Ordered Set). An orthomodular partially ordered set L, ≤, (.) , 0, 1 is a set L endowed with a partial order ≤ and a unary operation ( . ) (called orthocomplement), such that the following conditions are satisfied:

1. L has a least and a greatest element (respectively 0 and 1) and 0 = 1;

2. ∀ x, y ∈ L x ≤ y ⇒ y ≤ x ; 3. ∀ x ∈ L (x ) = x;
4. for any finite sequence {x i | i ∈ I} of elements of L such that i = j ⇒

x i ≤ x j , then i∈I x i exists in L;

5. ∀ x, y ∈ L x ≤ y ⇒ y = x ∨ (x ∧ y).
This latter condition is sometimes referred to as orthomodular law.

We will often denote L, ≤, (.) , 0, 1 simply by L, and assume that an omp

L i has underlying structure L i , ≤ i , (.) , 0 i , 1 i .
In general, ∧ and ∨ operations can be undefined for some pairs of elements of L. Two elements x and y in L are said to be orthogonal, noted by x ⊥ y, when x ≤ y . As a consequence of their basic properties, as listed at the end of Section 2.1 above, the set of regions of a cets is an omp. 0 and 1 are, respectively, ∅ and Q, the order is given by set containment and orthocomplement by set complement; moreover, two regions are orthogonal whenever their intersection is empty.

Since all the structures we consider are finite, an omp can be equivalently specified by its collection of atoms A(L) := {x ∈ L \ {0} | ∀y ∈ L : (y ≤ x) → (y = x or y = 0)} together with their orthogonality relation. Each element of L can then be retrieved as the join of a collection of pairwise orthogonal atoms. The reader is referred to [START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF] for the details, and to [START_REF] Foulis | Operational Statistics. I. Basic Concepts[END_REF][START_REF] Randall | Operational statistics. II. Manuals of operations and their logics[END_REF] for a full characterisation of this representation of omp's. A subset M of an omp L with the order relation and the orthocomplement operation inherited from A morphism f :

L is a sub-omp of L in case M is an omp ([16], definition 1.2.2). Definition 6 (OMP-morphism). ([16], finite case of definition 1.2.7) Let L 1 and L 2 be omp. A mapping f : L 1 → L 2 is a morphism of omp's if the following conditions are satisfied: 1. f (0) = 0; 2. ∀x ∈ L 1 f (x ) = f (x) ; 3. for any finite sequence {x i | i ∈ I} of mutually orthogonal elements in L 1 , f ( i∈I x i ) = i∈I f (x i ) ∅ b 1 b 2 b 3 b 4 b 5 b 1 b 2 b 3 b 4 b 5 Q b 3 b 4 b 5 b 2 b 1
L 1 → L 2 is an isomorphism if f is injective, maps L 1 onto L 2 and f -1 is a morphism. An isomorphism f : L 1 → M , when M is a sub-omp of L 2 , is called embedding.
Morphisms of omp's preserve compatibility, order and orthogonality [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF]. Two elements x, y of an omp L are said to be compatible if they admit a common orthogonal base. Formally, [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF], definition 1.2.1). This orthogonal base generates a Boolean algebra which contains both elements, and is a sub-omp of L. When considered as conditions of a net system, two compatible elements will induce sequential behaviour of their neighbouring events, imposed by the orthogonal base. Indeed, orthogonality is to be interpreted as mutual exclusion of the corresponding conditions, and in fact, when restricted to atoms, the compatibility relation reduces to orthogonality [START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF]. We say two elements are incompatible, denoted x $ y, when they are not compatible. Compatible elements x and y will be denoted by x $ y. An embedding that preserves and reflects compatibility is called strong ( [START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF], Definition 15).

x $ y ⇔ ∃a, b, c ∈ L : (a ⊥ b ⊥ c ⊥ a) and (x = a ∨ b and y = b ∨ c) ([
A two-valued state of an omp is a well-known concept in the literature on Quantum Logics ([16], Chapter 2). We will use here a specific definition of two-valued state, state for short, apt to the omp's of the regions of CETS.

Definition 7 (State of an OMP). ([2]

, theorem 41) Let A(L) be the set of atoms of an omp L. Let E be the set of the maximal cliques of ⊥ in L restricted to A(L); let C be the set of the maximal cliques of $ in L restricted to A(L); let

S = {s ∈ C | ∀e ∈ E |s ∩ e| = 1}, then the up-closure of s, ↑(s), in L is a state in L.
Since states of L are uniquely defined by means of up-closures of maximal cliques of $, we will frequently refer to these maximal cliques as the states of the omp L when clear from the context.

The collection of all states of an omp L will be denoted by S(L). It will also be useful to consider the collection of states which contain a given element. Given an element x ∈ L, let us define the notation S x := {s ∈ S(L) | x ∈ s}.

The following section is concerned with the representation of an omp as a collection of subsets of its states. The order in L will then be set inclusion, and the orthocomplement will simply be given by set complement. When this is the case, the omp is said to admit a concrete representation.

Properties of Regional OMP's

As explained in more detail in section 2.4 below, omps can be constructed from the regions of a cets. These omps will, from now on, be referred to as regional omps. They have been shown, in [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF], to satisfy the following properties. We shall not get in the details, and the reader is referred to [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF] for a full explanation of these, and the implications they provide in the interpretation of omps for system specification.

• Regularity: L is regular whenever every set of pairwise compatible elements is a compatible set [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF] • Richness: L is rich whenever it admits a concrete representation [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF] There is some inconsistency in the terminology referring to these concepts in the literature, and it has led to some confusion in the axiomatic definition of these properties. Indeed regularity is the term used in [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF], whereas [START_REF] Richard | The Structure and Interpretation of Quantum Mechanics[END_REF] refers to this property as coherence. The concept of coherence is however slightly more general in [START_REF] Foulis | Operational Statistics. I. Basic Concepts[END_REF][START_REF] Randall | Operational statistics. II. Manuals of operations and their logics[END_REF].

Richness is axiomatically defined in [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF], and assumed through [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF][START_REF] Bernardinello | Synthesis of transition systems from quantum logics[END_REF][START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF][START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF]] to be equivalent to the notion of primeness defined in [START_REF] Richard | The Structure and Interpretation of Quantum Mechanics[END_REF]. The two definitions however, differ enough to raise a doubt, as their equivalence might not seem straightforward. This equivalence is, in what follows, proven, for the sake of mathematical rigour.

An omp L is called unital when for every element x ∈ L \ {0} there is a state s ∈ S(L) such that x ∈ s.

Definition 8 (Rich OMP). [16] An omp L is rich iff ∀x, y ∈ L : S x ⊆ S y ⇒ x ≤ y
In Chapter 2 of [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF], a few equivalent characterisations of richness are shown. In particular we are interested in the following:

Theorem 1. [16] An omp L is rich iff it is unital and ∀x, y ∈ L : x $ y ⇒ ∃s ∈ S(L) :
x ∈ s and y ∈ s Rich omp's where shown in Chapter 2 of [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF] to also be unital.

Definition 9 (Prime OMP). [14] An omp L is prime iff ∀x, y ∈ L : x = y ⇒ ∃s ∈ S(L) : x ∈ s ⇔ y / ∈ s
We now prove two lemmas that will be used in the proof of Theorem 2 below:

Lemma 1. In an omp L we have that: ∀x, y ∈ L : x $ y ⇒ ∃s ∈ S(L) : x / ∈ s and y / ∈ s if and only if: ∀x, y ∈ L : x $ y ⇒ ∃s ∈ S(L) : x ∈ s and y ∈ s Proof. Suppose that L is an omp such that ∀x, y ∈ L : x $ y ⇒ ∃s ∈ S(L) : x ∈ s and y ∈ s and consider two elements x, y ∈ L : x $ y; then it is also true that x $ y . By hypothesis we have that ∃s ∈ S(L) : x ∈ s and y ∈ s and neither x nor y belong to s, which means that ∃s ∈ S(L) : x / ∈ s and y / ∈ s for every pair of incompatible elements in L.

Conversely, suppose that L is an omp such that ∀x, y ∈ L : x $ y ⇒ ∃s ∈ S(L) : x / ∈ s and y / ∈ s and consider two elements x, y ∈ L : x $ y; again it is true that x $ y . By hypothesis we have that ∃s ∈ S(L) : x / ∈ s and y / ∈ s and both x and y belong to s, which means that ∃s ∈ S(L) : x ∈ s and y ∈ s for every pair of incompatible elements in L. Lemma 2. A prime omp L is also unital.

Proof. For every element x ∈ L distinct from 0 we have that ∃s ∈ S(L) : x ∈ s ⇔ 0 / ∈ s and since 0 doesn't belong to any state of L then x belongs to s.

Theorem 2. An omp L is rich iff it is prime.

Proof. Suppose that L is rich and consider two elements x, y ∈ L : x = y. There are two cases:

x $ y : in this case there exists a Boolean subalgebra of L that contains both x and y. For both of them we calculate the set of atoms below them, A ↓ (x) = {a ∈ A(L) | a ≤ x} and A ↓ (y) = {a ∈ A(L) | a ≤ y}. Such sets differs in at least an element z, because otherwise it would mean, out of Axiom 4 in Definition 5, that x = y. If we consider a state s ∈ S(L) such that z ∈ s then s will contain only one element between x and y because z belongs to to only one of the two sets A ↓ (x) and A ↓ (y). The existence of s is ensured by the fact that L is unital, which means that ∃s ∈ S(L) : x ∈ s ⇔ y / ∈ s. x $ y : if x $ y then we also have that x $ y . From the hypothesis we have that ∃s ∈ S(L) : x ∈ s and y ∈ s which means that ∃s ∈ S(L) : x ∈ s ⇔ y / ∈ s. Now, suppose that L is prime, for Lemma 2 L is also unital. Consider two elements x, y ∈ L : x $ y. That means that x $ y and in particular x = y . By Definition 9 of prime omp we have that ∃s ∈ S(L) : x ∈ s ⇔ y / ∈ s which means that ∃s ∈ S(L) : x ∈ s ⇔ y ∈ s. Again, we consider two cases:

x ∈ s : then ∃s ∈ S(L) : x ∈ s and y ∈ s x / ∈ s : then ∃s ∈ S(L) : x / ∈ s and y / ∈ s which, as shown in Lemma 1, is equivalent to ∃s ∈ S(L) : x ∈ s and y ∈ s .

Saturated Transition System and the Stability Problem

As anticipated in section 2.3 above, we know that, given a cets A = (S, E, T ),

H(A) = (R(A), ⊆, ∅, S, ( .
) ) with ( . ) the set complement operation, is an omp. Moreover, this regional omp is regular and rich.

We know as well that, given a regular and rich omp L, J(L) = (S(L), E(L), T (L)) is a cets; where E(L) and T (L), the sets of, respectively, events and transitions, are constructed by exploiting symmetric differences between states as defined in [START_REF] Bernardinello | An algebraic model of observable properties in distributed systems[END_REF] and [START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF]. Formally, a state is a collection of elements of L, and each event e ∈ E(L) will be of the form e = [s, s ] = (s \ s , s \ s), as such, it is characterised by its sets of pre-conditions, and post-conditions seen as elements of L. The underlying graph of J(L) is complete, presenting a transition (s, [s, s ], s ) between each ordered pair of distinct states (s, s ), and is therefore called the saturated transition system of L.

In [START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF] it is shown how the construction of J(L) can be done by considering exclusively the atoms of L. In particular, states are determined by their atoms, and so are events, as symmetric differences of the atoms of states. A natural question concerns the full axiomatic definition of the regional omps. We know that, in the general case, given a rich and regular logic L, L embeds in H(J(L)). We consequently define as stable an omp L isomorphic to H(J(L)).

It is our long-term goal to fully characterise the class of stable omp's, and it is conjectured that it coincides with the class of regional omp's. In the next section, we will present a compositional approach which will allow us to prove the stability of a wide class of regional omps.

Composition of OMPs and their Stability

This section will start introducing a rather general composition operation for omp's. The result of this composition will not always be an omp but we show

0 I 1 I x x I 0 1 = φ 1 (0 I ) 1 1 = φ 1 (1 I ) a 1 b 1 c 1 = φ 1 (x) a 1 b 1 c 1 = φ 1 (x ) L 1 0 2 = φ 2 (0 I ) 1 2 = φ 2 (1 I ) a 2 b 2 φ 2 (x ) = c 2 a 2 b 2 φ 2 (x) = c 2 φ 1 φ 2 L 2 Figure 3: A V-formation
that it is the case in particular instances.

Composition of OMP's

We here present a composition operation for omp's, in the fashion of those presented in [START_REF] Bruns | Amalgamation of ortholattices[END_REF] for modular ortholattices and orthomodular lattices. We focus on the class of orthomodular posets, of which orthomodular lattices are a subclass. Note that, whereas the cases treated in [START_REF] Bruns | Amalgamation of ortholattices[END_REF] are classes of algebras, our class is not, since in omp's joins and meets are only partially defined.

Definition 10 (V-formation). A V-formation of omp's is a tuple (I, L 1 , L 2 , φ 1 , φ 2 ), such that I, L 1 , and L 2 are omp's, and φ i :

I → L i (i = 1, 2)
are morphisms of omp's.

A V-formation serves as specification for composing L 1 and L 2 on the common interface I. In categorical terms, it is simply a diagram in the category of omp's. Strictly speaking, the interface would only be φ -1 1 (L 1 ) ∩ φ -1 2 (L 2 ), so in order to simplify notation, we here consider V-formations in which φ i (i = 1, 2) are embeddings, and so for each i = 1, 2, φ i (I) is a sub-omp of L i isomorphic to I.

Given a V-formation of omp'embeddingss, we propose a straightforward composition operation. This construction is inspired by co-equalisers [START_REF] Goldblatt | Topoi: the categorial analysis of logic[END_REF], however, universality of the construction remains an open problem. As a matter of fact, the nature of the proposed composition interprets the interface I as a subomp of each component, so as to identify the two copies element-wise. In this sense, it requires the morphisms of the V-formation to be actual embeddings.

Definition 11 (Equivalence induced by a V-formation). Let V = (I, L 1 , L 2 , φ 1 , φ 2 ) be a V-formation of omp's. Consider L, the disjoint union of L 1 and L 2 such that φ i : I → L (i = 1, 2), with φ 1 (I) ∩ φ 2 (I) = ∅. The equivalence relation induced by V is the binary relation

∼ V := {(x, x) ∈ L2 | x ∈ L}∪{(φ 1 (x), φ 2 (x)) ∈ L2 | x ∈ I} ∪ {(φ 2 (x), φ 1 (x)) ∈ L2 | x ∈ I}.
It is straightforward to verify that ∼ V is an equivalence relation. It is reflexive, and symmetric by construction. If x, y, z ∈ L are such that x ∼ V y and y ∼ V z then it must be either z = x or z = y, so ∼ V is transitive. We will denote the equivalence class of an element x ∈ L by [x]. This equivalence relation satisfies these additional properties:

Proposition 1. Let V = (I, L 1 , L 2 , φ 1 , φ 2
) be a V-formation of omp's, and ∼ V as in Definition 11. Then:

1. ∀i ∈ {1, 2} : ∀x ∈ L i : [x] ∩ L i = {x} 2. ∀x, y ∈ L : x ∼ v y ⇔ x ∼ V y 3. [0 1 ] = [0 2 ] = {0 1 , 0 2 } and [1 1 ] = [1 2 ] = {1 1 , 1 2 } 4. ∀x, y ∈ L : (x ≤ y) ⇒ ¬(∃x ∈ [x], ỹ ∈ [y] : ỹ < x).
Proof.

1. By construction.

2. Let x, y ∈ L be two distinct elements such that x ∼ y. From the definition of ∼ V , it follows that ∃z ∈ I : φ 1 (z) = x and φ 2 (z) = y (up to swapping of x and y). Since φ 1 and φ 2 are omp-morphisms, it follows that φ 1 (z ) = x and φ 2 (z ) = y , hence x ∼ y .

3. It is a direct consequence of Axioms 1 and 2 in Definition 6.

4. Let x ≤ y then ∃i ∈ {1, 2} : x, y ∈ L i , and x ≤ i y. Analogously, ∃j ∈ {1, 2} : x, ỹ ∈ L j , and ỹ < j x. Clearly it must be that i = j. Now, x ∼ V x, and y ∼ V ỹ ⇒ ∃a, b ∈ I : φ i (a) = x, φ j (a) = x, φ i (a) = y, and φ j (b) = ỹ. Since φ i is an omp-embedding, it must reflect the order, yielding a ≤ b, but φ j preserves the order, and so x ≤ ỹ.

These results allow for endowing the quotient L/ ∼ V with a structure.

Definition 12 (I-pasting of OMP's). Consider the setting of Definition 11, and define: Then the I-pasting of L 1 and L 2 induced by V is the structure

1. L = L/ ∼ V , 2. 0 = [0 1 ] = [0 2 ] and 1 = [1 1 ] = [1 2 ], 3. ∀[x] ∈ L : [x] = [x ], 4
L 1 |I|L 2 = L, ≤ , (•) , 0, 1 .
It follows immediately from Proposition 1 (4.) that ≤ is an order relation. Proposition 1 (2.), states that ∼ V is congruence for complementation, and so (•) is well-defined on L. It is furthermore trivial to verify that 0 and 1 are respectively the minimal and maximal elements in L. Note that whenever x ⊥ y, then [x] ⊥ [y]. So defined, the composition of two omp's over an interface is simply obtained by identifying the elements whose pre-images through φ 1 and φ 2 coincide. In general, such an I-pasting will be an orthocomplemented partial order [START_REF] Richard | The Structure and Interpretation of Quantum Mechanics[END_REF]. This is, however, not sufficient to guarantee that it is in fact an omp.

The following example shows that it can fail to be.

Example 2. With reference to Figure 3. Let I = {0, 1, x, x } be an omp with 0 ≤ x, x ≤ 1. For i = 1, 2, let L i be Boolean algebras with three atoms each {a i , b i , c i }. Since φ i are omp-morphisms, φ i (0) = 0 i , and

φ i (1) = 1 i . Now let φ 1 (x) = c 1 , so that φ 1 (x ) = c 1 = a 1 ∨ b 1 ; and φ 2 (x ) = c 2 so that φ 2 (x) = c 2 = a 2 ∨ b 2 . In this case, ∼ is the reflexive and symmetric closure of {(0 1 , 0 2 ), (1 1 , 1 2 ), (c 1 , c 2 ), (c 1 , c 2 )}, let [x] denote its equivalence classes. In L 1 |I|L 2 , we have that [a 1 ] ≤ [c 1 ] = [c 2 ] ≤ [a 2 ]. Hence, [a 1 ] and [a 2 ] are or- thogonal, but they have no least upper bound. Indeed, [a 1 ], [a 2 ] ≤ [b 1 ], [b 2 ] and [b 1 ] ∧ [b 2 ] = [0].
In what follows, we will study cases in which this composition is actually an omp.

Extending a system with a sequential component

It is a known result [START_REF] Pták | Orthomodular Structures as Quantum Logics[END_REF] that whenever L 1 and L 2 are omp's, and I = {0, 1} is the trivial Boolean algebra, then L 1 |I|L 2 , as defined in the last subsection, is an omp. It was further shown in [START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF] that in this case, if L 1 and L 2 are stable, then so must be L 1 |I|L 2 . This composition operation corresponds to the parallel composition of the corresponding operand systems. Indeed, the two systems are simply considered as a whole, although they remain independent, they do not synchronise or exchange information. Since the result of L 1 |I|L 2 is stable, it can be itself an operand for a further composition, and so this composition can be iterated, in order to generate a wide class of stable systems. As an operation, it is associative and commutative.

We consider now the case in which the interface is a non-trivial Boolean algebra I = {0, 1, y, y } whose atoms are: A(I) = {y, y }. With such an interface, we impose that the corresponding saturated transition systems synchronise according to the specified embeddings. One of the components will be a finite Boolean algebra, and will be denoted B. Boolean algebras considered as omp's were shown to be stable in [START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF]. The proof of stability will require the notion of free atom. An atom is free in an omp if it belongs to just one maximal Boolean subalgebra. When seen as the region of a cets, a free atom is a minimal region belonging to a single regional partition. The considered embeddings will then identify a free atom of an omp with an atom of B. Theorem 3. Let L be an omp, and x ∈ A(L) be an atom. Let B be a finite Boolean algebra, and a ∈ A(B). Let I = {0, 1, y, y }, and define φ L : I → L, and φ B : I → B such that φ L (y) = x, and φ B (y) = a. Then L|I|B induced by the V-formation (I, L, B, φ L , φ B ) is an omp.

Proof. After Proposition 1, it suffices to show that orthogonal joins are well defined, and that the orthomodular law holds. First note that in this case, the only identifications are , for which the join must be [START_REF] Bernardinello | Synthesis of net systems[END_REF], and is well defined. Now, from the Definition 12 (4.) of ≺, it follows that for every pair of ordered elements [c] ≤ [f ], there must be one L ∈ {L, B} such that c ∈ L ∩ [c] and f ∈ L ∩ [f ], with c ≤ f . Now, this ensures, on one hand, that orthogonal joins (and meets) are well defined in the I -pasting, whenever they are well-defined on L and B.

[0] = {0 L , 0 B }, [1] = {1 L , 1 B }, [x] = {x,
Indeed if c ∈ L ∩ [c] and f ∈ L ∩ [f ], with c ≤ f holds for both L = L and L = B, φ L preserving order, it must be either c = 0 L or f = 1 L .
On the other hand, since

L is an omp, then c ≤ f implies that f = c∨(f ∧c ), hence [f ] = [c] ∨ ([f ] ∧ [c] ).
In the following, L|I|B will refer to this particular construction, and L will be assumed to be stable. Furthermore, we will suppose that φ L (y) = x is a free atom, and show that L|I|B is stable whenever L is.

We start defining J (L) = (Q L , E L , T L ) in the following way:

Q L = S y ∪ {s ∪ {v i }|s ∈ S y , φ B (y) = v i ∈ A(B)} E L = {[s, s ]|s, s ∈ Q L , s = s } T L = {(s, [s, s ], s )|s, s ∈ Q L , [s, s ] ∈ E L , s = s }. Lemma 3. J (L) is isomorphic to J(L|I|B). Furthermore, for every v i ∈ A(B) such that φ B (y) = v i , the subgraph of J (L) induced by S(y)∪S(v i ) is isomorphic to J(L).
Proof. Every state q ∈ Q L contains one, and only one, atom of B and since φ L (y) is a free atom, it has one, and only one, atom for every Boolean algebra of L. Hence the elements of Q L coincide with the elements of S(L|I|B). Since the construction of J (L) is completely determined by the set of states as in the construction of J(L|I|B), the two transition systems J(L|I|B) and J (L) are isomorphic.

We also observe that for every atom v i ∈ A(B) distinct from φ B (y) the elements in S(y) ∪ S(v i ) coincide with the elements of S(y) ∪ S(y ), which is the set of states of J(L). From this observation it is easy to see that there is an isomorphism between J(L) and any subgraph of J (L) induced by a set of states in the form S(y) ∪ S(v i ).

This last lemma will permit to consider J (L) instead of J(L|I|B). Lines represent transitions in both directions. Dashed lines have an incidence with respect to φ L (y) or [φ L (y)], whereas solid lines are independent from them.

Lemma 4. If a region H ∈ R(J (L)) contains a state s ∪ {v i } ∈ Q L and S vi ⊆ H then ∀s ∪ {v j } ∈ Q L : s ∪ {v j } ∈ H.
Proof. Since S vi ⊆ H there must be a state s ∪ {v i } / ∈ H, which means that the event [s, s ] is an event labeling a transition exiting H. Now suppose that there is a state s ∪ {v j } ∈ Q L that doesn't belong to H. This means that the transition from s ∪ {v j } to s ∪ {v j } ∈ Q L does not exit H, but such a transition is also labeled [s, s ], which is not possible since that would mean that H is not a region. Lemma 5. Every atomic region of J (L) is in the form S x for x ∈ A(L|I|B).

Proof. Assume that there is an atomic region H ∈ {S x | x ∈ A(L|I|B)}. Consider the subgraphs induced by S y ∪ S vi for all the φ B (y) = v i ∈ B and call H vi ⊂ H the sets H ∩(S y ∪S vi ). All the H vi are atomic regions in every subgraph of J (L) induced by S y ∪ S vi , since they are all isomorphic to J(L). The regions H vi can be atomic or not. If they are atomic, then they must coincide with an atomic region S x , with y = x ∈ L n . Hence, after Lemma 4, H ∈ {S x } x∈Ln+1 . If they are not atomic, then there are H v1 ⊂ H v1 ,...,H v k ⊂ H v k from which we can make the region i∈{1,..,k} H vi ⊂ H, hence H is not atomic. Proof. We wish to show that H(J(L|I|B)) L|I|B. With Lemma 3, we reduce it to showing that H(J (L))

L|I|B. Since H(J (L)) is a finite omp, it is characterised by the orthogonality relation among its atoms. Now, Lemma 5 states that each atom of H(J (L)) corresponds to an atom of L|I|B, and it was shown in [START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF] that every atom of L|I|B must be an atom of H(J (L)).

For each pair of elements x ∈ A(L), y ∈ A(B), there is a state in s ∈ J(L|I|B) such that [x] ∈ s and [y] ∈ s. Hence, the pasting must preserve incompatibility. Since the pasting also preserves orthogonality, we have that L|I|B, and H(J(L|I|B)) have same collection of atoms, with identical orthogonality relations. As it was shown in [START_REF] Bernardinello | On the decomposition of regional events in elementary systems[END_REF], this is sufficient to state that H(J(L|I|B)) L|I|B

Stability of Atom Refinement

The operation of refining an atom of an omp into two new atoms preserves stability.

Theorem 5. Let L be a stable omp. Let x ∈ A(L). Consider M a = (A(L) \ {x}) ∪ {y, z}, in which all orthogonal atoms of L remain orthogonal in M a , and all atoms orthogonal to x in L are orthogonal to both y and z in M a . Then the omp M generated by M a is stable.

Proof. We will consider only the atoms of L and M , and states as represented by maximal cliques of $ as in Definition 7. Let S x be the set of states of L not containing x. By construction of M a , for each state s ∈ S x of L there are two states of M , in S y and S z respectively. Furthermore, the states of S x all contain an atom orthogonal to x in L, and it will be orthogonal to both y and z in M . Thus, S x , S y and S z constitute a partition of the states of M .

Starting from the states of M as partitioned above, it is possible to define the following sets of events:

E x = {[s, s ] | s, s ∈ S x , s = s }, E y,z = {[s, s ]|s ∈ S y , s ∈ S z }, E y = {[s, s ]|s, s ∈ S y }, E z = {[s, s ]|s, s ∈ S z }, E x ,y = {[s, s ] | s ∈ S x , s ∈ S y } and E x ,z = {[s, s ] | s ∈ S x , s ∈ S z }.
Let A y be the transition system with the following sets of states and events: S x ∪ S y and E x ∪ E y ∪ E x ,y , let, symmetrically, A z be the transition system whose states are S x ∪ S z and whose events are E x ∪ E z ∪ E x ,z . We note that both R(A y ) and R(A z ) are isomorphic to the regions of the saturated system J(L) since in both cases of R(A y ) and R(A z ), atoms y and z replace uniformly x. Moreover, since states S y and S z are disjoint, it is possible to construct the cets A = A y ∪ A z endowed by all the new events in E y,z . We note that R(A) must contain R(J(M )) since cets A, having less events than J(M ), can have less constraints in the construction of its regions. We want to show that R(A) = R(J(M )). Let, by contradiction, r be a region in R(A) not belonging to R(J(M )).

If r ⊆ S x , then r ∈ R(J(M )) since all the labels in E x belong to both cets A and J(M ) and the new events in E x ,y and E x ,z are distinct copies of the original events E x ,x in J(L), so they do not create new regions. If r ⊆ S y and, symmetrically, for r ⊆ S z then r must be a region in R(J(M )) since all the labels in E y and E x ,y , and symmetrically E z and E x ,z are, by construction, isomorphic to the labels E x ,x in J(L) and all the new labels E y,z are exiting from or, respectively, entering in r. The only remaining case could be for r being a minimal region in R(A) and a non minimal region in R(J(M )) but this would be in contradiction with y and z being atoms in M .

i : I → B i (i = 1, 2), such that φ i (x) = c i . B 1 is a stable [b 3 ] [b 4 ] [b 5 ] [b 3 ] [b 4 ] [b 5 ]
Figure 5: Greechie diagrams of two omps. The omp depicted below, is obtained from the one depicted above, by refining the atoms as shown. Since the omp above is stable, so is the one below. The Greechie diagram of L|I|B 3 is depicted at the top of Figure 5. Thanks to Theorem 5, we can now split any atom of L|I|B 3 , obtaining, for example, the stable omp depicted at the bottom of Figure 5.

A free atom of the operands can be refined both after and before the composition operation, in this second way, only one of the two refining atoms will be used as interface with the appended sequential component, the other one remaining free for further composition. With this method, one can iterate the composition operation without worrying about exhaustion of available free atoms of the original system.

Example 5. Consider a system made of two sequential components each of which can get to a state for which they require the same resource. If each of these components can be in two additional states, the regional omp for this system is represented as L 1 in Figure 6. B 1 and B 2 represent the two sequential components, and x 1 , x 2 correspond to their mutually exclusive states. B c represents the state of the resource c, it can be in state x i , indicating that B i holds c, or in state y 1 , indicating that c is available. When c is available, no other component is involved in the coresponding local state of B c , and so y 1 is a free atom. L 1 is isomorphic to the omp at the top of Figure 5, and was shown to be stable in Example 4. We may want to make the resource c available for a third sequential component B 3 , so we can use Theorem 4 to compose L 1 and B 3 on y 1 , obtaining the stable L 2 = L 1 |I|B 3 a shown in Figure 6. However, in this new compound system, the resource must be held by one of the three components B i , as B c has no more free atoms. No additional component can be added to the system, to compete for c. Instead of performing the composition L 1 |I|B 3 directly, we can first make use of Theorem 5, and refine y 1 in L 1 into two new free atoms x 3 , y 2 , thus obtaining the omp L 3 of Figure 6. We can now compose it with B 3 on x 3 , thus obtaining L 4 = L 3 |I|B 3 , which is stable. One can see that y 2 remains a free atom, a local state representing that the resource is available. This process can be iterated, to obtain a system with n sequential components competing for the same resource.

L 1 x 1 x 2 y 1 B c B 1 B 2 L 2 x 1 x 2 y 1 B c B 1 B 2 B 3 L 3 x 1 x 2 x 3 y 2 B c B 1 B 2 y 1 L 4 x 1 x 2 x 3 y 2 B c B 1 B 2 B 3

Conclusion

In [START_REF] Bernardinello | On stability of regional orthomodular posets[END_REF], a collection of regional omps were shown to be stable. Furthermore, the parallel composition operation was shown to preserve stability. In the present work, we have formalised two additional operations which preserve stability. One corresponds to the refinement of a local state into two. The second corresponds to the extension of a system with a sequential component which synchronises with it over a local state which is not already a synchronisation. With these elements, we can define an algebra of system omps, such that all its elements are stable. However, not every regional omp can be obtained with the defined operations. For instance, a strong limitation is the restriction of the composition operation to interfacing on free atoms. Another limitation of this operation is that it does not allow for extending a system with a sequential component that would synchronise with the system at more than one state. A clear goal in our approach is to extend the composition operation so as to generate every regional omp. In this way, by showing that we can generate all regional omps with stability preserving operations, we would prove the conjecture that all regional omps are stable.

4 Figure 1 :

 41 Figure1: To the left, a CENS. To the right, its sequential case graph. The state q 1 corresponds to the marking {b 1 , b 4 }, at which e 1 and e 2 are concurrently enabled. To the right, the subset of states {q 1 , q 2 } is in evidence; it corresponds to the extension of b 1 , namely the set of markings containing b 1 .

Figure 2 :

 2 Figure 2: Two representations of the omp obtained by ordering the regions of the transition system in Figure 1. To the left, all regions are represented in a Hasse diagram. To the right, only the atoms are represented; each of the two solid lines is a maximal clique of orthogonality. This last representation is called Greechie diagram.

  . [x] ≺ [y] ⇔ ∃x ∈ [x] : ∃y ∈ [y] : x ≤ y in L, and 5. ≤ ⊆ L × L as the transitive closure of ≺.

Example 3 .

 3 With reference to Figure 2, b 1 , b 2 , b 4 , b 5 and b 3 are all atoms, however b 1 , b 2 , b 4 , b 5 are free, but b 3 is not. Indeed, b 3 belongs to two maximal Boolean algebras.

  a} and [x ] = {x , a }, all other equivalence classes being singletons. Since both x and a are atoms, we have that ≤=≺, in the setting of Definition 12. Furthermore, for each pair of orthogonal elements [c] ⊥ [d], there must be c ∈ L ∩ [c] and d ∈ L ∩ [d], where L ∈ {L, B} such that c ⊥ d in L . If this holds for both L = L and L = B, then the only possibility is [c] = [x] and [d] = [x ]

3 Figure 4 :

 34 Figure 4: Construction of J (L), where L is the omp of Figure 2, I is as in Theorem 3, and B is a Boolean algebra with three atoms: A(B) = {φ B (y), v 1 , v 2 }.Lines represent transitions in both directions. Dashed lines have an incidence with respect to φ L (y) or [φ L (y)], whereas solid lines are independent from them.

Theorem 4 .

 4 Let L be a stable omp, and x ∈ A(L) be a free atom. Let B be a finite Boolean algebra, and a ∈ A(B). Let I = {0, 1, y, y }, and define φ L : I → L, and φ B : I → B such that φ L (y) = x, and φ B (y) = a. Then L|I|B induced by the V-formation (I, L, B, φ L , φ B ) is stable.

Example 4 .

 4 Consider three Boolean algebras B 1 , B 2 and B 3 with three atoms each, A(B i ) = {a i , b i , c i } for i = 1, 2, 3. Let I = {0, x, x , 1} and consider the two omp-morphisms φ

logic, and c 1

 1 is clearly a free atom, so after Theorem 4,L = B 1 |I|B 2 is a stable omp. L is isomorphic to the omp in Figure 2, by considering b 3 = [φ 1 (x)] = [φ 2 (x)].Since L is stable, and b 5 is a free atom, we can compose it with B 3 , by means of the morphisms φ L and φ 3 , provided by φ L (x) = b 5 and φ 3 (x) = c 3 .

Figure 6 :

 6 Figure 6: Refining a free atom

Acknowledgements

We wish to thank Lucia Pomello, and Luca Bernardinello for the fruitful discussions. Thanks to the three anonymous reviewers for the useful remarks and suggestions. This work is partially supported by MIUR.


Roberto Ausas
email: rfausas@icmc.usp.br

Gustavo Buscaglia
email: gustavo.buscaglia@icmc.usp.br

Vitoriano Ruas
email: vitoriano.ruas@upmc.fr

Finite element method for incompressible viscous flow with immersed pressure jumps with applications to actuator disks and microfluidics

Keywords: Finite elements, Incompressible viscous flows, Navier-Stokes, Pressure jump, Immersed boundary, Surface tension

HAL is

Introduction

Simulation of fluid flow in several industrial problems requires taking into consideration singular forces at internal interfaces. An important example is the simulation of wind turbines, ship propellers and cooling fans by means of the actuator disk model (see e.g. [START_REF] Réthoré | Verification and validation of an actuator disc model[END_REF][START_REF] Harrison | Accuracy of the actuator disc-RANS approach for predicting the performance and wake of tidal turbines[END_REF][START_REF] Ammara | A Viscous Three-Dimensional Differential/Actuator-Disk Method for the Aerodynamic Analysis of Wind Farms[END_REF]). This model replaces the individual action of the blades of the turbine (or propeller, or fan) by a singular force distributed on an immersed imaginary disk that represents the whole rotor. Singular forces also arise, for example, in the simulation of surface tension phenomena in multiphase flows.

Actuator disk models, or immersed singular forces in general, are especially useful when combined with non-conforming meshes; i.e., meshes that do not conform to the geometrical location of the interfaces at which the forces arise. This allows one to use a single mesh to explore different geometrical configurations, thus reducing the pre-processing burden that can dominate the CFD design pipeline. Non-conforming meshes, however, require very robust approximation methods, that do not loose stability when the intersection of the interface with the grid cells happens to be oddly shaped. Indeed the accurate approximation of pressure and velocity exhibiting strong and weak discontinuities at immersed interfaces remains numerically challenging and has attracted intensive research in the last years. In the finite element literature modifications or enrichment of standard finite element spaces have been proposed to better accommodate the discontinuities and to obtain improved convergence rates. Among these, the XFEM method introduced in [START_REF] Pd | A finite element technique for multifluid incompressible flow using Eulerian grids[END_REF] is well known (see also [START_REF] Gross | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Reusken | Analysis of an extended pressure finite element space for two-phase incompressible flows[END_REF] for applications to two phase flow problems). Elementwise condensable enrichment can be found in the literature [START_REF] Coppola-Owen | Improving eulerian two-phase flow finite element approximation with discontinuous gradient pressure shape functions[END_REF][START_REF] Ausas | A new enrichment space for discontinuous pressures[END_REF][START_REF] Kamran | A locally extended finite element method for the simulation of multi-fluid flows using the particle level set method[END_REF] or locally modified spaces [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF] which do not alter the structure of the global system of equations to be solved. Variants of Nitsche's methods were also proposed to weakly enforce coupling conditions at immersed interfaces. Together with special function spaces, stabilization techniques are essential for robustness and accuracy. In [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method[END_REF][START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for stokes' problem[END_REF][START_REF] Hansbo | A cut finite element method for a stokes interface problem[END_REF] significant efforts have been devoted to tackle such issues. There, stabilization methods are proposed to attain improved accuracy and to avoid ill-conditioning of the linear systems resulting from the so called bad-cut elements. All these developments have mainly focused on multiphase flow models, and have not been assessed or optimized for general interfaces with singular forces, such as those that appear in actuator disk models. It is noticeable that our method avoids remeshing in order to conform to a moving interface where singular forces act, like in other known approaches that can be used to tackle such a problem.

In [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF], a stabilization of the pressure jump across an interface as a result of prescribed singular forces has been proposed. The method is close to the one studied in [START_REF] Buscaglia | Finite element methods for the Stokes system with interface pressure discontinuities[END_REF]. Both are based on a Petrov-Galerkin formulation of the Stokes problem of the type proposed by Franca, Hughes and Balestra [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF] or Douglas and Wang [START_REF] Douglas | An Absolutely Stabilized Finite Element Method for the Stokes Problem[END_REF], depending on a mesh-dependent parameter δ. The main difference between them is the addition, in [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF], of a pressure jump term multiplied by another mesh-dependent parameter ε > 0. Such a term penalizes errors in the mechanical equilibrium at the immersed interface, increasing the accuracy of the imposition of the singular force distribution. In this article we explore this term for problems with singular forces such as those arising from the actuator disk model. We combine this extra penalty with the Algebraic Subgrid Scale (ASGS) method of Codina [START_REF] Codina | A stabilized finite element method for generalized stationary incompressible flows[END_REF]. The ASGS variant is a popular generalization of the methods in [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations[END_REF][START_REF] Douglas | An Absolutely Stabilized Finite Element Method for the Stokes Problem[END_REF] that includes transient and inertial effects. The combined method is an ASGS method itself, now incorporating the aforementioned pressure jump term which is very easy to implement. As shown by numerical experiments, the proposed method is insensitive to bad-cut elements and yields solutions with enhanced accuracy at the interface. A specific formula is proposed for ε after extensive tests, leaving a method with no free parameters.

An outline of the paper is as follows. In Section 2 we formalize the mathematical setting, while in Section 3 we recall the finite element method of [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] and its combination with the ASGS method for the Navier-Stokes equations. We also give a priori criteria for the choice of method's parameters δ and ε. In section 4 we perform a series of numerical experiments considering both Stokes and Navier-Stokes regimes, with examples involving immersed actuator disks and a microfluidic application related to the thermocapillary Marangoni migration of a droplet. We conclude in Section 5 with a few remarks on the whole work.

Mathematical setting

We consider the flow of a viscous incompressible fluid with dynamic viscosity µ and density ρ, occupying two complementary sub-domains Ω 1 and Ω 2 of a bounded domain Ω ∈ R 2 with boundary Γ, separated by an interface Σ. Although in this article we restrict ourselves to the two dimensional case for simplicity, the extension to three dimensions is straighforward. Akin to [START_REF] Buscaglia | Finite element methods for the Stokes system with interface pressure discontinuities[END_REF], in this method basically a standard continuous piecewise linear representation of both velocity and pressure is adopted. However, in order to take into account discontinuities across Σ, the pressure interpolation is modified in the elements cut by it in the way specified in [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF] that we recall hereafter. In [START_REF] Buscaglia | Finite element methods for the Stokes system with interface pressure discontinuities[END_REF] and also in [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] when N = 2, Σ was assumed to be a Dini-smooth open arc (see e.g. [START_REF] Pommerenke | Boundary Behaviour of Conformal Maps[END_REF]) intersecting Γ at exactly two distinct points or a closed Dini-smooth Jordan curve [START_REF] Pommerenke | Boundary Behaviour of Conformal Maps[END_REF] completely immersed in Ω, or yet polygonal counterparts of both situations. Although we keep these assumptions here for the purpose of the formal results given in Section 3, they are not really fulfilled in all the numerical experiments performed in this work. On the other hand, akin to both papers, Ω is still assumed to be either a polygon or a domain such that Γ is of the piecewise C 1 class. For the same purpose we also assume that a non empty intersection of Σ and Γ does not contain a vertex of Γ in case Ω is a polygon, and the sinus σ of any angle between Σ and Γ is not zero. The problem to solve can be stated as follows: Let f be a given force field in [L 2 (Ω)] 2 , and ϕ be a given force distribution on Σ belonging to C 0 (Σ). Prescribing a velocity g ∈ [C 0 (Γ)] 2 satisfying Γ g • n = 0, where n is the unit outer normal vector on Γ, we wish to find a velocity u ∈ [H 1 (Ω)] 2 and a pressure p ∈ H 1 (Ω \ Σ) ∩ L 2 0 (Ω) satisfying the following Navier-Stokes system:

               ρ(∂ t u + u • ∇u) -µ∆u + ∇p = f in Ω \ Σ div u = 0 in Ω p = ϕ on Σ u = g on Γ, (1)
assorted with initial velocity conditions, where ∂ t (•) stands for partial differentiation with respect to time t and q represents the jump across Σ of the traces of a function q ∈ H 1 (Ω \ Σ) on both sides of Σ. As usual H 1 0 (Ω) is the subspace of H 1 (Ω) of those functions that vanish on Γ and L 2 0 (Ω) is the subspace of L 2 (Ω) consisting of functions f such that Ω f dx = 0.

In [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] we studied finite element approximations of (1) in the case where the inertia term ρu • ∇u is negligible as compared to the viscous term µ∆u. This corresponds to flow at a low Reynolds number, and in this case the Navier-Stokes system (1) reduces to a Stokes system, whose equivalent variational form can be stated as follows. Denoting by (•, •) the standard inner product of L 2 (Ω) in scalar, vector or tensor version, and by (•, •) D the standard inner product of L 2 (D) for D Ω, as long as µ is large enough so that inertia can be neglected and the flow is stationary, we may rewrite the underlying simplification of problem (1) in the following variational form of the Galerkin type:

               Find u ∈ [H 1 (Ω)] 2 and p ∈ L 2 0 (Ω) such that µ(∇u, ∇v) -(p, div v) = (f , v) + (ϕ, v • ν) Σ ∀v ∈ [H 1 0 (Ω)] 2 (div u, q) = 0 ∀q ∈ L 2 0 (Ω) u = g on Γ, ( 2
)
where ν is the unit normal vector on Σ oriented from Ω 1 towards Ω 2 , and the symbol • is used to represent the euclidean inner product of two vectors of R 2 .

Clearly enough, as a part of the solution of (2), the pressure p is a priori discontinuous on Σ. Observing that every v ∈ H 1 (Ω) belongs to L 4 (Ω) (cf. [START_REF] Adams | Sobolev Spaces[END_REF]), problem (2) can be easily modified into an equivalent variational form of (1), by adding the term ρ(u • ∇u, v) to the right hand side of the first (momentum) equation.

Finite element formulations

Preliminaries

Let P be a quasi-uniform family of partitions (cf. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) T h of Ω consisting of triangles with maximum edge length equal to h, and whose union is a polygonal approximation Ω h of Ω, in such a way that all of its vertices lie on Γ. We denote by Γ h the boundary of Ω h which coincides with Γ in case Ω is a polygon. The assumption that Ω be a convex domain in case it is not a polygon is further made in order to avoid some non essential complications related to external approximations. Of course this is an academic assumption, since our method works as well for non convex domains. It is natural to approximate problem (2) using a pressure finite element space consisting of discontinuous functions on Σ. One could do this by matching the mesh in such a way that Σ is approximated by the union of edges of a patch of elements in T h . However in doing so one would loose the possibility of working with structured or even uniform meshes, if Ω is rectangular for instance, since Σ may be of very irregular shape. Therefore our numerical approach is based on a modification of the usual continuous piecewise linear finite element space to represent the pressure, incorporating discontinuities over Σ (in fact, Σ h ) without requiring conforming meshes. The adopted space was originally proposed in [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF] and thoroughly assessed in [START_REF] Buscaglia | Variational formulations for surface tension, capillarity and wetting[END_REF].

Before recalling the construction of such a space and introducing the finite element formulation, let us include necessary geometrical definitions and assumptions.

For numerical purposes Σ will be approximated by a polygonal line Σ h constructed in connection with T h in the following manner (see also Figure 1). Let R h be the set of elements in T h whose closures have a non empty intersection with Σ, and S h be the subset of R h such that T ∈ S h means that there are two edges of T cut by Σ at least at one point different from their common vertex. Then Σ h is a polygonal approximation of Σ formed by the chords joining two points of T ∩ Σ on different edges of T ∈ S h that are not vertices of T . We assume Σ h to have an orientation, and denote by ν h its oriented normal. In case Σ h divides Ω h in two separate subdomains, these will be denoted as Ω h1 and Ω h2 . In this case ν h is assumed to point from Ω h1 towards Ω h2 .

With the above definitions of Σ h , we assume that whenever Σ is connected so is Σ h . Actually this assumption will hold true provided some conditions on the mesh listed in [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] are satisfied. For instance the case where a vertex or an edge of the mesh is exactly contained in Σ must be ruled out, by slightly perturbing the mesh (or Σ). We just recall here a basic condition concerning "h being small enough". Hypothesis H: The mesh size h is small enough for the maximum absolute value of the curvature [START_REF] Cartan | Leçons sur la géométrie des espaces de Riemann[END_REF] of Σ to be much smaller than h -1 . Moreover in case Ω is a polygon, h is sufficiently small for the distance from Σ measured along Γ of anyone of its vertices to be greater than h. Now if T ∈ S h , we subdivide the quadrilateral subset K of T lying on one side of Σ h , into two triangles T 1 and T 2 by means of the diagonal of K adjacent to its largest angle. We choose T i to be the triangle having 3 -i vertices coinciding with vertices of T , i = 1, 2 (cf. Figure 2). The pressure space associated with T h of discontinuous functions in the elements of S h , denoted here by Qh , was introduced in [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF]. For the sake of completeness we recall its definition as follows: Every function q ∈ Qh is necessarily continuous in Ω h \ Σ h ; in each triangle not belonging to S h , q is either constant or linear. Referring to Figure 2 for a generic triangle T ∈ S h , q |T3 is constant, q is linear in T i for i = 1, 2 and satisfies q(S 1 ) = q(S 5 ) and q(S 2 ) = q(S 4 ). An illustration of a set of basis functions for Qh in a triangle T ∈ S h is supplied in Figure 3.

q 3 q 2 q 1 Σ h S1 S3 S2 S4 S5 T 2 T 1 T 3 Figure 3
. Set of basis functions q i , i = 1, 2, 3 for the space of discontinuous pressures Qh on a triangle T ∈ S h .

Henceforth the intersection of Qh with L 2 0 (Ω h ) (cf. [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]) will be denoted by Q h . The velocity space V h is the standard space of continuous fields whose restriction to each triangle in T h is linear. We denote by V h0 the subspace of V h consisting of those fields that vanish identically on Γ h .

The proposed formulation for the Stokes case

In the sequel g h stands for the piecewise linear interpolate of g along Γ h at the vertices of T h lying on Γ. Moreover we denote by (•, •) m,D the standard inner product of H m (D) and by • m,D the associated norm. The subscript m is frequently omitted when m = 0.

Before recalling the solution method experimented in this work, we still have to define a piecewise gradient operator ∇ Σ h of a function q ∈ L 2 (Ω h ) such that q ∈ H 1 (Ω h \ Σ h ) as follows:

∇ Σ h q(x) = ∇ q| Ω h \Σ h (x), if x ∈ Ω h \ Σ h . (3)
Now given a strictly positive parameter δ and a non negative parameter ε, our finite-element approximation of ( 2) is as follows:

               Find u h ∈ V h and p h ∈ Q h such that ∀v ∈ V h0 and ∀q ∈ Q h , µ(∇u h , ∇v) Ω h -(p h , div v) Ω h = (f h , v) Ω h + (ϕ h , v • ν h ) Σ h ε( p h h , q h ) Σ h + δ(∇ Σ h p h , ∇ Σ h q) Ω h + (div u h , q) Ω h = ε(ϕ h , q h ) Σ h + δ(f h , ∇ Σ h q) Ω h u h = g h on Γ h . (4)
where η h is the jump across Σ h of a function η discontinuous on this manifold. In [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] the following well-posedness result was proven for problem (4): Proposition 3.1 Under Hypothesis H, if ε > 0 problem (4) has a unique solution.

If ε = 0, then (4) becomes [START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for stokes' problem[END_REF], which also has a unique solution, as shown below.

Proposition 3.2 Under Hypothesis H, given a strictly positive parameter δ, the following problem is well-posed:

               Find u h ∈ V h and p h ∈ Q h such that ∀v ∈ V h0 and ∀q ∈ Q h , µ(∇u h , ∇v) Ω h -(p h , div v) Ω h = (f h , v) Ω h + (ϕ h , v • ν h ) Σ h δ(∇ Σ h p h , ∇ Σ h q) Ω h + (div u h , q) Ω h = δ(f h , ∇ Σ h q) Ω h u h = g h on Γ h . (5)
Proof. Assume that ϕ, g and f vanish everywhere in their respective definition domains. Since ( 5) is a linear problem we have to prove that corresponding u h and p h vanish identically. Since u h ∈ V h0 let us take in the first and second equations of (5) v = u h and q = p h . Adding up the resulting relations we come up with,

µ ∇ u h 2 Ω h +δ ∇ Σ h p h 2 Ω h = 0. ( 6
)
Since u h = 0 on Γ, we have u h ≡ 0. On the other hand [START_REF] Hansbo | A cut finite element method for a stokes interface problem[END_REF] implies that p h is constant in each connected component of Ω h \ Σ h . If Σ h does not separate Ω h into two disjoint parts, p h is constant everywhere and thus p h ≡ 0. If Σ h separates Ω h into Ω h1 and Ω h2 , then p h = c i in Ω hi for i = 1, 2, where c 1 and c 2 are constants. Since p h ∈ Q h we have c 1 area(Ω h1 ) + c 2 area(Ω h2 ) = 0. Moreover from the first equation of [START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for stokes' problem[END_REF] we have (p h , div v) Ω h = 0 ∀v ∈ V h0 . Applying the divergence theorem in both Ω h1 and Ω h2 , taking into account that ∇ p h ≡ 0 in both sub-domains, this trivially yields

|c 1 -c 2 || Σ h v • ν h ds| = 0 ∀v ∈ V h0 .
Clearly enough there exists at least one field V h0 such that Σ h v • ν h ds = 0. Therefore c 1 = c 2 and thus p h ≡ 0.

We have thus recalled and adapted the jump-stabilized formulation of [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] for the case under study. Under suitable assumptions, a priori error estimates applying to problem (4) were proved in [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] for δ ∼ h 2 and ε ∼ h, as long as f , g, ϕ and the exact solution of (2) are sufficiently smooth. More precisely, it holds that

∇(u -u h ) Ω h + p -p h Ω h + p -p h Σ h ≤ C(u, p, f , ϕ) h, (7)
where C(u, p, f , g, ϕ) is a constant depending only on u, p, f , g and ϕ, and p is a suitable projection of p |Σ onto Σ h .

We now turn to the incorporation of inertial and transient effects by combining this formulation with the ASGS method for the Navier-Stokes equations.

Combined Algebraic Subgrid Scale and Pressure Jump Stabilization

In order to take into account the inertial term in the momentum equation with the method of the previous section it would be natural to adopt the following discrete formulation for δ > 0 and ε ≥ 0:

                     Find u h ∈ V h and p h ∈ Q h such that ∀v ∈ V h0 and ∀q ∈ Q h , µ(∇u h , ∇v) Ω h + ρ(u h • ∇u h , v) Ω h -(p h , div v) Ω h = (f h , v) Ω h + (ϕ h , v • ν h ) Σ h ε( p h h , q h ) Σ h + δ[ρ(u h • ∇u h , ∇ Σ h q) Ω h + (∇ Σ h p h , ∇ Σ h q) Ω h ] + (div u h , q) Ω h = ε(ϕ h , q h ) Σ h + δ(f h , ∇ Σ h q) Ω h u h = g h on Γ h . (8)
The formulation above works well in general, however, it does not consider stabilization of the convective terms in the momentum equation. Moreover, it is also desirable to consider non-stationary situations that are relevant in several applications. To this end, we propose the following combination of pressure-jump stabilization with the ASGS formulation proposed by Codina et al 2001 [START_REF] Codina | A stabilized finite element method for generalized stationary incompressible flows[END_REF]:

                       Find u h ∈ V h and p h ∈ Q h such that ∀v ∈ V h0 and ∀q ∈ Q h , (G u , v) Ω h + (2µ D(u h ), ∇ v) Ω h -(p h , div v) Ω h + T ∈T h δ T (G u + ∇ Σ h p h , u h • ∇v) T + + T ∈T h α T (div u h , div v) T -(ϕ h , v • ν h ) Σ h = 0 (q, div u h ) Ω h + T ∈T h δ T ρ (G u + ∇ Σ h p h , ∇ Σ h q) T + T ∈T h ε T ( p h -ϕ h , q h ) Σ h ∩T = 0 (9)
where we have defined

G u = ρ (∂ t u h + u h • ∇u h ) -f h and δ T = 4 µ ρ h 2 T + 2 u h h T -1 , (10)
α T = 2 µ + ρ u h h T , (11)
ε T = C ε δ T ρ h T ( 12
)
in which h T is the diameter of element T . Regarding this formulation the following observations are in order:

-Inertial and transient effects are now considered. The formulation is a natural extension of the singlephase ASGS formulation, which contains penalty terms for the momentum equation and for the continuity equation, respectively multiplied by δ T and α T . The proposed modification is a third penalty term, this time for the interface equilibrium equation p = ϕ. -Notice that some simplifications have been made in the discrete formulation by particularizing it to (a) piecewise linear u h , and (b) constant viscosity. If the viscosity exhibits a jump through Σ h , then the penalty term would also contain the jump of normal viscous stress across Σ h . -Elementwise penalty parameters δ T and α T are adopted, with the expressions suggested by Codina et al [START_REF] Codina | A stabilized finite element method for generalized stationary incompressible flows[END_REF]. Similarly, an elementwise ε T is implemented and given an expression with consistent units, with a dimensionless constant C ε whose ideal value found by means of our numerical experimentation is 2. -Instead of using ∇ u h , the symmetric gradient D(u h ) = 1 2 (∇ u h + ∇ T u h ) appears in the momentum residual equation. This is needed for problems with Neumann boundary conditions and/or variable viscosity.

-In problem (9) the time step index has been dropped in all unknown fields to simplify the notation. As for time discretization, this problem is solved at each step for the unknowns u h and p h in a coupled manner, using a Newton-Raphson iterative procedure together with the approximation

∂ t u h ≈ u h -u n h ∆t ( 13
)
where ∆t is the time step and u n h stands for the velocity field at the last computed step. Clearly, all formulations discussed in the previous section are particular cases of (9). For instance, in a stationary case, neglecting inertia and taking α T = 0, we recover formulation (4).

Numerical experiments

In this section we numerically assess the proposed discrete formulation with a few examples. We first consider the laminar flow in a horizontal channel in the presence of an actuator disk with several configurations and assess the convergence properties of the method. More complex configurations, considering inertial and transient effects, are also addressed. To show the versatility of the method, we finally solve a problem of interest in microfluidics, namely the capillary migration of a circular droplet with variable surface tension coefficient.

Actuator disks -Stokes regime

Let us consider the computational domain [0, L] × [0, H]. A horizontal singular force ϕ is applied along a straight vertical line defined by

Σ h = {(x 1 , x 2 ), x 1 = a, x 2 ∈ [b, H]} (14)
This is the simplest actuator disk model that can be considered. Periodic boundary conditions are assumed in the x 1 -direction and no-slip boundary conditions are taken on the bottom and top walls of the channel, i.e., u(

x 1 , x 2 = 0) = u(x 1 , x 2 = H) = 0 ( 15
)
The indeterminacy of the pressure is removed by imposing p(0, 0) = 0. An illustration of the computational domain is shown in Figure 4. Two cases are considered depending on the value adopted for b in ( 14), namely:

-Case I -b = 0: the interface Σ occupies the entire channel section. The normal force ϕ is taken equal to 1 along Σ. The exact solution is known in this academic case and corresponds to a Poiseuille flow with quadratic velocity profile. The pressure field is piecewise linear exhibiting a jump equal to the value of ϕ at the location of the interface x 1 = a. This problem has been solved several times before (see e.g. [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF][START_REF] Ausas | A new enrichment space for discontinuous pressures[END_REF]) and its solution reads

u 1 (x 1 , x 2 ) = 1 2 µ L (H -x 2 )x 2 , (16)
u 2 (x 1 , x 2 ) = 0 (17) p(x 1 , x 2 ) = - x 1 L + H(x 1 -a) (18)
where H(•) is the Heaviside function (i.e., H(s) = 1 if s > 0 and zero otherwise). A sequence of computational meshes is built of which the first one is shown in Figure 4. A characteristic mesh size h = 0.05 is assigned to this mesh. The other meshes in the sequence are obtained by subdivision of the triangles into four identical triangles. As one can observe, the interface is neither conforming nor parallel to the edges of the finite element partition. The problem is solved taking L = 3, H = 1, µ = 1 and a = 2. To illustrate the effect of the interface stabilization and the value adopted for the stabilization parameter (C ε = 2), we compare results to a case with C ε ≫ 1 and a case with C ε ≪ 1. Figure 5 shows a section of the pressure field p h along the horizontal line at x 2 = H/2 corresponding to the first mesh in the sequence in the three cases. By watching the insert, we clearly observe the formulation performs best for the proposed choice of stabilization parameter, whereas choosing the stabilization parameter being very large or very small has a deleterious effect in the pressure field in the proximity of the interface. Notice that this occurs not only in the sense being observed, i.e. the mean-square (L 2 ) sense, but also to a lesser extent in the pointwise (L ∞ ) sense.

We report in Table 1 the L 2 (Ω h ) and H 1 (Ω h ) error norms of pressure and velocity, respectively, and an estimation of the convergence rates. The pressure error is close to an O(h 1.5 ), which can be viewed as a super-convergence phenomenon, for both the unstabilized and the interface stabilized methods. However, the error is significantly smaller for the latter. -Case II -b = 0.5: this corresponds to an interface Σ h that occupies only half of the channel section.

C ε = 64 C ε = 2 C ε = 0.01
The same physical parameters and computational meshes of the previous case are used. The solution is not explicitly known in this case. The error is now measured by comparison of the solution for a given mesh to the solution obtained at the subsequent level of mesh refinement. As for the force to be applied over Σ h , in the actuator disk model, the ideia is to replace the effect of a real propeller with a permeable disk over which the forces acting on the fluid are distributed. In order to qualitatively mimic this force the following variation is chosen for ϕ (see e.g. [START_REF] Réthoré | Verification and validation of an actuator disc model[END_REF][START_REF] Mikkelsen | Actuator disc methods applied to wind turbines[END_REF])

ϕ = c (1 -x 2 ) tanh [10 (x 2 -0.5)] ( 19
)
where the factor c is set equal to 3.1344 such that the maximum value of ϕ approximately equals 1.

Pressure and velocity fields are illustrated in Figure 6. In Figure 7 the pressure jump and the applied force are plotted for the coarsest mesh level, comparing the stabilized and unstabilized cases. As one can infer from this figure, the interface stabilized formulation allows for a much better representation of the discontinuity in the pressure field on Σ h . Table 2 reports the L 2 (Ω h ) and H 1 (Ω h ) error norms for pressure and velocity respectively and the convergence rate in both cases. It should be noticed that the estimated rate of convergence is ca. 1.3, which is closer to the theoretical prediction of O(h), than the super-convergence rate of ca. 1.5 observed in the previous case. However, both are quite similar, which means that the use of the interface stabilization does not have any deleterious effect on the results considering these global norms. On the other hand, through the error measure of the jump, the advantage of the new formulation clearly appears as in the previous case.

Actuator disks -Inertial and transient effects

Stationary case

The problem setting considered is shown in Figure 8 with two inclined actuator disks positioned under a square obstacle. In this case we opt to provide all quantities and results in dimensional form. The computational domain has width L = 10 cm and height H = 20 cm and the obstacle is a square of size 2 cm (see Figure 8). Fluid properties correspond to that of air, with density ρ = 1.2 Kg•m -3 and viscosity µ = 2 × 10 -5 Pa•s. On the lateral sides of the domain no-slip boundary conditions are imposed. The bottom and top walls are respectively inlet and outlet boundaries, so the fluid is free to enter or leave the domain, however, at the outlet the horizontal component of the velocity field is set to zero. The domain is discretized with a mesh also shown in Figure 8 having 18,042 triangles and a characteristic mesh size h = 0.065 near the actuator disks. Over both disks a normal force given by

ϕ(s, t) =      C(t) s tanh [4 (0.5 -s)] if s < 0.5 C(t) (1 -s) tanh [4 (s -0.5)] if s ≥ 0.5 (20)
is applied, where for each disk s ∈ [0, 1] is a line parameter along the horizontal direction as indicated in figure 8 and C is used to control the magnitude of this force as a function of time. This way, one is able to produce some complex flow patterns around the obstacle and test the proposed method in a more challenging situation in the presence of both, inertial and transient effects. Notice that testing other geometrical configurations (disk's size, position and inclination) would very simple using this unique background mesh. We look for a steady state solution by taking a very large time step. (20) equal to 3.85 × 10 -4 Pa and the Reynolds number in this case is ca. 22.5 based on the obstacle length scale, that is 2 cm. The complex flow pattern in the present case is noteworthy, as compared to the previous simulations. Figure 10 shows the applied force and compares it to the numerical pressure jump along the interface for the stabilized and unstabilized formulations. These results certify again, the superior behavior of the former. It is important to note that the velocity field does not exhibit any significant difference with respect to the unstabilized formulation as indicated in Figure 11 where the magnitude of velocity u h 2 is plotted along the left interface and compared to the solution obtained with a fine grid having characteristic mesh size h = 0.02 near Σ h , that is, 3.25 times a smaller mesh size than the one of the orginal mesh showed in Figure 8. We have observed comparable results for other geometrical configurations of the actuator disks.

Transient case -Fixed disks

In the transient case the behavior is analogous. For illustration we progressively increased the amplitude C of the force (20) from the value 3.13 × 10 -5 to 5.95 × 10 -3 Pa for 0 ≤ t ≤ 1 s. From t = 1 s, the amplitude of ϕ remains constant at the value C(t = 1). Results are presented in Figure 12 where pressure and velocity contours are displayed at three different times, namely, t = 1, 3 and 5 s. The top part of Figure 13 shows the applied force and pressure jump for the stabilized and unstabilized formulations along the left actuator disk at two different times, namely, t = 0.5 and t = 4.5. Similar results are observed for the right disk and also at other times along the simulation. The bottom part of this figure shows the magnitude of the velocity field u h 2 along the interface for the values of t. Here again we observe the velocity along Σ h produced by the stabilized and unstabilized methods are very similar.

Transient case -Moving disk

We finally consider the geometrical setting shown in Figure 14. The computational domain is the square region [0, 10] × [0, 10] cm. A rotating disk of length 4 cm is placed at the center. Figure 14 shows the disk at several inclination angles, namely, π/16 (maximum inclination), 0 and -π/16 (minimum inclination). The angular velocity ω of the disk is prescribed accordingly by

x 2 x 1 x 2 = H = 20 cm u = 0 u = 0 Inlet Outlet, u 1 = 0 (4, 2) (2, 4) (6, 2) ϕ (8, 4)
x

ω(t) = 5π 2 8 cos(10 π t) (21)
The prescribed force is similar to previous cases except that now s ∈ [-0.5, 0.5] is a line parameter along the interface being equal to 0 at the disk's center. The amplitude C(t) is taken as

ϕ(s, t) =      -C(t) (s + 0.5) tanh (4s) if s < 0.5 C(t) (s -0.5) tanh (4s) if s ≥ 0.5 (22)
C(t) = a sign[ω(t)] ( 23
)
where a is a parameter to be chosen. Notice that with this choice for C(t) the force produced by the disk instantaneously changes sign in the same way ω(t) does. the bottom part shows results for a = 2 × 10 -3 Pa (the same value that corresponds to the contours showed in figure 15). For each plot, the mean value, standard deviation (std) and maximum values over time of these errors are indicated. A superior behavior of the stabilized formulation at almost all times is observed, except for a very few times when a = 2 × 10 -3 Pa.

From this numerical evidence we conclude that the proposed stabilization with the ASGS formulation allows to better capture pressure discontinuities on immersed interfaces in much more general situations than those originally considered in [START_REF] Buscaglia | Finite element methods for the Stokes system with interface pressure discontinuities[END_REF] and [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF].

Thermocapillary Marangoni migration of a droplet

To illustrate the versatility of the method we finally move to an application of interest in micro-fluidics which is the Marangoni migration of a spherical immiscible droplet induced by a varying surface tension coefficient. In practice, such a variation may be the result of heterogeneities in the temperature field or the concentration of surfactant agents. Inertial effects are neglected, which is a reasonable assumption due to the typically small size of the system. The problem to be solved reads

               -µ∆u + ∇ p = 0 in Ω \ Σ div u = 0 in Ω σ • ν = f Σ on Σ u(x → ∞) = 0, p(x → ∞) = p ∞ , (24)
When the surface tension coefficient γ depends on x, the singular force is given by

f Σ = -γ κν + ∇ S γ (25)
where κ is the mean curvature and ∇ S γ = P • ∇γ stands for the surface gradient of γ, P (:= Iν ⊗ ν) being the tangent projector onto Σ. The normal component of the singular force is thus ϕ = -γ κ. The force in [START_REF] Mikkelsen | Actuator disc methods applied to wind turbines[END_REF] has also a non-zero tangential component leading to a jump of the velocity gradient on Σ.

A linear variation of the surface tension coefficient is assumed, namely,

γ = γ 0 -γ 1 x 1 ( 26
)
where γ 0 and γ 1 are constants. Therefore, the tangential component of the singular force equals γ 1 [(1 -

ν 2 1 )ǐ -ν 1 ν 2 ],
where ǐ and  are unitary vectors along the x 1 and x 2 directions respectively. A spherical droplet of radius R is placed at the center of the domain. Since inertia is neglected the droplet instantaneously exhibits a rigid translation to regions of lower γ at constant speed given by Figure 17 shows the problem setting and relevant definitions. An analytical solution for this problem can be found in [START_REF] Balasubramaniam | An-Ti Chai, Thermocapillary migration of droplets: an exact solution for small Marangoni numbers[END_REF] and its numerical resolution by the finite element method was tackled in [START_REF] Buscaglia | Variational formulations for surface tension, capillarity and wetting[END_REF]. Setting v = u -U ǐ the velocity and the pressure fields in each region read

U = 2 15
γ 1 R µ (27) t = 5 s t = 3 s t = 1 s
             v + r = -U cos θ 1 - R 3 r 3 v + θ = U sin θ 1 + R 3 2 r 3 p + = p ∞ if r > R ( 28
)
10 0

             v - r = 3 2 U cos θ 1 - r 2 R 2 v - θ = 3 U sin θ r 2 R 2 - 1 2 p -= p ∞ + 2 R (γ 0 -γ 1 r cos θ) if r < R (29)
where the superscript + refers to the region outside the droplet and the superscript -refers to the region inside it and (r, θ) are the polar coordinates (see Figure 17). The velocity is continuous along Σ h , and its radial component for r = R is identically zero, thus preserving the droplet's shape. The pressure in turn exhibits a jump across Σ depending linearly on x 1

x 2

x 1 u = 0 u = 0 x 2 = 10 cm

x 1 = 10 cm ω(t) = 5π 2 8 cos(10πt) Σ(t) u = 0 u = 0 4 cm
= 2 R (γ 0 -γ 1 x 1 ) (30)
The singular force on Σ h is directly enforced considering an uniform mean curvature κ = 2/R along the interface. This is to avoid any numerical error related to the computation of κ that may pollute the results. The momentum equation in variational form thus reads

µ(∇ u h , ∇ v) Ω h -(p h , div v) Ω h = (-γ κ, v • ν) Σ + (∇ S γ, v) Σ ( 31
)
This approach is different from the one used in [START_REF] Buscaglia | Variational formulations for surface tension, capillarity and wetting[END_REF] in which the identity

Σ (-γ κ ν + ∇ S γ) • v dΣ = - Σ γ ∇ S • v dΣ = - Σ γ P : ∇ v dΣ (32)
is used (see [START_REF] Buscaglia | Variational formulations for surface tension, capillarity and wetting[END_REF] for more details). We solve this problem with an axisymmetric code. The computational domain is the rectangular region Ω = [0, L] × [0, H], where L = 12 R and H = 4R. We consider R = 0.25, µ = 1, γ 0 = 1.5, γ 1 = 1 and p ∞ = 0. The symmetry axis is considered to be the line x 2 = 0. On the rest of ∂Ω h , we prescribe the exact solution ( 28)-( 29) as Dirichlet conditions on u h and p h when assembling the discrete linear system. The computational domain is initially discretized with 6000 triangles to which we assign a mesh size h = 0.03334. From this, a sequence of meshes is built by uniform refinement. We display in Figure 18 streamlines for the velocity v (= u h -U ǐ) and contours of velocity and pressure.

Figure 19 shows the pressure jump along the interface using the proposed method. for the second mesh of the sequence and the surface tension force ϕ imposed. Figure 20 shows the L 2 (Ω h ) error norm of pressure and the H 1 (Ω h ) error norm of velocity. The error in the pressure field behaves clearly better in the stabilized case. The observed convergence rate is nearly an O(h 1.5 ). As for the velocity error, the results are almost identical. The suboptimal order observed in this case (∼ h 0.5 ) is a consequence of the discontinuity in the velocity gradient (cf. Figure 18) as a result of a non-zero tangential singular force on Σ. The discrete space for velocity is simply made up of standard piecewise linear continuous functions, with no special treatment in elements cut by the interface, as in the case of the pressure, for which the space proposed in [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF] is being employed. As already mentioned, a possibility to deal with such kinks is the use of enrichment functions near Σ h or the XFEM method (see e.g. [START_REF] Belytschko | Arbitrary discontinuities in finite elements[END_REF]). Statically condensable enrichments to tackle problems with weak discontinuities have also been proposed recently in [START_REF] Idelsohn | Elemental enriched spaces for the treatment of weak and strong discontinuous fields[END_REF], although, the authors only considered examples involving scalar unknowns. Applied force ϕ(θ)

Figure 19. Singular force and pressure jump along Σ for the thermocapillary migration of a droplet for the the second mesh of the sequence.

Conclusions

In this work we have proposed and extensively tested a finite element formulation for viscous incompressible flows in two dimensions, designed to simulate pressure embedded discontinuities along interfaces, owing, for instance, to actuator disks and surface tension effects. The new finite element method is a combination of the ASGS proposed by Codina [START_REF] Codina | A stabilized finite element method for generalized stationary incompressible flows[END_REF] combined with the pressure jump stabilization technique studied in [START_REF] Buscaglia | Finite element methods for the Stokes system with interface pressure discontinuities[END_REF] and [START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF] (there, from the strict theoretical point of view) as applied to the Stokes system. The pressure jump stabilization is easy to implement in any standard finite element code, since it is based on the standard P 1 -i.e. continuous piecewise linear -interpolation of both velocity and pressure, with a local modification of the pressure representation proposed in [START_REF] Ausas | An improved finite element space for discontinuous pressures[END_REF], allowing for discontinuities along interfaces arbitrarily located with respect to the mesh. This formulation contains the basic ingredients considered in [START_REF] Buscaglia | Finite element methods for the Stokes system with interface pressure discontinuities[END_REF][START_REF] Buscaglia | Approximation of the multi-dimensional Stokes system with embedded pressure discontinuities[END_REF], that is, a stabilization depending on a strictly positive parameter δ, in order to circumvent the Babuška-Brezzi condition violated by the P 1 -P 1 combination, but is much more general and allows to study more challenging problems such as those governed by the Navier-Stokes equations, including inertial and transient effects The pressure jump stabilization term depends on a positive real parameter ε which is basically proportional to the mesh step size in order to reinforce the relation stating that the pressure jump across the interface is either a prescribed force or a part of the solution process itself. It is noteworthy that in all cases that we have studied, the observed convergence rate of the pressure in the L 2 -norm is significantly greater than the theoretical O(h), h being the maximum mesh step size. This rate even attained in some cases the theoretical O(h 1.5 ) rate of the interpolation error using our discontinuous pressure space. The main focus of the numerical experiments has been the simulation of actuator disks. First, the convergence properties of the scheme have been assessed in an academic example and later on, considering a more realistic situation, the accuracy and robusteness of the method has been demonstrated in problems involving non-negligible inertial and transient effects. The method was also tested on a totally different field, by considering a problem of interest in microfluidics involving the thermocapillary Marangoni migration of a droplet and obtained an improved convergence rate for pressure with respect to the unstabilized case and no order erosion for velocity. The main conclusion of our numerical study is that the proposed method is always more effective in the computation of a pressure field accurately fulfilling the jump condition. By using the proposed formulation, more accurate numerical solutions can be expected without any substantial increase of computational effort and any deterioration of the velocity field.

Figure 1 .

 1 Figure 1. Local construction of Σ h .

Figure 2 .

 2 Figure 2. Sub-division of a triangle cut by Σ h into three triangles

Figure 4 .

 4 Figure 4. Setting for the actuator disk problem and the coarsest mesh used for discretization.

Figure 5 .

 5 Figure 5. Section of pressure field along the line x 2 = H/2 = 0.5 in Case I of the actuator disk problem at the coarsest mesh level considering the overstabilized case (Cε ≫ 1), the understabilized case (Cε ≪ 1) and the proposed value Cε = 2.

Figure 6 .

 6 Figure 6. Pressure and velocity field in Case II of the actuator disk problem for the second mesh level.

  Figure 9 shows pressure and velocity for the simulated flow field. These results correspond to a constant value of C in

Figure 7 .

 7 Figure 7. Singular force and pressure jump along Σ in Case II of the actuator disk problem for the interface stabilized and unstabilized case at the second refinement level.

Figure 8 .

 8 Figure 8. Problem setting to study the influence of two inclined actuator disks on the flow around a square obstacle taking into account inertial and transient effects.

Figure 9 .

 9 Figure 9. Pressure field (in Kg/ (cm•s 2 )) and velocity field (in cm/s) illustrating the influence of two inclined actuator disks inducing stationary flow around a square obstacle at a Reynolds number Re ≈ 22.5.

1 × 10 - 4 2 × 10 - 4 3 × 10 - 4 4 × 10 1 Figure 10 .

 110410410410110 Figure 10. Pressure jump and applied force along the left interface for the actuator disks inducing stationary flow around a square obstacle by considering inertia.

2 [cm • s - 1 ]Figure 11 .

 2111 Figure 11. Magnitude of velocity field u h 2 on the left interface for the actuator disks inducing stationary flow around a square obstacle by considering inertia.

Figure 12 .

 12 Figure 12. Pressure and velocity field at different times for the actuator disks inducing transient flow around a square obstacle by considering inertial effects.

1 Figure 13 .

 113 Figure 13. Pressure jump compared to the applied force (top) and magnitude of velocity field u h 2 (bottom) along the left interface at two different times for the actuator disks inducing transient flow around a square obstacle by considering inertial effects.

Figure 14 .

 14 Figure 14. Problem setting for the rotating actuator disk.

  p

20 tFigure 15 . 6 stda = 2 × 10 - 4 mean = 2 .77 × 10 - 4 max = 9 .48 × 10 - 4 std = 1 . 12 × 10 - 4 mean = 1 .39 × 10 - 4 max = 3 .42 × 10 - 4 std- 4 Figure 16 .

 2015621042104910411210411043104416 Figure 15. Pressure and velocity field at different times for the rotating disk problem taking a = 2 × 10 -3 Pa for amplitude of the force in Eq. (23).

Figure 17 .

 17 Figure 17. Problem setting for the thermocapillary (Marangoni) migration of a droplet.

Figure 18 .

 18 Figure 18. Contours of velocity v = u h -U ǐ (top) and pressure p (bottom) for the thermocapillary migration of a droplet.

= 2 . 4 ,

 24 ∼ h 0.5462 Unstabilized, ∼ h 0.5465

Figure 20 .

 20 Figure 20. L 2 (Ω h ) error of pressure (left) and H 1 (Ω h ) error of velocity (right) for the thermocapillary migration of a droplet for the unstabilized method (red) and stabilized method (blue). The estimated convergence rate in each case is also indicated.

Table 1

 1 Error norms for pressure and velocity in Case I (interface Σ occupying the entire channel section) of the actuator disk problem.

		Unstabilized	Stabilized
		p -p h 0,Ω h	u -u h 1,Ω h	p -p h 0,Ω h	u -u h 1,Ω h
	h/4	0.00043	0.00181	0.000134	0.00181
	h/2	0.00123	0.00383	0.000357	0.00362
	h	0.00353	0.00727	0.00096	0.00724
	Order	1.515	1.002	1.424	1.000

Table 2

 2 Error norms for pressure and velocity in Case II (interface Σ occupying half of the channel section) of the actuator disk problem.

		Unstabilized	Stabilized
		p -p h 0,Ω h	u -u h 1,Ω h	p -p h 0,Ω h	u -u h 1,Ω h
	h/8	-	-	-	-
	h/4	0.00223	0.00216	0.00219	0.00214
	h/2	0.00521	0.00424	0.00546	0.00417
	h	0.01341	0.00829	0.02131	0.00814
	Order	1.223	0.980	1.318	0.962

Acknowledgements

The authors gratefully acknowledge the financial support received from both FAPESP (grants 2014/19249-1, 2013/07375-0 -Cepid-CeMEAI) and CNPq (grants 307996/2008-5, 447607/2014-6 and 308728/2013-0).