Hyunsuk Hong

Hugues Chaté

Lei-Han Tang

Hyunggyu Park

Finite-size scaling, dynamic fluctuations, and hyperscaling relation in the Kuramoto model

Keywords: numbers: 89.75.-k, 05.45.Xt, 05.45.-a

We revisit the Kuramoto model to explore the finite-size scaling (FSS) of the order parameter and its dynamic fluctuations near the onset of the synchronization transition, paying particular attention to effects induced by the randomness of the intrinsic frequencies of oscillators. For a population of size N , we study two ways of sampling the intrinsic frequencies according to the same given unimodal distribution g(ω). In the 'random' case, frequencies are generated independently in accordance with g(ω), which gives rise to oscillator number fluctuation within any given frequency interval. In the 'regular' case, the N frequencies are generated in a deterministic manner that minimizes the oscillator number fluctuations, leading to quasi-uniformly spaced frequencies in the population. We find that the two samplings yield substantially different finite-size properties with clearly distinct scaling exponents. Moreover, the hyperscaling relation between the order parameter and its fluctuations is valid in the regular case, but is violated in the random case. In this last case, a self-consistent mean-field theory that completely ignores dynamic fluctuations correctly predicts the FSS exponent of the order parameter but not its critical amplitude.

I. INTRODUCTION

We revisit the Kuramoto model [1] to investigate the phase order parameter and its dynamic fluctuations near the onset of the synchronization transition. In particular, we focus on their critical finite-size scaling (FSS) behavior. Even though the Kuramoto model and its extended versions have been widely studied [2,[START_REF] Acebrón | [END_REF], some issues persist, in particular about the above points. In the few existing studies [4,5], the critical properties near the transition point have only been partially resolved, which prompted the present work.

FSS is particularly important for data analysis, since real systems are always finite, and particularly so in many experimental systems, biological or not [6]. Recently, the FSS of the synchronization order parameter in the Kuramoto model and several extended versions on complex networks has been investigated. An unusual FSS exponent ν = 5 2 has been found [7][8][9]. This value was obtained in the classic case, hereafter referred to as the random distribution, where the individual frequencies of oscillators are sampled independently according to a given unimodal distribution g(ω). In such a situation, oscillator frequencies can be arbitrarily close to each other even in a finite population. Interestingly, when the oscillator density fluctuation along the frequency axis is suppressed in a population with regular distribution, a very different FSS exponent ν = 5 4 is obtained. Here we present details of our numerical study that led to this previously reported value [10,11].

Dynamic fluctuations of the order parameter in the Kuramoto model is another long-standing problem [START_REF] Acebrón | [END_REF]12]. In the perturbation theory developed by Daido [4], it is claimed that they scale differently when approaching the synchronization threshold from above and from below, with the scaling exponent γ = 1 4 in the supercritical region and γ = 1 in the subcritical region. This result is quite surprising by itself, in the sense that conventional scaling theory predicts identical scaling behavior on either side of the transition. Note that Daido used the regular frequency distribution in numerical investigations as well as in his perturbation theory [4].

In this paper, we report extensive numerical results on the FSS of the order parameter and its dynamic fluctuation in the random and regular realizations of the intrinsic frequency distribution. The FSS exponent in the regular case is estimated numerically as ν 5 4 , half the value of its random counterpart. In contrast to Daido's predictions, we find that dynamic fluctuations scale with the same exponent γ = γ on both sides of the transition. Interestingly, our data supports γ γ 1 for the random frequency distribution and γ γ 1 4 for the regular case. Given the order parameter exponent β = 1 2 , the hyperscaling relation γ = ν -2β is thus violated in the random case, but obeyed in the regular case. We finally briefly discuss the implications of our observations on the ordering process and dynamic response in the Kuramoto model.

The paper is organized as follows. In Sec. II, we introduce the Kuramoto model and present the two types of arXiv:1503.06393v3 [cond-mat.stat-mech] 25 Aug 2015 frequency realizations that we consider in this paper. In Sec. III, we investigate the FSS of the order parameter and its dynamic fluctuations with the random frequency distribution and discuss the validity of the hyperscaling relation in the model. In Sec. IV, we repeat the investigation and analysis with the regular distribution. We summarize and discuss the results in Sec. V. Appendix A presents the mean-field approximation to the FSS in the case of the random frequency distribution.

II. THE KURAMOTO MODEL

The Kuramoto model [1] consists in a finite population of N globally-coupled phase oscillators governed by

dφ j dt = ω j - K N N k=1 sin(φ j -φ k ), j = 1, • • • , N. (2.1)
Here φ j is the phase of the jth oscillator, and ω j its intrinsic frequency. A given realization of the Kuramoto model corresponds to a particular choice of the N intrinsic frequencies {ω j }. In this paper, we shall assume that the number density of these frequencies on the real axis follows, on average, the normal distribution

g(ω) = (1/ √ 2π)e -ω 2 /2 . The coupling strength K is pos- itive (K > 0).
Phase synchronization is conveniently described by the complex order parameter defined as [1]

Z(t) ≡ ∆e iθ = 1 N N j=1 e iφj , (2.2)
where a nonzero (positive) value of ∆ in the asymptotic N → ∞ limit implies the emergence of phase synchronization, and θ is the phase of the global order parameter. With this definition, Eq. (2.1) is rewritten:

φj = ω j -K∆ sin(φ j -θ), (2.3)
where both ∆ and θ are dynamic variables defined by Eq. (2.2). In the limit N → ∞, the model exhibits a continuous symmetry-breaking transition at K = K c = 2/[πg(0)] = 8/π [1]. In this limit, the phase order parameter is strictly zero on the subcritical side (K < K c ), while on the supercritical side

(K > K c ), ∆ ∼ (K -K c ) β with β = 1 2 . (2.4)
This result holds for any unimodal and symmetric g(ω).

For a finite population of oscillators, ∆(t) exhibits both dynamic and sample-to-sample fluctuations whose scaling with N is the main focus of this work. The sampleto-sample fluctuations arise from the randomly drawn oscillator frequencies from the distribution g(ω). Consequently, the spacing between neighboring frequencies on the real axis is Poissonian in any given realization. It turns out that this quenched noise and the nonlinear dynamic fluctuations are intimately related to each other in the sense that, when the former is removed, a very different set of FSS exponents are obtained. The two types oscillator frequency distributions are discussed in Sec. III and IV, respectively.

III. RANDOM FREQUENCY DISTRIBUTION

In this section, we consider the realizations (samples) of frequency sets {ω j } by randomly and independently drawing frequencies from the Gaussian distribution g(ω). We measure the phase order parameter and its dynamic fluctuation for each sample and average over many samples. In practice, we numerically integrate Eq. (2.3), for system sizes up to N = 12 800, using the Heun's method [START_REF] See | Numerical Analysis[END_REF] with a discrete time step δt = 0.01, up to total 2 × 10 5 time steps. We start with random initial values of {φ j (0)} and average the data over time after the first 10 5 steps to avoid any transient behavior and also average over 10 2 -10 3 different samples of {ω j } and {φ j (0)}. At K = K c , we perform numerical integration for larger systems up to N = 204 800 with twice more realizations.

A. Finite-size scaling of the order parameter

In analogy with equilibrium critical phenomena, we expect the order parameter in the critical region to satisfy the FSS [START_REF] Um | [END_REF],

∆(K, N ) = N -β/ν f ( N 1/ν ) , (3.1)
where = K -K c and the scaling function f (x) is a monotonically increasing function with the following asymptotic properties,

f (x) ∼ x β , x 1; (-x) β-ν/2 , x -1. (3.2)
The exponent ν is known as the FSS exponent. For large N , this scaling ansatz reproduces Eq. (2.4) on the supercritical side, and yields ∆ ∼ N -β/ν at the transition and ∆ ∼ N -1 2 in the subcritical regime, respectively. In the usual mean-field treatment of the Kuramoto model, the entrained state is established by solving Eq. ( 2.3) at a constant ∆ and θ (dynamic fluctuations are ignored): Oscillators with |ω j | < K∆ get into an entrained (locked) state, reaching a fixed angle. On the other hand, the oscillators with |ω j | > K∆ perform periodic motion with modified frequencies. With this, the self-consistent equation for ∆ from Eq. (2.2), can be written as

∆ = K∆ -K∆ dωg(ω) 1 -(ω/K∆) 2 (3.3) in the N → ∞ limit. Eq. (3.3) has a nontrivial solution (∆ > 0) for K > K c = 2/[πg(0)].
For a finite-size system, the right-hand side of Eq. (3.3) is replaced with

∆ = Ψ(K∆) ≡ 1 N j,|ωj |<K∆ 1 -(ω j /K∆) 2 , (3.4) which contains a sample-dependent correction δΨ ≡ Ψ -[Ψ] ∝ √ N s /N [8].
Here N s is the number of (entrained) oscillators in the frequency interval (-K∆, K∆) and [•] denotes sample average. As was shown previously in Refs. [8,10] and explained in detail in the Appendix A, for K close to K c , Eq. (3.4) admits a scaling solution of the form (3.1) with the exponents

β = 1 2
, and ν = 5 2 .

(3.5)

The scaling function f (x), however, has significant sample-to-sample fluctuations.

The exponent value β = 1 2 is typical of various equilibrium mean-field models described by the so-called φ 4 theory. However, the FSS exponent value of ν = 5 2 is quite unusual and implies much stronger fluctuations than that in ordinary models with the conventional value ν = 2. This remarkable result has been confirmed by extensive numerical simulations [7,8] and the straightforward extension to the Kuramoto model on various sparse random networks in [9,11,[START_REF] Um | [END_REF].

B. Dynamic fluctuations

Now, we consider the dynamic (temporal) fluctuations of the order parameter by measuring the quantity [15]

χ(K, N ) ≡ N (∆ -∆ ) 2 = N ∆ 2 -∆ 2 (3.6)
where • denotes time average in the steady state of a given sample and [•] sample average, respectively.

For K < K c , we expect from simple power counting that χ ∼ O(1), because ∆ ∼ N -1 2 . For K > K c , the order parameter becomes O(1), but these O(1) terms cancel out each other exactly and the next order terms proportional to N -1 2 contribute. Again, χ becomes O(1).

However, at K = K c , strong critical FSS occurs since ∆ ∼ N -β/ν = N -1 5
, and one can expect that χ will diverge in the N → ∞ limit. Summing up, the critical behavior of χ near the transition may be expressed, in the N → ∞ limit, as

χ ∼ -γ , > 0; (-) -γ , < 0. (3.7)
The two exponents γ and γ describe the divergence of χ in the supercritical ( > 0) and subcritical ( < 0) regions, respectively. For homogeneous systems, scaling in the critical region is controlled by a single fixed point in the renormalization group (RG) sense, and hence γ = γ . The corresponding FSS takes the form

χ( , N ) = N γ/ν h( N 1/ν ) , (3.8)
irrespective of the side of the transition, and the scaling function h(x) has the limiting behavior

h(x) ∼    const, x = 0; x -γ , x 1; (-x) -γ , x -1.
(3.9)

In this context, the Kuramoto model is somewhat special in that each oscillator has its own intrinsic frequency and the effect of the mean-field is different on different oscillators. Therefore the usual RG argument with only a few relevant variables may not apply and one should entertain the possibility that γ = γ . In such a situation, the FSS exponent ν should also be different on the two sides of the transition, so that the ratio relation γ/ν = γ /ν be satisfied, in order to maintain the continuity of χ at = 0. Indeed, through perturbative calculation, Daido [4] obtained γ = 1 4 and γ = 1. We have carried out extensive numerical simulations to obtain accurate estimates of γ and γ as well as other FSS exponents. Figure 1 shows numerical data for χ versus for various system sizes N . As expected, finitesize effects become huge near the transition point. We observe two important features of the finite-size effects: First, the peak position at K = K max is always on the subcritical side ( < 0) and approaches the bulk critical point ( = 0) as N increases. Second, the peak height χ max grows rapidly with N . By analyzing the data for the peak position and height, we estimate γ and ν in the subcritical side, using the conventional FSS ansatz

χ max ∼ N γ /ν , δK max ∼ N -1/ν (3.10)
Even though our estimate of 1/ν is rather weak due to huge statistical uncertainties in locating K max , it includes our analytic value of 1/ν = 2 5 within the error bar. Taking a simple ratio of the above exponent values, we estimate γ = 1.2(3). If we take 1/ν = 2 5 , then γ = 1.00(3). We also measure the value of χ at the bulk critical point K = K c . As K c is located to the right of the peak position K max (N ) for all N in this case, the FSS behavior at K c may be regarded as a continuation of the scaling on the supercritical side. Again, the FSS ansatz leads to

χ c ≡ χ( = 0, N ) ∼ N γ/ν .
(3.12)

Figure 3 shows that the numerical data agree almost perfectly with the FSS ansatz with γ/ν = 0.39(2) .

This result confirms the ratio relation γ/ν = γ /ν , but does not provide estimates for γ and ν separately.

To estimate γ and γ directly from off-critical data in Eq. (3.7), we examine the local slopes of χ against | | in double logarithmic plots for very large N . The effective exponent γ eff is defined as

γ eff ( ) ≡ d ln χ d ln | | , (3.14)
and the critical exponents are obtained from their asymptotic values as γ = lim →0 + γ eff ( ) and γ = lim →0 -γ eff ( ). Figure 4 shows the inverse of the effective exponent γ -1 eff as a function of . We find that γ -1 eff approaches γ -1 linearly as goes to zero from below in the subcritical region, with its asymptotic value γ = 1.01(3). This is consistent with our previous numerical estimates and also with the theoretical prediction of γ = 1 by Daido [4]. Note that γ = 1 has been also obtained in a systematic 1/N expansion [5], where the random frequency distribution is prerequisite.

On the supercritical side, γ -1 eff varies in a nonlinear fashion far from the transition, but approaches γ -1 linearly near small with the asymptotic value γ = 0.98 (5). This suggests strongly that γ = γ = 1, and definitely excludes Daido's value γ = 1 4 . With γ = 1, our numerical estimate in Eq. (3.13) lead to ν = ν = 5 2 . Summing up all these results, we conclude that our numerical data strongly support simple scaling exponents:

β = 1 2 , ν = ν = 5 2
, and γ = γ = 1 .

(3.15) However, a proper analytic treatment for the dynamic fluctuations in the supercritical side is lacking, and further investigations are thus needed. One should note that the dynamic fluctuations of the order parameter are not necessarily proportional to the susceptibility (response of the order parameter to an infinitesimal external field) in general nonequilibrium steady states due to probable violation of the fluctuation-dissipation theorem. It is almost trivial to derive the susceptibility exponents in the Kuramoto model, which turn out to be γ sus = γ sus = 1 [16]. Yet this accordance with the dynamic fluctuation exponents in Eq. (3.15) may be merely coincidental.

C. Hyperscaling relation

In equilibrium critical phenomenon, the equation

γ = ν -2β , (3.16)
with ν = dν, is usually referred to as a hyperscaling relation. Here d is the dimension of the system and ν is the correlation length exponent. Equation (3.16) follows from a simple power counting of Eq. (3.6) at the transition. The left-hand side of Eq. (3.6) is χ c ∼ N γ/ν , while the right-hand side is proportional to N 1-2β/ν if the leading order terms do not exactly cancel out each other.

It is obvious that our result Eq. (3.15) does not satisfy this relation. Violating the hyperscaling relation is only possible when the leading order terms in ∆ 2 and ∆ 2 cancel out exactly and the sub-leading order terms show up, yielding γ = 1 instead of 3 2 . A plausible conjecture for leading and sub-leading terms for the large N expansion of the two quantities is then . These data are expected to show straight lines. The least-square fitting yields two straight lines with the same intercept. Lines between data are only guides to the eye.

∆ 2 = N -2β/ν c + dN -1 5 , ∆ 2 = N -2β/ν c + eN -1 5 , (3.17
with c = c . Indeed, our numerical data show in Fig. 5 that c ≈ c = 1.140(2), d = 0.10(4), and e = -0.04(4). Then, the dynamic fluctuation χ at the transition scales as

χ c = (d -e)N 1-2β/ν-1 5 ∼ N γ/ν , (3.18)
yielding

γ = ν -2β - 1 2 = 1. (3.19)
The cancellation of the leading order (c = c ) in Eq. (3.17) would be expected if dynamic fluctuations only produce subleading order corrections to the solution of the mean-field equation (3.4) at the transition. Interestingly, as we show in Appendix A, the meanfield approximation yields

[ ∆ 2 ] = [ ∆ 2 ] 0.93N -2
which is smaller in amplitude as compared to the numerically determined value c = 1.14. Therefore it appears that dynamic fluctuations renormalize the strength of the quenched noise δΨ, but do not change its scaling with N . It is interesting to compare the violation of hyperscaling in the Kuramoto model to known cases in equilibrium phase transitions. One well-known such case is the incompatibility between the correlation length exponent ν = 1 2 at the Gaussian fixed point and the mean-field exponents β = 1 2 and γ = 1 above the critical dimension d c = 4 [17]. Violation of hyperscaling is also reported in the three-dimensional random-field Ising model whose ordering transition is controlled by a zero-temperature fixed point [18][19][START_REF] Nattermann | Spin Glasses and Random Fields[END_REF][START_REF] Vink | [END_REF]. In both cases, the energy associated with the symmetry-breaking field is much stronger than the thermal energy on large length scales, hence the order parameter fluctuations, as measured by the exponent γ, become less than what is predicted by the hyperscaling relation. In this sense, the weaker dynamic fluctuations of the order parameter in the Kuramoto model may be interpreted in a similar way, despite the fact that we are now in a nonequilibrium situation.

IV. REGULAR FREQUENCY DISTRIBUTION

In this section, we consider the 'regular' frequency set {ω j } with minimal disorder among frequencies. This set can be generated following a deterministic procedure given by

j -0.5 N = ωj -∞ g(ω)dω, (4.1)
with j = 1, . . . , N and g(ω) = (1/ √ 2π)e -ω 2 /2 . The generated frequencies are quasi-uniformly spaced in accordance with g(ω) in the population. As this set is uniquely determined for each N , sample disorder may result from initial conditions only. Note that the regular frequency distribution is considered in numerical attempts by Daido [4].

Similar to the case of the random distribution, we numerically integrate Eq. (2.3) for system sizes up to N = 12 800 with a discrete time step δt = 0.01, up to total 10 7 time steps. We start with random initial values of {φ j (0)} and average the data over the latter half of time steps and also average over 20 -100 different samples of {φ j (0)}. At K = K c , we perform numerical integration for larger systems up to N = 25 600 with 10 -100 samples.

A. Finite-size scaling of the order parameter We measure the order parameter and its square during numerical integrations and average over time in the steady state, and then average over samples with different initial values for the phase variables. These averages are denoted by [ ∆ ] and ∆ 2 , respectively. Figure 6 shows their size dependence at the transition ( = 0), exhibiting the nice power laws

[ ∆ ] ∼ N -0.39 (2) and ∆ 2 ∼ N -0.78 (3) . Given the value β = 1 2 , we obtain ν ≈ 5 4 . As a consistency check, we investigate the FSS relation in Eq. (3.1) for = 0. We plot ∆N β/ν versus N 1/ν for various system sizes N in Fig. 7. (For a moment, we drop the average brackets as ∆ ≡ [ ∆ ] for convenience.) The data show a perfect collapse with the choice of ν = 5 4 . We also considered the regular Lorentzian distribution given by g(ω) = γ ω /[π(ω 2 +γ ω

2 )] with the half-width γ ω , where the oscillator frequency is chosen as ω j = γ ω tan[jπ/N -(N + 1)π/(2N )] for j = 1, . . . , N . We find a similar behavior again with ν ≈ 5 4 .

B. Dynamic fluctuations

In this subsection, we present numerical data for the dynamic fluctuations χ in Fig. 8. Compared to the random distribution case, χ is very much reduced in magnitude and peaks only in a much narrower region around the transition. The small amplitude makes it difficult to determine the peak height and position accurately. However, with the exact information of K c = 8/π, it is numerically feasible to investigate the critical increase of χ c at the transition as a function of N . The data shown in Fig. 9 fit well to a power law χ c ∼ N γ/ν with γ/ν = 0.22 (2) .

(4.4)

Using ν = 5 4 , we find γ = 0.27(3), which is close to 1 4 . This value is consistent with Daido's theoretical result of γ = 1 4 [4]. We also estimate γ and γ directly from off-critical data by examining the effective exponents γ eff for large N . Figure 10 shows (4γ eff ) -1 versus . Following a similar procedure as in the random case, we find γ = 0.24(2) and γ = 0.25 (1). We also fit the off-critical data directly in Fig. 8, which agree very well with χ ∼ | | -0.25 . All the numerical results strongly suggest scaling exponent values as

β = 1 2 , ν = ν = 5

C. Hyperscaling relation

The hyperscaling relation, γ = ν -2β, is seen to hold in the regular case. Furthermore, there is no detectable "disorder fluctuation" arising from different initial phase values, leading to [ ∆ ]

2 ≈ ∆ 2 .

V. SUMMARY

In this paper, we revisited the well-known Kuramoto model and investigated the FSS of the order parameter and its dynamic fluctuations near the onset of synchronization transition. We find that critical scaling behavior crucially depends on the presence of disorder among intrinsic frequencies of oscillators. In the random case with disorder, it is shown analytically and numerically that the FSS exponent ν = 5 2 and the dynamic fluctuation exponent γ ≈ 1, in both the supercritical and the subcritical region. In the regular case without disorder, we find ν ≈ 5 4 and γ ≈ 1 4 . In both cases, the exponent ν is different from its conventional value ν = 2 in usual mean-field theory of homogeneous systems with global coupling.

In the random case, disorder in the oscillator frequencies broadens the critical region and drives order parameter fluctuations in a sample-dependent manner. The hyperscaling relation that describes the interplay between static (i.e. averaged) and dynamic order parameter fluctuations is violated. Interestingly, a modified meanfield theory that completely ignores dynamic fluctuations yields correct critical FSS of the order parameter, although the amplitude is renormalized by dynamic fluctuations. This is one of the main findings in this work. In comparison, numerically determined exponents γ and ν in the regular case do satisfy the hyperscaling relation. Further work is required to explain the particular values observed in our study.

FIG. 1 :

 1 FIG. 1: (Color online) Dynamic fluctuation χ against = K -Kc for various system size N , with Kc = 8/π. The data have been obtained from the average over 100 -1000 samples with different sets of {ωj} and initial phases {φj(0)}. Statistical errors are represented by symbol sizes.

FIG. 2 :FIG. 3 :FIG. 4 :

 234 FIG. 2: (Color online) (a) Behavior of χmax is shown as a function of N in a log-log plot. The straight line obtained through the least-square fitting method displays χmax ∼ N 0.40 . (b) Behavior of δKmax(= |Kmax -Kc|) is shown as a function of N in a log-log plot. The straight line obtained through the least-square fitting method displays δKmax ∼ N -0.33 .

  )

5 FIG. 5 : 5 ∆ 2

 5552 FIG. 5: (Color online) (a) Critical decay of ∆ 2 and ∆ 2 at Kc(= 8/π) are plotted against N in a log-log plot. These two data are so close to each other that it is not possible to discern two data in this plot. The straight line is the fitting line of 1.14N -0.40 . (b) N 2 5 ∆ 2 and N 2 5 ∆ 2 versus N -1 5. These data are expected to show straight lines. The least-square fitting yields two straight lines with the same intercept. Lines between data are only guides to the eye.

FIG. 6 :FIG. 7 :

 67 FIG. 6: (Color online) Critical decay of the order parameter [ ∆ ] and ∆ 2 at Kc(= 8/π) is plotted as a function of N in a log-log plot, which are represented by (red) open boxes and (blue) open circles, respectively. The slopes of two straight lines are -0.39 and -0.78, respectively.

(4. 2 )

 2 Following the FSS theory in Eqs. (3.1) and (3.2), we estimate the decay exponent as β/ν = 0.39(2) .(4.3)

FIG. 8 :

 8 FIG.8:(Color online) Dynamic fluctuation χ is plotted as a function of (= K -Kc) for various system size N for the regular distribution, where Kc = 8/π. The solid (red) line represents a fitting line for data with largest system size.

4 , 4 . ( 4 . 5 )FIG. 9 :FIG. 10 :

 4445910 FIG. 9: (Color online) Critical increase of χc at Kc(= 8/π) is plotted as a function of N in a log-log plot for the regular The slope of the straight line is 0.22.

Acknowledgments

This research was supported by the NRF Grants No. 2012R1A1A2003678 (H.H) and 2013R1A1A2A10009722 (H.P), and by the NSFC Grant No. 11175013 (L.-H.T). We thank Korea Institute for Advanced Study for providing computing resources for this work.

Appendix A: Finite-size scaling in the mean-field approximation

In this appendix we consider in some detail the solution to the self-consistent mean-field equation (3.4). For this purpose, it is convenient to introduce a function

In terms of ψ, we have,

Since each term in the sum is a monotonically increasing function of z = K∆, independent of the sign of ω j , the sum is also a monotonically increasing function of z. In terms of z, Eq. (3.4) can be rewritten as,

When K is sufficiently small, only the trivial solution z = 0 is obtained. The nontrivial solution is obtained only for K > K c where K c depends on the specific choice of the frequencies {ω j }.

We now consider the statistical properties of the solution to Eq. (A3) when the frequencies {ω j } are drawn independently from a given distribution g(ω). In this case, ψ j ≡ ψ(ω j /z) are also independent random variables whose distribution can be obtained from g(ω). The quantity Ψ as defined by Eq. (A2) is simply the average of the random variables ψ j , j = 1, . . . , N . Its mean value is given by [ψ]. Let δΨ = Ψ -[ψ]. Its variance is given by

The variance of ψ can be easily computed from g(ω):

Invoking the central limit theorem, we conclude that the distribution of Ψ tends to a Gaussian at large N , with mean and variance as stated above.

For small z,

Equation (A3) can now be written in the form

where K -1 c = πg(0)/2, c = -πg (0)/16 > 0, and d = (4g(0)/3)

1 2 is a Gaussian random variable with zero mean and unit variance.

Solution to Eq. (A7) can be cast in the scaling form (3.1) with β = 1 2 and ν = 5 2 . Specifically, let

where x = N 2 5 (K -K c )/K is the scaled distance to the transition point.

We now consider solution to Eq. (A8) at K = K c or x = 0. Simple algebra yields

Here A = (8/3π)

. The nth moment of f (0) at the transition is given by (n > 0),

In the special case g(ω) = e -ω 2 /2 / √ 2π, A = (8π/3)


Denis Kuzzay
email: denis.kuzzay@cea.fr

Davide Faranda

Bérengère Dubrulle

Global vs local energy dissipation: the energy cycle of the turbulent von Kármán flow

In this paper, we investigate the relations between global and local energy transfers in a turbulent von Kármán flow. The goal is to understand how and where energy is dissipated in such a flow and to reconstruct the energy cycle in an experimental device where local as well as global quantities can be measured. In order to do so, we use PIV measurements and we model the Reynolds stress tensor to take subgrid scales into account. This procedure involves a free parameter that is calibrated using angular momentum balance. We then estimate the local and global mean injected and dissipated power for several types of impellers, for various Reynolds numbers and for various flow topologies. These PIV-estimates are then compared with direct injected power estimates provided by torque measurements at the impellers. The agreement between PIV-estimates and direct measurements depends on the flow topology. In symmetric situations, we are able to capture up to 90% of the actual global energy dissipation rate. However, our results become increasingly inaccurate as the shear layer responsible for most of the dissipation approaches one of the impellers, and cannot be resolved by our PIV set-up. Finally, we show that a very good agreement between PIV-estimates and direct measurements is obtained using a new method based on the work of Duchon and Robert [J. Duchon and R. Robert, Nonlinearity, 13, 249 (2000)] which generalizes the Kármán-Howarth equation to nonisotropic, nonhomogeneous flows. This method provides parameter-free estimates of the energy dissipation rate as long as the smallest resolved scale lies in the inertial range. These results are used to evidence a well-defined stationary energy cycle within the flow in which most of the energy is injected at the top and bottom impellers, and dissipated within the shear layer. The influence of the mean flow geometry and the Reynolds number on this energy cycle is studied for a wide range of parameters.

I. INTRODUCTION

Understanding how and where energy is dissipated in turbulent flows has been a great challenge for many years, and would have important implications in many areas such as fundamental research, aeronautics or industry. In the classical three-dimensional turbulence phenomenology, energy is injected at large scales by the forcing mechanism, transferred downscale at a constant rate following a self-similar cascade, and then dissipated into heat at the Kolmogorov length scale, where viscous effects become dominant. In ideal stationary, homogeneous and isotropic turbulence, the measurement of energy dissipation can therefore be achieved via 3 independent and equivalent means: i) by monitoring the injected energy; ii) by monitoring the dissipated heat; iii) by monitoring the cascade energy rate via multi-scale single points measurements of the velocity (via e.g. anemometers or array of hot wires). The first two measurements are global, the last one local, but given the homogeneity in space and time, they all provide the same information. In most realistic situations, however, the turbulence is anisotropic and/or inhomogeneous and/or non-stationary. In such cases, there is not any necessary equivalence between global energy injection, global energy dissipation and local energy dissipation. The study of these three quantities requires detailed knowledge of the forcing, the heat distribution and velocity over the whole domain.

In that respect, it is interesting to focus on intermediate situations, where the turbulence is generated by a well controlled mechanism, in a simplified geometry, as achieved for example in classical laboratory experiments such as Taylor-Couette or von Kármán set-ups. In such cases, it is easy to monitor the forcing and implement a cooling mechanism so as to achieve stationarity, where global energy injection and dissipation equilibrate on average. The local energy dissipation can then be computed from stationary energy budgets derived from Navier-Stokes equation, using local measurements of velocity obtained for example using the now classical Particle Image Velocimetry (PIV). This technique provides measurements of the instantaneous velocity field u i at several points of a plane (or of a volume) at the same time. From this, one may compute the Reynolds stress tensor S ij = 1 2 (∂ i u j + ∂ j u i ), and study the local dissipated power v = 2νS ij S ij . A detailed comparison between this local estimate and the global energy dissipation rate in a Taylor-Couette flow at various Reynolds number has recently been made by Tokgoz et al. [START_REF] Tokgoz | Spatial resolution and dissipation rate estimation in Taylor-Couette flow for tomographic PIV[END_REF] using tomographic PIV measurements. They observe that the local and global estimates coincide within 10% as long as the laminar flow and Taylor vortex flow regimes are fully resolved. However, as the Reynolds number increases, the dissipative scale decreases and becomes smaller than the finite resolution of the PIV (set by the camera resolution and the velocity reconstruction algorithm). In that case, the estimate of local energy dissipation based on the velocity gradients and on the viscosity v becomes increasingly inaccurate, and underestimate the global energy dissipation. In order to remedy this problem, it has been suggested to use techniques borrowed from Large Eddy Simulation (LES) [START_REF] Sheng | A large eddy piv method for turbulence dissipation rate estimation[END_REF] . This allows us to model the subgrid scales (SGS) in terms of the large scale velocity field resolved by the PIV and allows computation of all terms in the energy budget [START_REF] Saarenrinne | Turbulent kinetic energy dissipation rate estimation from PIV velocity vector fields[END_REF][START_REF] Tanaka | A correction method for measuring turbulence kinetic energy dissipation rate by PIV[END_REF][START_REF] Delafosse | Estimation of the turbulent kinetic energy dissipation rate from 2d-piv measurements in a vessel stirred by an axial mixel ttp impeller[END_REF][START_REF] Baldi | On the Direct Measurement of Turbulence Energy Dissipation in Stirred Vessels with PIV[END_REF] including the terms responsible for the scale-to-scale energy transfer. By the Richardson-Kolmogorov cascade picture, this allows an estimate of the local energy dissipation as long as the scale used in the computation lies in the inertial range.

In this paper, we test these methods in a turbulent flow generated by two contra-rotating impellers (von Kármán flow) for Reynolds numbers ranging from 10 3 up to more than 10 6 . At such Reynolds numbers, the dissipative scale ranges from a few millimeters to a few tens of microns. With fixed velocity of the impellers and cooling, the resulting flow is stationary, highly anisotropic and inhomogeneous [START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF], thereby providing a unique laboratory flow to test local energy dissipation procedure. This closed flow geometry permits direct estimates of the global energy injection by torque monitoring at the two impellers. Using a large scale Helium facility with calorimetric measurements [START_REF] Rousset | Superfluid high REynolds von Kármán experiment[END_REF], we were able to show that in a stationary state, the global injected power and the global dissipated heat coincide within a few percent at large Reynolds number, and for a wide range of operating conditions (including differential rotation of the impellers). Moreover, the simple cylindrical geometry allows for Stereoscopic PIV measurements of the velocity field over a vertical plane spanning the whole experiment, at a resolution of a few millimeters. Since measurements have only been made in a meridional plane, we do not have access to orthoradial derivatives. Hence, it is challenging to have an accurate estimate of the local energy dissipation rate despite this lack of information. To do this we first use a LES method using statistical axisymmetry and angular momentum budget to calibrate the model. Then, we test a generalization of the Kármán-Howarth formula derived by Duchon and Robert [START_REF] Duchon | Inertial energy dissipation for weak solutions of incompressible Euler and Navier-Stokes equations[END_REF] that provides a parameter-free estimate of the local energy dissipation for any nonhomogeneous, nonisotropic flow. Both estimates are then averaged over the whole volume for comparison with the global estimate of the energy dissipation based on torque monitoring. This is done for a wide range of control parameters, varying both the Reynolds number, the mean flow geometry and the flow asymmetry. We then use these measurements to evidence a stationary energy cycle within the flow where energy is injected at the top and bottom impellers, and dissipated within the shear layer. This paper is organized as follows: in Section II, we summarize the theoretical tools needed for the implementation of the LES technique and the Duchon-Robert (DR) formula in our analysis of the energy dissipation rate. In Section III, we review the von Kármán geometry and specialize these formulae to the case of cylindrical geometry. The flow diagnostics based on PIV measurements are derived and summarized. In Section IV, we apply these results to a set of measurements drawn from our data base of von Kármán flow. We first tune our LES model using angular momentum balance. Then, we compare PIV estimates of the global dissipated power with direct torque measurements using both the LES technique and the DR formula. Finally, in section V, we evidence the energy cycle of the von Kármán flow, and we discuss its evolution, as well the evolution of our diagnostics, as a function of the flow topology.

II. THEORETICAL BACKGROUND AND METHODOLOGY

In this section, we use a filtering approach to derive energy and angular momentum balance at a given scale from Navier-Stokes equations, and we give the expression of the different terms that we will use in our analysis of the injected and dissipated power. We consider Navier-Stokes equations

∂ t u i + u j ∂ j u i = -∂ i P + ν∂ j ∂ j u i (1) ∂ j u j = 0, (2)
where we use Einstein summation convention over repeated indices.

A. The filtering approach

Following the procedure in [START_REF] Germano | Turbulence -The filtering approach[END_REF], we define a coarse-grained velocity field at scale as

u i ( x) = d r G ( r)u i ( x + r), (3)
where G is a smooth filtering function, nonnegative, spatially localized and such that d r G( r) = 1. The function G is rescaled with as G ( r) = -3 G( r/ ). Coarse-graining the Navier-Stokes equations gives

∂ t u i + u j ∂ j u i = -∂ j τ ij -∂ i P + ν∂ j ∂ j u i (4)
∂ j u j = 0. (5)
In equation ( 4) we introduced τ ij = (u i u j ) -u i u j which is the stress tensor from the SGS.

In what follows, in order to simplify the formulae and for readability considerations, we will drop the index . Unless specified otherwise u i now denotes the i th component of the coarse-grained velocity field (same thing for P ).

B. The energy balance equation

We now take the scalar product of equation ( 4) with u i , and after a few lines of algebra we get the local energy balance equation

∂ t E + ∂ j (u j E) = -∂ j (u i τ ij ) + S ij τ ij -∂ i (u i P ) + ν(∂ j ∂ j E + ∂ i ∂ j (u i u j )) -2νS ij S ij , (6)
where E = uiui 2 is the large scale kinetic energy per unit mass, and

S ij = 1 2 (∂ i u j + ∂ j u i )
is the large scale strain rate tensor. The overline denotes the statistical average of a quantity. Let us now define a vector J j = u j (E + P ) + u i τ ijν(∂ j E + ∂ i (u i u j )). We can then rewrite (6) as

∂ t E + ∂ j J j = S ij τ ij -2νS ij S ij . (7)
This local equation is valid in any geometry, for any type of flow, for any filtering with the properties given above.

In section III D, we apply this formula to the specific axisymmetric von Kármán geometry to derive the local energy production and dissipation rate per unit mass.

C. The angular momentum balance equation

In a very similar fashion, it is also possible to derive an angular momentum balance equation at scale . We take the cross product of r j and equation ( 4), and after a few lines of algebra we get

∂ t L i + ∂ j (u j L i ) = -ijk ∂ k (r j P ) -ijk ∂ m (r j τ km ) + ν(∂ j ∂ j L i + ijk ∂ m ∂ k (r j u m )) -ν ijk ∂ j u k , ( 8
)
where ijk is the three dimensional, total antisymmetric Levi-Civita symbol so that L i = ijk r j u k is the i th component of the angular momentum per unit mass. Of course, it can be checked that it is also possible to derive (8) from ( 6) using ( 4) and the identity E = 1 2 r 2 ( L 2 + ( r. u) 2 ). As before, we can write (8) as an equation of conservation for each component of

L ∂ t L i + ∂ j T ij = 0, (9)
where

T ij = u j L i -ijk r k P + imk r m τ jk -ν(∂ j L i -ijk ∂ m (r k u m ) + ijk u k ).
In section IV A, we use this angular momentum balance as a constraint to calibrate the LES model, that we now describe.

D. LES method for balance equations

The computation of the different terms in [START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF] and ( 9) requires the knowledge of both the velocity field at scale and the Reynolds stress τ ij . In most practical situations, e.g. when the flow is turbulent and the velocity is measured through a PIV system, only the former is available, since we cannot resolve the dissipative scale. A traditional way to overcome this problem is to use a LES-PIV approach [START_REF] Tokgoz | Spatial resolution and dissipation rate estimation in Taylor-Couette flow for tomographic PIV[END_REF][START_REF] Sheng | A large eddy piv method for turbulence dissipation rate estimation[END_REF] to model τ ij in terms of the large scale velocity field [START_REF] Lesieur | Turbulence in Fluids[END_REF][START_REF] Eyink | Multi-scale gradient expansion of the turbulent stress tensor[END_REF]. Several models exist. In the present paper, we choose the gradient model [START_REF] Clark | Evaluation of subgrid-scale models using an accurately simulated turbulent flow[END_REF] where

τ ij = C∆ 2 r ∂ k u i ∂ k u j , ( 10
)
where C is a constant to be calibrated and ∆ r is the width of the filtering. This model ensures forward scatter and backscatter of energy between resolved scales and SGS.

E. Duchon-Robert energy balance equation

An alternative local energy balance equation has been derived by Duchon and Robert [START_REF] Duchon | Inertial energy dissipation for weak solutions of incompressible Euler and Navier-Stokes equations[END_REF] using Leray's weak solution formalism and Onsager's ideas [START_REF] Onsager | Statistical hydrodynamics[END_REF]. The latter amounts to consider a sequence of coarse-grained solutions of the Navier-Stokes equations Eq. ( 4) in the limit → 0 and derive the corresponding energy balance, that reads

∂ t E + ∂ j u j (E + P ) -ν∂ j E = -ν∂ j u i ∂ j u i -D(u i ), (11)
where D(u i ) is expressed in terms of velocity increments δ u( r) = u( x + r) -u( x) as

D( u) ≡ lim →0 D ( u) = lim →0 1 4 V d r ( ∇G )( r) • δ u( r) |δ u( r)| 2 , (12)
where the dependence of δ u and D in x is implied. As the Reynolds number tends to infinity, the scale can be chosen as small as one wants, and the quantity D(u i ) can be seen as the contribution to dissipation coming from a generalized cascade process (possibly linked with the formation of small-scale singularities). Since the result cannot depend on the filtering function G , Duchon and Robert specialized the expression to a radially symmetric filter to get the alternative expression, devoid of any free parameter

D( u) = - 3 16π lim →0 1 | χ|=1 dΣ( χ) |δ u( χ)| 2 δ u( χ) • χ , (13)
where dΣ denotes the area measure on the sphere. As noticed by Duchon and Robert, this expression coincides with the statistical mean rate of inertial energy dissipated per unit mass derived from the anisotropic version of the Kármán-Howarth equation. They therefore argue that the previous formula provides a local non-random form of the Kármán-Howarth equation, valid even for anisotropic, inhomogeneous flows. In the sequel, we apply this formula to our PIV measurements to test whether it can provide a parameter-free estimate of the global dissipation, as well as local instantaneous maps of the local energy dissipation.

III. APPLICATION TO A VON KÁRMÁN GEOMETRY

The goal of the present paper is to compare global estimates of dissipated power with torque measurements, check their coincidence, and study maps of local energy dissipation rate to get some insight of the detailed processes governing energy transfers in a von Kármán flow.

A. Experimental set-up

The von Kármán experiment has been extensively studied over the past years [START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF][START_REF] Marié | Transport de moment cinétique et de champ magnétique par un écoulement tourbillonnaire turbulent: influence de la rotation[END_REF][START_REF] Ravelet | Bifurcations globales hydrodynamiques et magnétohydrodynamiques dans un écoulement de von Kármán turbulent[END_REF][START_REF] Monchaux | Mécanique statistique et effet dynamo dans un écoulement de von Kármán turbulent[END_REF][START_REF] Saint | L'écoulement de von Kármán comme paradigme de la physique statistique hors équilibre[END_REF][START_REF] Cortet | Experimental Evidence of a Phase Transition in a Closed Turbulent Flow[END_REF][START_REF] Saint-Michel | A zero-mode mechanism for spontaneous symmetry breaking in a turbulent von Kármán flow[END_REF][START_REF] Saint-Michel | Probing quantum and classical turbulence analogy in von Kármán liquid helium, nitrogen, and water experiments[END_REF][START_REF] Ravelet | Toward an experimental von Kármán dynamo: Numerical studies for an optimized design[END_REF]. We give here a brief review of the main features of the set-up. Our von Kármán flow is generated in a vertical cylinder by two coaxial, contra-rotating impellers providing energy and momentum flux at the upper and lower end of the cylinder. The inner radius of the cylinder is R = 100mm and the distance between the inner face of the impellers is H = 180mm, which gives an aspect ratio of H/R = 1.8. The turbulence properties (anisotropy, fluctuations, dissipation) are influenced by the geometry of the impellers, i.e their non dimensional radius R t , the oriented angle α between the blades (see Fig. 1), the number n of blades and their heights h b [START_REF] Ravelet | Bifurcations globales hydrodynamiques et magnétohydrodynamiques dans un écoulement de von Kármán turbulent[END_REF]. In the present paper, we consider only impellers with h b /R = 0.2 and α = ±72 • . Those impellers are the so-called "TM60" (with 16 blades) and "TP87" impellers (with 8 blades), the characteristics of which are summarized in Table I. They are essentially similar, except for their number of blades and the material they are made of. A single impeller can be used to propel the fluid in two opposite directions, respectively associated to the concave or convex face of the blades going forward. This can be taken into account by a change of sign of the parameter α.

In the sequel, we denote (-) (resp. (+) ) an impeller used with the concave (resp. convex) face of its blades going forward. The impellers are driven by two independent motors which can rotate at frequencies up to 10Hz. The motor frequencies can be either set equal to get exact counter-rotating regime, or set to different values f 1 = f 2 . To change the viscosity, we have used either water or glycerol at different dilution rates.

B. Control parameters

In the sequel, we choose R and Ω -1 = (π(f 1 + f 2 )) -1 as unit of length and time. The von Kármán experiment is then characterized by two control parameters:

• the Reynolds number:

Re = π(f 1 + f 2 )R 2 ν -1 ,
where ν is the fluid kinematic viscosity, ranges from 10 2 to more than 10 6 so that we can span a full range of regimes, from the purely laminar to the fully turbulent one. • the rotation number:

θ = f 1 -f 2 f 1 + f 2 ,
measures the relative influence of global rotation over a typical turbulent shear frequency. Indeed, the exact counter-rotating regime corresponds to θ = 0. For a nonzero rotation number, our experimental system is similar, within lateral boundaries, to an exact counter-rotating experiment at frequency f = (f 1 + f 2 )/2, with an overall global rotation at frequency (f 1 -f 2 )/2 [START_REF] Marié | Transport de moment cinétique et de champ magnétique par un écoulement tourbillonnaire turbulent: influence de la rotation[END_REF][START_REF] Ravelet | Toward an experimental von Kármán dynamo: Numerical studies for an optimized design[END_REF]. In our experiments we vary θ from -1 to +1, exploring a regime of relatively weak rotation to shear ratio.

C. Measurements

The set-up allows for both global and local flow diagnostics. Torque (global) measurements at each impeller are performed with SCAIME technology and provide values over the kHz range of C 1 and C 2 , being respectively the torque applied to the bottom and top shafts. Following the procedure described in [START_REF] Marié | Transport de moment cinétique et de champ magnétique par un écoulement tourbillonnaire turbulent: influence de la rotation[END_REF], they are calibrated using measurements at different mean frequencies, so as to remove spurious contributions from genuine offsets or mechanical frictions. From this, we compute the nondimensional value K p1 and K p2 of the torque as:

K pi = C i /(ρR 5 Ω 2 )
, where ρ is the density of the working fluid.

Local measurements of the velocity field of the flow have been made by PIV techniques in the stationary regime. The typical size of the particles used is a few tens of micrometers and their density is 1.4. Two cameras take between 600 and a few thousand successive pictures of the flow at a 15Hz frequency. The resolution of our camera frame is 1600x1200 pixels, and the reconstruction is done using peak correlation performed over overlapping windows of size 16 to 32 pixels. As a result, we get measurements of velocity field on a grid of approximate size 60 2 in a vertical plane containing the axis of symmetry (Oz), in a cylindrical system of coordinate. The maximum spatial resolution we can reach for the velocity field with this set-up is therefore of the order of 200/60 ≈ 3mm, about 10 to 100 times larger than the dissipative scale. For more details about the experimental set-up or the measurement techniques see for instance [START_REF] Saint | L'écoulement de von Kármán comme paradigme de la physique statistique hors équilibre[END_REF]. We can therefore estimate the derivatives of the velocity field only along the r and z directions, but we do not have access to derivatives along θ. As a consequence, we will either set these derivatives to zero in the sequel, or take them into account using hypothesis of statistical axisymmetry. We have checked that the two procedure give essentially the same result. In any case, we use incompressibility to estimate ∂ θ u θ .

D. Diagnostics

Flow geometry

We are going to study different types of von Kármán flows (see Fig. 2). These types of flow happen depending on the forcing condition (+) or (-), and whether the system undergoes a spontaneous phase transition (bifurcation) or not [START_REF] Cortet | Experimental Evidence of a Phase Transition in a Closed Turbulent Flow[END_REF]. There is, then, four types of flows: one corresponding to the (+) forcing condition where a phase transition cannot be observed (Fig. 2(a)), one corresponding to the (-) forcing condition where a phase transition is not observed (Fig. 2(b)), and two more corresponding to the two states of the flow that can be observed once the flow has undergone its phase transition (Fig. 2(c)-(d)). The difference between these flows can be characterized through our PIV measurements by their mean velocity profile in the vertical plane of measurements (see Fig. 2). In the first two types, the mean flow is symmetric with respect to the equatorial plane z = 0 and there is a strong shear layer in the middle. In the bifurcated states, the flow is no longer symmetric with respect to the equatorial plane. It consists in a one-cell flow in the vertical direction, with a strong shear layer at the impeller that rotates in the direction opposite to the orthoradial mean flow. The two bifurcated states are symmetric to each other. The interest of considering these different mean flow geometries is that they are characterized by well resolved (resp. badly resolved) shear layer for the symmetric state (resp. bifurcated state). Since we expect an important fraction of energy dissipation to be localized where there is a strong velocity gradient, this difference may be a large source of error in local estimates based on PIV measurements, as we demonstrate later.

Global dissipation

Using torque measurements, we can get an accurate estimate of the global power injected into the flow. Indeed, since we study a stationary situation, this input must balance the rate of energy dissipation within the flow. This has been checked in a scale 4:1 version of our experiment (see [START_REF] Rousset | Superfluid high REynolds von Kármán experiment[END_REF]). From these measurements, we get the global mean dissipation rate as

D = 2π(C 1 f 1 + C 2 f 2 ).
From this, we compute diagnostics:

• the dimensionless mean dissipation rate:

D = D ρR 5 Ω 3 = K p1 (1 + θ) + K p2 (1 -θ), (14)
• the dimensionless mean dissipated power per unit mass:

= D R 3 V = DR πH , (15)
which should not be confused with the Levi-Civita symbol. V = πR 2 H is the volume of the experiment.

These quantities depend on the Reynolds number, the rotation number, the characteristics of the impellers and the mean flow geometry [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von K?rm?n closed flow[END_REF]. Examples of variation of as a function of Re are provided in Fig. 3(a) at θ = 0, with TP87 and TM60 impellers, for the different flow geometries illustrated in Fig. 2. At low Reynolds numbers, the dissipation rate decays as Re -1 , until Re ≈ 300 where the turbulence sets in. The dissipation rate then levels off at a value which depends upon the flow geometry: it is the lowest for the symmetric (+) flow, then increases for the (-) symmetric flow and is the largest for the two bifurcated (-) flows. In Fig. 3(b), the mean dissipation rate is shown as a function of θ for TP87 (±) impellers at Re = 3 × 10 5 . In the (+) case (pink disks), it may be seen that the minimum of is obtained at θ = 0. When θ is varied from, say, 0 to 0.5, increases and reaches a value twice as large as its minimum. In the (-) case, there is a discontinuity of the energy dissipation at θ = 0 due to the global symmetry breaking. We see that the symmetric state (blue rhombi) dissipates less energy than the two bifurcated states (red triangles). Another striking feature of Fig. 3(b), is the coexistence of two branches of dissipation, corresponding to the two bifurcated states, for a certain range of values of θ. This coexistence vanishes around θ = ±0.2 (this value depends on the Reynolds number, see [START_REF] Saint-Michel | Probing quantum and classical turbulence analogy in von Kármán liquid helium, nitrogen, and water experiments[END_REF]).

Local production and dissipation rates

Since the von Kármán flow is neither isotropic, nor homogeneous, it is interesting to study more locally the production and dissipation of energy. Our set-up being statistically axisymmetric, it is sufficient to perform the study in the vertical plane including the rotation axis, where we perform our velocity measurements. We thus consider the main stationary terms of the local energy balance Eq. ( 7), and divide them by the experimental volume V = πR 2 H to get three local quantities:

• the local production rate per unit mass:

Γ = - 1 V ∂ j J j . ( 16
)
Integrating over the whole volume, one may then get from this the corresponding total production rate per unit mass Γ T :

Γ T = - 1 V V ∂ i J i dV = - 1 V S J i dS i , (17)
where the last equality comes from Green-Ostrogradsky formula.

• the positive local dissipation (transfer) rate per unit mass:

Π = 1 V -S ij τ ij + 2νS ij S ij . ( 18
)
The first term represents the energy dissipation due to energy transfers towards SGS while the second one is the laminar energy dissipation. • the positive local singularity dissipation (transfer) rate per unit mass:

Π DR = 1 V D( u). (19)
As before, we denote by Π T and Π DRT the total dissipation rate per unit mass, obtained by volume integration of Π and Π DR .

With our stereoscopic PIV measurements, we have access to all terms to compute these quantities, except for the pressure P in J i , and the terms involving derivatives with respect to the azimuthal angle. In the sequel, we will present maps of the local quantities without the terms we do not have access to. In other words, we neglect all terms involving pressure and azimuthal derivatives. Moreover, the computation of the transfer terms requires the calibration of the constant C involved in the gradient model (see Eq. ( 10)). In section IV A, we describe a calibration procedure using the angular momentum balance equation.

IV. DATA PROCESSING

A. Calibration of the gradient model using angular momentum balance

The first step is then to find the value of the C constant in order to be able to estimate the contribution of the terms containing τ ij in ( 16) and [START_REF] Saint | L'écoulement de von Kármán comme paradigme de la physique statistique hors équilibre[END_REF]. For this, we follow the work of [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF] who has shown, by using high-resolution Laser Doppler Velocimetry (LDV) measurements, that in a symmetric situation (θ = 0), the vertical flux of angular momentum is a constant equal to the torque injected by the impellers. This result stems from the z-component of the angular momentum balance Eq. ( 9) integrated over a volume V (z) that describes a cylinder extending from the bottom impeller, to an altitude h(z) = H/2 + z (see Fig. 4). In this case, the global z-component of the angular momentum balance equation ( 8) reads

K p1 = -Re -1 S(z) (u θ + r∂ r u θ )dS - Σ(z) r∂ z u θ dΣ + S(z) r(u r u θ + τ rθ )dS + Σ(z) r(u θ u z + τ θz )dΣ ≡ -K v (z) + Σ(z) r(u θ u z + τ θz )dΣ, (20)
where Σ(z) is a surface at altitude h(z) from the bottom impeller, and S(z) is the lateral boundary (see Fig. 4). As discussed in [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF], this equation states that the angular momentum transmitted by the motor to the fluid (measured by the dimensionless torque at the bottom K p1 ) is either evacuated through the lateral boundary, or transported to the upper layers of the fluid to be received by the upper motor as a drag. It was further shown in [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF] that in a symmetric situation (θ = 0), where the torques at the bottom and at the top are equal, the K v (z) term is negligible, so that there is a constant flux of angular momentum from the bottom to the top. For each θ = 0 case, we thus compute the quantity

Φ Lz (z) = Σz r(u θ u z + τ θz )dS, (21)
and adjust the constant C, so that Φ Lz (z) = K p1 = K p2 ≡ K p for any z. In Fig. 4, we present the results of our computations for TM60(+) impellers. The adjustment of the constant C is done statistically, by assuming that it does not depend on the Reynolds number (as it is supposed to depend only on the velocity reconstruction algorithm).

We then fix ∆ 2 r = ∆ x ∆ z , where ∆ x = 2R/58 and ∆ z = (H -2h b )/63 being the spatial resolutions of our PIV and we adjust C using 38 symmetrical flows at different Reynolds numbers so that their statistical mean provides a constant value of Φ(z)/K p = 1 within the (statistical) standard deviation. The optimal value of C has to be taken equal to C ≈ 4. This is around 50 times more than the conventional choice made in [START_REF] Eyink | Multi-scale gradient expansion of the turbulent stress tensor[END_REF][START_REF] Leonard | Energy cascade in large-eddy simulations of turbulent fluid flows[END_REF] for a Gaussian filter. We have identified several factors that may explain this difference: our data are not filtered in a Gaussian way; since we ignored azimuthal derivatives, the constant must be higher to compensate; the turbulence is neither isotropic nor homogenous.

Dividing the flux into the contribution due to the time-average flow and the contribution due to the fluctuations, we also recover that near the impellers the flux is only due to the mean flow while at the center it is almost only due to fluctuations [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF]. We also checked that the diffusive terms in [START_REF] Saint-Michel | A zero-mode mechanism for spontaneous symmetry breaking in a turbulent von Kármán flow[END_REF], going as Re -1 , are negligible at high Reynolds number. Likewise, the flux through S(z) is also small, meaning that almost all of the flux goes in the vertical direction.

B. Implementation of the Duchon-Robert formula

As we have seen in section II E, D ( u) provides the local dissipation for a given velocity field u, at a given length scale . Since this formula is a generalization of the Kármán-Howarth equation to any kind of flow, it provides an estimate of the subgrid energy transfer as long as the considered scale is in the inertial range. Moreover, the considered scale must be sufficiently large with respect to the PIV smallest resolved scale, so as to guarantee statistical convergence through sufficient average in the scale space (on the sphere of radius ). To check these two points, we present on Fig. 5 D ( u) in two plots, averaged in the radial (resp. vertical) direction, as a function of z (resp. x) and , for an experiment using TM60(+) impellers, at very large Reynolds number (Re ≈ 8 × 10 5 ). We see that this quantity is close to zero at large scales ( > 0.4), but that it increases at small scales in the domain |x| < 0.4, |z| < 0.7, i.e. at the location of the median shear layer. For between 0.1 and 0.15 (i.e approximately 4 to 5 times the smallest scale resolved by our PIV set-up), there is the start of a saturation, indicating the beginning of the inertial range. While the extent of the inertial range is likely to vary (increase) with the Reynolds number, its largest scale is likely to be independent of the Reynolds number, as long as the flow is turbulent. Indeed, as discussed in [START_REF] Thalabard | A statistical mechanics framework for the large-scale structure of turbulent von Kármán flows[END_REF], the geometry of the largest scales of the von Kármán flow appear fairly independent of the Reynolds number, except around Re = 10 5 where they may experiment abrupt changes due to the equivalent of a phase transition [START_REF] Cortet | Experimental Evidence of a Phase Transition in a Closed Turbulent Flow[END_REF]. Since in the Kolomogorov picture energy cascades from large to small scales, it is reasonable to assume that the beginning of the inertial range is solely determined by the large scale topology, thereby becoming independent of the Reynolds number (except maybe around Re = 10 5 ). To check this, we report on Fig. 6(a) the comparison between the total dissipated power Π DRT at = 0.1 (blue symbols) and direct measurements of the injected power (black symbols). We observe that our estimates are in good agreement with direct measurements, especially for symmetric flows. In Fig. 6(b), the results of our computations to estimate the global dissipated power using the LES method are displayed. These results will be discussed in more detail later, but we can already observe that we obtain a good agreement with direct measurements. We can conclude from all our computations that the method using the DR formula may be seen as an interesting alternative to the widespread LES-PIV method, since it relies on very few arbitrary hypotheses. We explore its performances with other flow geometries in the sequel. FIG. 6: Comparison between direct measurements of energy injection obtained using torque measurements (black symbols) and PIV estimates at various Reynolds numbers, for TM60 and TP87 impellers for the four different mean state geometries displayed in Fig. 2: disks: (+) symmetric. Rhombi: (-) symmetric. Up triangles: (-) shear layer sent downwards. Down triangle: (-) shear layer sent upwards. a) energy dissipation ΠDRT using the DR method (blue symbol). b) energy injection ΓT (red symbols) and dissipation ΠT (blue symbols) using the LES-PIV method. The dotted line represents the "laminar fit" = 37Re -1 . The estimates are computed based on 2 to 15 realizations of the experiment where at least 600 instantaneous velocity snapshots have been taken for each of them. The symbols represent the mean of our computations while the error bars represent the standard deviation.

V. RESULTS

A. Symmetric case θ = 0 a. Energy production and transport To understand the flux of energy in the von Kármán flow, it is interesting to focus first on the symmetric case (θ = 0), specifically on the case when both impellers rotate in the (+) sense with the same frequency, in a stationary regime. In such regime, the dissipated power equals the injected power. Through Green-Ostrogradski theorem, the total energy production Γ T is equal to the entering flux of J i at the boundaries. In this symmetric case, most of the flux is provided by the component J z , and appears to provide a fairly good estimate of the injected power at large enough Reynolds numbers (Re >= 10 5 ). Indeed, we show on Fig. 6(b) the results of our computations for Γ T (empty red symbols). We see that for Re >= 10 5 , these estimates coincide within 20%. However, at lower Reynolds numbers (Re ≈ 4 × 10 3 ), the PIV estimates only capture 55% of the actual injected power. This may be due to the fact that as the Reynolds number is decreased, an increasing part of the injected power is either through pressure effects, azimuthal variations or viscous boundary layers that are not resolved by our measurements. Given that the J z contribution dominates the total energy production, it is interesting to focus on the spatial variation of this quantity. Fig. 7(a As can be seen, at the center of the cylinder, i.e. x/R ∈ [-0.8 0.8], big structures reflect the advection of energy towards the impellers, through the mean Ekman pumping by the impellers. In contrast, at the walls, smaller structures of opposite sign are observed, reflecting the injection of energy within the flow. These two kinds of structures are mirrored by the recirculation cells that we observed on Fig. 2. The local structure of J z can be further used to get information about energy transport in the flow. Indeed, performing an integration of equation ( 7) over the heightvarying volume V (z) (see Fig. 4) and setting ∂ t E = 0, we can get an equation for the energy transport in the flow as

S J i dS i = -K p1 (θ + 1) -K v (z) + Sz J z dS, (22)
where K v is the contribution due to the lateral boundaries at height z. Ignoring this contribution, we can write

S J i dS i = -K p1 + Φ E (z), (23)
where Φ E (z) = Sz J z dS. This quantity is displayed on Fig. 8 . We observe that near the impellers, the total flux of energy (filled black squares) equals the energy injected by the impellers. The change of sign comes from the fact that the impeller at the bottom injects energy in the +z direction 8: Plot of the normalized vertical flux of energy as a function of z. The dark grey area represents convective transport due to the mean flow while the light grey area represents convective transport due to fluctuations. Filled squares represent the sum of the two contributions so that the total convective transport is a decreasing function of z. Each point has been obtained by taking the mean of several computations obtained from several flows, and the error bars represent the statistical standard deviation.

-0.5 0 0.5 -3 -2 -1 0 1 2 3 Φ E (z) K p z FIG.
whereas the upper impeller injects energy in the -z direction. At the center of the experiment the flux is zero, meaning that on average there is not any energy going from one half of the cylinder to the other through the shear layer. Finally, it is interesting to see that, as before, the mean flow plays an important role in the transport of energy near the impellers, whereas near the shear layer energy is carried by fluctuations. A slight difference is that near the impellers, as the mean flow sends energy towards the center of the cylinder, fluctuations create a flux which goes in the opposite direction and tries to send energy back to the impellers. Altogether, the local map of injected power is shown on Fig. 9(a). We see that near the impellers, the divergence term brings energy into the system as is expected. We also see that energy leaves the center of the recirculation cells and tends to be advected towards the center of the impeller through Ekman pumping. b. Energy dissipation We now turn to the total energy dissipation estimated through the LES method Π T . Its values at various Reynolds numbers in the symmetric case are reported in Fig. 6(b). We observe that it is in very good agreement with direct measurements at Re ≈ 4 × 10 3 . At such Reynolds number, the dissipative scale is of the order of 1mm, close to the PIV resolution. This result is therefore in agreement with the observation of Tokgoz et al. [START_REF] Tokgoz | Spatial resolution and dissipation rate estimation in Taylor-Couette flow for tomographic PIV[END_REF] obtained for a Taylor-Couette flow. At larger Reynolds numbers, the estimates using the LES method are in good agreement with respect to direct measurements since we are able to capture up to 90% of the actual energy dissipation. The map of the local dissipation with this method is provided in Fig. 9(b), and appears fairly uniform across the vessel, with larger intensity along the vertical axis around r = 0 and near the lateral boundaries around z = ±0.3. This map can be compared with the one obtained using the DR formula D ( u) at = 0.1 Fig. 9(d). One observes the same localized structures of energy dissipation near the walls at x = ±0.6, which are symmetrically distributed with respect to the mid-plane. These structures are statistically significant, since they are not observed on plots of the instantaneous local dissipation D ( u) (see Fig. 9(c)). They may, therefore, trace the intense vortices of the shear layer. In addition, one observes a clear localization of energy injection (red areas) at the tip of the impellers, with an energy dissipation in the middle part of the cells. Regarding instantaneous maps, it is interesting to note that it is characterized by intense, localized events, that can reach 10 to 20 times the maximum average local dissipation. Whether these intense events are connected to localized quasi singular structures is an interesting open question that we leave for future work.

Overall, these dissipation maps provide clear evidence that the maximum energy dissipation lies within the shear layer and that the DR formula provides a better estimate of the energy dissipation than the LES method in the symmetric case. c. Description of the energy cycle Altogether, our results regarding energy production, transport and dissipation can be summarized into a simple picture of the "energy cycle" for the von Kármán flow, which is sketched in Fig. 10: the energy is advected to the impellers via the Ekman pumping, the flow is then accelerated inside the impellers and expelled at the walls, providing an energy injection towards the mixing layer. In that mixing layer, turbulent fluctuations dissipate an energy equal to . The fluid is then pumped again into the impellers for further reinjection, closing the energy cycle and providing a stationary situation with energy injection and dissipation equal to .

In the sequel, we study the influence the flow topology and the forcing conditions onto this energy cycle. Since in this picture most of the energy is dissipated within the middle shear layer, there is good hope that we can capture its main contribution by the PIV measurements, provided that the shear layer is not too close to the impeller since in that case, the PIV measurements cannot resolve the flow. We check in section V B-C that this is indeed the case, and provide detailed informations about the localization of energy dissipation to complete the picture of the energy cycle.

FIG. 10: Energy cycle in the von Kármán experiment: the energy is advected to the impellers via the Ekman pumping, the flow is then accelerated inside the impellers and expelled at the walls, providing an energy injection towards the mixing layer.

In that mixing layer, turbulent fluctuations dissipate an energy equal to . The fluid is then pumped again into the impellers for further reinjection, closing the energy cycle and providing a stationary situation with energy dissipation .

B. Influence of the flow topology

Let us now consider the case when both impellers rotate in the (-) sense at θ = 0. For this type of forcing, at sufficiently high Reynolds numbers, there is coexistence of three different flow geometries for TP87 type impellers.

Considering the case where the flow has not undergone any phase transition, Fig. 6 shows that for this type of forcing we measure a dissipation almost three times bigger than what was observed with the previous forcing condition. In this symmetric case, our estimates of the injected and dissipated power are within 20% of the measured value using the LES-PIV method, whereas we reach 98% of the actual dissipation rate of energy with the DR formula (see Fig. Here, the remarks are essentially the same as in Fig. 9: near the impellers, the divergence term brings energy into the system while energy leaves the center of the recirculation cells to accumulate at the center of the cylinder. The dissipation term does not change much either and is approximately constant throughout the plane of measurements. A noticeable difference however, is that the color scale is bigger for this configuration. We thus recover the fact that the flow dissipates more energy than for a forcing in the (+) sense and, as a consequence, the torque of the impellers has to be higher. The situation changes for the bifurcated (-) states. The shear layer is now very close to the upper or lower impeller, depending on the state "chosen" by the system. The mean velocity field for these two states is essentially the same, differing only by the transformation z → -z. The measured injected power is also the same in these two states. However, in terms of PIV-estimated injected or dissipated power, a clear asymmetry occurs. Indeed, when the shear layer is sent downwards (up-pointing triangles), the LES-PIV method provides good estimates of both the injected and the dissipated power (see Fig. 6(b)), while the DR-method provides about 65% of the dissipated power. In the case where the shear layer is sent upwards (down pointing triangles), however, both the LES-PIV method and the DR method totally fail to reproduce the measured dissipated energy (giving the wrong sign). Since there is not any asymmetry observed in the mean flow, this difference must be attributed only to fluctuations and tiny asymmetries of the experimental set-up (laser sheet location, focalization of cameras, ..). This points out the the importance of resolving the shear layer in the PIV-estimates of the injection or dissipation energy rate. We thus conclude that neither the LES-PIV, nor the DR method are appropriate in cases when the shear layer is at the location of the impellers, where it cannot be resolved by our PIV set-up.

Local maps of injected and dissipated power in the bifurcated states are plotted on Fig. 12. They correspond to the flows (c) and (d) displayed on Fig. 2 where the forcing is in the (-) sense. Maps (a) (c) and (e) correspond to a flow where the shear layer has been sent near the lower impeller. As a consequence, we observe that energy enters from the top and is advected downwards to the shear layer. The map of the dissipation term for this flow looks quite like the production term. We see that most of the dissipation we are able to capture happens at the center of the cylinder but, as we said before, we miss all the dissipation that happens near the lower impeller where there is the shear layer. At first sight, maps (b) (d) and (f) are just symmetric of maps (a) (c) and (e) with respect to z → -z. However, by looking closely, tiny differences are observed, especially near the walls and impellers. These differences explain the asymmetry between the two bifurcated states shown on Fig. 6.

C. Influence of the forcing asymmetry θ

We now investigate the case where the impellers are not exactly counter-rotating, so that the parameter θ varies from -0.5 to 0.5. In the (-) case, the geometry observed at |θ| > 0 is always bifurcated (see Fig. 3(b)), with the shear layer located at the top or bottom impeller. In the (+) case, the transition is more gradual: the shear layer is increasingly shifted upwards (resp. downwards) as θ goes from 0 to 1 (resp. -1), allowing for finer tests of the accuracy of our PIV-estimates as a function of the flow geometry. On Fig. 13 are plotted our estimates for the total injected (red symbols) and dissipated power (blue symbols) using both the LES-PIV method (Fig. 13(b)) and the DR formula (Fig. 13(a)). For each θ, our computations are based on sets of at least 600 instantaneous velocity snapshots. Because of the symmetry z → -z, we expect all estimates to be symmetric with respect to θ → -θ, provided the statistics are well converged and that the shear layer is sufficiently resolved. We see that this is indeed the case for (+) forcing condition, where the shear layer always lies in between the two impellers. However, we observe that both methods become inaccurate when |θ| is too high i.e |θ| > 0.3.

For the (-) forcing condition, the estimates give very good agreement for the symmetric state. However, when the phase transition occurs, one observes the same asymmetry as in the θ = 0 case between the two possible states: when the shear layer is sent downwards we get a good agreement between measurements and PIV estimates with the LES method. However, the DR method systematically underestimates the dissipation. In the case where with the shear layer is sent upwards, the estimates are really bad with both methods.

Corresponding maps of the injected and dissipated power are provided in Fig. 14 at θ = -0.1 in the (+) and (-) sense. They obey qualitatively the same behavior as in the symmetric case, so that the energy cycle description is qualitatively the same, with the shear layer location being moved as θ varies. Right: (-). Top line represents the estimated injected energy, using LES-PIV method. Middle line represents the estimated dissipated energy using the LES-PIV method. Bottom line represents the estimated dissipated energy using the DR method.

Areas where energy accumulates are represented in red, and those where energy leaves are represented in blue.

VI. CONCLUSION

In this paper, we investigate the energy cycle of a turbulent von Kármán flow. In this kind of set-up, it is possible to control the mean energy flux inside the flow in order to get a statistically stationary regime. It is also possible to monitor the torque and frequency of the impellers. In our case, we use Particle Image Velocimetry (PIV) measurements to study local and global energy transfers inside the flow. This kind of measurements gives us access to the effective velocity field on a grid that does not resolve the dissipative scale. Therefore, we use a Large Eddy Simulation (LES) approximation to model the influence of unresolved scales. This procedure involves a free parameter which has to be calibrated for our set-up. This is achieved by imposing angular momentum balance at the smallest resolved scale [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF]. After deriving an energy balance equation at a fixed scale , we proceed to estimate four quantities from our PIV measurements: the local and global mean power injected by the impellers and the local and global mean dissipated power. This computation is performed for various Reynolds numbers and for various flow topologies. These PIV-estimates are then compared with direct injected power estimates provided by torque measurements at the impellers. The agreement between PIV-estimates and direct measurements depends on the flow topology. In symmetric situations, we capture up to 90% of the actual energy dissipation. However, our results become increasingly inaccurate as the shear layer responsible for most of the dissipation approaches one of the impellers, and cannot be resolved by our PIV set-up. At the same time, we explore a new method for PIV-estimates of the energy dissipation, based on the work of Duchon and Robert [START_REF] Duchon | Inertial energy dissipation for weak solutions of incompressible Euler and Navier-Stokes equations[END_REF] that generalizes the Kármán-Howarth equation to nonisotropic, non homogeneous flows. This method provides parameter-free estimates of the energy dissipation as long as the smallest resolved scale lies in the inertial range and the shear layer is resolved by the PIV set-up. With this method, we obtain a very good agreement between PIV estimates and direct measurements, and we are able to capture up to 98% of the actual dissipated power in symmetric situations. However, this method also gives increasingly inaccurate results as the mixing layer approaches one of the impellers. These results are used to evidence a well-defined stationary energy cycle within the flow in which the energy is injected by the top and bottom impellers towards the shear layer. There, turbulent fluctuations dissipate energy, and the flow is then pumped towards the impellers, closing the energy cycle.

FIG. 1 :

 1 FIG. 1: Illustration of a TP87 type impeller, with 8 blades. The convention used to name the two different forcing conditions is represented: (+) when the convex face of the blade goes forward, (-) for the other. The angle α characterizes the curvature of the blades and is equal to α = 72 • in the case of TP87 impellers. TM60 impellers look essentially the same, with 16 blades instead of 8.

5 FIG. 2 :

 52 FIG.2:The four types of flow characterized by their mean velocity profiles. Arrows represent the velocity field in the radial plane while colours represent the orthogonal component to that plane. a) The flow is forced under the (+) condition, no phase transition can be observed. b), c) and d) show the three different states of a flow forced under the (-) condition. b) no phase transition has occurred, c) (resp. d)) the shear layer has been sent downwards (resp. upwards). In each of these flows Re ≈ 3 × 10 5 , and the two TP87 impellers rotate in opposite directions at θ = 0.

FIG. 3 :

 3 FIG.3: Plots of the dimensionless dissipated power per unit mass: a) As a function of Re for TP87 and TM60 impellers for the four geometries depicted in Fig.2at θ = 0. Dotted line: Fit ≈ 37Re -1 . Pink disks: symmetric state (+). Blue rhombi: symmetric state (-). Red triangles: bifurcated states (-). b) As a function of θ for TP87 impellers, Re = 3 × 10 5 . Pink disks: symmetric state (+). Blue rhombi: symmetric state (-). Red triangles: bifurcated states (-).

FIG. 4 :

 4 FIG. 4: Left: Control volume for the angular momentum budget. Right: Plot of the normalized vertical flux of vertical angular momentum as a function of z. The dark grey area represents convective transport due to the mean flow while the light grey area represents convective transport due to fluctuations. Filled squares represent the sum of the two contributions so that the total convective transport is almost constant with z and equals to Kp. Each point has been obtained by taking the mean of several computations obtained from several flows at different Re. The error bars represent the standard deviation.

FIG. 5 :

 5 FIG. 5: Vertical (a) and radial (b) average of the local dissipation D ( u) for TM60(+) at Re = 8 × 10 5 as a function of x, z and .

FIG. 7 :

 7 FIG. 7: Typical maps of J z : a) mean energy transfer through the system, b) contribution of the mean flow, c) contribution of the fluctuations

FIG. 9 :

 9 FIG. 9: Top: Maps of the energy production rate (a) and energy transfer rate (b) for the same flow as in Fig. 2(a), using LES-PIV estimates. Bottom: Instantaneous (c) and time averaged (d) maps of the local dissipation D u) at = 0.1 for the same flow. Areas where energy accumulates are represented in red, while those where energy leaves are represented in blue.

  6(a)). Local maps of injected and dissipated power are plotted on Fig.11. They correspond to the flows (b) displayed on Fig.2where the forcing is in the (-) sense. Here Fig.11, (a) and (b) represent LES-PIV estimates while Fig.11, (c) and (d) represents instantaneous and time averaged maps using the DR method.

FIG. 11 :

 11 FIG. 11: Top: Maps of the energy production rate (a) and energy transfer rate (b) for symmetric (-) geometry corresponding to the flow in Fig. 2(b), using LES-PIV estimates. Bottom: Instantaneous (c) and time averaged (d) maps of the local dissipation D ( u) at = 0.1 for the same flow. Areas where energy accumulates are represented in red, while those where energy leaves are represented in blue.

FIG. 12 :

 12 FIG.12: Maps of the energy production rate and energy transfer rate for the two bifurcated geometries (-): Left: for the flow represented in Fig.2(c). Right: for the flow represented in Fig.2(d). Top line represents the estimated injected energy, using LES-PIV method. Middle line represents the estimated dissipated energy using the LES-PIV method. Bottom line represents the estimated dissipated energy using the DR method. Areas where energy accumulates are represented in red, and those where energy leaves are represented in blue.

FIG. 13 :

 13 FIG.13: Comparison between direct measurements of energy injection obtained through torque measurements (black symbols) and PIV based estimates at various θ, for TP87 impellers for the four different mean state geometries of Fig.2: disks (+). Rhombi: (-) symmetric. Up triangles: (-) shear layer sent downwards. Down triangle: (-) shear layer sent upwards. a) energy dissipation using DR method (blue symbols). b) energy injection ΓT (red symbols) and dissipation ΠT (blue symbols) using the LES-PIV method. The estimates are computed based on 1 to 2 realizations of the same experiment where at least 600 instantaneous velocity snapshots have been taken for each of them. The symbols represent the mean of our computations while the error bars represent the standard deviation.

FIG. 14 :

 14 FIG.14:Maps of the energy production rate and energy transfer rate at θ = -0.1 for the two rotation directions: Left: (+). Right: (-). Top line represents the estimated injected energy, using LES-PIV method. Middle line represents the estimated dissipated energy using the LES-PIV method. Bottom line represents the estimated dissipated energy using the DR method. Areas where energy accumulates are represented in red, and those where energy leaves are represented in blue.

TABLE I :

 I Parameter space explored in this paper.