Numerical study of impeller-driven von Karman flows via a volume penalization method

The von Karman flow apparatus produces a highly turbulent flow inside a cylinder vessel driven by two counter-rotating impellers. Over more than two decades, this experiment has become a very classic turbulence tool, studied for a wide range of physical systems by many groups, with incompressible flow, compressible flow, for magnetohydrodynamics and dynamo studies inside liquid metal, for particle tracking purposes, and recently with turbulent super-fluid helium. We present a direct numerical simulation (DNS) version the von Karman flow, forced by two rotating impellers. The cylinder geometry and the rotating objects are modelled via a penalization method and implemented in a massive parallel pseudo-spectral Navier-Stokes solver. We choose a special configuration (TM28) of the impellers to be able to compare with set of water experiments well documented. But our good comparison results implied, that our numerical modelling could also be applied to many physical systems and configurations driven by the von Karman flow. The decomposition into poloidal, toroidal components and the mean velocity fields from our simulations are in agreement with experimental results. We analyzed also the flow structure close to the impeller blades and found different vortex topologies.


Introduction
A von Kármán experiment is a cylindrical vessel in which a flow is generated by the rotation of two impellers at the extremities of the vessel [1]. Studying turbulence by means of von Kármán fluid experiments has a strong tradition in last two decade. Several teams set up such experiments with different designs in incompressible [2,3] and compressible flows [4,5,6,7]. This type of experiments was also one of the first setups used to study Lagrangian statistics of turbulent flows [8,9,10,11], by tracking solid particles or bubbles. Recently, a helium super-fluid experiment has been built with the classic von Kármán configuration [12] to reach even higher Reynold number, and to study the interaction of the super and classic fluid.
In the last decade, in order to gain a better understanding of the underlying processes of magnetohydrodynamic and dynamo effect, many experimental groups have investigated experiments using liquid sodium [13,14,15,16]. A very successful experiment used the von Kármán apparatus with Sodium liquid (VKS) hosted in Cadarache which was able to reproduce dynamo action in a turbulent flow [17,18,19,20,21,22]. Before starting with sodium experiments, prototypes filled with water were set up, which are smaller than the final VKS machine by a factor of two. They compared and optimized different impeller designs to seek the highest kinematic dynamo growth rate [23,24].
In this paper we will focus on the purely hydrodynamic properties of the von Kármán flow driven by rotating impellers and keep the investigation of the magnetized dynamics, or other experimental application as a future goal.
In this work we perform direct numerical simulations of a VKS flow, thus a three-dimensional impeller-driven turbulent flow in cylindrical geometry. The geometry of rotating impellers assembled of several basic geometric objects is modeled via an immersed boundary technique (IBM) and implemented in a massive parallel pseudospectral Navier-Stokes solver. High resolution simulations allow for the development of a fully developed turbulent flow. We first compare our numerical results with the water campaigns of the Saclay group [25,26,27,28,29] to validate our numerical modelization approach. We will also study the flow dynamics in the vicinity of impellers, a region not easily accessible experimentally.

Basic equations
We consider the incompressible Navier-Stokes equation with the velocity field u(x, t), pressure p, viscosity ν. The velocity field additionally fulfills the incompressibility condition ∇ · u = 0. At the boundaries we impose a no-slip condition such that the velocity of the fluid equals the velocity of the boundary itself u| boundary = V penalized . It is zero on the fixed outer cylinder and equals the solid rotation velocity on the disks and the blade structures. The cylindrical wall and the moving impellers are implemented by a penalization method (see 2.3).

The Fourier-spectral scheme
To solve the equation system a standard pseudo-spectral method is applied using the 2 /3 rule for dealiasing and resolutions of 256 3 and 512 3 grid points. The time derivative on the left-hand side of the Navier-Stokes equation is discretized via a strongly stable third order Runge-Kutta method [30]. Incompressible turbulent flows have been intensively studied in a periodic box, a classical mathematical framework for theories [31] as well as for numerical simulations of isotropic and homogeneous turbulence [32,33]. In this geometry the pseudo-spectral numerical method is the most precise global numerical method for a fixed mesh size and the success of this method is essentially due to the efficiency of the Fast Fourier Transform. In the present work, we are using an immersed boundary method to impose no-slip boundary conditions inside the simulation domain.
In this framework we lose the spectral precision near the boundaries, but we are able to design any geometry of static or moving structure while keeping the usability of a pseudo-spectral code. The implementation of the penalization boundaries and the moving impellers increase the CPU time by 2 times. The most time-consuming part of this implementation is the recalculation of the moving boundaries in each step. For future MHD simulations we expect a factor below 2, as the penalization then comsumes less time compared to the solution of the basic equations. Those methods are sufficiently accurate to reproduce standard benchmarks [34] and benchmarks depending crucially on the boundary layer solution [35]. The used method will now be explained in more detail.

Penalization method
To simulate flows within a solid cylindrical boundary and moving impellers with complex shape, a penalization method is applied. For points inside the boundaries a "pseudo" forcing term is added to the right-hand side of the Navier-Stokes equation (1), which adjusts the velocity exactly to the prescribed value in and on the wall or the impellers. The method we used to calculate the force, which is first introduced in [36,37], is called direct forcing and allows to calculate the force directly from the Navier-Stokes equation without the necessity of further auxiliary parameters. To increase the precision of the boundary layers, we used a predictor for the pressure gradient [38]. As we deal with complex geometries the boundaries of the objects do not coincide with the rectangular grid. This makes it necessary to interpolate at the boundary, for which we take the volume fraction V b occupied by the solid object into account. Regarding the volume of each cell V c = n 3 ∆x∆y∆z, the force is weighted with the factor V b /Vc. Practically, this is performed via a Monte-Carlo method using 50 random points within each cell to calculate the volume [37]. The details of this penalization method were described and tested in another fluid context [35].
The solid rotating velocity of the impellers is simply given by the angular velocity Ω and the distance from the axis r as V boundary = Ω × r. The angular velocity of the impellers can be set independently for each impeller. We create an embedded cylinder inside the periodic box, using our penalization method, with a radius R c = 3.0. Though the periodicity is kept along the z axis to decrease Gibbs oscillations, the velocity at the end of the cylinder is close to zero due to the symmetry of the system. The height of the cylinder is 2π, giving an aspect ratio of 2π/R c . In experimental setups this height varies from 2 up to 3 cylinder radii. The interior of this cylinder represents ∼ 71% of the total computing domain.
For one set of simulations we choose a very similar configuration of the curved disk-blades to the setup called "TM28" [23] with a distance between the disks of 1.8R c , a radius of the rotating disk of R d = 0.9R c , a height of the eight blades of 0.2R c and a curvature radius of the blades of C = 0.5R c . The angle of the expelled flow at the end of the blades is given by α = arcsin ( R d /2C) ≈ 1.11976 rad ∼ 64.15 deg. The simulated straight blade configuration is close to the "TM70" and "TM80" configuration [26,28,24], with eight blades per disk. A difference is that our disk radius is R d = 0.9R c instead of (R T M 70 = 0.75R c , R T M 80 = 0.95R c ). We call this configuration straight configuration. Here, C = ±∞ and α = 0.
In the following we analyze these two blade geometries (TM28 curved and straight blades). The disks are always counter-rotating with the same rotation rate. Depending on the rotation direction we have two different blade curvatures (called + or −) (see figure 1).

Figure 1: Simulated impellers with 8 curved blades (left) and 8 straight blades (right):
A positive curvature, denoted (+), means that the convex side of the blade points in the turning direction of the impeller, while for a negative value, denoted (-), it is the concave side.

Non dimensional numbers
We have chosen six simulations to highlight our results. We vary the viscosity ν controlling the dissipative term, the rotation rate Ω of the impellers, and the curvature of the blades. Of course other parameters such as the geometry of the blade, the number of blades, the difference in the rotation rate of the disks might change the topology of the flow and the quantitative results. From our numerical data we compute a set of meaningful non-dimensional numbers (see table 1). The Reynolds number accessible with direct numerical simulations with 512 3 grid points is as usual a few orders of magnitude lower than experimental Reynolds numbers. Nevertheless, at this Reynolds number the flow is turbulent. Our measured fluctuation level δ = u 2 / u 2 defined in [39] is similar to that obtained from the water experiments for synchronized rotating disks and with an annulus deviator configuration in the "TM60" and "TM73" setup. With asynchronous disk rotations or without annulus, this level of fluctuation could be higher (above 2.0) [39].
The efficiency E f = Umax /ΩR d of the impellers, which measures how much energy is injected into the fluid, is defined as the maximum velocity of the fluid in the bulk induced by the impeller motion. The variation of the E f as a function of the expelled flow angle α for different experimental configurations is confronted with our numerics in figure  4.b. With a Reynolds number three orders of magnitude lower than the experiments, we nevertheless have a good agreement with the experimental efficiency. The efficiency number decreases with the same slope found in the water experiments. We noticed that our straight blade efficiency is almost identical with the "TM70" configuration efficiency. For the positive curved blades our numerical data is closer to the "TM60" configuration (16 blades) than to the "TM28" (8 blades). The ratio of the root mean square velocity and the maximum velocity of the impellers (U rms /V max ) is pretty close to the TM28 experiment for our higher rotation rate simulation (1.6Ω = 2.4). Another non-dimensional number is the ratio of rotation period of the impellers T Ω = 1/f = 2π Ω and the large eddy turn over time T nl . Our simulations are in the same disk rotation regime as the experiments (see the last line of table 1) showing that there is a bit more than two large eddy turn over times during one turn of the impellers.   [23,26,39] and specially the configuration "TM28". We collected and defined several quantities or non-dimensional numbers: ν the kinematic viscosity of water or in our simulations, Ω the rotation rate, U rms = 2E(t), U max is the (spatial) maximum of the (time-averaged) mean velocity field in the bulk in the range −0.8R c < z < 0.8R c , where z=0 is the center of the cylinder. L = 2π E(k) E(k)/k is the integral length scale computed with E(K), the isotropic spectral density of the kinetic energy, and T nl = L/U rms is the eddy turn over time. We used the experimental Reynolds number definition of VKE-VKS experiments R exp = ΩR d Rc /ν. To compare with numerical works, we define another Reynolds number R num = UrmsL /ν. The efficiency E f = Umax /ΩRc of the impellers states how much energy is injected into the fluid. δ = u 2 / u 2 is the fluctuation level defined in [39]. We present also the ratios U rms /V max and T Ω /T nl = 2πUrms ΩL . Our simulations are during more than 20T Ω (turns of the impellers), which is the duration of the time average computation. of the von Kármán flow. For comparison, a snapshot of the enstrophy (figure 2.c) shows interacting vortex filaments in the bulk region which is characteristic to a turbulent flow. The vorticity is produced and thus very high near the blades. This observation will be analyzed in detail in section 4.

Poloidal and toroidal components
For further analysis of the simulated flows we performed a decomposition into poloidal and toroidal components of the form with unique potential functions Ψ and Φ and the unit vector in z direction e z . To be precise, in the periodic box, we normally need to add a component F (z) depending only on z (see [40] and (course 2, C. A Jones) [41]). Although our penalized cylinder is periodic along the z axis there is no mean flow crossing the box as the gap between the cylinder and the disks is small. We checked that F is zero (up to the numerical digit precision). This decomposition has also been performed with the experimental data with the TM28 impellers, which allows a comparison of experimental and simulated data. Supposing axial symmetry around the z axis, we can compute the poloidal and toroidal two dimensional fields ∇ × [Ψ(r, z)e z ] respectively ∇ × ∇ × [Φ(r, z)e z ] (see figure 3) and compare them to figure 3 of [23]. We have visually a good agreement, the poloidal flow consists of two large scale vortices in the regions z < 0 and z > 0, where the toroidal components respectively point in opposite directions. Impellers with different curvature (negative curvature and straight blades) produce the same kind of flow structure. From the images it is difficult to distinguish between the different blades configurations. We therefore present in table 2 the mean and maximum velocity of the poloidal/toroidal components, all of them rescaled by their respective maximum velocity of the impellers V max = ΩR d . Those values can be compared to the experimental ones [23]. Our velocities of the (+) configuration are 10% − 20% lower than the experimental data of "TM28". This could be explained by the fact that our efficiency coefficient is lower than the "TM28" configuration (see figure 4 right), implying a smaller velocity in the central region of the vessel. We also compare the ratio of the poloidal over the toroidal velocity versus the expulsion angle of the blades α with different water experiment results [24,28,23]. This ratio was used to seek the dynamo onset as a control parameter. Like the efficiency the ratio Γ mean for the (+) configuration is closer to the "TM60" than the "TM28" setup (figure 4 left). However we stress that our ratios have also a positive    Table 2: Quantities from experiments and simulations : the maximum and the mean of the poloidal and toroidal velocity and the respective ratio Γ max = u pol,max /utor,max and Γ mean = u pol,mean /utor,mean. All the velocities are normalized by the maximum velocity of the impellers V max = ΩR d . A quantification of the poloidal and toroidal components is done by extracting the maximum and mean values in the bulk, the region −0.8R c < z < 0.8R c . In addition, the average torque T on the impellers has been computed for the three simulations with lower resolution and normalized to the value obtained in the run (+).
slope. Even if there is not a perfect agreement our ratios are quite close to the different experimental measurements.

Impact of viscosity and rotation speed
In order to study impact of the viscosity and the disk rotation rate on the mean flow structure we performed additional simulations with curved impeller blades and positive turning direction. In the first of these runs we increase the resolution to 512 3 grid points (Run (+) 512 in   Figure 4: From different water experiment setups given in data tables (TM7x [24], TM8x [28], TM60 and TM28 [23]) and from our numerical results, we plot (left) the ratio of the poloidal and the toroidal mean velocity (Γ mean ) and (right) the efficiency E f = Umax /ΩRc of the impellers both versus the flow expulsion angle α at the end of the blades.
In the second run (run (+) 1.6Ω in tables 1&2), the angular velocity of the impellers is increased by a factor of 1.6, while all other parameters as well as the resolution are unchanged. The Reynolds number then increases to Re = 3888. The listed values remain constant, only the ratio of the rotation period over the eddy turn over time slightly decreases. In the last run the viscosity is lowered by factor 2 (run (+) ν/2 in tables 1&2) and the resolution is increased to 512 3 grid points. In this case the Reynolds number increases up to Re = 4860. Regarding all quantities obtained in the simulations, there is evidently no major influence of rotation speed and viscosity on the mean flow profile. Only the fluctuation rate δ slightly increases. The mean flow quantities which we present do not change with the rotation rate or the viscosity. This implies that the corresponding simulations are already in an asymptotic regime, where the Reynolds number has only little effect on the large scale structure of the flow. Indeed, the range of Reynolds number numerically acheived in our work (2500 − 4800), is at the edge the inertial regime of water experimentals results ( [29] see their Fig. 5 and Fig. 7). In this experimental campain, the Reynolds number has been increased progressivelly by growing the rotation frequency of the disks, to explore different phases : viscous, transition and inertial regimes.

Vortex generation behind blades
Besides the general flow structure imposed by the impellers, our main interest lies in the near blade flow in the frame of reference of the rotating impeller. Here, we present and analyze the structures obtained in the three different configurations ((+), (−), straight) with n = 256. We construct the mean flow in the rotating frame of reference by first averaging the flow each time the blades pass the same positions, which is every eighth of a turning period. As we expect to find the same mean flow at each blade we also take the spatial average for rotations around 90 • . As the running time of our simulation corresponds to 20 impeller turns we were able to average over 640 realizations. The mean flow in the rotating frame is obtained from this average by subtracting a solid rotation u * = u − Ω × r.
To get clear pictures of the flow we consider two different planes: one parallel to the disk plane, with a view from the top and the other perpendicular to the disk (the position of the perpendicular plane is indicated by a red line in the top view (left column of figure  5)). In the top view, we show the velocity projected onto the plane while its amplitude is represented by a color map. This perspective shows how the flow is sucked in from above to the center, moving between the blades and how it is finally expelled from the disk outwards. Just by looking at the left column of figure 5, we can easily distinguish between the different configurations, specially the (+) and the (−), where the expelled flow in the rotating frame has the same rotational direction as the disk rotation. The horizontal component represents the major part of the magnitude of the full velocity along the blade. There is clearly an acceleration from the center to the expulsion area. The most rapid velocity is generally along the pushing blade wall. A small sucking area is located just behind at the end of the blade. In the perpendicular plane we decide to visualize the velocity by streamlines to highlight the topology of perpendicular flows. Note that the projected velocity is not solenoidal, thus streamlines may have an end point. In all three simulations, (right column of figure  5) the streamline plots show vortex rolls emerging directly behind the moving blade, as vortices are ripped off at the blade's edge. Those vortices appear to take most of the space between the blades. Clearly, the negative configuration (-) has a different topology than the two others. This observation agrees with the horizontal velocity in the top views. The cut along the red line allows the visualization of the mean flow vortex at different radii for different cells. In the positive configuration (+), it can be deduced that a cone vortex is produced along the pushing blade in each inter-blade cell. Those vortices are also present in the straight configuration which is, however, less clear for the negative configuration (−).

Numerical and experimental comparisons
By means of direct numerical simulations using a penalization technique we are able to reproduce the large scale structure of experimental von Kármán flows produced by moving impellers. Our good agreement with the experimental results could be explained by the fact that the mean flow geometry and other global quantities are converging  (straight)). The velocity is averaged in time, projected on the plane, and computed in the rotating frame (u * = u − Ω × r). In all images, the color map represents the magnitude of the velocity. The transverse view planes are perpendicular to the red line shown on the top view planes. Not all of the perpendicular planes associated with the red lines is plotted (only 70% of the red line of the left side). The position of the blades helps to relate the top and side views. Note that this projection of the velocity in the considered plane is not solenoidal, thus streamlines could end at the boundary, where the projected velocity tends to zero. Streamlines of the 3D flow do not enter in the solid object, but instead slide along the boundary. This behaviour is not captured by the 2D projection. rapidly even at low Reynolds number. The natural next step could be to seek the small scale properties of the von Kármán flow, like the studies on filaments [3,42], effect of large scales on small scales [7] or the energy injection [43,44]. Of course, direct numerical studies with confined flows in the full von Kármán geometry remain a challenge asking for a big increase of spatial resolution. Our modelization using the penalization could be easily used to explore different experimental setups -at least for the large scale properties.

Vortex and Dynamos
We found a characteristic outwards spiraling vortex between the blades. Note that the influence of the vortex generated around the blades is suspected to have a strong impact on the dynamo effect [45]. Some numerical results assuming a dynamo mechanism concentrated around the disk-blade structure [46,47] have found that the magnetic mode has a dipole structure according to the experimental results. The geometry of the experimental magnetic mode cannot be explained by a mean flow dynamo only. Recent numerical studies using FLUENT (k − ǫ)-RANS simulations [48] computed the α-tensor produced by the vortex dynamics, given a switch between α 2 and α − Ω dynamo types. Despite these vortex dynamics, without soft iron impellers the dynamo threshold was not achieved. The material of the impellers plays a crucial role for the efficiency of the dynamo action [49,50,51]. The material properties of sodium make in-situ diagnostics very difficult. Direct numerical simulations provide a unique tool to assess spatially and temporally resolved variables. With our dynamics of the blade vortex and correct treatment of the magnetic properties of the solid impellers, all the physical ingredients should be present to address those different questions, to seek the interaction between the conducting fluid and the ferromagnetic impellers. We could also check the different dynamo onset prediction or measurement of the different configurations produce by the variation of the blade geometries or material properties [51].

Long term dynamics
When the water experiments are running during a long time, the von Kármán mean flow can change to different topology solutions, breaking symmetries [27,29,52]. In the experiment, the typical time scale to record this multi-stability is around 10 5 to 10 6 hydrodynamic large eddy turnover times. In our present simulations, we are completely out of reach to record such dynamics: we computed 20 disk turns, which represent around 45 eddy turnover times only. It will thus be very interesting, even if less popular than the highest resolution limit, to perform long numerical runs with moderate grid points to reach and study the long time physical behaviors or to improve statistical data.