L Cappanera

P Debue

H Faller

D Kuzzay

E-W Saw

C Nore

J.-L Guermond

F Daviaud

C Wiertel-Gasquet

B Dubrulle

Turbulence in realistic geometries with moving boundaries: when simulations meet experiments

Keywords: Penalty method, Large Eddy Simulation model, Turbulent flow, Bifurcation, Transition to turbulence

Considering the current advances in experimental capabilities in fluid mechanics and the advances in computing power and numerical methods in computational fluid mechanics, a question that naturally arises is whether the two sets of techniques are approaching a level of sophistication sufficiently high to deliver results on turbulent flows in realistic geometries that are comparable. The purpose of this paper is to give elements of answers to this question by considering the so-called von Kármán flow where the fluid in a cylindrical container is driven by two counter-rotating impellers. We compare in the mentioned flow the torque and the flow topology obtained by experiments, direct numerical simulations (DNS), and large eddy simulations (LES) at various Reynolds numbers ranging from R e = O(10 2 ) to R e = O(10 5 ). In addition to validating the proposed LES model, the level of agreement that is observed between the numerical and the experimental data shows that the degree of accuracy of each of these techniques is reaching a threshold beyond which it is possible to use each of them with high confidence to explore and better understand turbulence in complex flows at R e = O(10 5 ) and beyond.

Introduction

The equations modeling the motion of incompressible homogeneous fluids are well-established. Denoting by u the three-dimensional velocity field, p the pressure, ρ the (constant) density, f the forcing, and ν the kinematic viscosity, the incompressible Navier-Stokes equations take the following form:

∂ t u + u•∇u = - 1 ρ ∇p + ν∆u + f , (1)
∇•u = 0. ( 2
)
It is known since the pioneering work of Reynolds that after proper rescaling of the equations by L and U , some characteristic length and velocity scales, the Navier-Stokes equations only depend on one parameter: the Reynolds number R e := LU/ν. The above system of balance equations is simple to formulate, but this apparent simplicity is deceiving since the Navier-Stokes system is the source of the notoriously hard problem commonly referred to as turbulence [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF].

Over the years, two competing strategies have been used to advance the understanding of turbulence: numerical simulations and laboratory experiments. The constant progresses made in computer technology, computing techniques, and imaging softwares have pushed the limits of applicability of both numerical simulations and laboratory experiments. For instance, the availability of larger parallel computers and progresses made in parallel linear algebra libraries now allow Direct Numerical Simulation (DNS) and turbulence models, like Large Eddy Simulation (LES), to handle larger Reynolds numbers and to perform longer time statistics. Recent algorithmic progresses made on penalty methods now allow to use Fourier techniques to simulate complex geometries, thereby unlocking the power of highly performant spectral methods. Similarly, technological progresses over the years have steadily increased experimental capabilities. For instance, the use of CCD cameras combined with the continuous increase of the performance of computers have made digital particle image velocimetry a very accurate tool, which is now widely available. Particle image velocimetry (PIV) gives access to all the velocity components in laboratory experiments.

A natural question that comes to mind is whether we are at a point where numerical techniques and laboratory experiments have reached a common ground where both techniques can be used at the same time to investigate the same complex flow, at the same values of the control parameter, and be either compared to each other or used to complement each other in order to advance the theory of turbulence. One goal of this paper is to give elements of answer to this question. The other two goals of the paper are to validate the performance of two numerical techniques: (i) an entropy viscosity LES technique and (ii) a pseudo-penalty technique to simulate moving boundaries.

Comparisons between LES (and/or DNS) simulations and experimental results have already been conducted on similar setups [START_REF] Zadghaffari | Large-eddy simulation of turbulent flow in a stirred tank driven by a Rushton turbine[END_REF][START_REF] Posa | Large-eddy simulations in mixed-flow pumps using an immersed-boundary method[END_REF][START_REF] Benard | Large-Eddy Simulation of wind turbines wakes including geometrical effects[END_REF], but the experimental data used therein were obtained from the literature. Instead here we produce both experimental and numerical results in the very same configuration.

In the present paper we focus our attention on the socalled von Kármán flow which consists of a fluid in a cylindrical container driven by a pair of counter-rotating impellers. The turbulence generated in this setup is neither isotropic, nor homogeneous, and is therefore an example of turbulence in a complex geometry that is rarely discussed in the turbulence literature. To account for the moving impellers driving the fluid in the numerical simulations, we adapt the pseudo-penalty method by Pasquetti et al. [START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF] to a mixed finite element/Fourier approximation setting. We also adopt an entropy-viscosity-based LES technique to handle large Reynolds numbers. We demonstrate in this paper that in this setting the combination of the pseudo-penalty method and the entropy-viscositybased LES model allows us to reach levels of accuracy that are comparable to laboratory experiments. By carefully comparing global and local indicators of turbulence at similar values of R e , we show that the numerical simulations and the laboratory experiments agree with each other for Reynolds numbers ranging from R e = O(10 2 ) up to R e = O(10 5 ). Moreover, the numerical simulations reproduce different bifurcated states of the turbulent flows (i.e., turbulent states that break symmetries and coexist at the same high Reynolds number) that are experimentally observed.

The paper is organized as follows. In §2 we describe the von Kármán flow together with the relevant control parameters. The numerical methods used to approximate the Navier-Stokes equations with moving solid obstacles are described in §3. We also briefly describe in this section the entropy viscosity stabilization model that allows us to approximate large Reynolds number flows with coarse meshes. Section §4 introduces the experimental setup and additional parameters that characterize von Kármán flows. Comparisons between the experimental results and the numerical simulations are done in §5. One highlight of this section is Figure 15 where we compare experimental measurements of the global energy dissipation in the considered setup with numerical simulations in the range R e ∈ [10, 10 5 ]. Concluding remarks are reported in §7.

Setup description

The von Kármán (VK) flow is produced by the stirring of a fluid contained in a cylindrical tank. This setup has been extensively studied using various fluids such as glycerol, water, air, helium gas, superfluid helium, etc. [START_REF] Labbé | Study of the von Kármán flow between coaxial corotating disks[END_REF][START_REF] Ravelet | Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation[END_REF][START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF][START_REF] Monchaux | Statistical mechanics and dynamo action in a turbulent von Kármán flow[END_REF][START_REF] Saint-Michel | Influence of Reynolds number and forcing type in a turbulent von Kármán flow[END_REF][START_REF] Rousset | Superfluid high reynolds von kármán experiment[END_REF], since it is a canonical configuration for the investigation of turbulence in a confined geometry. The radius and height of the cylindrical vessel are denoted R and H v , respectively. The stirring is done with two counterrotating impellers that are separated by a distance denoted H. A schematic representation of the experimental setup is shown in Figure 1 and a sketch of the impellers is shown in Figure 2. The impellers may rotate at different signed frequencies f 1 and f 2 . In the present paper, we focus our attention on the cases where the two impellers are counterrotating at the same frequency, i.e.,

f i := |f 1 | = |f 2 | and f 1 = -f 2 .
Looking at Figure 1, we adopt the following conventions: the vertical axis is oriented upward, and we use the right-hand rule to define the sign of rotations about the vertical axis. For instance, for the situation shown in Figure 1 the bottom impeller rotates in the positive direction and the top impeller rotates in the negative direction; that is 0 ≤ f 2 = -f 1 . This operating condition is referred to in the rest of the paper with the adjective CONTRA. Similarly when -f 1 = f 2 ≤ 0, we say that the configuration is ANTI. Notice that the CONTRA and the ANTI operating conditions produce different flows since the blades are curved as shown in Figure 2. The convention adopted in Figure 2 is that the supporting disk is in the background, the blades are in the foreground, and the shaft is behind the supporting disk. In this figure, the bottom impeller works in the CONTRA condition when it rotates counter-clockwise, and it works in the ANTI condition when it rotates clockwise. The ANTI flow regime is characterized by a larger level of fluctuations and larger energy dissipation than the CONTRA regime; as will be shown below, it is also subject to a spontaneous symmetry breaking bifurcation of the mean flow topology.

In the entire paper, all the velocity fields are non dimensionalized using the typical forcing velocity V 0 := 2πRf i based on the radius of the cylinder and the rotation frequency of the impellers. The control parameter of the VK flow is the Reynolds number, defined as:

R e := 2πf i R 2 ν -1 . (3)
Unless explicitly specified otherwise, all the lengths are now expressed in units of the cylinder radius R.

Numerical approximation method

The simulations corresponding to the VK experimental setup have been done using a code henceforth referred to as SFEMaNS (for Spectral/Finite Elements code for Maxwell and Navier-Stokes equations). This section describes the algorithm we have adopted to represent the moving solid

Top impeller

Bottom impeller

Frame of work: the code SFEMaNS

SFEMaNS uses a hybrid spatial discretization combining spectral and finite elements. In a nutshell the approximation in space is done by using a Fourier decomposition in the azimuthal direction and the continuous Hood-Taylor Lagrange element P 1 -P 2 (linear approximation for pressure and quadratic approximation for the velocity in the meridian section). All the discrete functions A are written in the generic form: [START_REF] Benard | Large-Eddy Simulation of wind turbines wakes including geometrical effects[END_REF] with (r, θ, z) the cylindrical coordinates, t the time and M the number of Fourier modes considered. The functions A m,cos h and A m,sin h belong to a finite element space (piecewise linear approximation for pressure and piecewise quadratic approximation for the velocity). The approximation in time is done by using a pressure-correction method described in Guermond and Shen [START_REF] Guermond | On the error estimates for the rotational pressure-correction projection methods[END_REF]. The moving counter-rotating impellers are accounted for by using a pseudo-penalty technique described in Pasquetti et al. [START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF].

A(r, θ, z, t) = A 0,cos h (r, z, t) + M m=1 A m,cos h (r, z, t) cos(mθ) + M m=1 A m,sin h (r, z, t) sin(mθ),
The full algorithm is detailed in §3.2. Modulo the computations of nonlinear terms with the fast Fourier transform, the linear problems at each time step for each Fourier mode in the meridian section are uncoupled and are thereby parallelized by using the message passing interface. The solution of each linear problem in the meridian section is further parallelized by using graph partitioning techniques from the METIS library (Karypis and Kumar [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF]) and subroutines from the portable extensible toolkit for scientific computation library (PETSc) (Balay et al. [START_REF] Balay | PETSc users manual[END_REF]), for the linear algebra. SFEMaNS has been thoroughly validated on numerous manufactured solutions and against other hydrodynamic codes (Giesecke et al. [START_REF] Giesecke | Influence of highpermeability discs in an axisymmetric model of the Cadarache dynamo experiment[END_REF], Hollerbach et al. [START_REF] Hollerbach | Parity-breaking flows in precessing spherical containers[END_REF], Jackson et al. [START_REF] Jackson | A spherical shell numerical dynamo benchmark with pseudo-vacuum magnetic boundary conditions[END_REF], Marti et al. [START_REF] Marti | Full sphere hydrodynamic and dynamo benchmarks[END_REF]).

Algorithm to enforce moving domains

Since the impellers move with opposite angular velocities, it is impossible to find a frame of reference where the fluid domain is time-independent. This problem is addressed by combining into a single computational domain the counter-rotating impellers, Ω solid (t), and the fluid domain, Ω fluid (t). Using the cylindrical coordinate system (r, θ, z) about the vertical axis with the convention that the vertical axis is oriented upwards, we then define the indicator function of the fluid domain, χ, as follows:

χ(r, θ, z, t) = 1 if (r, θ, z) ∈ Ω fluid (t) 0 if (r, θ, z) ∈ Ω solid (t). (5)
As shown in Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], the problem can be reformulated in the global computational domain by adding a penalty term of the form (1 -χ)(u -u obst )/α on the left-hand side of the momentum equation [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] where α is a user-dependent penalty parameter and u obst the velocity of the disks and the blades given by:

u obst (r, θ, z) = -sgn(f 2 )re θ if z > 0, sgn(f 2 )re θ if z ≤ 0, (6)
with sgn(f 2 ) equal to 1 if f 2 > 0, 0 if f 2 = 0, and -1 otherwise. A significant difficulty with this formulation when working with the Fourier approximation is that the term χu involves FFTs and thereby cannot be made implicit.

An elegant solution to this problem has been proposed in Pasquetti et al. [START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF]. It is showed therein that the problem can also be reformulated by solving the following balance equations:

(1 -χ) u -u obst α -1 Re ∆u + ∇p = χ(-∂ t u -(∇×u)×u + f ), (7a) ∇•u = 0, ( 7b
)
where u is the velocity field, p is the pressure field, f is the forcing term introduced in equation ( 1) and α is a penalty parameter. The key idea in Pasquetti et al. [START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF] is the observation that replacing α by τ , with τ the time step used to discretize the time derivative ∂ t u, gives a scheme that is stable and does not involve the term χu on the left-hand side.

The system [START_REF] Cappanera | Nonlinear stabilization of magnetohydrodynamic equations and applications to multiphase flows[END_REF] is approximated in time by using a pressure-correction method. For any time-dependent function v(t), we denote by v n the approximation of v at time t n = nτ . The velocity is updated by using the following time-stepping scheme:

3u n+1 2τ - 1 R e ∆u n+1 = -∇p n + (1 -χ n+1 ) 3u n+1 obst 2τ + χ n+1 4u n -u n-1 2τ -∇( 4ψ n -ψ n-1 3 ) + χ n+1 -(∇×u * ,n+1 )×u * ,n+1 + f n+1 , (8)
where u * ,n+1 = 2u n -u n-1 . Then, the pressure increment ψ n+1 is obtained by solving the following Poisson problem:

∆ψ n+1 = 3 2τ ∇•u n+1 . (9)
Finally, the pressure is updated as follows:

p n+1 = p n + ψ n+1 - 1 R e ∇•u n+1 . ( 10
)
The adaptation of the pseudo-penalty technique to pressurecorrection methods relies on the observation that the pressure increment ψ also needs to be penalized by the function χ in (8) to get a stable scheme. A proof of the stability of the algorithm with u obs = 0 is established in Cappanera [START_REF] Cappanera | Nonlinear stabilization of magnetohydrodynamic equations and applications to multiphase flows[END_REF].

Notice that the velocity and the pressure are solutions of the Navier-Stokes equations in the region where χ = 1, i.e., in the fluid domain Ω fluid (t). In the region where χ = 0, i.e., in Ω solid (t), the momentum equation reduces to

2 ) near the solid-fluid interface limits the global convergence rate in time to 1 2 , but when the Reynolds number is large enough the convergence rate in time becomes 3 2 as the time step τ scales like R e -1/2 , see section 1. So the higher the kinetic Reynolds number the smaller the term τ Re , i.e., the more accurate the method. Another method was used in [START_REF] Balaras | Modeling complex boundaries using an external force field on fixed cartesian grids in large-eddy simulations[END_REF][START_REF] Posa | Large-eddy simulations in mixed-flow pumps using an immersed-boundary method[END_REF][START_REF] Nicolaou | A robust direct-forcing immersed boundary method with enhanced stability for moving body problems in curvilinear coordinates[END_REF] based on the direct-forcing approach proposed by [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for threedimensional complex flow simulations[END_REF] where boundary body forces allow the imposition of boundary conditions on interfaces not coinciding with the computational grid.

Manufactured tests

To illustrate the convergence properties of the algorithm, we consider a set of manufactured solutions (χ, u, p) with a relative small Reynolds number R e = 100. The domain of computation Ω is set to {(r, θ, z) | 0.2 ≤ r ≤ 1; 0 ≤ θ ≤ 2π; -1 ≤ z ≤ 1} and the solutions considered are defined as follows:

                       χ(r, θ, z, t) = 1 r≥0.5 , u r (r, θ, z, t) = (2r -1) 2 sin(z + t)1 r≥0.5 , u θ (r, θ, z, t) = 0, u z (r, θ, z, t) = (2 - 1 r )(6r -1) cos(z + t)1 r≥0.5
+ (r -0.5) sin(2θ)1 r≥0.5 , p(r, θ, z, t) = r 2 z 3 cos(t) + r cos(θ),

The source term f in the momentum equation ( 8) is computed accordingly. It involves the first five Fourier modes, meaning M = 4 with M defined in equation (4).

To investigate the convergence properties of the algorithm, we perform two sets of tests. First we focus on the time convergence of the method by analyzing results obtained with a fixed mesh size h equal to 5×10 -3 . Then we study the global convergence of the algorithm by setting τ = 0.4h 2 and by performing tests on five different grids of mesh size h ∈ {0.05, 0.025, 0.0125, 0.005, 0.0025}. The tests are performed using M = 4, meaning that only the Fourier modes m ∈ {0, 1, 2, 3, 4} are computed.

Table 1 displays the L 2 -norm of the error on the velocity and the L 2 -norm of the error on the pressure with a fixed mesh size and time step τ ∈ {10 -3 , 5×10 -4 , 2.5×10 -4 , 10 -4 , 5×10 -5 }. The observed convergence rate is larger than or equal to 1 2 which is expected since the Reynolds number is not large (R e = 100). The results shown in table 2 are performed with τ = 0.4h 2 . As the error is dominated by the time error, which is of order 1 2 , we expect the global convergence rate to be equal to 1. Indeed, the displayed L 2 -norms of the error on velocity and L 2norm of the error on the pressure are compatible with the rate O(h). We note that the method has also been validated in Cappanera [START_REF] Cappanera | Nonlinear stabilization of magnetohydrodynamic equations and applications to multiphase flows[END_REF] against classical test cases such as the flow past a sphere at various Reynolds numbers. time step τ u (L 2 -error) rate p (L 2 -error) rate 10 -3

LES Scheme: Entropy viscosity stabilization

To avoid the accumulation of energy at the grid scale at high Reynolds numbers when the grid is not fine enough to resolve the Kolmogorov scale, we use a LES-like technique called entropy viscosity. This method, developed in Guermond et al. [START_REF] Guermond | Subgrid stabilized projection method for 2D unsteady flows at high Reynolds number[END_REF][START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF][START_REF] Guermond | From suitable weak solutions to entropy viscosity[END_REF], consists of adding a local artificial viscosity made proportional to the residual of the kinetic energy balance. This artificial viscosity is added on the right-hand side of (7a) in the form ∇•(ν E ∇u). This induces a nonlinear diffusion proportional to the local energy imbalance that in turn allows the unresolved scales to be better accounted for. The method has its roots in the notion of suitable weak solutions introduced by Scheffer [START_REF] Scheffer | Nearly one-dimensional singularities of solutions to the Navier-Stokes inequality[END_REF] and which has been shown by Caffarelli et al. [START_REF] Caffarelli | Partial regularity of suitable weak solutions of the Navier-Stokes equations[END_REF] to be the only reasonable notion of solution currently available for the 3D Navier-Stokes equations.

We now give some technical details on the computation of the entropy viscosity. Since the approximation mixes finite elements and Fourier approximation, we construct a three-dimensional mesh by considering the tensor product of the finite element mesh in the meridian section with the uniform azimuthal one-dimensional mesh induced by the Fourier approximation. Denoting by M the number of complex azimuthal Fourier modes, the mesh size in the azimuthal direction at the radius r is 2πr/(2M -1). For each two-dimensional finite element cells K, we denote by h K = min(min x∈K 2πr 2M -1 , diam(K)). Assuming that n ≥ 2, we define the residual of the momentum equation as follows:

Res n NS = u n -u n-2 2τ + (u n-1 • ∇)u n-1 - 1 R e ∆u n-1 + ∇p n-1 -f n-1 . ( 12
)
This residual is then computed at each time step and over every mesh cell in the real space. The local artificial viscosity is defined on each cell K by:

ν n R|K = h 2 K Res n NS • u n L ∞ (D K ) u n 2 L ∞ (D K ) . ( 13
)
where D K is the patch composed of the cells sharing one face with the cell K in the real space. The quantity ν n R|K is expected to be as small as the consistency error in smooth regions and to be large in the regions where the Navier-Stokes equations are not well resolved. To be able to run with CFL numbers of order O(1), we finally define the entropy viscosity as follows:

ν n E|K = min c max h K u n L ∞ (D K ) , c e ν n R|K , (14)
where c max = 1 8 and c e is a tunable constant O(1). In the following we set c e = 1. Thus defined, and given that we use P 2 polynomials to approximate the velocity, the entropy viscosity scales like O(h 3 K ) in smooth regions and scales like O(h K ) in regions with very large gradients.

This LES technique has been validated with c max = 1 and c e = 1 in Cappanera et al. [START_REF] Cappanera | Two spinning ways for precession dynamo[END_REF] for flows in precessing cylinders and in Wang et al. [START_REF] Wang | An entropy-viscosity large eddy simulation study of turbulent flow in a flexible pipe[END_REF] for turbulent flows in a flexible pipe (notice that the parameter α, defined therein in equation (2.13), is equal to c e /2). We have used this method in Nore et al. [START_REF] Nore | Numerical simulation of the von Kármán sodium dynamo experiment[END_REF] to perform high Reynolds number computations in a magnetohydrodynamics version of the von Kármán experiment.

Tools and flow description

We describe in this section the different diagnostic tools, the von Kármán flow, the experimental setup and the numerical parameters that we use.

Diagnostic tools 4.1.1. Time averages

Since we are going to regularly invoke time averages, for any time dependent quantity A : [0, T ] → R n we define:

A = 1 T T 0 A(τ )dτ . (15)
Moreover, given an experimental time series (A k ) 1≤k≤N measured at a fixed point in the inertial frame of reference of the laboratory, we define the inertial discrete time average as follows:

A • = 1 N N k=1 A k . ( 16
)
Notice that when the flow is time-dependent this definition implies averaging the characteristics of the quantity A not only in time but also with respect to the azimuth. Using the cylindrical coordinates (r, θ, z), this corresponds to extracting the time average of the angular Fourier mode m = 0 of A(r, θ, z, t).

Energy dissipation and torque

The power that is injected inside the VK flow and that is eventually transformed into heat by viscous effects can be estimated from the two torques C 1 and C 2 that are exerted on the top and the bottom impellers, respectively; the injected power in question is given by the expression

P inj = 2πf i (|C 1 | + |C 2 |).
In the statistically stationary regime, the mean input power equals the mean dissipation power. A useful diagnostic of the dissipation power can then be derived using a non-dimensional number P * as follows:

P * = P inj 2ρR 5 (2πf i ) 3 , (17)
where we recall that ρ is the density of the fluid. A nondimensionalized torque can also be defined as

K p = |C 1 | + |C 2 | 2ρR 5 (2πf i ) 2 . ( 18
)
Notice that K p = P * . In the experiment, the torque measurements at each impeller are performed using either the engines or torque meters. The torque applied to the top shaft is denoted C 1 , and the torque applied to the bottom shaft is denoted C 2 . Following the procedure described in Marié [START_REF] Marié | Angular momentum and magnetic field transport in a turbulent swirling flow: influence of rotation[END_REF], the torques are calibrated using measurements at different mean frequencies, so as to remove spurious contributions from genuine offsets or mechanical frictions. The torque measurements give the power dissipation.

Numerically, we can compute the dimensionless torque K p as follows:

K p = 1 2 Ω solid |(r×f s ) • e z |dΩ, (19)
where f s is the non-dimensional body force that induces the solid rotation of the impellers. The force f s can be written as ∇•(pI -R e -1 ∇u), where u and p are any smooth extensions of the velocity and the pressure in the solid. Notice that the term f s and the source term f in the Navier-Stokes equations (1) are unrelated. In the sequel f is set to zero. The fluid is driven by the movement of the counter-rotating impellers, represented by the pseudopenalty method, and not by an exterior forcing f . Using the notation from ( 5)-( 6), we deduce from the expression of the discrete momentum balance (8) that the torque at time t n+1 is given by

K p = 3 4 Ω r(1 -χ)sgn(f 2 z) u n+1 -u obst τ •e θ dΩ, ( 20
)
where sgn is the sign function introduced in (6).

The von Kármán flow 4.2.1. Symmetries

When the two impellers counter-rotate at the same frequency, i.e., when f 1 = -f 2 , the VK setup is symmetric with respect to any rotation of angle π about any axis in the equatorial plane that crosses the rotation axis. This type of symmetry is henceforth referred to as R πsymmetry (Nore et al. [START_REF] Nore | The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF]), and acts as follows in cylindrical coordinates:

R π   u r u θ u z   (r, θ, z) ≡   u r -u θ -u z   (r, -θ, -z) (21)
In the absence of blades on the supporting disks, the setup is also axisymmetric. At very low Reynolds numbers, the instantaneous velocity field obeys these symmetries. At large Reynolds numbers, the time averaged velocity is statistically axisymmetric but can experience breaking of the R π -symmetry in certain conditions as we will report below.

Mean flow topology

The topology of the time-averaged and azimuthallyaveraged mean flow (in short mean flow) is simple and depends on whether the forcing conditions are CONTRA or ANTI:

1. The CONTRA mean flow is divided into two toroidal recirculation cells separated by an azimuthal shear layer. In a vertical plane containing the axis of rotation, the corresponding mean velocity field is symmetric under R π (i.e., u r (r, -z) = u r (r, z) ,

u θ (r, -z) = -u θ (r, z) , u z (r, -z) = -u z (r, z) )
and there is a strong shear layer in the middle. 2. The ANTI mean flow has two possible geometries depending on the Reynolds number and the shape of the impellers. (1) The flow can have the same geometric symmetries as the CONTRA mean flow.

(2) The flow can have a bifurcated geometry resulting from the merging of the two toroidal counterrotating recirculation cells into a single cell (see Ravelet et al. [START_REF] Ravelet | Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation[END_REF] and Ravelet et al. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF]). In this case, the mean velocity field in any vertical plane containing the axis of rotation is no longer invariant under R π . The mean flow is then mainly composed of one cell in the vertical direction with a strong shear layer at the impeller that rotates in the direction opposite to the orthoradial mean flow. This bifurcated state only exists for Reynolds numbers that are large enough (see Ravelet et al. [START_REF] Ravelet | Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation[END_REF]) and for impellers that are fitted with blades that are sufficiently curved, which is the case of the TM87 impellers studied in the present paper. This turbulent bifurcation results in multistability between the two turbulent flow states, with possible complex dynamics between them (see challenge to reproduce numerically this configuration.

2.5×10 -3 2.5×10 -3 1.25×10 -3 1.25×10 -3 10 -4 h min 2×10 -2 5×10 -3 5×10 -3 5×10 -3 2.5×10 -3 h max 2×10 -2 10 -2 2×10 -2 2×10 -2

Experimental implementation 4.3.1. Experimental setup

The VK experimental setup used for the present study has been thoroughly described in Ravelet et al. [START_REF] Ravelet | Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation[END_REF][START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF], Monchaux [START_REF] Monchaux | Statistical mechanics and dynamo action in a turbulent von Kármán flow[END_REF], Saint-Michel et al. [START_REF] Saint-Michel | Influence of Reynolds number and forcing type in a turbulent von Kármán flow[END_REF]. The fluid is confined inside a cylinder of radius R = 100 mm, and put in motion by two rotating impellers of radius 92.5 mm (see Figure 1 and Figure 2). We recall that the lengths are nondimensionalized with respect to the radius of the cylinder, R. The aspect ratio of the experiment is defined as the (non-dimensional) distance between the inner faces of the two disks supporting the blades, H = 1.8. The turbulence properties (anisotropy, fluctuations, dissipation) are influenced by the geometry of the impellers, their non dimensional radius R t , the oriented angle α (see Figure 2) and height h b of the blades, and the number n of blades (Ravelet [31]). In the present paper, we consider mainly TM87 impellers, with n = 8 blades, h b = 0.2, |α| = 72 • and R t = 0.925. Mixtures of water and glycerol with different dilution rates are used to change the viscosity of the fluid. The temperature is maintained constant by means of heat exchangers. A summary of the experimental cases presented in this paper is given in Table 3.

Both the torques and the mean flow topology are very sensitive to geometric parameters like the height and the curvature of the blades, the radius of the impellers, the aspect ratio, as discussed in Marié [START_REF] Marié | Angular momentum and magnetic field transport in a turbulent swirling flow: influence of rotation[END_REF], Ravelet et al. [START_REF] Ravelet | Towards an experimental von Kármán dynamo : numerical studies for an optimized design[END_REF], Ravelet [START_REF] Ravelet | Bifurcations globales hydrodynamiques et magnétohydrodynamiques dans un écoulement de von Kármán turbulent[END_REF], Monchaux [START_REF] Monchaux | Statistical mechanics and dynamo action in a turbulent von Kármán flow[END_REF], Burnishev and Steinberg [START_REF] Burnishev | Torque and pressure fluctuations in turbulent von Karman swirling flow between two counter-rotating disks[END_REF]. Therefore to be able to reproduce numerically the experiments, we have observed that it is essential that the geometry of the setup be reproduced as accurately as possible in SFEMaNS.

Multi-scale velocity measurements

The velocity measurements used in the present paper result from two types of techniques: (i) Stereoscopic Particle Image Velocimetry (SPIV) and (ii) Laser Doppler Velocimetry measurements (LDV).

The SPIV system provides the three components of the velocity field in a meridian plane on a grid of typical size 90×70 points.The optical device can be adapted so that the horizontal and the vertical distance between two measurement points can be controlled in the range [2.4×10 -3 , 2.4×10 -2 ]. The meridian section that can be explored is {(r, z) ∈ [0, 1]×[-0.75, 0.75]}. The time series are composed of about 3000 to 30000 frames regularly sampled at 15 Hz. The time average of time series of SPIV measurements is denoted . • , see [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF]. Because of the constraints on the measurements technique and the unsteadiness of the velocity pattern, such time-average automatically produces an azimuthal average, i.e., time averaging projects the measurements onto the angular Fourier mode m = 0.

We have also performed a few Laser Doppler Velocimetry measurements. These measurements provide the mean azimuthal and mean vertical components of the velocity in a meridian plane on a grid composed of 11×17 points located at 0 ≤ r ≤ 10∆r with ∆r = 0.1 and -8∆z ≤ z ≤ 8∆z with ∆z = 0.1125.

Summary of the numerical parameters

The numerical parameters that have been used in the various simulations reported in this paper are listed in Table 4. The computations are done only in the TM87 configuration.

The spatial resolution in the meridian plane and in azimuth of DNS simulations is set to match the Reynolds number, i.e., the computational grid is refined as R e increases. The meshes are usually coarser for LES runs than for DNS runs. Fine meshes are used to simulate the ANTI bifurcated flow since very thin shear layers are created in this case (see Figure 13 (b)-(e) for 0.6 ≤ z ≤ 0.8). The mesh sizes h min and h max in Table 4 are the typical distance between two grid points on the velocity mesh (i.e., the P 2 mesh). Between 128 to 512 real Fourier modes are typically used. The shape of the impellers and a computational grid are displayed in Figure 3. The parallelization is done with one complex Fourier mode per processor, and the meridian plane is further divided among the processors by using a domain decomposition technique, the graph partitioning being done by METIS. The linear algebra in the meridian section is handled by PETSc and the fast Fourier transforms are done with FFTW3. Depending on the spatial resolution and the Reynolds number, one rotation period (one turn) requires between 2 to 63 wall-clock hours on the cluster IBM x3750-M4 from GENCI-IDRIS. Each run does between 15 to 60 turns. The cumulated computing time for the runs presented in this article is about 7×10 5 CPU hours on one processor.

Comparisons Experiment vs. Numerics

In this section we compare experimental and numerical velocity profiles in a meridian section of the cylindrical container at various Reynolds numbers and in different operating conditions. The comparisons are done on snapshots and on the time-averaged velocity fields. In order to have good quantitative comparisons, we always use the same colorbars for the experimental data and the numerical results. At the end of the section we also compare torque measurements with computations over a wide range of Reynolds numbers.

Flow topology at low Reynolds numbers

We start by investigating the VK flow at low Reynolds numbers; i.e., R e ≤ 10 3 .

Low Reynolds numbers: steady-state regime

At very low Reynolds numbers (ie R e < 5×10 2 ), the velocity field is stationary. Moreover, the CONTRA and the ANTI operating conditions give very similar results. Therefore we do not show the two cases but rather focus on the CONTRA operating mode. We show in Figure 4 the three components of the stationary and axisymmetric velocity field at R e = 10 2 using the cylindrical coordinate representation. We compare in panels (a) and (b) the experimental measurements of the radial component of the velocity with the computational results. The comparisons for the azimuthal component are shown in panels (c) and (d). The comparisons for the axial component are shown in panels (e) and (f). The white zones in the experimental fields correspond to areas where measurements were not possible due to the presence of the blades. The resolution of the LDV measurements being much lower than the numerical results, we have interpolated the numerical data on large pixels corresponding to the experimental resolution. This process allows for a better comparison between the experiments and the numerical simulations. Figure 4 shows that the numerical data and the experimental results are very similar with slight differences close to the border of the experimental acquisition zone.

Low Reynolds numbers: time-dependent regime

The flow becomes time-dependent at R e = 500. We compare in Figure 5 the experimental and the numerical time-averaged velocity fields in the CONTRA configuration at R e = 10 3 using the same pixelization procedure as described above. Notice that this time we use the averaging operator • • defined in [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF] to average the experimental data in time. For the numerical simulations, the averaging is done in azimuth and in time; that is, we show the time average of the Fourier mode m = 0. Figure 6 shows radial profiles of averaged azimuthal and vertical velocity components at z ∈ {-0.4, 0, 0.4}. For the LDV experimental data, measurements are performed in one meridian section. Since for all z ∈ Z := {-0.4, 0, 0.4} the quantities u θ (r = 0, z), u θ (r = 1, z) and u z (r = 1, z) must be equal to zero, we estimate the experimental error on u θ at z ∈ Z as follows: ∆ θ (z) := max{|u We now compare the CONTRA and the ANTI regimes at R e = 10 3 . In Figure 7 we compare the numerical results obtained in the CONTRA operating configuration with the numerical results obtained in the ANTI operating configuration. The flow patterns look different. Close to the blades, the radial centrifugal component of the velocity in the CONTRA case is stronger than in the ANTI configuration, but in the equatorial shear layer the radial velocity component is more intense and focused in the ANTI regime. Note also that the azimuthal component of the velocity is overall stronger in the ANTI than in the CON-TRA regime. These features persist at larger Reynolds numbers for the averaged velocities.

Flow topology at high Reynolds numbers

We now describe the flow topology at high Reynolds numbers. The experimental data reported below have been obtained at R e = 3×10 5 and the LES computations have been done at R e = 10 5 . Therefore the comparisons are performed at sligthly different Reynolds numbers. At these Reynolds numbers, the CONTRA configuration gives a solution that is highly turbulent but on average the flow is organized into two cells that are invariant under the R π transformation. This structure is very robust with respect to the preparation of the flow. This is not the case for the ANTI configuration. The flow is also highly turbulent in this operating mode, but, depending on the preparation of the flow, one observes two possible timeaveraged states. One state is composed of two cells as in the CONTRA operating conditions, but the other one is composed of one cell only. We have observed this bifurcation to occur in the ANTI configuration at R e ∼ 10 4 (see §5.3 and Figure 15 for a detailed discussion).

CONTRA operating mode

We start by comparing the results for the CONTRA operating mode. We show in Figure 8 the time and azimuthally averaged velocity field at R e ∼ 10 5 . The simulation has been done by using the final snapshot of a statistically converged simulation at R e = 10 4 as initial data. About 20 turns have been performed at R e = 10 5 .

Since the simple time averages over the 20 rotation periods of the LES computations are not long enough to be fully converged statistically, we show in Figure 8 to the boundary). Figure 9 shows averaged radial profiles of the azimuthal and the vertical velocity components at z ∈ {-0.4, 0, 0.4}. For the SPIV experimental data, measurements are performed over an entire meridian plane, i.e., measurements are done simultaneously on two diametrically opposite meridian sections. The quantities displayed are u exp z (r, z) = 1 2 (u z (r, 0, z) + u z (r, π, z)) and u exp θ (r, z) = 1 2 (u θ (r, 0, z) -u θ (r, π, z)). The error bars for each component are the variance of the temporal sig- nals. This estimation of the experimental error is also used in Figure 12 and Figure 14. The numerical profiles fit the error bar region of the experimental data. This suggests that the difference in the Reynolds numbers (10 5 for the numerical simulations and 3×10 5 for the experiments) has a moderate impact on the average profiles. Furthermore these profiles are similar to the ones obtained at the much smaller R e = 10 3 : the local maximum of |u θ | is around r = 0.8 at z = ±0.4 and the change of sign of u z occurs at r ≈ 0.75 for z = ±0.4.

We compare in Figure 10 experimental and computational snapshots of the Cartesian components of the velocity field in a vertical plane passing through the axis. These figures clearly show that the flow is highly turbulent. Of course, these instantaneous snapshots are not identical but they share similar amplitudes and structures.

ANTI operating mode: symmetric solution

We now show in Figure 11 the time-averaged experimental and numerical velocity fields for the ANTI configuration.

Here the experimental flow at R e = 3×10 5 is obtained by progressively increasing the angular frequency of the impellers and making sure that the two angular velocities are all the time exactly opposite. The preparation of the numerical simulations at R e = 10 5 is also done by progressively increasing the Reynolds number and by always enforcing the two angular velocities to be exactly opposite as explained in §5.2.1. This process leads to a highly turbulent flow that is organized on average into two cells that are R π symmetric. The organization into two cells is clearly visible in panels 11(c) and 11(f). Notice, though, that the CONTRA and ANTI mean flow fields are very different. The differences are particularly noticeable when comparing the radial component of the velocity in the panels 8(a) and 8(d) with that in the panels 11(a) and 11(d). Also, by comparing the panels 8(b) and 8(e) with the panels 11(b) and 11(e), we observe that the ANTI configuration produces stronger azimuthal components of the velocity than the CONTRA configuration. In the ANTI configuration the large values of the azimuthal component are concentrated near the lateral wall. Notice also that the radial and axial components of the velocity have smaller values than in the CONTRA configuration.

Figure 12 compares the experimental and numerical profiles for the symmetric ANTI configuration at high Reynolds numbers. Again the agreement is good since the numerical profiles are inside the error bar region of the exper-

ANTI operating mode: bifurcated solution

We now focus on the bifurcated flow for the ANTI configuration. This solution is obtained by preparing the flow in a nonsymmetric way. In the experiment, the angular frequencies of the top and the bottom impellers are not increased simultaneously. For some time one of the impellers rotates faster than the other one. Then, eventually, exact counter-rotation is prescribed and maintained. This produces a symmetry breaking that allows the flow to explore another solution branch that is not invariant under R π . For the numerical simulations we proceed as follows. We use a snapshot of the ANTI configuration at R e = 10 4 as initial data. We perform 15 rotation periods at R e = 5×10 4 in the ANTI configuration, but we reduce the angular velocity of the top impeller by setting f 1 = -1 2 f 2 with the Reynolds number defined with respect to the angular velocity of the bottom impeller only. Then we do 10 more rotation periods with f 1 = 0, f 2 being unchanged. After 10 rotation periods, the structure of the mean flow shifts from two recirculation cells to a single recirculation cell. The angular velocity of the top impeller is then increased so as to match the angular velocity of the bottom impeller: 5 rotation periods are done with f 1 = -f 2 . Finally the Reynolds number is increased to R e = 10 5 and 12 periods are done at this Reynolds number.

We show in Figure 13 the averaged velocity fields in a meridian section at R e = 10 5 for the LES computations and at R e = 3×10 5 for the experiments. We clearly see that the mean flow is composed of one large recirculation cell only, which is very different from the symmetric case. Here again, the numerical results agree reasonably well with the experimental observation despite the SPIV's being blind to the impellers and to the boundary layer at r = 1.

This statement is supported by inspecting the radial profiles shown on Figure 14. The profiles of the azimuthal component of the velocity from the numerical and from the experimental data do not coincide precisely but they have similar shapes; we observe in particular that u θ is mostly negative at the three vertical heights z ∈ {-0.4, 0, 0.4}. The agreement on the vertical component of the velocity is significantly better; the vertical profiles are almost invariant with respect to z. The behavior of u θ and u z indicates that only one cell is present and that it is the bottom impeller that imposes its sense of rotation.

Torque vs. R e

We now compare the measurements and the numerical computations of the non-dimensional torque K p defined in [START_REF] Hollerbach | Parity-breaking flows in precessing spherical containers[END_REF]. (Recall that K p also measures the dissipation power as discussed in §4.1.2.) All the results are reported in one single graph shown in Figure 15. We show there the experimental measurements and the numerical estimations of K p as a function of the Reynolds numbers over the range R e ∈ [10, 10 6 ]. The grey symbols, crosses and stars correspond to measurements. The circles correspond to numerical simulations; empty circles are for DNS simulations and yellow-filled circles are for LES simulations. We observe three plateaus at large Reynolds numbers. The blue line (bottom horizontal line) corresponds to data from the CONTRA configuration, the red line (middle horizontal line) corresponds to data from the symmetric flow in the ANTI configuration, and the green line (top horizontal line) corresponds to data from the bifurcated flow in the ANTI configuration.

For R e ≤ 4×10 2 the CONTRA configuration (blue symbols) and the ANTI configuration (red symbols) give the same torque, both in the numerical simulations and in the experiments. The two curves split at R e ∼ 4×10 2 . The torque for the CONTRA configuration seems to decrease monotonously with respect to the Reynolds number over the entire range of Reynolds numbers and to converge to an asymptotic value K ∞ cont ∼ 0.05 at very large values of R e . After reaching a minimum in the range R e ∈ [5×10 2 , 5×10 3 ], the torque for the ANTI configuration grows again after R e = 10 3 and seems to converge towards an asymptotic value K 2,∞ anti ∼ 0.14 at very large values of the Reynolds number. These asymptotic values of the torque coincide with experimental measurements in liquid sodium (R e ∼ 10 7 ) and in superfluid Helium, see St-Michel et al. [START_REF] St-Michel | Evidence for Forcing-Dependent Steady States in a Turbulent Swirling Flow[END_REF] and Dubrulle [START_REF] Dubrulle | Beyond Kolmogorov cascades[END_REF].

The bifurcation in the ANTI configuration discussed in §5.2.2 and §5.2.3 occurs around R e ≈ 10 4 . The bifurcated ANTI solution is composed of one recirculation cell rotating somewhat in phase with one of the impellers (Ravelet et al. [START_REF] Ravelet | Towards an experimental von Kármán dynamo : numerical studies for an optimized design[END_REF][START_REF] Ravelet | Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation[END_REF]). The bifurcated ANTI flow dissipates far more energy than its symmetric counterpart. The probable cause for this higher energy dissipation rate is that the velocity undergoes very large shears in the vicinity of the top impeller as can be seen in Figure 13(e). Ac-tually this simulation requires a finer meridian grid and a larger number of azimuthal modes than the symmetric ANTI flow. The mesh refinement in the meridian section is done locally in the vicinity of the top impeller (see table 4). The torque applied to the bifurcated ANTI flow seems to converge to the asymptotic value K 1,∞ anti ∼ 0.56 for large values of R e .

Overall the experimental measurements and the numerical estimations coincide up to 10 percents over the entire range of Reynolds numbers explored and for the two forcing conditions, including in the bifurcated case. The fact that we have been able to reproduce the bifurcated branch and to estimate accurately the torque (i.e., the dissipation power) at R e = 10 5 is quite remarkable considering that we are using a LES model in this range. This means that the LES model dissipates the energy properly. The method stabilizes the computation without introducing excessive numerical dissipation; i.e., the energy is allowed to cascade freely and is dissipated at the smallest mesh scale at the correct rate. To the best of our knowledge, it is the first time that numerical simulations reproduce such bifurcated branches of turbulent flows.

Energy Spectra

This section investigates the energy spectra of the numerical simulations done at R e = 10 5 using the LES model. The results are interpreted in the context of the previous studies [START_REF] Dubrulle | Beyond Kolmogorov cascades[END_REF][START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF].

Spatial spectrum

In 3D periodic numerical simulations it is common to compute the spatial energy spectrum as:

E(k) = u(x, t) • u(x + r, t)e ik•r dr x,t . (22)
Since for isotropic flows E(k) only depends on k , it is therefore natural to consider:

E(k) := E(k) k =k . (23)
We use the LES simulations at R e = 10 5 to estimate E(k).

The value of E(k) is computed from snapshots of the numerical velocity field extracted from the cube (-1 2 , 1 2 ) 3 located at the center of the tank. The spectra are averaged over the snapshots. The energy spectrum E(k) is supposed to scale as k -5 3 in homogeneous isotropic turbulence. For the family of von Kármán flows considered in the paper, Dubrulle [START_REF] Dubrulle | Beyond Kolmogorov cascades[END_REF] suggests the universal scaling function E(k)/( πH K p and η = (R e 3 ) -1/4 . Figure 16 shows E(k)/(

2 3 η 5
3 ) for the three flow configurations. The low wavenumbers correspond to large scales and the large wavenumbers correspond to the inertial range. We observe that in the inertial range the spectra are in agreement with the k -5/3 law (see dotted line on Figure 16). Note that the ANTI bifurcated flow has much less energy than the CONTRA and symmetric ANTI flows in the low wavenumber region. Recall that the ANTI bifurcated flow is essentially composed of one recirculation cell with a thin region where the velocity gradients are very large, see the region 0.6 ≤ z ≤ 0.8 in Figure 13. The slopes of the lines fitting best these curves in the inertial zone are -1.7 for the CONTRA case (for 3 10 -3 ≤ kη ≤ 4 10 -2 ), -1.6 for the symmetric ANTI (for 3 10 -3 ≤ kη ≤ 4 10 -2 ), and -1.7 for the bifurcated ANTI case (for 2 10 -3 ≤ kη ≤ 1 10 -1 ).

Temporal spectra

In addition to spatial spectra, one can also compute temporal spectra from time series of velocity measurements at a fixed point. Following Ravelet et al. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF] we select the following point x 0 (r = 0.9, θ = 0, z = 0) in the computational domain. The power spectral density (PSD) is then defined as:

E(f ) =
u θ (x 0 , t)u θ (x 0 , t + s)e -2iπf s ds . [START_REF] Marti | Full sphere hydrodynamic and dynamo benchmarks[END_REF] For the ANTI configuration and R e ≤ 6.5×10 3 , Ravelet et al. [START_REF] Ravelet | Supercritical transition to turbulence in an inertially driven von Kármán closed flow[END_REF] observed the power-law f -1 in the low-frequency regime (below the impeller frequency f i ) and the powerlaw f -5 3 in the inertial range (above the impeller frequency f i ). Note that the power-law f -1 is difficult to observe because it requires very long statistics. Figure 17 presents the PSD extracted from the numerical simulations at R e = 10 5 for the three flow configurations. Depending on the configuration, the time integration window in [START_REF] Marti | Full sphere hydrodynamic and dynamo benchmarks[END_REF] ranges from 10 to 20 rotation periods. We observe in Figure 17 three regions in the frequency domain: (i) There is a low frequency range corresponding to slow motions of the shear-layer in the CONTRA and symmetric ANTI cases; (ii) There is the intermediate inertial range. In this range the spectra behave like f -5/3 ; (iii) There is the large frequency range corresponding to the dissipation zone. Slope fitting in the inertial range gives -1.7 for 1 ≤ f /f i ≤ 50 in the CONTRA case, -1.5 for 1 ≤ f /f i ≤ 50 in the symmetric ANTI case, and -1.7 for 6.5 ≤ f /f i ≤ 125 in the bifurcated ANTI case. These slopes are in agreement with the Kolmogorov exponent. Figure 17: Power spectral density E(f ) of u θ (x 0 , t) for different forcing conditions at Re = 10 5 : blue is for the symmetric CONTRA branch, red for the symmetric ANTI branch, and green for the bifurcated ANTI branch. f is the analysis frequency and f i is the impeller rotation frequency. The dashed line shows the behavior f -5 3 in the inertial range.

Conclusion

In this paper, we have compared experimental data and numerical simulations (DNS & LES) for a complex turbulent system which is sometimes humorously referred to in the literature as the "French washing machine." We have successfully reproduced numerically the different flow types observed experimentally in the CONTRA configuration and in the ANTI configuration. We have been able to reproduce the bifurcation observed at high Reynolds numbers in the ANTI regime: there is a turbulent symmetric solution whose temporal average consists of two recirculating zones that are invariant under the R π transformation, and there is a bifurcated solution whose temporal average consists of one recirculation cell only. We have focused our comparisons on two criteria: the torque value and the flow topology at various Reynolds numbers. The qualitative comparisons of the flows show good agreements between the numerical simulations and the experiments, except near the edges of the experimental measurement area. The measurements and the numerical computations of the torque coincide over a wide range of Reynolds numbers for the three flow types. Obtaining similar outputs for criteria that are so dependent on forcing conditions makes the comparison a successful one. Moreover spatial and temporal spectra computed at R e = 10 5 for the three configuration flows show an inertial range compatible with the -5/3 Kolmogorov exponent characterizing fully turbulent flows.

This study also proves the effectiveness of the proposed entropy viscosity stabilization as a LES model. The simulations reported in the paper validate the performance of the proposed model on a highly anisotropic turbulent flow in a complex geometry that is time-dependent with counter rotating motions that make the geometry complex in any referential. The entropy viscosity model involves two parameters c max and c e (see section 3.4), but we have observed that the following choice (c max , c e ) = (1/8, 1) is robust in the sense that it performs well independently of the problem. In contrast to other LES models, like the model of Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations, part i: the basic experiment[END_REF] which involves a parameter that needs to be tuned depending of the problem (Li et al. [START_REF] Li | Subgrid-scale modeling of helicity and energy dissipation in helical turbulence[END_REF], Delafosse et al. [START_REF] Delafosse | LES and URANS simulations of hydrodynamics in mixing tank: comparison to PIV experiments[END_REF]), the robustness of the entropy viscosity stabilization with respect to (c max , c e ) facilitates its use on a wide range of problems.

Figure 1 :

 1 Figure 1: Experimental setup with H the distance between the inner faces of the disks, R the cylinder radius, f i the frequency of each disk. Impellers with 8 blades are called TM87. The arrows show the rotation direction for the CONTRA configuration.

Figure 2 :

 2 Figure 2: TM87 impellers (8 blades) with radius Rt and angle |α| = 72 • . The supporting disks are in the background, the blades are in the foreground, the shafts are behind the supporting disks.

8

 8

Figure 3 :

 3 Figure 3: Shape of the impellers and the computational grid in (r,z) plane with 128 Fourier modes, hmax = 2×10 -2 and h min = 5×10 -3 .

Figure 4 :

 4 Figure 4: Stationary and axisymmetric experimental and DNS velocity field in the CONTRA configuration at Re = 10 2 : (a)-(d) radial component ur, (b)-(e) azimuthal component u θ , (c)-(f) vertical component uz.

Figure 5 :

 5 Figure 5: Time-averaged experimental and DNS velocity field (m = 0) in the CONTRA configuration at Re = 10 3 : (a)-(d) radial component ur, (b)-(e) azimuthal component u θ , (c)-(f) vertical component uz.

Figure 6 :

 6 Figure 6: Radial profiles of u θ and uz for time-averaged experimental and DNS velocity field (m = 0) in the CONTRA configuration at Re = 10 3 : u exp • in solid line with errorbars, u m=0 in dashed line at various z as indicated. The color convention applies to all the profiles in other figures.

Figure 7 :

 7 Figure 7: Time-averaged DNS velocity field (Fourier mode m = 0) in CONTRA and ANTI operating configurations at Re = 10 3 : (a)-(d) radial component ur, (b)-(e) azimuthal component u θ , (c)-(f) vertical component uz.

  only the time average of the Fourier mode m = 0 of the velocity field. This corresponds to the experimental data mean operator • • which only detects the Fourier mode m = 0 anyway. Notice that this time the image resolutions are now different; one sees the blades in the numerical simulations but these are not visible in the experiments. The blades start at r = 0.1 as shown in Figure 2. Also the SPIV allows us to reach only r = 0.95. This explains why we see a boundary layer in the LES simulation close to r = 1 in panel 8(d) which is not captured by the SPIV measurements in the panel 8(a) (white pixels indicate that there are no experimental data close

Figure 8 :Figure 9 :

 89 Figure 8: Time and azimuthal averaged velocity field in the CON-TRA configuration. Re = 3×10 5 for the experiments; Re = 10 5 for the computations: (a)-(d) radial component ur, (b)-(e) azimuthal component u θ , (c)-(f) vertical component uz.

Figure 10 :

 10 Figure 10: Instantaneous velocity fields in the CONTRA configuration in a vertical plane passing through the axis. Top row: experiments at Re = 3×10 5 ; bottom row: numerical simulations at Re = 10 5 ; same colorbars. Left column: ux; central column: uy; right column: uz. Horizontal axis: x; vertical axis: z.

Figure 11 :Figure 12 :

 1112 Figure 11: Time and azimuthal averaged velocity field in the symmetric ANTI configuration. Experiments at Re = 3×10 5 ; numerical simulation at Re = 10 5 : (a)-(d) radial component ur, (b)-(e) azimuthal component u θ , (c)-(f) vertical component uz.

Figure 13 :

 13 Figure 13: Time and azimuthal averaged velocity field in the bifurcated ANTI configuration. Experiments at Re = 3×10 5 ; numerical simulations at Re = 10 5 : (a)-(d) radial component ur, (b)-(e) azimuthal component u θ , (c)-(f) vertical component uz.

Figure 14 :

 14 Figure 14: Radial profiles of u θ and uz for the time-averaged experimental and LES velocity fields (m = 0) in the bifurcated ANTI configuration at Re = 3×10 5 for the experiment and Re = 10 5 for the computation: u exp • in solid line with errorbars, u m=0 in dashed line. Color labels are defined in Figure 6.

  f (kη) where = 2R

Figure 16 :

 16 Figure 16: Energy spectrum E(k) as a function of the wavenumber k for different forcing conditions at Re = 10 5 : blue is for the symmetric CONTRA branch, red for the symmetric ANTI branch, and green for the bifurcated ANTI branch. Dotted line corresponds to k -5 3 .

Table 1 :

 1 L 2 -norm of the errors on the velocity and outer pressure at time t = 0.1 and rates of convergence. The mesh size h is set to 5×10 -3 in P 2 .

		6.16E-3	-	8.55E-3	-
	5×10 -4	4.40E-3	0.49	6.10E-3	0.49
	2.5×10 -4	3.14E-3	0.49	4.36E-3	0.48
	10 -4	1.90E-3	0.55	2.64E-3	0.55
	5×10 -5	1.12E-3	0.76	1.51E-3	0.81
	mesh size h u (L 2 -error) rate p (L 2 -error) rate
	0.05	1.91E-2	-	3.05E-2	-
	0.025	1.23E-2	0.63	1.72E-2	0.83
	0.0125	5.67E-3	1.12	8.10E-3	1.09
	0.005	2.51E-3	0.89	3.45E-3	1.23
	0.0025	5.31E-4	2.24	6.14E-4	2.49

Table 2 :

 2

L 2 -norm of the errors on the velocity and outer pressure at time t = 0.1 and rates of convergence. The time step is set to 0.4h 2 with h the mesh size in P 2 .

Table 3 :

 3 Experimental parameters: kinetic Reynolds number Re, flow type, kinematic viscosity ν, rotation frequency f

	R e	10 2	10 3	3×10 5	3×10 5	3×10 5
	Flow type CONTRA CONTRA CONTRA ANTI symmetric ANTI bifurcated
	ν(m 2 s -1 )	4.3×10 -4	4.3×10 -4	1.0×10 -6	1.0×10 -6	1.0×10 -6
	f (Hz)	0.68	6.8	5.0	5.0	5.0
	R e	10 2	10 3	10 5	10 5	10 5
	Flow type CONTRA CONTRA CONTRA ANTI symmetric ANTI bifurcated
	Model	DNS	DNS	LES	LES	LES
	τ

Thalabard et al. [41]

, Saint-Michel et al.

[START_REF] Saint-Michel | A zero-mode mechanism for spontaneous symmetry breaking in a turbulent von Kármán flow[END_REF]

, St-Michel et al.

[START_REF] St-Michel | Evidence for Forcing-Dependent Steady States in a Turbulent Swirling Flow[END_REF]

). It is therefore a genuine

Table 4 :

 4 Numerical parameters for the computations: kinetic Reynolds number Re, flow type, numerical model DNS or LES, timestep τ , mesh size in the blade region h min , mesh size at the outer boundary hmax (the meridian mesh is non-uniform), number of grid points in the P 2 meridian mesh n df , number of real Fourier modes, number of processors.

	10 -2

  Non-dimensional torque Kp as a function of Re for different forcing conditions: blue is for the symmetric CONTRA branch, red for the symmetric ANTI branch, and green for the bifurcated ANTI branch. Circles are numerical simulations with TM87 (empty circles are for DNS and yellow-filled circles are for LES). The other symbols correspond to other experimental results: stars are TM87 experiments using a mix of glycerol and water and diamonds are TM87 experiments with water. Grey dots are TM60 experiments (same impellers but with 16 blades rather than 8) with a mix of glycerol and water from Ravelet et al.[START_REF] Ravelet | Towards an experimental von Kármán dynamo : numerical studies for an optimized design[END_REF]. Dashed lines correspond to asymptotic values measured in liquid sodium (Re ∼ 10 7 ) and in superfluid Helium 4 at 2.3 • K obtained from the SHREK experiment.

		DNS LES Experiments
	10 0
	10 -1
	10 2	10 4	10 6
	Figure 15:

Acknowledgments

The HPC resources for SFEMaNS were provided by GENCI-IDRIS (grant 2018-0254) in France and by the Texas A&M University Brazos HPC cluster. J.-L. Guermond acknowledges support from University Paris Sud, the National Science Foundation under grants DMS 1620058, DMS 1619892, the Air Force Office of Scientific Research, USAF, under grant/contract number FA9550-18-1-0397 and the Army Research Office, under grant number W911NF-15-1-0517. L. Cappanera is thankful to Texas A&M University, LIMSI, and CNRS for their financial support. Part of this work has been supported by ANR EXPLOIT, grant agreement no. ANR-16-CE06-0006-01.


Mario Sánchez Aguilar

Uffe Linda Marie Ahl

Thomas Jankvist

Ola Helenius

Mario Sánchez Aguilar

Uffe Thomas Jankvist

A Delphi-inspired study on the concepts of large and small-scale projects in mathematics education research

Keywords: Implementation research, large-scale project, small-scale project, conceptualization

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Expressions such as "large-scale study" or "scaling-up" frequently appear in discussions and conceptualizations about the implementation of educational innovations (e.g., [START_REF] Century | Implementation research: Finding common ground on what, how, why, where, and who[END_REF][START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF], including the field of mathematics education (e.g., [START_REF] Ahl | Recent discussions on implementation research at CERME12[END_REF][START_REF] Heck | Scaling up innovative learning in mathematics: Exploring the effect of different professional development approaches on teacher knowledge, beliefs, and instructional practice[END_REF][START_REF] Roesken-Winter | Effective CPD on a large scale: Examining the development of multipliers[END_REF]. Although the notion of "scale" is widely used in educational research, its meaning tends often to be handled tacitly, and there is diversity in how the term "scale" is used [START_REF] Morel | The multiple meanings of scale: Implications for researchers and practitioners[END_REF]. The field of mathematics education research is not exempt from this lack of conceptual clarity in relation to the use of the notion of scale. With some exceptions (e.g., [START_REF] Krainer | Reflections on the increasing relevance of large-scale professional development[END_REF], notions such as "scaling-up", "at scale", or "large-scale" are used implicitly, giving rise to a variety of interpretations of these terms. In pursuit of conceptual clarity, we have undertaken a Delphi-inspired study that explores how the notions of "large-scale project" and "smallscale project" are understood in mathematics education research. The results of our study constitute a first attempt at how to characterize a "large-scale project" and a "small-scale project" in mathematics education. The following research question guides the study: In relation to implementation research, how are the notions of a "large-scale project" and a "small-scale project" interpreted in mathematics education research?

Although it may be impossible to answer this question in a definite manner, at least in this paper, an approximation is offered through an analysis of the opinions of experienced researchers. To analyze the experts' opinions, we adopt a classic conceptualization of the notion of "scale" proposed by [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF], which is informed by research on reform implementation. [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF] conceptualization of "scale" emerges as a response to the dominant perspective in educational research that tends to conceptualize studies at scale in a unidimensional way, that is, "involving solely or predominantly the expansion of numbers of schools reached by a given reform effort" (p. 3). Coburn argues that an adequate conceptualization of scale should go beyond merely spreading and considering all implementation challenges of reforms documented by decades of implementation research. Hence, Coburn proposes a multidimensional conceptualization of scale that not only focuses on the issue of spread but also considers other critical elements in implementing and adopting innovations, such as substantial changes in classrooms, endurance over time, and changes in authority. This conceptualization of scale consists of four interrelated dimensions: depth, sustainability, spread, and shift in reform ownership.

Conceptual framework

As part of the data analysis for this study-which will be explained in more detail in the method section-indicators were defined for each of these four dimensions. Next, the dimensions and the indicators defined for the data analysis are explained. The indicators are illustrated with mathematics education literature related to the notion of scale.

Depth. This dimension focuses on the quality and impact of reforms within individual classrooms or schools. In other words, it assesses whether there is a profound and consequential change that goes beyond surface structures and procedures and that has the potential to modify teachers' beliefs, social interaction norms, and underlying pedagogical principles.

Indicators. The depth indicators considered are utterances that refer to some kind of deep change within the classroom or the school system. For example:

• "To what extent can resources for mathematics classroom teaching, assessment, and professional development lead to changes in day-to-day classroom teaching?" (Maass et al., 2019b, p. 307) Sustainability. This dimension refers to the possibility of reforms to persist over time despite political, financial, or organizational changes. Whether the reforms can endure and continue to impact even after the initial implementation can be measured in terms of sustainability. The dimension involves the critical conditions that must be present in schools and school districts to support and sustain classroom change over time.

Indicators. The indicators associated with this dimension are utterances that refer to the sustainability of the implementation of an innovation or to the conditions that must be present in a specific context so that the changes produced by an innovation or reform are lasting and sustained. For instance:

• "What comes after the professional development program is over? How can we sustain what was gained and how can the impact be scaled up?" (Tirosh et al., 2015, p. 153) Spread. In this dimension, the focus is on spreading reforms or innovations to a more significant number of classrooms and schools. However, this dimension not only considers the possibility of spreading solely in terms of outward expansion to more and more schools; it also considers the spread of a reform or innovation within a classroom, a school, or a district.

Indicators. The spread indicators are those utterances that refer to an outward expansion of reform or innovation but also those that refer to a 'spreading within'. For instance:

• "This extends to around 40,000 teachers spread over the country's 290 municipalities and nearly 6,000 schools. The government has set aside SEK 650,000,000 (EUR 75,000,000) for the 4-year programme. This is one of several large-scale initiatives to improve mathematics teaching in Sweden" (Boesen et al., 2015, p. 129) Shift in reform ownership. This dimension concerns the transfer of control and ownership of reforms from outside organizations to teachers, schools, or districts. It assesses whether the reforms are becoming embedded within the organizational culture and whether they are being led and sustained by those they directly impact.

Indicators. The indicators of this dimension are utterances that reference the presence of structures or mechanisms favoring ongoing learning about the reform or innovation, which provide funding and support for the implementation of the reform, or that reference how internal actors take responsibility for a continued spread of the reform. An example of such an indicator is:

• "Effective scaling requires teachers to obtain ownership of the reform and ultimately transfer the new approaches into an "internal reform" in their own school" (Potari et al., 2019, p. 419)

Method

We claim that the study reported in this article is a "Delphi-inspired" study because although it follows some of the methodological procedures of the Delphi technique, it does not strictly adhere to these. In the following, we describe the Delphi technique and the procedures from this technique in more detail. We also describe how the data analysis informed by [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF] was developed.

The Delphi technique and the procedures implemented

The Delphi technique is a method used to gather expert opinions or predictions on a specific topic. Its use as a research technique is not new in mathematics education (e.g., [START_REF] Manizade | Using Delphi methodology to design assessments of teachers' pedagogical content knowledge[END_REF]. It is a systematic, interactive forecasting process that utilizes a panel of experts to provide input through several rounds of questionnaires. The purpose of the Delphi technique is to achieve a consensus or to converge towards a common understanding, based on the experts' knowledge and experience [START_REF] Barrett | What are Delphi studies?[END_REF].

The Delphi technique typically involves the following methodological procedures:

1. Identifying and selecting a panel of experts in the relevant field. 2. Sending an initial questionnaire to the experts to gather their initial opinions and predictions.

3. Compiling and summarizing the results from the first round of questionnaires. 4. Sending a second questionnaire to the experts, incorporating the results of the first round and asking them to revise their opinions in light of the results. 5. Repeating this process until a consensus is reached or there is no further change in the opinions of the experts.

The current study followed Delphi technique's procedures 1, 2, and 3.

In the case of procedure 1, based on our experience and knowledge of the educational field and its literature, we identified and selected a group of 14 educational researchers from different regions of the world who had experience in leading or investigating the large-scale implementation of educational innovations and/or reforms. All educational researchers except for two are mathematics educators. However, these two exceptions are experts in scale studies, whose concepts are widely used in implementation research, including mathematics education research. All 14 experts were invited to participate in developing the second procedure of the Delphi technique.

For procedure 2, the 14 researchers were invited via email to answer two scale-related questions. All invited researchers agreed to participate in the study. They come from Europe, Latin America, and North America. Some decided to participate anonymously, while others consented to us using their names. It was explained to the researchers that the questions posed and their answers were part of the development of a review on implementation projects that included small-scale and large-scale projects [START_REF] Ahl | A theoretical framework for analyzing factors of influence for the success or failure of implementation projects[END_REF]. However, deciding if an implementation project is small-scale or largescale is difficult, and some may fall between the two categories. With this context as background, researchers were asked the following two questions:

• How do you define a large-scale project or a project at scale? • Are "large scale" and "at scale" synonyms?

As part of procedure 3, the researchers' email responses were collected in a shared document for analysis. Although the questions posed to the researchers involved various scale-related notions, in this paper, we only report on the interpretations related to the notions of large-scale and small-scale projects.

Data analysis

The analysis of the empirical data (the experts' responses) was guided by the four dimensions that constitute the conceptualization of scale proposed by [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF]. First, indicators for each of these dimensions were defined (see the "Conceptual framework" section). Afterward, the presence of the indicators in the explicit definitions or characterizations of the concepts of "large-scale project" and "small-scale project" stated by the experts were sought out. The absence or presence of these indicators in the empirical data were considered as markers of the dimensions of scale that are considered in the interpretations of a "large-scale project" and a "small-scale project" in mathematics education research.

Results

All of the surveyed experts explicitly referred to either the notion of a large-scale project or smallscale project or both-this, of course, was motivated by the fact that they were asked to define largescale projects in one of the questions posed to them. The responses of the experts can be divided into two groups.

The first group of responses state that large-scale and small-scale projects are fuzzy concepts, difficult to define, or without a clear definition. Seven experts referred to these conceptual difficulties (e.g., "There is no absolute border between small and large scale", Expert #3; "The only thing we could say about it was that different people treat it differently", Expert #5; "There is no clear definition", Expert #11; "The difference between small-scale and large-scale projects is quite challenging", Expert #13).

The second group of responses refers to those that explicitly characterize these concepts. Nine experts provided these types of responses, which were analyzed in search of the indicators of the dimensions of the [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF] framework. The results show that all explicit references to the notions of largescale and small-scale contain spread indicators. For example, "Projects that focus on instructional improvement at scale can be quite small scale (i.e., require limited funding and involve just a few people).", Expert #1; "Of course, no one would deny that mathematics education research studies involving thousands of students or teachers are large-scale studies", Expert #6; "We consider our PD project as a large-scale project as we offer PD for facilitators who themselves conduct PD for teachers", Expert #11. However, a couple of the researchers interviewed warned us that the issue of size in large-scale projects is relative. For instance, Expert #3 states: "Above 100 teachers is not small scale any more, but what is really large might also depend on the size of the country".

Indicators of depth and sustainability were only identified in the responses of one of the interviewed mathematics education researchers. In one of his answers, the expert refers to how large-scale projects should "permeate the system", which is interpreted as a reference to a significant change in the system in which an innovation is implemented: "'Large scale' means 'big' (often costly as well) but something may be big without necessarily permeating a system.", Expert #14. When referring to large-scale projects, this researcher mentions "stability", which is interpreted as a reference to the potential of the implemented projects to remain stable over time, i.e., sustainability: "There are teaching support programs called 'math circles'. There are many such organizations-see https://mathcircles.org/. The number of organizations, and their stability, says that they are large scale", Expert #14.

Nevertheless, in four formulations of large-scale or small-scale projects, we identified features that are not included in [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF] conceptualization of scale. These features refer to differences in methodological conditions and tools and notions such as representativeness and generalization: "There is no doubt for me, that whatever are the distinctions made, studies at different scales differentiate in terms of methodological conditions and constraints, which has impact on their possible outcomes", Expert #6; "One criterion could be if you can do quantitative, statistical analysis on the group? Versus small scale if only qualitative analysis makes sense?", Expert #12; "Large scale is rather vague, I may suggest; is it all the population or large enough so that then a generalization is accepted? Actually, it is the question of generalization which is important and must be the criterion", Expert #8.

Discussion

The responses from seven out of the 14 experts indicate that it is challenging to distinguish between the notions of large-scale and small-scale projects. This fact emphasizes the relevance of developing a study that looks for conceptual clarity on these notions.

Regarding the question that guides this study, "In relation to implementation research, how are the notions of a 'large-scale project' and a 'small-scale project' interpreted in mathematics education research?", a first approximation to the answer would be that these notions are generally interpreted from a perspective of scale that emphasizes the spread dimension. The dominant perspective consists of quantitatively interpreting these notions, prioritizing the number of teachers, schools, or districts that are reached by implementing a reform or innovation. Although there is evidence of mathematics educators' interpretations that include dimensions of depth and sustainability, this perspective is not dominant.

The prevalence of this scale interpretation among mathematics educators, which tends to emphasize the spread dimension, is corroborated when contrasting their conceptualizations of scale with those of the two educational researchers (non-mathematics educators) who participated in the study. The educational researchers' conceptualizations of scale also include the dimension of spread, but the notion as a whole is portrayed as an ultimate state fully realizing all targeted scaling dimensions:

There is a difference between a scaling project and a project that is "at scale." Saying "at scale" suggests that the project has accomplished the goal of having fully realized all targeted dimensions of scaling., Expert #4

We define "spread" as the process of an innovation becoming more widespread. Similar terms include diffusion and scaling. We define "scale" as the end-state of spread-that is how something is used once it has become widespread. So, for us, scale is a goal or end-state., Expert #9

Moreover, some experts' large-scale and small-scale formulations include features not encompassed in [START_REF] Coburn | Rethinking scale: Moving beyond numbers to deep and lasting change[END_REF] conceptualization of scale, such as differences in methodological conditions and tools and the notions of representativeness and generalization. These characterizations point to the need for a conceptualization of scale that not only focuses on aspects of impact, reach, stability, and independence of the implemented project, but which also considers the differences in approaches and tools that enable the study of its effects, as well as the potential for generalization to other populations and contexts.

The study reported in this paper has continued to evolve, complementing the interpretations identified in this Delphi-inspired study with specialized literature on implementing innovations at scale in mathematics education (see Aguilar et al., 2023 for details). In particular, we have reviewed 22 articles from two special issues related to implementing and scaling innovations, published in Educational Studies in Mathematics [START_REF] Maass | Different ways to large scale implementation of innovative teaching approaches[END_REF] and ZDM -Mathematics Education [START_REF] Roesken-Winter | Effective CPD on a large scale: Examining the development of multipliers[END_REF]. The articles were read and analyzed, extracting data on scale definitions and explaining what was referred to when talking about small scale, large scale, scaling, and at scale. In the elaboration of this paper, we suggest considering factors such as contact, material, and organizational factors to differentiate between small-and large-scale implementation projects. Additionally, we introduce the concept of medium-scale implementation projects, which can bridge the gap between small-and large-scale projects. We also provide examples of scenarios in which medium-scale projects can be conceptualized.

The insights gained from the elaboration of this research can inform future implementation research in mathematics education by providing a characterization for understanding and discussing scale and scaling-and this not only in terms of quantitative measures but also regarding the qualitative features of the innovations in question.

Acknowledgment

This paper is part of grant 2022-04811 under the Swedish Research Council.