Thomas Bonald
email: thomas.bonald@telecom-paristech.fr

Davide Cuda
email: davide.cuda@telecom-paristech.fr

Rate-Optimal Scheduling Schemes for Asynchronous Input-Queued Packet Switches

The performance of input-queued packet switches critically depends on the scheduling scheme that connects the input ports to the output ports. We show that, when packets are switched asynchronously, simple scheduling schemes where contention is solved locally at each input or output can achieve rate optimality, without any speed-up of the internal transmission rate.

INTRODUCTION

We consider an input-queued packet switch with virtual output queueing and no speed-up (i.e., the internal transmission rate of the switching fabric is equal to the line card rates). Current switches typically operate synchronously, the connections between input ports and output ports changing at the ticks of a common clock. Since packets have variable sizes, they must be segmented into fixed-size data units, which are switched independently and reassembled at the switch output.

The growth of the line card rates questions the viability of such an architecture [START_REF] Bianco | Asynchronous vs synchronous input-queued switches[END_REF]. In this paper, we consider an asynchronous switch, which does not require any packet segmentation. The connections between input ports and output ports are changed at arbitrary times. We describe some simple, distributed scheduling schemes that are provably rate-optimal, i.e., stabilize the queue lengths whenever possible. The contentions at each input or output are solved locally, based on the lengths of the corresponding virtual output queues (VOQ).

This stability result is remarkable in that synchronous switches are known to require complex scheduling schemes like maximum-weight matching to achieve rate optimality [START_REF] Tassiulas | Stability properties of constrained queueing systems and scheduling policies for maximum throughput in multihop radio networks[END_REF][START_REF] Mckeown | Achieving 100% throughput in an input-queued switch[END_REF][START_REF] Dai | The throughput of data switches with and without speedup[END_REF], except for the specific case of 2 × 2 switches [START_REF] Brzezinski | Greedy weighted matching for scheduling the input-queued switch[END_REF]. In fact, practical schedulers like PIM [START_REF] Anderson | High-speed switch scheduling for local-area networks[END_REF] or iSLIP [START_REF] Mckeown | The iSLIP scheduling algorithm for input-queued switches[END_REF] only guarantee maximal matchings (i.e., both the input and the output associated with a non-empty VOQ are matched) and, as such, are suboptimal [START_REF] Dimakis | Sufficient conditions for stability of longest-queue-first scheduling: Second-order properties using fluid limits[END_REF]. It turns out that, in the limit of large queue lengths, maximum matchings (i.e., the largest number of input-output connections associated with non-empty VOQs) are imposed by the switch asynchrony, which stabilizes the system.

Our result differs from those derived for wireless networks, that either rely on random access algorithms to approximate maximum-weight scheduling [START_REF] Jiang | Distributed random access algorithm: Scheduling and congestion control[END_REF] or apply simple, greedy algorithms under restrictive assumptions on the underlying interference graph [START_REF] Dimakis | Sufficient conditions for stability of longest-queue-first scheduling: Second-order properties using fluid limits[END_REF][START_REF] Feuillet | Random capture algorithms: Fluid limits and stability[END_REF][START_REF] Maguluri | The stability of longest-queue-first scheduling with variable packet sizes[END_REF]. Note in particular that the stability result of Maguluri, Hajek and Srikant [START_REF] Maguluri | The stability of longest-queue-first scheduling with variable packet sizes[END_REF] on Longest Queue First requires the local pooling condition of [START_REF] Dimakis | Sufficient conditions for stability of longest-queue-first scheduling: Second-order properties using fluid limits[END_REF], while that of Feuillet, Proutière and Robert [START_REF] Feuillet | Random capture algorithms: Fluid limits and stability[END_REF] on random capture algorithms assumes multipartite interference graphs. In both cases, the condition on the underlying conflict graph of an N × N switch imposes that N = 2, as in [START_REF] Brzezinski | Greedy weighted matching for scheduling the input-queued switch[END_REF]. We here use the bipartite structure of the graph connecting the inputs to the outputs to prove rate optimality for all N .

SCHEDULING SCHEME

Consider an N × N input-queued switch. Denote by I the set of inputs and by J the set of outputs. A feasible schedule is a matching on the bipartite graph G = (V, E) where V = I ∪ J and E = I × J.

Packets arriving at input i and destined to output j are directed to VOQ ij. If the VOQ is empty at the packet arrival and both the input and the output are unmatched, a new edge is added to the matching and the packet is directly transmitted to the output. Otherwise, the packet is queued. We denote by Qij the length of VOQ ij.

It remains to determine the change of matching at the end of a packet transmission, say from VOQ ij. Edge ij is then removed from the matching so that both input i and output j may select a new edge. By convention, we assume that the input takes the first decision. Specifically, for some α > 0:

• input i selects, if any, unmatched output j with probability proportional to Q α ij ;

• if there is no such output or if j = j, output j selects, if any, unmatched input i with probability proportional to Q α i j . The parameter α controls the sensitivity of the scheduling scheme to the queue lengths, from a purely random decision when α → 0 to a greedy algorithm based on Longest Queue First when α → ∞. Note that contentions are solved locally at input i then at output j. There are only two steps, unlike usual scheduling schemes in synchronous switches that are iterative and may involve, in theory, up to N steps [START_REF] Anderson | High-speed switch scheduling for local-area networks[END_REF][START_REF] Mckeown | The iSLIP scheduling algorithm for input-queued switches[END_REF].

TRAFFIC MODEL

We assume that packets arrive at VOQ ij according to a Poisson process with intensity λij and have exponential sizes with parameter µij (these assumptions are discussed in Section 6 below). We denote by ρij = λij/µij the corresponding load. Traffic is said to be admissible if:

∀i ∈ I, j∈J ρij < 1 and ∀j ∈ J, i∈I ρij < 1.
(

The vector Q(t) of VOQ lengths and the matching M (t) at time t define a Markov process. The scheduling scheme is said to be rate-optimal if this Markov process is ergodic for any admissible traffic.

EXAMPLE OF THE 3 × 3 SWITCH

To illustrate the stability result, we start with the simple example of a 3 × 3 switch with unit mean packet sizes and bidiagonal traffic pattern, i.e., λii = λ i|i+1|3 = λ for all i = 1, 2, 3. The load of each input or output is ρ = 2λ. Figure 1 shows the evolution of the length of an arbitrary VOQ at load ρ = 0.95, starting from the initial state Qij(0) = 100 for all VOQs ij receiving traffic and a matching M (0) of maximum size. Here α = 1 so that the input and output selection probabilities of the scheduling scheme are proportional to the queue lengths. Packet sizes are assumed to be constant and equal to the size of switched data units in the synchronous case, with a centralized scheduling scheme that iteratively selects active edges with probabilities proportional to the queue lengths. We observe that, while the VOQs eventually empty in the asynchronous case, they grow continuously in the synchronous case. The reason is that, in the latter case, the scheduling scheme regularly selects suboptimal matchings of size m = 2 (refer to Figure 2). Specifically, such matchings are selected with probability 1/3 when the lengths of the VOQs are equal and then limit the total service rate to 2/3 of the maximum service capacity. We obtain a maximum load of (1/3) × (2/3) + (2/3) × 1 = 8/9, as observed in [START_REF] Dimakis | Sufficient conditions for stability of longest-queue-first scheduling: Second-order properties using fluid limits[END_REF]. Since ρ > 8/9, the system is unstable.

In the asynchronous case, matchings of size m = 2 are unstable equilibriums: at the end of a packet transmission, the matching size increases to m = 3 whenever the input selects the other available output. In contrast, maximum matchings are stable equilibriums since packets are transmitted asynchronously: at the end of a packet transmission, the input cannot select another output, so that the same matching is used until one of the corresponding VOQs empties. In this sense, the scheduling scheme behaves as a random capture algorithm [START_REF] Feuillet | Random capture algorithms: Fluid limits and stability[END_REF].

RATE OPTIMALITY

We now state and prove the main result of the paper, namely the rate optimality of the scheduling scheme described in Section 2 for asynchronous packet switches. We start with two monotonicity properties of the matching size. We denote by E(t) = {(i, j) ∈ E, Qij(t) > 0} the set of active edges at time t.

Lemma 1. The matching size |M (t)| is a non-decreasing function of time t over any time interval where E(t) is constant.

Proof. Since E(t) is constant, ouput j can still be selected by input i at the end of the transmission of a packet from VOQ ij. 2

Lemma 2. If the matching M (t) is not maximum, the probability that |M (t+1)| > |M (t)| is lower bounded by some constant which only depends on the ratios Qij(t)/Q ik (t), for all (i, j), (i, k) ∈ E(t).

Proof. Let λ = (i,j)∈E λij, µ = (i,j)∈E µij and θ = min

(i,j)∈E(t) µij µ Qij(t) α k:(i,k)∈E(t) Q ik (t) α .
If the matching M (t) is not maximum, there exists some augmenting path P such that the matching

M (t)∆P ≡ (M (t) \ P ) ∪ (P \ M (t))
has size |M (t)| + 1, see, e.g., [START_REF] Lovász | Matching theory[END_REF]. The new matching is obtained by flipping the edge states (matched or unmatched) along path P . Since the length of any augmenting path is at most 2N , this can be done in one time unit with probability at least:

((1 -e -ν 2N )θ) 2N × e -(λ+µ)
with ν = min (i,j)∈E µij, where the first factor accounts for the flipping of the edge states along path P , each over a time interval of duration 1/(2N ), starting from the extreme edge of path P with matched input, and the second factor guarantees that no other event occurs. 2

Theorem 1. The Markov process (Q(t), M (t)) is ergodic under condition (1).

Proof. We prove the result by fluid limits, see, e.g. [START_REF] Dai | On positive Harris recurrence of multiclass queueing networks: a unified approach via fluid limit models[END_REF]. Consider a sequence of initial states

{(Q (n) (0), M (n) (0))} n≥1 such that: lim n→∞ Q (n) (0) n = Q(0), with | Q(0)| = 1.
Standard arguments based on the strong law of large numbers show the existence of an absolutely continuous, piecewise deterministic, stochastic process Q(t) on R N ×N + , referred to as the fluid limit, such that, for any such sequence of initial states, there is some subsequence indexed by n k such that, almost surely and uniformly on compact sets:

Q (n k ) (n k t) n k → Q(t) when k → ∞.
The corresponding fluid model is stable if there exists some finite time t0 such that Q(t) = 0 a.s. for all t ≥ t0. This is a sufficient condition for the ergodicity of the Markov process (Q(t), M (t)), see [START_REF] Dai | On positive Harris recurrence of multiclass queueing networks: a unified approach via fluid limit models[END_REF][START_REF] Robert | Stochastic networks and queues[END_REF]. Denote by Ē(t) = {(i, j) ∈ E, Qij(t) > 0} the set of active edges in the fluid model. Let [u, v) be some time interval such that Ē(t) is constant and not empty. We denote this set by Ē. By Lemmas 1 and 2, the fluid volumes are served by maximum matchings of the graph Ḡ = (V, Ē). Let m be the maximum matching size of Ḡ and define the workload of the fluid model at time t as:

W (t) = i∈I,j∈J 1 µij Qij(t).
Note that Qij(t) = 0 for all (i, j) ∈ Ē and t ∈ [u, v) so that:

∀t ∈ [u, v), dW (t) dt = (i,j)∈ Ē 1 µij d Qij(t) dt , =   (i,j)∈ Ē ρij   -m. (2)
Now by König's theorem [START_REF] Lovász | Matching theory[END_REF], the maximum matching size m is the minimum size of a vertex cover of Ḡ (set of vertices such that all edges of Ē are incident to at least one element of this set). Let C be a vertex cover of Ḡ of size m. We have:

(i,j)∈ Ē ρij ≤ i∈ C∩I j∈J ρij + j∈ C∩J i∈I ρij ≤ mδ, (3)
with

δ = max max i∈I j∈J ρij, max j∈J i∈I ρij .
The proof then follows from ( 2) and (3) on noting that δ < 1 and W (t) = 0 if and only if Q(t) = 0. 2

DISCUSSION

The assumption of Poisson packet arrivals is reasonable if traffic is generated by a large number of flows, each with a low peak rate relative to the line card rates. In contrast, flows with high peak rates may be limited by the line card rates and thus generate bursty packet arrivals with statedependent intensity due to the closed-loop control of TCP. The stability result easily extends to the case where each such flow maintains a constant number of packets in the queue, mimicking the behavior of TCP, until the completion of the data transfer. Such flows may, for instance, be assumed to arrive according to a Poisson process.

It proves more difficult to relax the assumption on the packet size distribution. In practice, this distribution is discrete with high weights on the minimum (40B) and maximum (around 1500B) sizes, which may synchronize packet departures. To keep rate optimality, it is then necessary to maintain asynchrony in the scheduling scheme by removing and adding edges successively, in some arbitrary order, as if packets would leave at different times. Note that the randomness in packet arrivals prevents the switch from the cyclic, pathological behavior described in [START_REF] Ganjali | Cell switching versus packet switching in input-queued switches[END_REF].

In future work, we plan to study the packet delay performance of the proposed scheduling schemes, under various traffic patterns.

Figure 1 :

 1 Figure 1: The stabilization effect of asynchrony.

(a) m = 2 (b) m = 3 Figure 2 :

 232 Figure 2: Maximal matchings of size m in the 3 × 3 switch with bidiagonal traffic pattern.


Sidi Ould Souvi

Kamal Sharkas

Julien Toulouse

Double-hybrid density-functional theory with meta-generalized-gradient approximations

We extend the previously proposed one-parameter double-hybrid density-functional theory [K.

I. INTRODUCTION

The double-hybrid (DH) approximations [1] have become ones of the most accurate approximations for electronic-structure calculations within densityfunctional theory (DFT). They consist in mixing Hartree-Fock (HF) exchange with a semilocal exchange density functional and second-order Møller-Plesset (MP2) correlation with a semilocal correlation density functional:

E DH xc = a x E HF x + (1 -a x )E x [n] +(1 -a c )E c [n] + a c E MP2 c , (1)
where the first three terms are calculated in a selfconsistent hybrid Kohn-Sham (KS) calculation, and the last MP2 term is usually evaluated with the previously obtained orbitals and added a posteriori (see, however, Ref. 2 for a double-hybrid scheme with orbitals optimized in the presence of the MP2 correlation term). The two empirical parameters a x and a c are usually determined by fitting to a thermochemistry database. A variety of such double hybrids have been constructed with different a x and a c parameters and various density functionals [1,3,4,[4][5][6][7][8]. Double-hybrid approximations with more parameters have also been proposed [9][10][11][12][13][14][15]. The so-called multicoefficient correlation methods combining HF, DFT and MP2 energies can also be considered to be a form of double-hybrid approximation [16][17][18][19].

Recently, Sharkas et al. [20] provided a rigorous theoretical justification for double hybrids based on the adiabatic connection formalism and which lead to a densityscaled one-parameter double-hybrid (DS1DH) approxi- * Present address: Institut de Radioprotection et Sûreté Nucléaire, PSN-RES/SAG/LETR, Cadarache, 13115 Saint-Paul-lès-Durance, France † Electronic address: sidi.souvi@irsn.fr ‡ Electronic address: julien.toulouse@upmc.fr mation

E DS1DH,λ xc = λE HF x + (1 -λ)E x [n] +E c [n] -λ 2 E c [n 1/λ ] + λ 2 E MP2 c , (2)
where E c [n 1/λ ] is the usual correlation energy functional evaluated at the scaled (squeezed) density n 1/λ (r) = (1/λ) 3 n(r/λ). In this class of double hybrids only one independent empirical parameter λ is needed instead of the two parameters a x and a c . The connection with the original double-hybrid approximations can be made by neglecting the density scaling

E c [n 1/λ ] ≈ E c [n], (3)
which leads to the one-parameter double-hybrid (1DH) approximation [20]

E 1DH,λ xc = λE HF x + (1 -λ)E x [n] +(1 -λ 2 )E c [n] + λ 2 E MP2 c . ( 4
)
corresponding to the standard double-hybrid approximation of Eq. ( 1) with parameters a x = λ and a c = λ 2 . DS1DH and 1DH approximations have been constructed using the Perdew-Burke-Ernzerhof (PBE) [21] and the Becke-Lee-Yang-Parr (BLYP) [22,23] exchangecorrelation density functionals, and it was found that when neglecting density scaling the accuracy on atomization energies largely deteriorates for PBE but in fact improves for BLYP (see, also, Ref. 24 for a 1DH approximation based on the modified Perdew-Wang [25] exchange functional and the Perdew-Wang-91 [START_REF] Perdew | Electronic Structure of Solids '91[END_REF] correlation functional). By appealing to the high-density limit of the correlation functional, Toulouse et al. [START_REF] Toulouse | [END_REF] argued that a more sensible approximation to the density-scaled correlation functional is

E c [n 1/λ ] ≈ (1 -λ)E MP2 c + λE c [n], (5)
leading to the linearly scaled one-parameter doublehybrid (LS1DH) approximation

E LS1DH,λ xc = λE HF x + (1 -λ)E x [n] +(1 -λ 3 )E c [n] + λ 3 E MP2 c , (6)
which works reasonably well for PBE and corresponds to the form used for the PBE0-DH [28] and PBE0-2 [29] double hybrids. Finally, note that Fromager and coworkers [30,31] explored the theoretical basis of the twoparameter double-hybrid approximations and proposed new double-hybrid schemes.

The vast majority of all these double hybrids have been applied with generalized-gradient approximation (GGA) density functionals. Only recently, Goerigk and Grimme [7] proposed two-parameter double-hybrid schemes based on meta-GGA density functionals. In these schemes, the opposite-spin MP2 correlation energy is combined with a refitted Tao-Perdew-Staroverov-Scuseria (TPSS) [32] meta-GGA exchange-correlation functional, or with a refitted Perdew-Wang [START_REF] Perdew | Electronic Structure of Solids '91[END_REF] GGA exchange functional and a refitted Becke-95 (B95) [33] meta-GGA correlation functional. Kozuch and Martin [13,15] also tested numerous spin-component-scaled double hybrids, with four or more optimized parameters, including the B95 meta-GGA correlation functional and the TPSS, B98 [34], BMK [35], τ HCTH [36] meta-GGA exchange-correlation functionals.

In this work, we reexamine the theoretical basis of double-hybrid approximations using meta-GGA functionals, and we construct one-parameter double hybrids using the TPSS meta-GGA functional. While the 1DH and LS1DH schemes can be readily applied with a meta-GGA functional, the DS1DH scheme requires an extension of the density scaling to the non-interacting kinetic energy density τ (r) that is used in the TPSS functional. We then assess the accuracy of these one-parameter meta-GGA double hybrids on test sets of atomization energies and reaction barrier heights.

II. DOUBLE-HYBRID META-GGA APPROXIMATIONS

In addition to the explicit dependence on the density n(r), its gradient ∇n(r) and possibly its Laplacian ∇ 2 n(r), a meta-GGA density functional also generally depends on the non-interacting positive kinetic energy density τ (r) which can be seen as an implicit functional of the density (for simplicity, we only write the equations for the spin-unpolarized case; the extension to the general spin-polarized case is straightforward):

τ [n](r) = Φ[n]|τ (r)|Φ[n] = 1 2 occ i |∇ r φ i [n](r)| 2 , ( 7
)
where the positive kinetic energy density operator τ (r) = (1/2) σ=↑,↓ ∇ r ψ † σ (r) • ∇ r ψσ (r) has been introduced (with the creation and annihilation field operators ψ † σ (r) and ψσ (r)), Φ[n] is the (KS) single-determinant wave function minimizing the kinetic energy Φ| T |Φ and giving the density n, and φ i [n](r) are the associated (KS) spin-orbitals.

Defining the scaled non-interacting kinetic energy density τ γ as the non-interacting kinetic energy density corresponding the scaled density n γ (r) = γ 3 n(γr) (where γ is an arbitrary positive scaling factor), one can show

τ γ [n](r) ≡ τ [n γ ](r) = 1 2 occ i |∇ r φ i [n γ ](r)| 2 = γ 3 2 occ i |∇ r φ i [n](γr)| 2 = γ 5 2 occ i |∇ γr φ i [n](γr)| 2 = γ 5 τ [n](γr), (8)
where the scaling relation

φ i [n γ ](r) = γ 3/2 φ i [n](γr)
has been used. The scaling relation on τ [n](r) is consistent with the well-known quadratic scaling of the non-interacting kinetic energy [37]. Furthermore, this scaling relation can easily be verified for one-electron systems with the von Weizsäcker kinetic energy density, τ W (r) = (1/8) |∇ r n(r)| 2 /n(r), and for the uniform electron gas,

T s [n γ ] = γ 2 T s [n] where T s [n] = τ [n](r)dr
τ unif = (3/10)(3π 2 ) 2/3 n 5/3 .
Working with τ [n](r) as an implicit functional of the density in a self-consistent KS calculation requires to use the optimized effective potential approach to calculate the functional derivative of the exchange-correlation energy E xc [n] with respect to the density n [38], which is computationally impractical. The standard practice [39][40][41][42][43][44][45] is to consider the exchange-correlation energy as an explicit functional of both n and τ , E xc [n, τ ]. Following this practice, we define a density-scaled one-parameter hybrid (DS1H) approximation with a τ -dependent functional as

E DS1H,λ = min Φ Φ| T + Vext + λ Ŵee |Φ + Ēλ H [n Φ ] + Ēλ xc [n Φ , τ Φ ] , (9)
where Φ is a single-determinant wave function, T is the kinetic energy operator, Vext is the external (e.g., electron-nucleus) potential operator, Ŵee is the Coulomb electron-electron interaction operator, Ēλ H [n Φ ] and Ēλ xc [n Φ , τ Φ ] are the complement Hartree and exchange-correlation functionals evaluated at the density and kinetic energy density of Φ, n Φ (r) = Φ|n(r)|Φ and τ Φ (r) = Φ|τ (r)|Φ . The complement Hartree and exchange functionals are linear with respect to λ:

Ēλ H [n] = (1 -λ)E H [n], (10)
Ēλ x [n, τ ] = (1 -λ)E x [n, τ ], (11)
where E H [n] and E x [n, τ ] are the usual KS Hartree and exchange functionals. The complement correlation functional is obtained via the extension of uniform coordinate scaling of the density [46][47][48][49] to the kinetic energy density

Ēλ c [n, τ ] = E c [n, τ ] -E λ c [n, τ ] = E c [n, τ ] -λ 2 E c [n 1/λ , τ 1/λ ], (12)
where E c [n, τ ] is the usual KS correlation functional, E λ c [n, τ ] is the correlation functional corresponding to the interaction λ Ŵee , n 1/λ (r) = (1/λ) 3 n(r/λ) is the scaled density and τ 1/λ (r) = (1/λ) 5 τ (r/λ) is the scaled kinetic energy density.

The minimizing single-determinant wave function Φ λ in Eq. ( 9) is calculated by the self-consistent eigenvalue equation:

T + Vext + λ V HF Hx [Φ λ ] + V λ H [n Φ λ ] + V λ xc [n Φ λ , τ Φ λ ] |Φ λ = E λ 0 |Φ λ , (13)
where V HF Hx is the nonlocal HF potential operator, V λ H is the complement local Hartree potential operator, and V λ xc is the complement exchange-correlation potential operator

V λ xc [n, τ ] = δ Ēλ xc [n, τ ] δn(r) n(r)dr + δ Ēλ xc [n, τ ] δτ (r) τ (r)dr, (14)
where n(r) = σ=↑,↓ ψ † σ (r) ψσ (r) is the density operator and the second term in Eq. ( 14) corresponds to a nonmultiplicative "potential" operator.

The nonlinear Rayleigh-Schrödinger perturbation theory of Refs. 20, 50-52 can readily be extended to start with the DS1H reference of Eq. ( 9) with a τ -dependent functional (details are given in the supplementary material [53]). For the second-order energy correction, due to Brillouin's theorem, only double excitations contribute and consequently the nonlinear terms of the perturbation theory vanish since they involve expectation values of the one-electron operators n(r) and τ (r) between determinants differing by two spin orbitals. The secondorder energy correction to be added to the DS1H energy has thus a standard MP2 form

E λ,(2) = λ 2 i<j a<b | ij||ab | 2 ε i + ε j -ε a -ε b = λ 2 E MP2 c , (15)
where i, j and a, b refer to occupied and virtual DS1H spin-orbitals, respectively, with associated orbital eigenvalues ε k , and ij||ab are the antisymmetrized twoelectron integrals. The DS1DH exchange-correlation energy for a τ -dependent functional is thus (dropping from now on the explicit dependence on Φ λ )

E DS1DH,λ xc = λE HF x + (1 -λ)E x [n, τ ] +E c [n, τ ] -λ 2 E c [n 1/λ , τ 1/λ ] + λ 2 E MP2 c . ( 16
)
Neglecting the scaling in the correlation functional,

E c [n 1/λ , τ 1/λ ] ≈ E c [n, τ ], gives the 1DH exchange- correlation energy E 1DH,λ xc = λE HF x + (1 -λ)E x [n, τ ] +(1 -λ 2 )E c [n, τ ] + λ 2 E MP2 c , (17)
and using the approximate scaling

E c [n 1/λ , τ 1/λ ] ≈ (1 -λ)E MP2 c + λE c [n, τ ] gives the LS1DH exchange- correlation energy E LS1DH,λ xc = λE HF x + (1 -λ)E x [n, τ ] +(1 -λ 3 )E c [n, τ ] + λ 3 E MP2 c . (18)
We apply these double-hybrid schemes with the TPSS exchange-correlation functional and refer to them as DS1DH-TPSS, 1DH-TPSS and LS1DH-TPSS.

III. COMPUTATIONAL DETAILS

Calculations have been performed with a development version of the MOLPRO program [START_REF] Werner | Molpro, version 2012.1, a package of ab initio programs[END_REF], in which the DS1DH-TPSS, 1DH-TPSS and LS1DH-TPSS approximations have been implemented. The scaling relations for the scaled exchange-correlation energy and its derivatives for a general meta-GGA functional are given in the appendix. The empirical parameter λ is optimized on the AE6 and BH6 test sets [START_REF] Lynch | [END_REF]. The AE6 set is a small representative benchmark set of six atomization energies consisting of SiH 4 , S 2 , SiO, C 3 H 4 (propyne), C 2 H 2 O 2 (glyoxal), and C 4 H 8 (cyclobutane). The BH6 set is a small representative benchmark set of forward and reverse hydrogen barrier heights of three reactions, OH + CH 4 → CH 3 + H 2 O, H + OH → O + H 2 , and H + H 2 S → HS + H 2 . All the calculations for the AE6 and BH6 sets were performed at the optimized QCISD/MG3 geometries [START_REF]The geometries are available in the Minnesota Databases for Chemistry and Solid[END_REF] using the Dunning cc-pVQZ basis set [START_REF] Dunning | [END_REF]59]. The performance of the best double hybrid is then checked on the larger benchmark set of 49 atomization energies of Ref. 54 (G2-1 test set [60,61] except for the six molecules containing Li, Be, and Na) at MP2(full)/6-31G* geometries using the Dunning cc-pVQZ basis set. Core electrons are kept frozen in all our MP2 calculations. Spinrestricted calculations are performed for all the closedshell systems, and spin-unrestricted calculations for all the open-shell systems.

IV. RESULTS AND DISCUSSION

In Figure 1, we plot the mean absolute errors (MAEs) for the AE6 and BH6 test sets as functions of the parameter λ for the DS1DH-TPSS, 1DH-TPSS and LS1DH-TPSS approximations. The MAEs and mean errors (MEs) of the double hybrids based on the TPSS functional at the optimal values of λ minimize the MAEs on the AE6 and BH6 sets are also reported in Table I, and compared to those obtained with standard BLYP, PBE, TPSS and MP2, as well as with other double-hybrid approximations based on the BLYP and PBE functionals.

For λ = 0, all these double hybrids reduce to a standard KS calculation with the TPSS functional, while for λ = 1 they all reduce to a standard MP2 calculation. For the AE6 set, DS1DH-TPSS gives by far the smallest MAE with 0.91 kcal/mol at the optimal value of λ = 0.70. LS1DH-TPSS gives a larger MAE of 3.28 kcal/mol for an optimal value of λ = 0.85, and 1DH-TPSS gives a yet larger MAE of 5.74 kcal/mol for an optimal value of λ = 0.05, providing virtually no improvement over the TPSS KS calculation at λ = 0. For the BH6 set, the three double-hybrid approximations are very similar for the entire range of λ, indicating that the scaling in the correlation functional contribution is not as crucial for barrier heights as for atomization energies. The MAE minima are 0.59 kcal/mol at λ = 0.75, 0.74 kcal/mol at λ = 0.78 and 0.96 kcal/mol at λ = 0.78 for DS1DH-TPSS, 1DH-TPSS, and LS1DH-TPSS, respectively.

Neglecting the scaling of the density and of the kinetic energy density in the TPSS correlation functional,

E c [n 1/λ , τ 1/λ ] ≈ E c [n, τ ],
i.e. going from DS1DH-TPSS to 1DH-TPSS, largely deteriorates the accuracy of atomization energies. A similar deterioration is obtained when neglecting the scaling of the density in the PBE correlation functional, whereas a large improvement of atomization energies is observed when neglecting the scaling of the density in the LYP correlation functional (see Table I and Ref. 20). These different behaviors between the PBE and TPSS correlation functionals on a one hand and the LYP correlation functional of the second hand may be related to the observation on atoms with nondegenerate KS systems that PBE and TPSS are more accurate than LYP for the high-density limit (or weakinteraction limit), E

c = lim λ→0 E c [n 1/λ ] [62-65]. (2)
Contrary to most other double hybrids, the DS1DH-TPSS double-hybrid approximation gives very close optimal values of λ on the AE6 and BH6 sets, i.e. λ = 0.70 TABLE II: Atomization energies (in kcal/mol) of the 49 molecules of the set of Ref. 54 (G2-1 test set except for the six molecules containing Li, Be, and Na). The calculated values were obtained using the functional TPSS and the double hybrids DS1DH-TPSS (with λ = 0.725) and B2-PLYP with the cc-pVQZ basis set and MP2(full)/6-31G* geometries. The zero-point energies are removed in the reference values. For each method, the value with the largest error is indicated in boldface. and λ = 0.75. For general applications, we propose to use the value of λ = 0.725 for this double hybrid when using the cc-pVQZ basis set.

Finally, in Table II, we compare DS1DH-TPSS (with the optimal parameter λ = 0.725) with TPSS and the standard double hybrid B2-PLYP [1] on the larger set of 49 atomization energies of Ref. 54. With MAEs of 3.9 and 2.3 kcal/mol for TPSS and DS1DH-TPSS, respectively, it is clear that the improvement in accuracy brought by DS1DH-TPSS over TPSS observed on the small AE6 set remains (although smaller) for this larger set. DS1DH-TPSS is however slightly less accurate on average on this set than B2-PLYP (MAE of 1.6 kcal/mol).

V. CONCLUSIONS

We have constructed one-parameter double-hybrid approximations using the TPSS meta-GGA exchangecorrelation functional and tested them on test sets of atomization energies and reaction barrier heights. We have shown that neglecting the scaling of the density and of the kinetic energy density in the correlation functional largely deteriorates the accuracy on atomization energies, in contrast to what was previously found for double hybrids based on the BLYP functional. We thus propose the density-scaled double-hybrid DS1DH-TPSS approximation with a fraction of HF exchange of λ = 0.725 as a viable meta-GGA double hybrid for thermochemistry calculations, improving over both standard KS TPSS and MP2 calculations. We hope that this work will lead to more investigations of meta-GGA double hybrids with minimal empiricism. Possible extensions of this work include introducing meta-GGA functionals in multiconfigurational hybrids [66] or in Coulomb-attenuated double hybrids [67].

We give the expressions for the scaled meta-GGA correlation functional E

FIG. 1 :

 1 FIG.1: MAEs for the AE6 (left) and BH6 (right) test sets as functions of the parameter λ for the DS1DH, 1DH and LS1DH approximations with the TPSS exchange-correlation density functional. All calculations were carried out with the cc-pVQZ basis set.

a

  Data from Ref. 20. b This is a local minimum. The global minimum is for λ = 1.0, i.e. MP2. c Data from Ref. 27.

  λ c [n, τ ] = λ 2 E c [n 1/λ , τ 1/λ ]and its derivatives (see, also, Refs.20,66). Starting from a standard meta-GGA density functional, depending on the density n(r), the square of the density gradient |∇ r n(r)|2 , the Laplacian of the density ∇ 2 r n(r) and/or the non-interacting kinetic energy density τ (r)E c [n, τ ] = e c n(r), |∇ r n(r)| 2 , ∇ 2 r n(r), τ (r) dr,(A1)the corresponding scaled functional is written asE λ c [n, τ ] = e λ c n(r), |∇ r n(r)| 2 , ∇ 2 r n(r), τ (r) dr,(A2)where the energy density is obtained by the scaling relation e λ c = λ 5 e c relations apply for spin-dependent functionals E c [n ↑ , n ↓ , τ ↑ , τ ↓ ].

TABLE I :

 I MAEs and MEs (in kcal/mol) on the AE6 and BH6 test sets for several methods. For the double-hybrid DS1DH, 1DH, and LS1DH approximations, the results are for the optimal values of λ which minimize the MAEs of the AE6 and BH6 sets, separately. All calculations were carried out with the cc-pVQZ basis set.

			AE6		BH6
	Method	λ	MAE ME	λ	MAE ME
	BLYP a PBE a		6.52 -1.18 15.5 12.4		8.10 -8.10 9.61 -9.61
	TPSS MP2 a		5.79 3.78 6.86 4.17		8.41 -8.41 3.32 3.11
	DS1DH-BLYP a λ = 0.80 1DH-BLYP a λ = 0.55 B2-PLYP a	4.73 -2.52 1.46 0.07 1.39 -1.09	λ = 0.65 0.60 λ = 0.75 0.80 -0.18 0.24 2.21 -2.21
	DS1DH-PBE a λ = 0.65	3.78 1.30	λ = 0.80 1.32	0.48
	1DH-PBE a

Acknowledgments

We thank Andreas Savin (UPMC/CNRS, Paris) for stimulating discussions.


Ahlem Khlass
email: ahlem.khlass@orange.com

Salah Eddine Elayoubi
email: salaheddine.elayoubi@orange.com

Thomas Bonald
email: thomas.bonald@telecom-paristech.fr

Capacity Gains from Multipoint Single Frequency Transmission in HSPA+

Keywords: HSPDA, Multipoint Transmission, HS-SFN, Sector Coordination, Cellular Networks, Flow Level Modeling, Queuing Theory

High Speed-Single Frequency Network (HS-SFN) is one of the possible multipoint transmission techniques proposed in the 3GPP standard for High Speed Downlink Packet Access (HSDPA) in order to improve network performance, especially at the cell edge. It allows neighboring cells to transmit simultaneously the same data stream to a User Equipment (UE) in the Handover region (HO). In order to evaluate the user-level performance of this technique, we develop a method based on network simulation coupled with Markov chain analysis. It shows that when the HS-SFN technique is performed in the HO region between adjacent cells, the user data rates increase significantly. However, this is true only when involved cells are partially loaded, which is not always the case. We thus propose an optimized approach that adapts the coordination area based on the average offered traffic observed in the network. Network performance is then improved at any load.

I. INTRODUCTION

Since its initial release in 1999, the WCDMA system [START_REF] Holma | WCDMA for UMTS: Radio Access for Third Generation Mobile Communications[END_REF] has been upgraded to offer substantial improvements in data rates and capacity. The first step was to improve the downlink by introducing the High-Speed Downlink Packet Access (HSDPA) system in release 5. Subsequently, the uplink has been upgraded by defining the High-Speed Uplink Packet Access (HSUPA) system in release 6 [START_REF] Dahlman | 3G Evolution[END_REF]. Introduced in release 7, the Enhanced HSPA system (also known as HSPA+) provides advanced features such as multiple input-multiple output (2x2 MIMO) technologies and higher order modulation (64 QAM). Recently, new features have been added such as the Multi-carrier HSDPA (MC-HSDPA) where two or more carriers are involved in serving a user equipment (UE). More recently, within the framework of Release 11, the concept of HSDPA Multipoint Transmission has been introduced [START_REF][END_REF]. The difference between multi-carrier and multipoint HSDPA is that the former allows multiple carriers, belonging to the same sector, to serve simultaneously the UE, while the latter coordinates different sectors to serve cell edge UEs.

There are three types of multipoint transmission schemes:

• In the Single Point data transmission scheme, the UE is allowed to be attached to two sectors but can be served only by one of them at each time slot, depending on its instantaneous channel.

• In the Multiflow data transmission, all sectors can simultaneously schedule independent transport blocks to the UE. In this scheme, coordinated sectors may or may not belong to the same NodeB but, in the latter case, the two NodeBs must exchange some signaling information.

In [START_REF] Zhang | System-level analysis and evaluation of sf-dc transmit mode in hspa system[END_REF], the authors propose a system-level evaluation of this scheme where two cells schedule independent data streams to the same UE in the same frequency. Results show that it performs better than the classical HSDPA transmission in terms of throughput and has the ability to balance load in asymmetric loading environment. Furthermore, the effect of coordination on UEs which are not concerned by the joint scheduling is shown to be negligible. • In the Single Frequency Network (SFN) data transmission, the same transport block is scheduled simultaneously by a pair of sectors belonging to the same NodeB.

As its name suggests, transmissions are performed synchronously so that the signals are combined and seen as a single multipath signal at the UE receiver [START_REF] Puchko | Phase adjustment in hs-sfn for hsdpa[END_REF]. We focus in this paper on this coordination scheme as its performance evaluation has been rarely addressed in the literature.

These HSDPA Multipoint techniques enable to enhance network performance by improving user experience, especially at cell edge, and efficiently balancing load [START_REF]Hspa+ advanced smart networks: Multipoint transmission[END_REF]. It turns out that the user dynamics are rarely taken into consideration when evaluating the performance gains of such techniques. The aim of this work is to evaluate the performance gains of HS-SFN in HSDPA at flow level, accounting for flow arrivals and departures. To do that, we develop a method based on network simulation coupled with Markov chain analysis. We consider that within a site, each pair of sectors serves users in the Handover (HO) region by simultaneous transmissions of the same data stream in the same carrier frequency. Performance is evaluated in terms of throughput in average and at cell edge.

The main contributions of this work are the following:

1) We propose an original hybrid performance evaluation approach based on a realistic radio simulator coupled with a Markovian analysis, where the former assesses the impact of coordination on the peak throughputs while the latter catches the behavior of the resource allocation mechanism.

2) We take into account the flow level dynamics in the evaluation and assess the impact of traffic on the coordination gain. 3) We propose an optimization of the coordination scheme that chooses the appropriate coordination zone. This paper is organized as follows. In section II, the concept of sector coordination using HS-SFN for HSDPA and the resource allocation algorithm are introduced. The performance analysis is presented in section III. Simulation results are provided in section IV. In section V, we propose an optimized coordination. Finally, conclusions and further work are drawn in section VI.

II. HS-SFN FOR HSDPA NETWORK A. Concept of sector coordination using HS-SFN

In classical HSDPA, each UE is associated to the socalled serving base station (BS) that offers the best signal. When signals from non-serving BSs exceed the signal from the serving one by a certain threshold, a handover (HO) mechanism is triggered in order to select the best BS to switch to. This threshold is denoted by δ handover and a user in this region verifies:

Received signal from the serving BS Received signal from the non-serving BS < δ handover

This region can be represented also by the fraction α of the total surface covered by the BS.

As proposed by the 3GPP standard [START_REF][END_REF], once a UE falls in this region, the non-serving BS with the best signal coordinates with the serving one in order to schedule the same transport block to this UE in the same frequency and using the same scrambling code, hence signals are combined and the receiver sees just one multipath signal. Thus, the Signalto-Noise-plus-Interference-Ratio (SINR) is improved leading to a better user experience at the cell edge. In fact, focusing on a tri-sectored site as shown in figure 1 where sectors are loaded symmetrically, the SINR measured at location u using the single point transmission technique (without coordination) is as follows:

SIN R u = P 0,s (u)

3 i=1 i =s P 0,i (u) + l 3 i=1 P l,i (u) + n ,
where P l,i (u) denotes the received power from sector S i for i ∈ N = {1, 2, 3} in site l (sector s of site 0 is supposed to be the serving sector of user u) and n is the thermal noise. However, when the HS-SFN technique is performed at the same location (in the HO region between adjacent sectors), the SINR becomes:

SIN R u = P 0,s (u) + P 0,c (u)

3 i=1 i =s,c P 0,i (u) + l 3 i=1 P l,i (u) + n .
We notice that the signal received from sector c is converted from interfering signal to useful signal. The SINR is thus largely improved as the received signal is getting stronger.

B. Resource Allocation Algorithm

As it is shown in figure 1, we can divide the tri-sectored site into six zones: We consider a centralized scheduler per site. Resources are divided according to a Proportional Fair resource allocation [START_REF] Viswanath | Opportunistic beamforming using dumb antennas[END_REF], [START_REF] Khlass | Resource allocation in coordinated cellular networks[END_REF] so that BS 1 , BS 2 and BS 3 can serve separately and independently users in zone Z 1 , Z 2 and Z 3 , respectively. But, when serving users in Z 12 for instance, BS 1 and BS 2 join their resources and block users from Z 1 and Z 2 , similarly for users in zones Z 23 and Z 31 . Since we neglect the fast-fading variations of the channel, the resource allocation algorithm can be described as follows:

At each timeslot,

• Select a user at random,

• Remove all other users sharing the same resources with it from the list of eligible users, • Iterate until all sectors are occupied (i.e. all resources are used) or the list of eligible mobiles is empty; • Serve the selected mobiles. We can thus express the fraction of time allocated to each type of users at a given time. Let K = {1, 12, 2, 23, 3, 31} be the set of zones in the site. We denote by y k the number of users in each zone Z k and y their corresponding vector. The Total number of users in the site is given by ȳ. Then, the fraction of time allocated for instance to users in zone Z 12 for instance is given by

φ 12 (y) = y 12 ȳ + y 3 ȳ y 12 y 1 + y 2 + y 12 .
In fact, a user in zone Z 12 is served if it is selected at first with a fraction of y 12 /ȳ or after selecting a user in zone Z 3 with a fraction of y 3 /ȳ (blocking by this way users in Z 23 and Z 13 ). Similarly, we get the fraction of time allocated to users in Z 1 :

φ 1 (y) = y 1 + y 23 ȳ + y 2 ȳ y 1 + y 3 y 1 + y 3 + y 31 + y 3 ȳ y 1 + y 2 y 1 + y 2 + y 12 ,
and so on for the other zones.

III. PERFORMANCE ANALYSIS

In this section, we first present the achievable rates obtained by the technique of HS-SFN in a HSDPA network. Then, these rates are used in the Markov process that describes the system dynamics. Finally, performance is evaluated in terms of average throughput and cell edge throughput for different load situations.

A. Achievable Rates

We consider a homogeneous, trisectored, HSDPA network in a hexagonal grid topology. Two interfering rings (i.e., a total of n = 19 cells) are considered when calculating intercell interference. We divide the cell depicted in figure 1 into elementary surfaces defined by an angle a relative to the main direction of the antenna (i.e. a is an integer ∈ [0, 120])) and a distance from the site center r ∈ [50..350].

Based on the analysis of section II.A, we measure the SINR of a single user at different positions (a, r) in two cases: i) Classical transmission (i.e. without sector coordination) and ii) Multipoint transmission using the HS-SFN technique (i.e. with sector coordination). Once this SINR is obtained, link level curves are used that associate to each SINR the corresponding throughput by taking into account the fast fading as described in [START_REF] Saadani | Realistic performance of hsdpa evolution 64-qam in macro-cell environment[END_REF]. We thus obtain two matrices of peak full buffer throughputs, R(a, r) and R (a, r), without and with coordination, respectively. Figure 2 illustrates the classical HSDPA throughput over a half-sector (PA3 channel, dense urban area with an inter-site distance of 500 meters), while figure 3 illustrates the case where α = 30% of the cell surface is subject of coordination. It can be observed that the peak throughput is largely improved when coordination is applied.

B. Markov Analysis

In the previous section, we showed the impact of coordination on the throughput when the user is alone in the cell. In this section, we perform a Markovian analysis that takes into account the impact of the centralized scheduling and the flow level dynamics on the performance.

Let C = {1, 2, ..., M } be the set of elementary surfaces such as each element i ∈ C corresponds to a combination of (a, r) and belongs to one of the six zones defined previously. We assume the traffic demand is uniformly distributed in the site. Users arrive to the site at random following a Poisson process of intensity λ and leave the network once they finish their data transfer. Precisely, users in position j arrive following a 3D representation of the peak throughput using single point transmission (hexagonal, tri-sectorized network). Poisson process of intensity λ i = p i λ. Flows are characterized by their size which is independent and identically distributed with mean σ. The mean service rates without and with coordination are given by µ i and µ i respectively such as

µ i = R i /σ; µ i = R i /σ.
The system can be described as a Markov process. Each state is denoted by x = {x i } i∈C and the transition rates between two states x and x are given as follows

q(x, x ) =          λ i for x = x + e i xi y k φ k (y)µ i for x = x -e i with i ⊂ Z k , k = {1, 2, 3} xi y k φ k (y)µ i for x = x -e i with i ⊂ Z k , k = {12, 23, 31} 0 else
where e i denotes an M -dimensional vector where the i'th component is equal to one and all the others are equal to zero.

We can obtain, by simulation, the steady-state probabilities π(x), leading to the average numbers of users in class i [START_REF] Bonald | Network Performance Analysis[END_REF] [11] such as Using Little's law, we obtain that the mean delay users in i as

E(x i ) = x x i π(x)
E(T i ) = E(x i ) λ i
The throughput, defined as the ratio of the mean flow size to the mean flow delay, is given by

γ i = σ E(x i )

C. Simulation Results

In this part, we evaluate the HS-SFN technique by injecting the full buffer throughputs of Figures 2 and3 in the Markovian flow level analysis. We define two performance measures: The average flow throughput over the cell and the cell edge throughput as the throughput achieved by the worst 5% of users. Figures 4 and5 represent the average and the cell edge throughputs obtained for different traffic values λ. We notice an important increase in both average and cell edge throughputs in low load situations. However, for high load situations, we see a deterioration of the throughput, followed by a decrease of the stability region when using the HS-SFN coordination technique. This is due to the fact that, even if sectors are loaded, they have to schedule users from other sectors leading to less scheduling opportunities for their own users. This motivates us to look for an optimized coordination area where to perform the HS-SFN for HSDPA without sacrificing the system performance.

IV. OPTIMIZED COORDINATION AREA

As seen so far, the performances of HS-SFN are mainly related to the average traffic observed in the cell. Allowing an α = 30% of the total area of a sector to enable Multipoint transmission is beneficial until some point where the single point transmission (without cell coordination) performs better. We thus study the impact of coordination area on the system performance. We start by illustrating, in figure 6, the average throughput γ(α, λ) for all possible values of α ranging from 0 % to 100 % and for different traffic values λ. We can observe that the extent of the coordination area has a large impact on performance and that this impact differs depending on the offered traffic. Figure 7 further illustrates this impact for some values of α ∈ {0, 15, 30, 50, 75, 100}. It can be observed that a value of α can be better than another in some load situations, and worse in others. More precisely, for low load situation, the best performances are obtained with large coordination areas while lower coordination is preferred for high loads.

We thus propose an optimized coordination scheme where the coordination zone is broadened or narrowed depending on the average traffic. For that, an information about the handled traffic need to constantly available in each sector. Hence, the optimal coordination value, denoted by α * , is the solution of the optimization problem described as

α * = arg max α γ(α, λ)
This optimal value is illustrated in figure 8 where α * is plotted for different traffic situations. We notice that the value of α * decreases with the load, until reaching the value of 0 near the network stability condition where the classical HSDPA operation is the preferred one. Using these results, each sector has to configure its coordination area in respect to the average traffic load. Figures 9 and 10 show the gain achieved when using the optimized coordination scheme compared to the common coordination scheme (with α = 30%)in average and at cell edge, respectively. We see that with a constant coordination threshold, we may have negative impact on network performances, while the optimized scheme brings almost always a positive gain.

V. CONCLUSION

In this paper, we analyze the flow-level performances of HS-SFN for HSDPA as proposed by the 3GPP standard. We develop a method based on network simulation coupled with Markov chain analysis. Results show an important increase in the average throughput at a price of a significant decrease in the stability region of the system. We hence propose an optimized coordination scheme that adapts the coordination area with respect to the load observed in the network and show that it achieves good performances in all traffic situations.

As future work, we aim at proposing self optimizing algorithms that allow adapting the coordination threshold based on radio conditions and traffic measurements. We aim also at studying other cell coordination algorithms (Single/Dual frequency Multiflow data transmission) in HSPA+ and also Coordinated Multipoint schemes in LTE-Advanced networks.

Fig. 1 .

 1 Fig. 1. A tri-sectored site.

•

  Non-coordination zones where users are served exclusively by a unique BS noted Z 1 , Z 2 and Z 3 . • Coordination zones where users are served simultaneously by two BSs noted Z 12 , Z 23 and Z 31 . For instance, zone Z 1 ⊂ S 1 includes users served exclusively by BS 1 and zone Z 12 ⊂ S 1 ∩ S 2 includes users served simultaneously by BS 1 and BS 2 .

  Fig. 2.3D representation of the peak throughput using single point transmission (hexagonal, tri-sectorized network).

Fig. 3 .

 3 Fig. 3. 3D representation of the peak throughput using Multipoint transmission (α = 30%)

Fig. 5 .

 5 Fig. 4. Average throughput VS traffic

Fig. 6 .Fig. 7 .

 67 Fig. 6. Average throughput for all possible values of α

Fig. 8 .Fig. 9 .

 89 Fig. 8. Optimal coordination area.

Fig. 10 .

 10 Fig.10. Coordination gain on cell edge throughput for the optimized case and the case where α = 30%.


Amélie Danlos
email: amelie.danlos@cnam.fr

Florent Ravelet

Olivier Coutier-Delgosha

Farid Bakir

Cavitation regime detection through Proper Orthogonal Decomposition: dynamics analysis of the sheet cavity on a grooved convergent-divergent nozzle

Keywords: Proper Orthogonal Decomposition, Sheet cavity, Cloud cavitation, Passive control, Distributed organized roughness, Venturi-type section

come

Introduction

Many methods are used to delay the cavitation inception, responsible for many issues like erosion, noise and vibrations. When pressure becomes inferior to the vapor pressure of the fluid, cavitation appears in the case of a flow over an obstacle, like hydrofoils, convergent-divergent nozzles or on blades in turbines and propellers [START_REF] Mejri | Comparison of computational results obtained from a homogeneous cavitation model with experimental investigations of three inducers[END_REF]). In some industrial configurations, it is difficult to avoid the cavitation inception. However, in some devices, cavitation can be an advantage for drag reduction or speed augmentation, especially in marine vehicles, by means of supercavitation [START_REF] Ceccio | Friction drag reduction of external flows with bubble and gas injection[END_REF]). Moreover it is important to have a better knowledge of the two-phase flow dynamics once cavitation has appeared in order to reduce its negative effects and/or to increase its positive influence. A convergent-divergent step like a Venturi type section is a simple geometry that leads to a rich sheet cavitation dynamics [START_REF] Barre | Experiments and modeling of cavitating flows in venturi: attached sheet cavitation[END_REF]; [START_REF] Campos-Amezcua | Numerical analysis of unsteady cavitating flow in an axial inducer[END_REF]).

The development of the instability of a sheet cavity on a Venturi-type section has been widely studied (Stutz and Reboud (1997a,b); [START_REF] Callenaere | The cavitation instability induced by the development of a reentrant jet[END_REF]; [START_REF] Dular | Relationship between cavitation structures and cavitation damage[END_REF]). Sheet cavity grows from the Venturi throat until a re-entrant jet appears in the sheet cavity closure, and flows upstream, near the wall, below the cavity. The reentrant jet then separates the sheet cavity when it reaches the liquid-gas interface of the flow. A large vapour cloud is shed and convected downstream, while the sheet cavity length is substantially reduced. Then, a cloud of vapor collapses in a higher pressure zone and sustains the re-entrant jet because of cloud implosion, which later produces a new shedding of vapor cloud [START_REF] Callenaere | The cavitation instability induced by the development of a reentrant jet[END_REF]; De Lange and De Bruin (1997)), that is characterized by the shedding frequency f s .

Studies about the sheet cavity cycle observed on Venturi profiles have resulted in the definition of two types of oscillatory behaviors [START_REF] Sayyaadi | Instability of the cavitating flow in a venturi reactor[END_REF]). A cloud cavitation regime appears when the sheet cavity is large enough and generates periodical shedding which generates vapor clouds before collapsing. When the sheet cavity length is shorter because of a higher pressure level and/or a Reynolds number smaller than a critical value [START_REF] Keil | Cloud cavitation and cavitation erosion in convergent divergent nozzle[END_REF][START_REF] Keil | On the transition from sheet to cloud cavitation[END_REF]), no cloud cavitation shedding occurs. In this regime, the sheet cavity oscillates but the cavity closure position does not vary and no large scale cavitation structure shedding is observed. These two cavitation regimes are characterized by different values of the Strouhal number St L = Lf s /v throat where L is the cavity length (defined in section 4.1) and v throat is the flow velocity measured at the Venturi throat. For a sheet cavity regime, St L 0.1 whereas a cloud cavitation regime is character-ized by 0.2 ≤ St L ≤ 0.4 [START_REF] Coutier-Delgosha | Internal structure and dynamics of sheet cavitation[END_REF]; [START_REF] Dular | Scale effect on unsteady cloud cavitation[END_REF]). The above thresholds are however not universal and may depend on the geometry or on the Reynolds number.

Furthermore, the cloud cavitation regime is more aggressive and leads to severe damages of solid structures. If a passive control method allows to reduce the sheet cavity length, it may then be possible to control the cavitation regime, in order to promote nearly stable regime instead of periodical unsteady behaviors. This can be done by modifying foil surface roughness. This technique has been used to delay cavitation inception, using local protuberances with different geometries or distributed irregularities [START_REF] Arndt | Rough surface effects on cavitation inception[END_REF]; [START_REF] Arndt | Cavitation in fluid machinery and hydraulic structures[END_REF]). These modifications act on the turbulent boundary layer, which drives the onset of sheet cavitation. The role of surface condition has also been studied in a fully developed cavitating flow. Most of works about random distributed, transversal distributed, or local roughness show that it is possible to decrease sheet cavity length, in order to increase oscillation frequency of the cavity or even to change cavitation regime [START_REF] Kawanami | Mechanism and control of cloud cavitation[END_REF]; [START_REF] Stutz | Influence of roughness on the two phase flow structure of sheet cavitation[END_REF]; [START_REF] Coutier-Delgosha | Effect of wall roughness on the dynamics of unsteady cavitation[END_REF]; [START_REF] Yongjian | Effects of grooves on cavitation around the body of revolution[END_REF]). But, to the best of our knowledge, there is no study about the effect of organized longitudinal roughness like grooves on the cavitation regime. This passive control method, applied on the suction side surface of a Venturi-type section with convergent and divergent angles of respectively 18 • and 8 • is discussed in this work.

In state-of-the-art studies, the different cavitation regimes are detected by observing the evolution of the sheet cavity development on acquired images. The lack of any rigorous, robust and efficient method to detect the change of cavitation regime handicaps studies about cavitation modeling or about passive control. In order to obtain a thorough analysis of the passive control effect on the sheet cavity dynamics, Proper Orthogonal Decomposition (POD) calculation [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF]) is applied on images derived from fast visualization of the sheet cavity development. This method has already been applied on cavitating flows to characterize several differences between cavitation models (Utturkar et al. (2005a,b)) using multiscale differences of dynamics provided by each model. Other studies have shown the interest of the use of this type of fluid mechanics analysis tools for cavitating flows to better understand the sheet cavity cycle, with wavelets [START_REF] Kjeldsen | Joint time frequency analysis techniques: a study of transitional dynamics in sheet/cloud cavitation[END_REF]) or POD and Dynamic Mode Decomposition (DMD) [START_REF] Prothin | Traitement d'image utilisant la pod et la dmd pour l'étude du développement de la cavitation sur un naca0015[END_REF]). After a description of the POD method and the advantages of its application on the sheet cavity regime study, results obtained with this rigorous analysis of the cloud cavitation shedding passive control with longitudinal grooves are presented to show the efficiency of this new cavitation regime detection technique.

Proper Orthogonal Decomposition

Proper Orthogonal Decomposition (POD) is a powerful tool to identify spatial and temporal coherent structures in a flow, providing a basis for a modal decomposition of a set of functions [START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF]). In other domains, this procedure is named Karhunen-Loeve decomposition, Principal Components Analysis, Singular Systems Analysis or Singular Values Decomposition and finds applications in image processing, signal analysis, data compression, oceanography or fluid mechanics. [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] was the first to transpose this technique in turbulence studies and to apply this post-processing on velocity fields [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF]). Large scale structures, which contribute for a large part to the total kinetic energy of a velocity field, may then be highlighted and selected to reconstruct flow dynamics filtering stochastic events [START_REF] Patte-Rouland | Proper orthogonal decomposition used for determination of the convection velocity of the initial zone of the annular jet. aerodynamic study and control of instabilities[END_REF]). By definition, POD is an optimal technique in an energetical sense that only a few eigenfunctions are necessary to recover flow dynamics. This linear process is based on spectrum theory of compact, self-adjoint operators. Eigenfunctions ϕ are chosen in order to maximize the average projection of a data field u on ϕ:

max ϕ∈L 2 ([0,1]) | (u, ϕ) | 2 ϕ 2 (1)
where |.| is the modulus, . may be a spatial or a temporal averaging operation and . is the L 2 -norm defined as:

f = 1 0 |f (x) | 2 dx 1 2
Maximization of | (u, ϕ) | 2 with the constraint ϕ 2 = 1 leads up to the corresponding functional:

J [ϕ] = (u, ϕ) 2 -λ ϕ 2 -1 , (2)
which obeys for all variations ϕ + δψ to:

d dδ J [ϕ + δψ] | δ=0 = 0. (3)
This condition reduces to the integral equation:

1 0 u (x) u * (x ) ϕ (x ) dx = λϕ (x) . (4)
The eigenfunctions {ϕ j } produce the optimal basis derived from Eq. (4) using the averaged autocorrelation function: R (x, x ) = u (x) u * (x ) .

(5)

Using Hilbert-Schmidt theory, we can assume that there is an infinity of orthogonal eigenfunctions associated to eigenvalues obtained by the diagonal decomposition of the averaged autocorrelation function. As the averaged autocorrelation R (x, x ) is non-negative, eigenvalues λ j are ordered so that:

λ j ≥ λ j+1 ≥ 0. ( 6
)
Each data field u may then be reconstructed from the modal decomposition based on the eigenfunctions ϕ j :

u (x) = ∞ j=1 a j ϕ j (x) (7)
where a j is a reconstruction coefficient defined as a j a * k = δ jk λ j . If studied data fields are velocity fields, eigenvalues represent twice the average kinetic energy in each eigenmode ϕ j . In this way, first modes, which are associated to the largest eigenvalues, refer to the most energetic structures of the flow field. POD is a method to extract disturbances in the flow at different energetic scales in the case of a turbulent velocity field application.

Here, POD is used to study sheet cavity dynamics by calculating the decomposition of sequences of images. In order to improve the calculation time, the snapshots method has been chosen [START_REF] Sirovich | Method of snapshots[END_REF]). Indeed, in the classic method proposed by [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF], . is a temporal average. R (x, x ) is then the spatial correlation tensor. In this case, when the N fields u present M space points each with n c vector components, the size of the eigenvalues problem (4) is N P OD = M n c . But if N << M , the snapshot methods allows to reduce the calculation time by using . as a spatial average, and the correlation tensor R is replaced by the temporal correlation tensor C (t, t ). In the snapshot method, the size of the eigenvalues problem becomes N P OD = N .

Experimental set-up

Experiments were conducted in the closed loop test rig of the DynFluid laboratory water tunnel. Two storage tanks, with a capacity of 4 m 3 each, provide water in the rig. The water flow is driven by a centrifugal pump. The test-section used for experiments is a 120× 100× 800 mm 3 volume (Figure 1). A Venturi profile is placed in the test section with convergent and divergent angles respectively 18 • and 8 • because this configuration presents an unstable cavitation regime which is the most inconvenient case for industrial applications [START_REF] Aeschlimann | Unsteady cavitation analysis using phase averaging and conditional approaches in a 2d venturi flow[END_REF]). At the Venturi throat, the minimum height of the flow is H throat = 67 mm. All positions in the test section are expressed with non-dimensional values x * = x/H throat , ties are acquired and analyzed. In order to visualize sheet cavity, 4000 images are acquired in each configuration. A CamRecord 600 Optronics camera with a Zeiss Makroplanar T* objective lens is used. Parameters used for visualization measurements are presented in Table 1. The flow is illuminated in a continuous mode by a Super Long Life Ultra Bright (SLLUB) White Led Backlight from Phlox on a 200 mm × 200 mm light output area. Its minimal luminance is 3000 cd.m -2 and its uniformity is 99.54%. The bottom wall of the test section, in the divergent part of the Venturi, is made of 7 mm thick, 242 mm long, and 120 mm wide interchangeable brass plates, to cover the entire Venturi surface downstream of the Venturi throat as it is presented in figure 3. The junction between the plate and the Venturi basis is done upstream of the Venturi throat to minimize disturbances on the flow. Plates are screwed on the basis and sealing between both components is ensured with vacuum grease. Table 2 presents the different grooved plates used for experiments. Plate A will be considered as reference, with a smooth surface on the suction side of the Venturi. Plates B to E are grooved surfaces made with a d = 1 mm or 2 mm diameter bur, with a depth h and the width of the ridge is e. The distance between the middle of two adjacent grooves is equal to λ = e+2 dh -h 2 for plates B to C and λ = e + d for plates D and E. Grooves open to the throat, so roughness affects the flow already from the onset of the sheet cavity.

d(mm) h(mm) e(mm) λ(mm) N Symbols for graphics A 0 0 0 0 0 × B 1 0.25 0.1 0.97 124 C 2 1 0.1 2.1 57 + D 2 2 1 3 40 ♦ E 2 2 0.1 2.1 57 2
Table 2: Characteristics of the studied plates used on the Venturi suction side

Detection of the cavitation regime with POD calculation in the reference case of a smooth Venturi

In order to highlight the efficiency of the POD method to detect the cavitation regime, the first step is to figure out the standard cavitation regime detection method. This method is described to analyze cavitation regimes obtained with the reference smooth case using post-processes cavitation images of a reference smooth plate. Then, these results are compared to POD analysis of the reference smooth case. The aim of this study is to explain the robustness of this efficient method.

Standard method of cavitation images post-processing

The cavitation regime of a sheet cavity developed on a Venturi profile can be determined by studying the Strouhal number associated with the variations of the closure position of the sheet cavity. This position is calculated with side view images of the sheet cavity, acquired with a camera. Mean values and standard deviation of the grey levels of 4000 images in the (yOz) plane (side views of the attached sheet cavity as shown in Fig. 1) are calculated.

The first step is to normalize instantaneous images I by a reference image I 0 (acquired with a non cavitating flow at rest). This normalization consists of the calculation of (I -I 0 ) /I 0 . The closure of the sheet cavity is localized with the maximum value of the standard deviation of the grey levels calculated with the images stack [START_REF] Dular | Relationship between cavitation structures and cavitation damage[END_REF]), as illustrated in figure 4. If we look at the profile of the root mean square of normalized images RMS(I) (figure 4c), we can see a first peak in about y * = 0.7. This value of RMS(I) corresponds to the vibration of the top part of the sheet cavity, close to the Venturi throat. At the beginning of the sheet cavity, grey levels oscillate between very low values (near 0 for the liquid phase) and very high values (near 1 for the vapor phase) [START_REF] Dular | Relationship between cavitation structures and cavitation damage[END_REF]). Filters have been applied in the image processing to reduce this peak, in order to detect the maximum of RMS(I), corresponding to the closure position of the sheet cavity. The mean cavity length L is then calculated from the Venturi throat to the closure line. The uncertainty of L, evaluated with the image processing by studying the effects of the threshold level for the binarization and filters applied on images is inferior to 2% (evaluated with the image processing by studying the effects of the threshold level for the binarization and filters applied on images). Considering the smooth plate A, mean cavity lengths obtained for different cavitation numbers σ ranging from 1 to 1.9 with an uncertainty of 1% are presented in figure 5(a). When σ decreases, the cavity length increases [START_REF] Callenaere | The cavitation instability induced by the development of a reentrant jet[END_REF], [START_REF] Dular | Scale effect on unsteady cloud cavitation[END_REF]). In order to study the different regimes of the sheet cavity, frequencies of sheet cavity oscillations have been determined using Fast Fourier Transform (FFT). A study of the grey levels at the sheet cavity closure position allows the calculation of the power spectral density associated with this phenomenon to determine the characteristic frequency of the sheet cavity. The uncertainty on the sheet cavity frequency is 0.2%, according to convergence and images processing tests (number of images used, choice of grey level threshold...). Strouhal numbers associated to these sheet cavity frequencies are then calculated as St L = Lf L /v throat , with 6% of uncertainty. Figure 5

POD post-processing of cavitation images

POD is applied on normalized images of the side view of the sheet cavity. A convergence test allows to determine the number of snapshots necessary to validate the decomposition. Figure 6 presents the contribution of the first modes j in the reconstruction of instantaneous images calculated as λ j / N j=1 λ j , from the eigenvalues λ j obtained by the decomposition with N snapshots. This analysis shows that 600 snapshots are sufficient for the decomposition to obtain converged results.

As the optimal basis is given by the eigenfunctions ϕ j of Eq. ( 4), if u (x) is a turbulent velocity field, the eigenvalues λ j represent twice the average kinetic energy in each mode ϕ j . In this work, u (x) is a grey level field and then, λ j is an information about the "weight" of modes in the reconstruction.

Figure 7 shows the first modes ϕ j (regardless the reconstruction coefficient a j ) obtained with POD of images in the smooth case study, for σ = 1.17 and σ = 1.60. These modes represent two different cavitation regimes according to the classic method of images post-processing (Fig. 5). As the sheet cavity length is small for σ = 1.60, it is difficult to discuss this result, but we can see that contributions of each mode from modes 1 to 3 are smaller compared to the case with σ = 1.17 (contributions are inferior to 0.3% for σ = 1.6 and superior to 0.6% for σ = 1.17). In the case of a sheet cavity regime, without a cloud cavitation shedding, the first mode 0 is almost sufficient to reconstruct the sheet cavity dynamics.

To obtain the temporal behaviour of the sheet cavity and the cloud cavitation, we study phase portraits of the temporal amplitude a j , corresponding to the reconstruction coefficient of the mode ϕ j . Reconstruction coefficients a 2 ) has a frequency twice the frequency of the sine wave a 3 . The mode 3 is relative to the cloud cavitation in the case of σ = 1.17. This result is very different from the sheet cavity regime presented for σ = 1.60 on figure 8(b'). These phase portraits are then an efficient tool to detect the cavitation regime, using a robust mathematics method.

5. Analysis of the longitudinal grooves effects on the cavity sheet developed on the divergent side of the Venturi 5.1. Effects of the grooves on the sheet cavity closure position In order to discuss the ability of longitudinal grooves to change the cavitation regime of the sheet cavity, the sheet cavity length is first determined in different operating conditions (grooves geometry or pressure conditions characterized by the cavitation number σ). Mean cavity lengths obtained for the different grooves geometries (using all plates presented in Table 2), and for cavitation number σ ranging from 1 to 1.9 are presented in Figure 9. As we have seen previously, the mean cavity length decreases when the cavitation number increases and all grooved plates follow the same trend. The geometric parameters of different plates have effects on the mean cavity length for a fixed value of σ. For large grooves diameter (d = 2 mm), the grooves depth depends is a crucial parameter: the comparison of results obtained for plates C and E, which have same diameters, same spacing (number of grooves is N = 57 for both cases) allows to measure the depth effect. For the plate C with h = 1 mm (h * = h/H throat 0.016), the cavity mean length is larger than for other plates, for all 1 σ 1.9. On the other hand, plate E has a larger grooves depth h = 2 mm (h * 0.031) and presents smaller mean cavity lengths than the reference plate A. As we can see in literature, the cavity regime seems to depend on the mean length of the sheet cavity [START_REF] Callenaere | The cavitation instability induced by the development of a reentrant jet[END_REF]), then it should be possible to modify the cavity regime by changing grooves geometry on the suction side of the Venturi-type section. Strouhal numbers, calculated for each grooved plate, are presented on Figure 10. For grooved plates, the range of cavitation numbers for each regime (defined with 0.2 St L 0.4 for the unstable regime and with 0.05 St L 0.15 for the stable regime) is similar to the smooth case, except for plates D and E. For these two plates, the unstable regime is suppressed for the entire studied range 1 σ 1.8 (see Media 1,2 proposed with Figure 11). Only one point for plate E gives a Strouhal number St > 0.5. This Strouhal number is calculated with the frequency f L = 44 Hz, which is an harmonic of the sheet cavity oscillation frequency f L = 9 Hz. Then, this particular point corresponds also to a sheet cavity regime.

Figure 11: Visualisation of (a-j) the cloud cavitation shedding for the smooth plate A and (a'-j') the sheet cavity pulsation for the plate E, when σ =1.17 (∆t = 2 ms between two images and colorbars represent normalized grey level). See Media 1-2

Dynamics analysis of the sheet cavity by POD

POD is applied on 600 normalized images for the different grooved plates. Figure 12 shows the evolution of λ j according to modes, for different plates, with σ = 1.17 in order to discuss the contribution of each mode in the reconstruction. Eigenvalues for the first mode 0 are 84.94% for the smooth case A, 84.74% for plate C, 74.27% for plate B and 92.95% for plate E. For the other modes 1 to 50, these eigenvalues are lower for the plate E than for the other plates. This result highlights that for σ = 1.17, the cavitation regime is different for plate E. As this regime is characterized only by a small variation of the closure position of the cavity sheet, the mode 0 is the most important in the reconstruction, in comparison with the other plates. First modes ϕ j are presented for plates A, B, C and E on figures 13 and 14. Results for plate D are not presented on these figures because they are very similar to plate E. As these figures present only modes, regardless the reconstruction coefficient a j , the results interpretation does not allow to compare the order of magnitude of each mode in these different cases. We focus then on the position of vapor and liquid zones which are represented by different grey levels in the instantaneous images used for the decomposition, and on the global shape of modes. We can see that these first modes are quite similar for plates A and B, even if the sheet cavity is shorter for the grooved plate B. The modes 0 to 3 represents 90% of the reconstruction of the instantaneous images for these two plates. For the plate C which is characterized by a longer sheet cavity for the same cavitation number, modes are different. Only modes 0 to 2 are necessary to obtain 90% of the images reconstruction in this case. These modes are not very different from the first two cases, but from the mode 4, the morphology of modes are different. Results obtained for the grooved plate E show that the first mode 0 is sufficient to obtain 90% of the reconstruction. We can see that downstream of the closure position of the sheet cavity, there is no structure representing bubbles. This result confirms the previous analysis: no cloud cavitation is convected with plate E for the entire range of studied cavitation number. The analysis of eigenmodes can then identify the type of cavitation regime, using the spatial behavior of the sheet cavity dynamics.

We then study phase portraits of the temporal amplitude a j , corresponding to the reconstruction coefficient of the mode ϕ j , as presented in figures 15 and 16. Reconstruction coefficients a 2 (mode 2) according to a 1 (mode 1) on figure 15 provide phase portraits similar to Lissajous figures [START_REF] Oberleithner | Threedimensional coherent structures in a swirling jet undergoing vortex breakdown: stability analysis and empirical mode construction[END_REF]) for plates A and B but not for plate E. This result shows again the efficiency of POD to detect the cavitation regime and then to measure the effect of a sheet cavity passive control.

Figure 16 presents phase portraits with the reconstruction coefficients a 3 for mode 3 according to a 1 for mode 1. The mode 3 is relative to the cloud cavitation in this case (as it is explained in the section 4.2 with the analysis of the Lissajou's figures). We can see on figure 16 that the smooth plate A and the grooved plate B present a cloud cavitation regime. For the grooved plate C, the result is less characteristic of this regime than the smooth plate, in contrary to the result obtained previously with the analysis of the gray level variations. This information shows that the grooves depth or the wavelength of grooves acts on the sheet cavity dynamics (λ > 2 mm and h > 1 mm for plates C, D and E). For plate E, the phase portrait a 3 = f (a 1 ) shows the sheet cavity regime. We can see on Figures 17 (a) and (c) the temporal evolution of the reconstruction coefficient a 1 corresponding to the mode 1, which is linked to the cloud cavitation, for plates A and E respectively. Results show the difference between these two plates: in the smooth case, the temporal

Conclusion

Side view images of the sheet cavity dynamics on a converging-diverging profile have been analyzed with a POD method, in order to detect the effects of grooved surfaces, with different geometries. Studying the reconstruction coefficients a j calculated with the eigenvalues λ j , we can identify the sheet cavity regime with no cloud cavitation convected, which allows a reduction of surface erosion. Analyzing the energy content of the second mode 1 allows to detect this regime: in the case of a sheet cavity regime the energy content of mode 1 is lower than 1% (Fig- ures 7(b) and13(d)), while this value is superior to 1% for a cloud cavitation regime (Figures 7(a), 13(a), 13(b) and 13(c)). This threshold value is not sufficient to obtain a robust method. However, the study of phase portraits allows to differentiate regimes with the observation of Lissajous figures. Then, this tool is efficient to characterize the cavitation regime but also to calculate the cloud cavitation shedding, avoiding the sheet cavity length calculation.

Proper Orthogonal Decomposition allows to study grooves effects. Results show that grooves with a large depth (h > 1 mm) and a large wavelength (λ > 2 mm) change the sheet cavity dynamics. For example, plates D and E (h = 2 mm and λ 2.1 mm) suppress the cloud cavitation: for the entire studied range of cavitation number 1 σ 1.8 we obtain a sheet cavity regime, characterized by a Strouhal number 0.05 St L 0.15. For the other plates, we can have a cloud cavitation regime, for 1 σ 1.3, characterized by 0.2 St L 0.4. As we can see with results obtained with different grooved plates, the cavitation regime is related to cavity length, but does not depend only on it (for examples plates D and E). Avoiding cloud cavitation, which is the most aggressive phenomenon, can reduce erosion of surfaces. It is then interseted to modify the boundary layer with longitudinal grooves (with chosen values of depth and grooves wavelength). But this analysis is possible only with an efficient tool to identify the two cavitation regimes (stable or unstable regimes), based on a robust method, as presented in this study. In order to complete the analysis of the sheet cavity dynamics, it could be interesting to deepen these results with an analysis of top views of the sheet cavities on these grooved plates.

Figure 1 :

 1 Figure 1: Test section of the experimental device using a Venturi profile (H throat =67 mm is the reference length and y * ref = y ref /H throat = -3.7).

Figure 2 :

 2 Figure 2: Characterization of the flow in the test section inlet in y * -4.5: (a) non-dimensional longitudinal velocity profile v * y = vy /v ref and (b) the profile of the turbulence intensity RMS (vy ) <vy > , where < vy > is the time-averaged velocity.

Figure 3 :

 3 Figure 3: Venturi components: (a) Venturi basis with C1 to C7 pressure sensors on the basis surface (untapped in this study), (b) junction between the Venturi basis and a plate and (c) zoom of a grooved plate with the definition of the grooves geometric parameters (d and h are respectively the grooves diameter and depth, e is the width of the ridge and λ is the grooves wavelength).

Figure 4 :

 4 Figure 4: Sheet cavity on the smooth plate A for σ =1.18: (a) normalized average image, (b) root mean square of normalized images (the line represents the non-dimensional sheet cavity length L * = L/H throat 1.62) and (c) a profile of the root mean square of normalized images RMS(I ) plotted in z*=0.

Figure 5 :

 5 Figure 5: Standard post-processing images applied on the reference smooth plate A: (a) non-dimensional sheet cavity length L * = L/H throat and (b) Strouhal numbers St L = Lf L v throat according to the cavitation number σ

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Convergence of the POD calculation for the plate C, with σ = 1.37: energy of modes (a) 0, (b) 1 and (c) 2 according to the number of images used for the POD

Figure 9 :

 9 Figure 9: Non-dimensional sheet cavity length L * = L/H throat according to the cavitation number σ: ×: plate A, : plate B, +: plate C, ♦: plate D, 2: plate E.

Figure 10 :

 10 Figure 10: Strouhal number St L = Lf L v throat according to the cavitation number σ: ×: plate A, : plate B, +: plate C, ♦: plate D, 2: plate E.

Figure 12 :

 12 Figure 12: Contribution of POD modes in the reconstruction calculated with the eigenvaluesλ i Σ N i=1 λ i, for σ = 1.17 with ×: plate A, :

Figure 15 :

 15 Figure 15: Phase portrait of POD reconstruction coefficients for modes 1 and 2, with σ = 1.17: (a) plate A, (b) plate B, (c) plate C and (d) plate E.

Figure 13 :Figure 14 :

 1314 Figure 13: POD modes and their contribution in the reconstruction of instantaneous images (in percentages) for modes 0 to 3, with σ = 1.17: (a) plate A (L * = 1.52), (b) plate B (L * = 1.14), (c) plate C (L * = 1.75) and (d) plate E (L * = 1.23).

Figure 16 :

 16 Figure 16: Phase portrait of POD reconstruction coefficients for modes 1 and 3, with σ = 1.17: (a) plate A, (b) plate B, (c) plate C and (d) plate E.

Figure 17 :

 17 Figure 17: Reconstruction coefficient a 1 of the mode 1: (a) temporal evolution for plate A, (b) frequencies spectrum for plate A, (c) temporal evolution for plate E and (d) frequencies spectrum for plate E.

Table 1 :

 1

	(a)	(b)	(c)

Parameters of visualizations acquired with a CamRecord 600 Optronics camera and a Zeiss Makro-planar T*

  (b) presents values of Strouhal numbers St L for each studied plate, with respect to the cavitation number σ. Two different sheet cavity regimes can be identified. The first regime is called unstable regime and is defined by 0.2

				St L	0.4 which
	corresponds to a periodic cloud cavitation shedding, with
	large variations of the sheet cavity length. In the case of
	a smooth plate A, this cloud cavitation regime is encoun-
	tered for 1	σ	1.4. The second regime, called stable
	regime, is characterized by a small oscillation of the cavity
	closure, with an almost constant sheet cavity length and
	no convected cloud cavitation. Strouhal numbers in the
	stable regime are ranging from 0.05	St L	0.15. This
	regime is obtained for 1.4 σ 2 in the case of a smooth
	plate.