Nicolas Boutry

Thierry Géraud

Laurent Najman
email: l.najman@esiee.fr

On making nD images well-composed by a self-dual local interpolation

Keywords: Digital topology, gray-level images, well-composed sets, wellcomposed images

Natural and synthetic discrete images are generally not wellcomposed, leading to many topological issues: connectivities in binary images are not equivalent, the Jordan Separation theorem is not true anymore, and so on. Conversely, making images well-composed solves those problems and then gives access to many powerful tools already known in mathematical morphology as the Tree of Shapes which is of our principal interest. In this paper, we present two main results: a characterization of 3D well-composed gray-valued images; and a counter-example showing that no local self-dual interpolation satisfying a classical set of properties makes well-composed images with one subdivision in 3D, as soon as we choose the mean operator to interpolate in 1D. Then, we briefly discuss various constraints that could be interesting to change to make the problem solvable in nD.

Introduction

Natural and synthetic images are usually not well-composed. This fact raises many topological issues. As an example, the Jordan Separation theorem, stating that a simple closed curve in R 2 separates the space in only two components is not true anymore for binary 2D discrete images [START_REF] Latecki | Well-composed sets[END_REF]. To solve this problem, we have to juggle with two complementary connectivities: 4 for the background and 8 for the foreground, or the inverse. 2D well-composed binary images have the fundamental property that 4-and 8-connectivities are equivalent. Hence, such topological issues vanish. In the same manner, well-composed nD images, with n > 2, have 2n-and (3 n -1)-connectivities equivalent [START_REF] Rosenfeld | Connectivity in digital pictures[END_REF]. Other advantages of well-composed images are for example the preservation of topological properties by a rigid transform [START_REF] Ngo | Topology-preserving rigid transformation of 2d digital images[END_REF], simplification of thinning algorithms [START_REF] Marchadier | Thinning grayscale well-composed images[END_REF] and simplification of graph structures resulting from skeleton algorithms [START_REF] Latecki | Well-composed sets[END_REF]. Also, and it is our most important goal, one can compute the Tree of Shapes [START_REF] Najman | Discrete set-valued continuity and interpolation[END_REF][START_REF] Géraud | Self-duality and discrete topology: Links between the morphological tree of shapes and well-composed gray-level images[END_REF] of a wellcomposed image with a quasi-linear algorithm [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF]. An introduction to the Tree of Shapes in the continuous case can be found in [START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF].

Section 2 recalls the definitions of 2D and 3D well-composed sets and grayvalued images, and introduces a characterization of 3D gray-valued well-composed images. Because we do not want to deteriorate the initial signal, we use an interpolation that produces a well-composed image. Furthermore, in order to treat in the same manner bright components on dark background and dark components over bright background, this interpolation process will be self-dual. We present in Section 3 a general scheme that recursively defines local interpolations satisfying a classical set of properties with one subdivision. We show that such interpolations, with the added property of being self-dual, fail in 3D (and then further) to make well-composed images. We conclude in Section 4 with some perspectives that could work in nD even if n > 2 (in local and non-local ways).

A characterization of 3D well-composed gray-valued images

2D WC Sets and Gray-Valued Images

Let us begin by the definitions of a block of Z n . We will then be able to recall the definition and the characterization of 2D well-composed sets and images.

A block in nD associated to z ∈ Z n is the set S z defined such that

S z = z ′ ∈ Z n ||z -z ′ || ∞ ≤ 1 and ∀i ∈ [1, n], z ′ i ≥ z i (where z i represents the i th coordinate of z). Moreover, we call blocks of D ⊆ Z n any element of the set S z ∃z ∈ D, S z ⊆ D . Definition 1 (2D WC Sets [5]) A set X is weakly well-composed if any 8- component of X is a 4-component. X is well-composed if both X and its com- plement X c = Z 2 \ X are weakly well-composed.

Proposition 1 (Local Connectivity and No Critical Configurations [5])

A set X ⊆ Z 2 is well-composed iff it is locally 4-connected. Also, a set X is wellcomposed if none of the critical configurations 1 0 0 1 or 0 1 1 0 appears in X.

Definition 2 (Cuts in nD) For any λ ∈ R and any gray-valued map u :

D ⊆ Z n → R, we denote by [u > λ] and [u < λ] the sets [u > λ] = M ∈ D u(M ) > λ and [u < λ] = M ∈ D u(M ) < λ .
We call them respectively upper strict cuts and lower strict cuts [START_REF] Géraud | A quasi-linear algorithm to compute the tree of shapes of n-D images[END_REF].

We remark that an image u :

D ⊆ Z 2 → R with a finite domain D can only be well-composed if D is itself well-composed (since [u < max(u) + 1] = D).
Definition 3 (Gray-valued WC 2D Maps [START_REF] Latecki | Well-composed sets[END_REF]) A gray-level map u : Z 2 → R is well-composed iff for every λ ∈ R, the strict cuts [u > λ] and [u < λ] result in well-composed sets.

We recall that the interval value of the couple (x, y) ∈ R 2 is defined as intvl(x, y) = [min(x, y), max(x, y)].

Proposition 2 (Characterization of 2D WC maps [START_REF] Latecki | Well-composed sets[END_REF]) A gray-level map u : Z 2 → R is well-composed iff for every 2D block S such that u S = a b c d , the interval values satisfy intvl(a, d) ∩ intvl(b, c) = ∅. As we will see, for n = 3, the equivalence between local connectivity and wellcomposedness is no longer true. This led Latecki [START_REF] Latecki | 3D well-composed pictures[END_REF] to introduce the continuous analog.

3D WC Sets and Gray-Valued Maps

Definition 4 (CA and bdCA [START_REF] Latecki | 3D well-composed pictures[END_REF]) The continuous analog CA(z) of a point z ∈ Z 3 is the closed unit cube centered at this point with faces parallel to the coordinate planes, and the continuous analog of a set X ⊆ Z 3 is defined as CA(X) = CA(x) x ∈ X . The (face) boundary of the continuous analog CA(X) of a set X ⊆ Z 3 is noted bdCA(X) and is defined as the union of the set of closed faces each of which is the common face of a cube in CA(X) and a cube not in CA(X).

Definition 5 (Well-composedness in 3D [START_REF] Latecki | 3D well-composed pictures[END_REF]) A 3D set X ⊆ Z 3 is wellcomposed iff bdCA(X) is a 2D manifold, i.e., a topological space which is locally Euclidian. [START_REF] Latecki | 3D well-composed pictures[END_REF]) A set X ⊆ Z 3 is well-composed iff the following critical configurations of cubes of type 1 or of type 2 (modulo reflections and rotations) do not occur in CA(X) or in CA(X c ) (see Figure 1).

Proposition 3 (No Critical Configurations

We remark that if a set X ⊆ Z 3 is well-composed, then X is locally 6connected. The converse is false (see Figure 2).

Definition 6 (WC Gray-valued Maps)

We say that a 3D real-valued map u :

D ⊆ Z 3 → R is well-composed if its strict cuts [u > λ] and [u < λ], ∀λ ∈ R, are well-composed.
(P1) (P3) (P2) (P4) (P5) (P10) (P9) (P7) (P8) (P6)
Fig. 4. The ten characteristical relations of wellcomposedness of a gray-valued image u restricted to a 3D block S.

To characterize 3D gray-level well-composed images, we first give two lemmas concerning the detection of the critical configurations of respectively type 1 and type 2.

Lemma 1

The strict cuts [u > λ] and [u < λ], λ ∈ R, of a gray-valued image u defined on a block S, such as depicted in Figure 3, do not contain any critical configurations of type 1 iff the six following properties hold:

intvl(a, d) intvl(b, c) = ∅ (P 1), intvl(e, h) intvl(g, f ) = ∅ (P 2) intvl(a, f ) intvl(b, e) = ∅ (P 3), intvl(c, h) intvl(g, d) = ∅ (P 4) intvl(a, g) intvl(e, c) = ∅ (P 5), intvl(b, h) intvl(f, d) = ∅ (P 6)
Proof : Let us assume that one of these properties (P i ), i ∈ [START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF][START_REF] Latecki | Well-composed sets[END_REF], is false. Let us say it is the case of (P 1). Then two cases are possible: either max(a, d) < min(b, c), and that means that there exists λ = (max(a, d) + min(b, c))/2 such that [u < λ] contains the critical configuration {a, d} (of type 1), or min(a, d) > max(b, c), and there exists λ = (min(a, d)+max(b, c))/2 such that one more time [u > λ] contains the critical configuration {a, d}. The reasoning is the same for all the other properties. Conversely, let us assume that there exists λ ∈ R such that either [u > λ] or [u < λ] contains a critical configuration of type 1. That means immediately that one of the 6 properties P i , i ∈ [START_REF] Caselles | Geometric Description of Images as Topographic Maps[END_REF][START_REF] Latecki | Well-composed sets[END_REF], corresponding to each of the six faces of the block S, is false (see Figure 4 for the faces and the corners concerned by the properties).

⊓ ⊔ Recall that the span of a set of values Proof : Let us assume that one of these properties (P i ), i ∈ [7, 10], is false. Let us say it is the case of (P 7). Then two cases are possible:

E ⊆ R is span(E) = [inf(E), sup(E)].
-either max(a, h) < min(b, c, d, e, f, g). Then there exists λ = (max(a, h) + min(b, c, d, e, f, g))/2 such that [u < λ] contains the critical configuration {a, h} (of type 2), -or min(a, h) > max(b, c, d, e, f, g). Then there exists λ = (max(b, c, d, e, f, g) + min(a, h))/2 such that (again) [u > λ] contains the critical configuration {a, h}. The reasoning is the same for all the other properties. Conversely, let us assume that there exists λ ∈ R such that either [u > λ] or [u < λ] contains a critical configuration of type 2. That means immediately that one of the 4 properties P i , i ∈ [START_REF] Levillain | Writing reusable digital topology algorithms in a generic image processing framework[END_REF][START_REF] Ngo | Topology-preserving rigid transformation of 2d digital images[END_REF], corresponding to each of the four diagonals of the block S, is false (see Figure 4). ⊓ ⊔ We are now ready to state the main theorem of this section, characterizing the well-composedness on a 3D gray-valued image.

Theorem 1 (Characterization of well-composedness in 3D) Let us suppose that D is a hyperrectangle in Z 3 . A gray-valued 3D image u : D → R is well-composed on D iff on any block S ⊆ D, u S satisfies the properties (i+½,j+½) (

(P i ), i ∈ [1, 10].
(i+½,j) (i+1,j+½) (i,j+½) (i+½,j+1) Using interpolations with one subdivision does not deteriorate the initial signal. The size of the original image is multiplied by a factor of 2 n , where n is the dimension of the space of the image. Figure 5 illustrates this subdivision process.

Subdivision of a Domain and Z

n as a poset

The subdivision of a block allows us to provide an order to the elements. Using this order, the subdivided space is a poset.

Let z be a point in Z n , and S z its associated block. We define the subdivision of S z as s(S z ) = {z ′ ∈ Z subdivision of a domain D ⊆ Z n is the union of the subdivisions of the blocks contained in D, i.e., s(D) = S⊆D s(S).

Definition 7 (Order of a point of Z 2 n ) Assume e i is a fixed basis of Z n . We note 1 2 (z) = {i ∈ [1, n] z i ∈ Z 2 \ Z}. The sets ❊ k , for k ∈ [0, n], are defined such that ❊ k = {z ∈ Z
(z) = ∪ i∈ 1 2 (z) {z -ei 2 , z + ei 2 }. The parents of z ∈ Z 2
n of order 0 are P 0 (z) = {z} and of order k > 0 are defined recursively by P k (z) = ∪ p∈P(z) P k-1 (p). Definition 9 (G(z) and ❆(z)) Let z be a point of

Z 2 n . The ancestors of z ∈ Z 2 n are ❆(z) = P ♦(z) (z). We set G(z) = ∪ k∈[0,♦(z)] P k (z). Notice that ❆(z) ⊆ Z n and that any point z ∈ ❊ k , k ∈ [1, n], has its parents in ❊ k-1 . Hence {❊ k } k∈[0,n] is a (hierarchical) partition of Z 2 n , and ( Z 2 n , P)
is a poset (see Figure 6).

Definition 10 (Opposites) Let z be a point of Z 2 n . The (set of ) opposites of z ∈ Z 2 n is the family of pairs of points opp(z) = ∪ i∈ 1 2 (z) {z -ei 2 , z + ei 2 } .

Interpolations with one subdivision

Let us recall that the convex hull convhull(Z) of a set of m points Z = {z 1 , . . . , z m } ⊆ Z n is:

convhull(Z) = m i=1 α i z i m i=1 α i = 1 and ∀i ∈ [1, m], α i ≥ 0
Definition 11 (Subdivision of edges, faces, and cubes)

Let E = {z 1 , z 2 } be an edge in Z n . The subdivision of E is s(E) = {z ∈ Z 2 n z ∈ convhull(E)}.
The subdivision of a face

F = {z 1 , z 2 , z 3 , z 4 } is s(F) = {z ∈ Z 2 n z ∈ convhull(F)}.
The subdivision of a cube

C = {z 1 , . . . , z 8 } is s(C) = {z ∈ Z 2 n z ∈ convhull(C)}.

A set of properties that an interpolation has to satisfy

An interpolation of a map u :

D ⊆ Z n → R to a map I(u) : s(D) ⊆ Z 2 n → R
is a transformation such that I(u) S = u S for any block S ⊆ D. Let u : D ⊆ Z 3 → R be any 3D gray-valued image. We say that an interpolation I : u → I(u) is self-dual iff I(-u) = -I(u). A self-dual interpolation does not overemphasize bright components at the expense of the dark ones, or conversely.

An interpolation I : u → I(u) in 3D is said ordered if the new values are inserted firstly at the centers of the subdivided edges, secondly at the centers of the subdivided faces, and finally at the centers of the subdivided cubes.

An ordered interpolation is said in between iff it puts the values at a point z in between the values of its parents P(z).

Finally, we say that an interpolation is well-composed iff the image I(u) resulting from the interpolation of u is well-composed for any given image u.

We are interested in interpolations I with the following properties.

(P) ⇔            I is invariant by translations, π
2 's rotations and axial symmetries

I is ordered I is in-between I is self-dual I is well-composed 3.

The scheme of local interpolations verifying P

A local interpolation I is an interpolation such as for any block S ⊆ D, I(u) on s(S) is computed only from its nearest neighbours belonging to ❊ 0 (we see an image as a graph). For convenience, we will write u ′ instead of I(u) for local interpolations in the sequel. verifying P can be characterized by a set of functions {f k } k∈ [1,n] such that:

∀ z ∈ Z 2 n , u ′ (z) = u(z) if z ∈ ❊ 0 f k (u ❆(z) ) if z ∈ ❊ k , k ∈ [1, n]
We denote such an interpolation I f1,...,fn .

Proof : The interpolation process on the subdivided edges depends only on the values of u at the vertices of the original edges due to the locality of the method. Furthermore it has to be invariant by axial symmetries and rotations.

Hence, there is a unique function f 1 characterizing the interpolation on the subdivided edges. The reasoning is the same on the faces and the cubes respectively for f 2 and f 3 .

⊓ ⊔ Notice that it is an implication and not an equivalence: an interpolation verifying this scheme does not verify all the properties in P.

I 0 (u, z) (def ) = {z -,z + }∈opp(z) intvl(u ′ (z -), u ′ (z + ))
Then, an ordered local interpolation I is said in-between iff u ′ (z) ∈ I 0 (u, z) for any image u : D → R and z ∈ s(D) \ ❊ 0 .

Definition 13 (I W C and I sol ) Let u : D → R be a gray-valued image, z be a point of s(D) \ ❊ 0 , and I be a given local interpolation. We define the set I W C (u, z) associated to I such as for any z ∈ ❊ 1 , I W C (u, z) = R and for any z ∈ ❊ k with k ≥ 2:

I W C (u, z) = { v ∈ R | u ′ (z) = v ⇒ u ′ G(z) is well-composed } Last, let us denote I sol (u, z) = I 0 (u, z) ∩ I W C (u, z).
The following scheme is necessary to satisfy P (but not sufficient).

Theorem 2 Any local interpolation I satisfying P is such that:

∀z ∈ Z 2 n , u ′ (z) = u(z) if z ∈ ❊ 0 f k (u ❆(z) ) ∈ I sol (u, z) if z ∈ ❊ k , k ∈ [1, n]
Notice that such a local interpolation I is ordered, in-between, well-composed, but not necessarily self-dual.

Determining f 1 for self-dual local interpolations

Let us begin with the study of f 1 , i.e., the function setting the values at the centers of the subdivided edges. This function has to be self-dual, symmetrical, and in-between. We choose one of the most common function satisfying these constraints: the mean operator f Lemma 4 Let z be a point in ❊ 2 . Modulo reflections and symmetries, an αconfiguration implies that u ❆(z) is not well-composed, whereas a U -or Z-configu- ration implies that u ❆(z) is well-composed. Let us begin with the computation of I 0 (u, z) for z ∈ ❊ 2 . From the values already set in u ′ on P(z) ⊆ ❊ 1 by f 1 during the recursive process (see Fig-

1 : R 2 → R : (v 1 , v 2 ) → f 1 (v 1 , v 2 ) = (v 1 + v 2 )/2.

Fig. 1 .

 1 Fig. 1. Illustration of the bdCA of a set containing a critical configurations of type 1 (left), and of type 2 (right).

Fig. 2 .

 2 Fig.2. A set locally 6-connected but not well-composed.

Fig. 3 .

 3 Fig. 3. The restriction u S of u to a 3D block S.

Lemma 2

 2 The strict cuts [u > λ] and [u < λ], λ ∈ R, of a gray-valued image u defined on a block S such as depicted in Figure3, do not contain any critical configuration of type 2 iff the four following properties are true: intvl(a, h) span{b, c, d, e, f, g} = ∅ (P 7) intvl(b, g) span{a, c, d, e, f, h} = ∅ (P 8) intvl(c, f ) span{a, b, d, e, g, h} = ∅ (P 9) intvl(d, e) span{a, b, c, f, g, h} = ∅ (P 10)

Fig. 5 .

 5 Fig. 5. Illustration of the subdivision process on a block S.

Fig. 6 .

 6 Fig. 6. s(S) ⊆ Z 2 n as a

2 n 1 2 2 n 2 n) 2 n 2 n

 212222 (z) = k} (where E denotes the cardinal of the set E), and represent a partition of Z . We call order of a point z the value k such that z ∈ ❊ k and we note it ♦(z).Definition 8 (Parents in Z Let z be a point of Z . The set of the parents of z ∈ Z , noted P(z), is defined by P

Fig. 7 .

 7 Fig.7. From left to right: an image, and its interpolations with the median, the mean/median, the min and the max.

3. 5

 5 I 0 , I W C , and I sol for local interpolations Let us introduce some useful sets to express recursively the local interpolations satisfying the properties P. Definition 12 (I 0 and definition of a local in-between interpolation) Let u : D → R be a gray-valued map, let z be a point of s(D) \ ❊ 0 , and let I be a given local interpolation. We define the set I 0 (u, z) associated to I by:

Fig. 8 .

 8 Fig. 8. The 3 possible configurations in 2D (modulo reflections and rotations).

2 Fig. 9 .

 29 Fig. 9. u ′ G(z) for z ∈ ❊2 for any self-dual local interpolation after the application of f1 (with m any value ∈ R).

3. 7

 7 Equations of f 2 for self-dual local interpolations Concerning f 2 , i.e., the function which sets the values of u ′ at the centers of the subdivided faces, let us compute I 0 (u, z) and I W C (u, z) for any given z ∈ ❊ 2 to deduce I sol (u, z). Their values depend on the configurations of u ❆(z) . Let us assume that u ❆(z) = a b c d . Then a total of 4! = 24 increasing orders are possible for these 4 values. Modulo reflections and axial symmetries, we obtain a total of 3 possible configurations: the α-configurations correspond to the relation a ≤ d < b ≤ c, the U -configurations to a ≤ b ≤ d ≤ c, and the Z-configurations to a ≤ b ≤ c ≤ d (see Figure 8).

Fig. 10 .

 10 Fig. 10. The Hasse diagrams for the α-and the U -configurations (left) and for the Z-configuration (right).

n ||zz ′ || ∞ ≤ 1 and ∀i ∈ [1, n], z ′ i ≥ z i }. The

ure 9), we can compute I 0 (u, z) using the Hasse diagram 3 for each configuration (see Figure 10). We obtain finally that I 0 (u, z) = intvl( a+c 2 , b+d 2 ) for the three configurations, with one important property: the median value of u ❆(z) always belongs to I 0 (u, z).

Let us follow with the computation of I W C (u, z), where u ′ G(z) (see Figure 9) satisfies four conditions:

In the case of the α-configuration, (2) ⇒ m ≤ b+d 2 and (4) ⇒ m ≥ b+d 2 . That implies that m = b+d 2 , which also satisfies ( 1) and (3). Consequently,

{med{u ❆(z) }}, and because I W C (u, z) ⊆ I 0 (u, z), I sol (u, z) = {med{u ❆(z) } in the not well-composed case.

In the cases of the U -and the Z-configurations, we obtain that

] ⊇ I 0 (u, z), so we conclude that I sol (u, z) = I 0 (u, z).

Theorem 3 Given an image u : D → R, any local interpolation I f1,f2,f3 satis-

otherwise.

Let z be a point in s(D)∩❊ 2 . Amongst the applications f 2 satisfying P, there exists (at least) the median method (see Figure 7), consisting in setting the value of u ′ (z) at med{u ❆(z) } (in this case f 2 is an operator and not only a function), the mean/median method of Latecki [START_REF] Latecki | Well-composed sets[END_REF] consisting in setting the value u ′ (z) at mean{u ❆(z) } in the well-composed case and to med{u ❆(z) } otherwise, and also the min/max method, consisting in setting the value u ′ (z) at 1 2 (min{u ❆(z) } + max{u ❆(z) }) in the well-composed case and to med{u ❆(z) } otherwise.

Equations of f 3 for local self-dual interpolations

Theorem 4 No local interpolation satisfies P for n ≥ 3 with one subdivision as soon as we chose the mean operator to interpolate in 1D.

Proof : Let z be the center of a subdivided cube. We have u ′ ❆(z) as in the Figure 11 (on the left). We apply the first interpolating function f 1 , i.e., we set 3 Recall that a Hasse diagram is used to represent finite partially ordered sets with the biggest elements at the top the values of u ′ at the centers of the subdivided edges at the mean of the values on the vertices. Then we apply the second interpolating function f 2 , which fixes the values of u ′ at the centers of the subdivided faces at the median of the values of u ′ at the four corresponding corners (because u is well-composed on none of the faces of the cube). Finally, referring to the properties that a function u ′ has to satisfy to be well-composed (see theorem 1), f 3 must also satisfy the constraints c ≥ 3 and c ≤ 1 (both are constraints of type 2) that are incompatible. So, no local interpolation of this sort can satisfy the set of constraints P as soon as n ≥ 3. ⊓ ⊔

Conclusion

In this paper, we have presented a characterization of well-composedness for 3D gray-valued images. We proved that no local interpolation satisfying P with one subdivision is able to make 3D well-composed images as soon as we choose the mean operator as interpolation in 1D.

Although our formalism is developped in the continuous domain (the interpolations take their values in R), it is in fact a discrete setting. Indeed, the image u ′ can easily be computed in Z as soon as the space image of u is also Z. We just have to multiply the values of the original image u by a factor k ∈ Z where k depends on the interpolation we use (e.g., k = 2 for the median method and k = 4 for the mean/median method in 2D). Another way to deal with images having values in Z/k is to use a generic image processing library [START_REF] Levillain | Writing reusable digital topology algorithms in a generic image processing framework[END_REF].

Future research should tackle the two following directions. The first direction is to use an alternative to f 1 such as (a, b) → med(a, b, c) (where c is the center of the space of the image u). The second direction is to use a non-local approach, e.g., a front propagation algorithm. In that case, we do not have to use any systematic operator f 1 anymore, nor to use an ordered interpolation. First results in this second direction are promising.


Maxime Vicentini
email: maxime.vicentini@onera.fr

Renaud Lecourt

Olivier Rouzaud

Virginel Bodoc

Olivier Simonin

CONGRÈS FRANCOPHONE DE TECHNIQUES LASER (CFTL2014) MARSEILLE, 15-19 SEPTEMBRE 2014 1 REGIMES DE COMBUSTION DANS LES FOYERS AERONAUTIQUES : MISE EN PLACE EXPERIMENTALE DE L'ECOULEMENT DIPHASIQUE

Keywords: régimes de combustion, diphasique, moyenne de phase, PIV, imagerie rapide, PDA

Introduction

Contexte.

De nos jours, la combustion diphasiqueinjection de carburant liquide dans un foyerest très largement répandue dans de nombreuses applications (turboréacteurs, moteurs-fusées, moteurs automobiles, fours industriels). Lorsque la théorie de l'évaporation des gouttes de carburant a été élaborée, s'est naturellement posée la question de savoir, si, dans un brouillard, les gouttes se consumaient individuellement ou en paquets. Chiu [START_REF] Chiu | Group combustion of liquid fuel sprays[END_REF] et Chigier [START_REF] Chigier | Group combustion models and laser diagnostic methods in sprays : a review[END_REF] ont introduit l'idée qu'il puisse exister différentes structures de λlammes au sein d'un nuage de carburant présentant une distribution spatiale hétérogène de gouttes. La compréhension phénoménologique des régimes de combustion de brouillard de gouttes pourrait être essentielle en vue d'optimiser les systèmes d'injection et l'architecture des chambres de combustion. Par la suite, de récents travaux ont montré l'inλluence de certains réμimes de combustion sur la production d'oxyde d'azote [START_REF] Koch | Investigation of the effect of incomplete droplet evaporation prevaporization on NOx emissions in LDI combustion systems[END_REF].

Objectifs de la Présente Etude.

En vue d'améliorer la modélisation de la combustion diphasique dans les solveurs CFD actuels, l'ONΔRA souhaite mettre en évidence, expérimentalement, l'existence ou non des différents régimes de combustion. Il s'aμirait de proposer une cartographie de ces régimes au moyen de grandeurs physiques et/ou de nombres adimensionnels. Pour y parvenir, un nouveau moyen d'essais, le montaμe PROMÉTHÉE (figure 1), spécifiquement dédié à cette problématique, a été conçu. L'objet de cette présente étude est la caractérisation de l'écoulement diphasique en conditions non-réactives à l'aide de techniques laser (PIV, imagerie rapide, PDA). L'un des points importants de ce travail consiste en la mise en place et la validation du conditionnement de mesures aλin d'eλλectuer des traitements statistiques de type « moyenne de phase ». La principale différence entre le banc d'essais présenté en [START_REF] Koch | Investigation of the effect of incomplete droplet evaporation prevaporization on NOx emissions in LDI combustion systems[END_REF] et ce montage réside dans la « bidimensionnalité » de l'écoulement. Δn eλλet, pour être mieux adaptée aux diagnostics optiques et pour s'aλλranchir des phénomènes complexes associés aux flammes stabilisées par un écoulement swirlé, la configuration a été dessinée de manière à être la plus « bidimensionnelle » possible. Toutefois, les résultats restent transposables à un foyer aéronautique réaliste car la configuration reproduit les principales caractéristiques des écoulements rencontrés dans ces foyers : carburant similaire au kérosène, brouillard de gouttes polydisperse, gamme de richesse, zone de recirculation, échelles de turbulence, confinement. Le premier travail a consisté à définir la forme de l'obstacle servant d'accroche-flamme.

Investigation expérimentale

Concept et

Conditions d'Essais

Dans la littérature, il existe un grand nombre de travaux au sujet d'écoulements turbulents dans lesquels un obstacle génère des lâchers de structures tourbillonnaires qui sont convectés vers l'aval [START_REF] Strzelecki | Etude du détachement tourbillonnaire[END_REF]. Ces études sont souvent, elles-mêmes, rattachées à la problématique des débitmètres à vortex. Par le passé, diλλérentes μéométries ont été examinées. Il s'aμissait d'optimiser la précision de ces appareils et d'estimer le mieux possible la corrélation entre la fréquence connue et le débit recherché. De plus, l'objectiλ était de garantir leur fiabilité sur la plus large plage de nombre de Reynolds en ayant un phénomène de détachement tourbillonnaire aussi régulier que possible.

En général, il a été observé que les obstacles présentant des arêtes vives favorisaient le développement de structures plus stables, plus cohérentes [START_REF] Strzelecki | Etude du détachement tourbillonnaire[END_REF]. En revanche, l'obstacle de type cylindre, par exemple, génère des instabilités supplémentaires, au niveau de la couche limite et dans la couche de cisaillement. C'est pourquoi cette μéométrie n'a pas été retenue.

Les trois formes retenues sont les obstacles rectangulaire, triangulaire et trapézoïdal. L'accrocheflamme rectangulaire a été souvent employé dans les études portant sur les statoréacteurs. L'obstacle triangulaire a été retenu car il s'apparente à l'accroche-flamme « V-Gutter » utilisé notamment pour la postcombustion des avions supersoniques. Le trapèze, à faible angle de découpe (4°), fait également l'objet de cette étude paramétrique car c'est une géométrie courante dans les débitmètres à vortex.

En plus de la géométrie, différents rapports de blocage ont été testés (figure 3).

Le rapport de blocage est défini par la formule (eq. 1) :

β = (eq.1)
où Ao, A sont respectivement la surface maître-couple de l'obstacle et la section de passage de la veine d'essais. Le spectre de l'accroche-flamme triangulaire est très différent du précédent. Cette fois, un pic de fréquence, autour d'un nombre de Strouhal égal à 0,29, ressort plus nettement. En revanche, cette géométrie produit encore un bruit excessif pour des fréquences supérieures à 100 Hz.

Pour la géométrie du trapèze, les deux spectres sont très similaires. Il existe, dans ce cas, un très bon rapport signal sur bruit. La géométrie trapézoïdale est donc retenue pour cette étude. Le trapèze ayant le rapport de blocage le plus important (β = 42%) est toutefois préféré. Celui-ci devrait rendre l'accrochaμe d'une λlamme plus aisé en générant une zone de recirculation plus allongée que le petit trapèze.

Algorithme pour le traitement conditionnel.

L'application simultanée, stricto sensu, de plusieurs techniques laser est une tâche difficile, voire impossible encombrement, accès optique, … . C'est pourquoi, l'un des points importants de ce travail consiste en la mise en place et la validation du conditionnement de mesures par rapport à un signal de référence. Δn eλλet, à partir d'une analyse fréquentielle, des moyennes de phase peuvent être réalisées aλin de caractériser l'instationnarité de l'écoulement. Δn outre, il permet de superposer / synchroniser, a posteriori, différentes mesures optiques.

Dans cette étude, deux signaux permettent de conditionner les mesures. Le

Traitement du Flot Optique

Le traitement d'imaμes est réalisé μrâce à l'alμorithme FOLKI-SPIV développé à l'ONΔRA [START_REF] Sant | [END_REF]. Cet algorithme permet une estimation précise d'un flot optique dense en minimisant la différence quadratique entre deux images consécutives. Sa technique consiste en un recalage de fenêtres de type Lucas-Kanade itératif et multi-résolution. Une autre caractéristique importante est sa structure hautement parallèle permettant de profiter pleinement de la puissance et de la rapidité des processeurs graphiques (GPU) [START_REF] Sant | [END_REF].

Pour cette étude les paramètres principaux sont les suivants :

nombre d'itération = -nombre de niveau pyramidal = 3 -taille de fenêtre d'interroμation (Diw) = 31 px -interpolation bilinéaire

Plans de Mesures

On s'intéresse ici uniquement aux plans de mesures médians à la veine d'essais. Deux plans de mesures, de largeur 100 mm chacun, sont étudiés à l'aval de l'obstacle. Les plans ont un taux de recouvrement de plus de 15% pour permettre une superposition des champs lors du post-traitement. A noter que par convention, la composante longitudinale U de la vitesse est selon l'abscisse X, tandis que la composante transversale V est selon l'ordonnée Y.

Histogramme de Vitesses

Pour juger de la qualité des mesures PIV, un histogramme des vitesses est obtenu sur 1880 images, en aval de l'accroche-λlamme, à l'abscisse X = 63 mm (figure 8). L'écoulement en amont de l'obstacle est caractérisé à partir d'un lot de 1000 images. La figure 9 représente le profil de vitesse longitudinale en fonction de la hauteur de la veine d'essais. Une légère dissymétrie peut être observée entre le haut et le bas de la veine d'essais (inférieure à 2 % en amplitude).

Cette dissymétrie disparaît en aval de l'obstacle. L'accélération de l'écoulement, provoquée par un effet de blocage important (β = 42%), rend le proλil de vitesse en aval de l'obstacle peu sensible aux conditions amont [START_REF] Strzelecki | Etude du détachement tourbillonnaire[END_REF]. Les traitements statistiques de type « moyenne de phase » consistent à moyenner toutes les images appartenant à la classe de phase ϕ ± ϕ. Le choix du paramètre ϕ vient du compromis entre la précision d'attribution des phases et le nombre d'imaμes à moyenner.

Dans cette partie, on compare les profils de vitesse obtenus avec diλλérentes valeurs de ϕ en maintenant le nombre d'imaμes moyennées à 100 images. la figure 10, le profil de vitesse longitudinale correspond à la phase 45° ± ϕ, en aval de l'accroche-flamme, à l'abscisse X = 63 mm. Sur la figure 13-a, l'écoulement, présente bien une zone de recirculation centrale dans le sillage de l'obstacle vitesses négatives ou nulles). Sa longueur varie au cours du battement. Sa longueur moyenne vaut 1,14 fois la hauteur du trapèze, soit 57 mm. Sur la 13-b, le critère 2 de Graftieaux [START_REF] Graftieaux | Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows[END_REF] permet d'identiλier les structures tourbillonnaires aux grandes échelles. Conformément aux prédictions, des tourbillons contrarotatifs sont générés de manière périodique. Durant leur phase de développement, on observe une dilatation puis une accélération des vortex vers l'aval, jusqu'à une vitesse de convection proche de la vitesse débitante. Des visualisations tomographiques par imagerie rapide ont été eλλectuées sur l'écoulement. Il s'aμit d'observer, d'un point de vue qualitatiλ, le comportement de la phase gazeuse et de la phase dispersée. Le post-traitement permet de réaliser des moyennes de phase et de superposer ces visualisations (figure 15). Le signal de luminance est utilisé pour le conditionnement des mesures.

Pour cette nouvelle étude, la tranche laser a été générée par un laser pulsé à haute cadence Quantronix Darwin-527 (classe IV) de type Nd :YLF. La lonμueur d'onde du laser est de nm. Une caméra rapide Phantom V341, équipée d'un objectif de focale 90 mm, est placée perpendiculairement à la nappe laser. Le système de synchronisation est réglé pour une acquisition de 500 images par seconde. La résolution de caméra est de 1024x768 px² pour une profondeur d'intensité lumineuse de 12 bit/pix. Comme pour la PIV, le champ de visualisation correspond au sillaμe de l'obstacle sur une lonμueur de 185 mm. Pour les essais en phase gazeuse, l'ensemencement de ΓΔHS est réalisé au culot de l'obstacle via le canal d'injection. Les moyennes de phases ont été réalisées à partir de 110 images pour la phase gazeuse. La taille et la position des tourbillons visibles grâce au traceur sont en adéquation avec les résultats obtenus en PIV, (figure 13-b).

Pour les essais en diphasique, le carburant liquide est injecté via l'injecteur « flat spray ». La carburation de la chambre de visualisation à froid (293 K) produisant un encrassement rapide des faces latérales a empêché l'enreμistrement de lonμues séries d'imaμes. En conséquence, les moyennes de phase n'ont été réalisées que sur des lots de 20 images pour la phase dispersée, ce qui est insuffisant. pour interpréter les résultats avec exactitude. Le seuillage sur l'intensité lumineuse des enregistrements de la phase dispersée fait ressortir différentes organisations quant à la distribution spatiale des gouttes. Les plus petites gouttes semblent être concentrées à l'intérieur des tourbillons. Les plus grosses sont, en général, projetées sur les parois supérieures et inférieures de la veine d'essais. Les gouttes de tailles moyennes semblent être ségrégées de façon préférentielle en périphérie des tourbillons comme le montre la figure 15. Pour confirmer ces résultats, des essais complémentaires telles que des visualisations directes et de l'Interλérométrie Phase Γoppler, à chaud (453 K) sont prévus.

Conclusion

Un nouveau moyen d'essais, le montaμe PROMÉTHÉΔ a été spécialement conçu pour l'étude des régimes de combustion. Une étude paramétrique a permis de retenir une μéométrie d'accrocheflamme. La forme trapézoïdale avec un rapport de blocage de 42 % semble optimale pour obtenir un détachement tourbillonnaire régulier et intense. L'un des points importants de ce travail a consisté en la mise en place et la validation du conditionnement de mesures par rapport à un signal de référence pertinent. Un algorithme de traitement conditionnel a été développé puis appliqué sur des mesures. L'alμorithme permet de réaliser des traitements statistiques de type « moyenne de phase ». Ces moyennes de phase ont permis de caractériser l'instationnarité de l'écoulement sur un cycle complet. En outre, elles ont montré qu'il était de synchroniser a posteriori différentes mesures. Ainsi, des visualisations de la phase gazeuse ont été superposées à des visualisations portant sur la phase dispersée. Des essais complémentaires, à chaud puis en combustion sont prévus dans la suite de ce travail en vue de montrer l'inλluence de la distribution spatiale de gouttes sur la structure des flammes.

Références :

  figure 1 : Présentation du prototype PROMÉTHÉE

  tableau A-Conditions opératoires du montage PROMÉTHÉEChoix de la Géométrie d'Accroche-Flamme.

figure 3 :

 3 figure 3 : Différents accroche-flammes étudiés Les obstacles ont été percés sur leur face supérieure et inférieure de trous de pression (1,7 mm de diamètre). Les trous sont reliés à des canaux à l'intérieur de l'obstacle. Ces canaux permettent l'établissement d'une pression hydrostatique. Celle-ci remonte vers un capteur de pression différentielle de type capteur à membrane Validyne P305D avec une sensibilité ± 0,86 kPa. Lors du régime oscillatoire de détachement tourbillonnaire, une différence de pression est mesurée par le capteur. Le signal résultant est alternatif et idéalement périodique. Pour chacun des accroche-flammes, des analyses spectrales sont réalisées sur le signal de pression. Dans un premier temps, un analyseur de spectre temps réel B&K permet de vérifier que l'écoulement est stabilisé. Δnsuite, le siμnal est enregistré, pendant une durée de 750 s, via un boîtier NI BNC-209 (22 voies), sur une carte d'acquisition PCI-MOI-16E-1. L'analyseur de spectre détecte le détachement tourbillonnaire entre 20 et 50 Hz pour les accroche-flammes retenus. La fréquence d'échantillonnaμe étant fixée à 2000 Hz, le critère de Nyquist-Shannon est donc largement respecté. Sur les figures 4, 5, 6, 7 sont données les densités spectrales normalisées par leur amplitude maximale. Pour chaque spectre, il existe un sommet de fréquence. Ce sommet peut-être associé à un nombre de Strouhal (Sth). Physiquement, le nombre de Strouhal représente le rapport du temps d'advection et du temps caractéristique de l'instationnarité. Il est défini comme suit :

figure 4 :figure 5 :figure 7 :

 457 figure 4 : Densité spectrale de puissance normalisée -Rectangle, β = 25 % -Sth = 0,11

figure 8 :

 8 figure 8 : Densité de probabilité des vecteurs en fonction des déplacements de particules -X = 63 mm L'histoμramme est une λonction μaussienne dont l'espérance mathématique correspond à un déplacement moyen de 3 px. Un critère courant est de vérifier que le déplacement des particules reste inférieur ou égal au quart de la taille de fenêtre d'interroμation Γiw [6]. Dans le cas présent, le critère de déplacement est toujours respecté. Il reste adapté à des zones dans lesquelles les vitesses atteignent 15 m/s.L'enveloppe de l'histogramme présente des petites oscillations avec une période de 1 px. Ce phénomène est typique d'un blocage de pic (peaklocking), récurrent dans ce genre de mesures. Ici, le phénomène existe mais n'est pas prononcé car il n'y a pas d'absence totale de mesure au voisinaμe des valeurs proches du demi-pixel. La qualité des données est donc satisfaisante.

Figure 11 :

 11 Figure 11: figure 9 : Profil U, amont stationnaire L'intensité de turbulence est calculée comme le rapport de la vitesse longitudinale fluctuante (0,27 m/s) sur la vitesse débitante moyenne (6 m/s). Par conséquent, l'intensité turbulente vaut , %. Il s'aμit donc d'une turbulence moyenne.

  figure 10: Profil U, 45 ± ϕ, 100 images, aval, X = 63 mm On constate une différence de vitesse de part et d'autre du sillaμe de l'obstacle. Pour la phase °± ϕ, cette différence produit une zone de cisaillement accompaμnée d'un enroulement

  figure 14 : Vélocimétrie (couleur), granulométrie (taille)

figure 15 :

 15 figure 15 : Phase gazeuse (couleur), phase dispersée (blanc)

en Ecoulement Non-Réactif.

		Le montage PROMÉTHÉE est alimenté par
	deux liμnes distinctes, l'air et le carburant.
		Le débit d'air est contrôlé par un col sonique
	à aiguille. Celui-ci est muni de capteurs analogiques
	calibrés (pression, température, déplacement aiguille).
	Ils permettent, au moyen d'une corrélation, d'estimer
	le débit d'air à % près. Tous les essais présentés ici
	correspondent à une vitesse débitante de 6 m/s dans
	la	chambre	de	visualisation	de	forme
	parallélépipédique carrée 120x120 mm².
		Pour les essais en non-réactiλ, l'éthanol	%
	est retenu comme fluide de similitude au kérosène.
	Ses propriétés physico-chimiques sont relativement
	identiques dans les conditions d'essais déλinies ;
	notamment la tension de surface (2,3 mN/m) et la
	viscosité dynamique (1,2 mPa.s) qui jouent des rôles
	importants dans les phénomènes d'atomisation. Le
	débit massique de carburant est mesuré par un
	débitmètre à effet Coriolis étalonné de la marque
	Micro Motion. Sa précision est d'environ , % de la
	mesure. L'injecteur est alimenté en carburant par un
	réservoir pressurisé à l'air sous bar. Les paramètres
	retenus sont résumés dans le tableau A.

de Vélocimétrie par Imagerie des particules (PIV) et leur validation.

  premier signal est enregistré par le capteur Valydine pendant les essais, en parallèle des mesures optiques. Ce signal correspond aux fluctuations de pression causées par les lâchers tourbillonnaires de part et d'autre de l'obstacle. Le deuxième est un signal obtenu lors du post-traitement des imaμes. Il s'aμit, dans ce cas, des fluctuations de luminance causées par des variations locales en concentration du traceur d'ensemencement ou des gouttes de carburant. Quel que soit le signal brut, le même algorithme est utilisé. Une transformée de Fourier est appliquée au signal temporel, auquel a été retranchée sa valeur moyenne (signal sans composante continue . Au pic d'amplitude de la FFT correspond sa fréquence principale ; par exemple la fréquence des allées de Von Karman. Le signal contenant éμalement d'autres λréquences, il est nécessaire de les éliminer aλin d'obtenir un signal de référence exploitable. Dès lors, le choix d'un λiltre de posttraitement s'impose. La donnée la plus importante étant la phase, il λaut choisir un λiltre n'introduisant aucun décalage temporel. Une autre spécificité attendue est d'avoir un λiltre ayant un μain aussi constant que possible dans sa bande passante pour minimiser les pertes d'inλormations. Pour ces raisons, un filtre passe-bande Butterworth d'ordre est préféré. Le filtre est centré sur la fréquence principale et limité à une bande de ± 5% autour de celle-ci. Le but de ce traitement est d'obtenir un signal filtré où les passages à zéro sont plus faciles à déterminer. Les passages à zéro sont estimés lors de la phase montante (passage du signe négatif au signe positif) par une interpolation linéaire du premier ordre. Ce sont ces passages à zéro qui définissent le temps zéro des cycles. L'intervalle entre deux de ces passages à zéro correspond à la période instantanée d'un cycle. Lorsqu'une imaμe de l'écoulement est prise, un signal TTL est enregistré par la carte d'acquisition. La détection du front montant permet d'associer chaque imaμe à son temps d'apparition. Une comparaison entre le temps d'apparition et deux passages à zéro consécutifs, encadrant l'arrivée de cet événement, permet d'y associer une valeur de phase ϕ comprise entre 0 et 360°. Une fois ce classement en phase réalisé, il est aisé de procéder à des traitements statistiques de type « moyenne de phase », i.e. toutes les images appartenant à l'intervalle ϕ ± ϕ sont moyennées pixel à pixel pour qu'il en résulte une seule imaμe. Le paramètre ϕ est essentiel. Il pilote, d'une part, la précision d'attribution de phase et, d'autre part, le nombre d'imaμes utiles pour effectuer des moyennes.

	une épaisseur d'environ	mm. La nappe laser
	traverse la manchette dans son plan médian. Une
	caméra CCD de résolution 2048x2048 px² est placée
	perpendiculairement à la tranche laser. Un objectif
	CANON d'une λocale de 60 mm est utilisé. Après
	étalonnage sur une mire, la résolution est estimée à
	17,6 px/mm. De même que pour les lasers, la caméra
	est reliée au système de synchronisation.
	Le dernier élément est le générateur de
	gouttelettes. Il s'aμit d'un atomiseur DEHS/air
	permettant d'ensemencer l'écoulement porteur d'un
	brouillard de fines gouttelettes <	m de diamètre .
	L'aérosol λormé se comporte alors comme un traceur.
			Mesures Sur la configuration choisie, des mesures
			PIV 2C (Vélocimétrie par Images de Particules à deux
			composantes) ont été effectuées. Le post-traitement
			des mesures permet d'accéder aux champs et aux
			profils de vitesses (moyennes temporelles, moyennes
			de phase, partie fluctuante) de l'écoulement. Le
			signal du capteur de pression est utilisé pour le
			conditionnement des mesures.
			Description du Dispositif PIV 2C
			Le dispositif de mesures PIV 2C est
			constitué d'une source de lumière, d'un μénérateur
			de nappe, d'un système d'acquisition d'imaμes et de
			synchronisation de données LaVision, et d'un
			ensemenceur. La source de lumière provient de deux
			lasers de classe IV pulsés Nd :Yag (2,5 Hz, =
			nm). Ces lasers fournissent une énergie de 250 mJ par
			pulse pendant une durée d'environ ns. La période
			de temps entre deux impulsions laser t est réglée à
			s. Le système optique utilisé pour générer la
			tranche laser est constitué d'un jeu de lentilles
			cylindriques convergentes et de miroirs. La nappe
			laser résultante possède une largeur de 100 mm et


Francois Roure

Magdalena Scheck-Wenderoth

Liviu Matenco

Kristaq Muska

Shaqir Nazai

Dynamics and active processes : the albanian natural laboratory and analogues

come