Merci À Wenwu

Liu Guodong

Pengyu, Zhihui Nan Tao

Fang, Hepeng, Yimei, Rongrong, Lijun Jean Marc

Lin Qi

Keywords: Écoulements tournants, Surface libre, Simulation numérique, Expérimentation, Polygons tournants, Instabilités. Experimental deformable free surface Rotating flows, Free surface, Numerical simulation, Experiments, Rotating polygons, Instabilities

Cette thèse était effectuée au sein de deux laboratoires: LIMSI et FAST avec le financement de

List of Notations

x h(r, t) measurements of surface deformation by Gocator along a fixed diameter h s (r, t) measurements by Gocator along a fixed diameter in established regime H s (r, f ) time Fourier transform of h s (r, t) h (r, θ) instantaneous surface height map δh (r, θ) azimuthal deviation from azimuthally-averaged surface profile Magnitude (or Mag) magnitude of one-side Fourier transform of signal, normalised by the length of signal N PSD calculated practically by Mag 2 /(2∆ω) ∆ω the increment of pulsation given by 2πf s /N , with f s the sampling frequency

Chapter 1 Introduction

Background and motivations

Rotating flows sometimes exhibit surprising instabilities which can give birth to interesting patterns but the physical mechanisms to destabilize the flow are sometime difficult to understand. Take the example of a simple experiment, in which a fixed cylindrical recipient is filled with a liquid, say water, and a disk placed at the bottom is driven into rotation (as shown in figure 1.1a). One may guess that the surface will be deformed in an axisymmetric way. Actually, when the rotating speed exceeds a certain value, the axisymmetry will break and a stable rotating polygonal occurs, as shown in figure 1. 1(b). This spectacular rotating polygonal pattern reminds us of the long life hexagon at the north pole of Saturn, which was first discovered by the Voyager spacecraft [Godfrey, 1988] and then photographed by Cassini after 2009. A recent work [START_REF] Fletcher | A hexagon in saturn's northern stratosphere surrounding the emerging summertime polar vortex[END_REF] shows that this structure may extend up to the stratosphere. Since observations and measurements are difficult in such a cosmic environment, some astrophysical phenomena remain unexplained. Laboratory tests on simplified flows are thus always interesting to help understand the mechanisms driving such phenomena, which are often a combination of several physical processes. The rotating flow is also of great importance on Earth. For example, it is very common in atmosphere and ocean motion. Mangkhut. On this photograph, one can clearly distinguish the eye of the cyclone which is a warm descending zone within the central depression area bordered by a strong ascending whirling flow. Such vortex flows can also be observed in water motion. For example, when two flows with different velocities meet, one or several whirlpools may be created. As shown in figure 1.3(b) on the Naruto whirlpool, the rotating flow is associated to a depression area at the center so that the water surface can become strongly deformed. The interaction between the fluid surface and the flow motion can be complex, and a better understanding of these large-scale rotating flows is most useful in the fields of meteorology and geophysics. Besides, we can also find numerous applications of rotating flows in industry. To list several, a very common application is in turbomachinery. Turbomachines such as compressor and turbine transfer energy between a rotor and a fluid and they are widely used in transport and power generation. Even in our daily life, we are also making use of rotating flows. As a simple example, tea and coffee drinkers may all experience the stirring of their beverage with a spoon in order to accelerate sugar dissolution. To this rotation are associated other motions, for instance from bottom to top at the center and vice-versa along the lateral side.

Useful and complex as it is, the rotating flow arouses a great interest among the scientists in the field of fluid mechanics. Recently, experimental techniques and numerical solvers have been greatly developed and, together with theoretical modeling, are likely to boost the investigations on rotating flows in any configuration. In the following, we first outline the most frequently encountered academic configurations used in the laboratory, and then describe the one used to generate the rotating polygons mentioned earlier.

Rotating flow in the laboratory

A rotating flow can be easily created by simply rotating a plane disk in a still fluid. This configuration was studied by [von Kármán, 1921] who derived a self-similar solution for this three-dimensional boundary layer. This model still holds for configurations where the fluid itself rigidly rotates far from the disk [Batchelor, 1951, Zandbergen andDijkstra, 1987].

These pioneering works have shown great interests of rotating flow induced by a disk for understanding the three dimensional boundary layers, which are common problems in geophysics and aeronautics engineering. For instance, a comprehensive understanding of the instabilities growing in the von Kármán boundary layer has developed, in order to conduct flow control [Pier, 2003, Pier, 2007, Imayama et al., 2014].

Rotating flow in confined cavity

When two rotating disks are enclosed by a fixed cylindrical wall, some additional effects come up due to the confined geometry * . According to the boundary conditions associated to the rotation speeds of the top and/or bottom disk, configurations can be generally classified into the three following categories. (a) Rotor-stator configuration: one disk is rotating while the other does not move. This configuration was investigated experimentally by [START_REF] Gauthier | Axisymmetric propagating vortices in the flow between a stationary and a rotating disk enclosed by a cylinder[END_REF], Schouveiler et al., 2001], and numerically by [START_REF] Serre | Annular and spiral patterns in flows between rotating and stationary discs[END_REF]. It displays instability patterns as axisymmetric vortices and positive spirals, due to the boundary layer instabilities. A review about this kind of rotating flow was given by [START_REF] Launder | Laminar, transitional, and turbulent flows in rotor-stator cavities[END_REF]. This flow is controlled by two parameters: the container height-to-radius ratio and the Reynolds number based on the cylinder height or radius [START_REF] Marques | Mode interactions in an enclosed swirling flow: a double hopf bifurcation between azimuthal wavenumbers 0 and 2[END_REF], Lopez, 2006].

(b) Co-rotating disk configuration: the two disks rotate in the same direction at independent speeds. As a result, the aspect ratio, the ratio of the two angular speeds and two Reynolds numbers based on the disk radii and angular speeds are usually the characteristic non-dimensional parameters. [START_REF] Gauthier | Instabilities in the flow between co-and counter-rotating disks[END_REF] found that the co-rotation case is very similar to the rotor-stator case for a low height-to-radius aspect ratio. [Lopez and Marques, 2004] specifically investigated the case in which the two disks rotate with the same speeds. The authors found two distinct instability mechanism depending on the aspect ratio.

(c) Counter-rotation: the two disks rotate in opposite directions at independent angular speeds. [START_REF] Gauthier | Instabilities in the flow between co-and counter-rotating disks[END_REF] observed experimentally specific negative spirals compared to the co-rotating case at the same low aspect ratio. Later, [START_REF] Moisy | Experimental and numerical study of the shear layer instability between two counterrotating disks[END_REF]) investigated numerically and experimentally the shear-layer instability in a much larger range of aspect and angular speed ratios. The authors confirmed and explained the previously observed negative spirals. [START_REF] Nore | The 1:2 mode interaction in exactly counter-rotating von Kármán swirling flow[END_REF], Nore et al., 2004, Nore et al., 2005] investigated the flow instabilities for a height-to-radius ratio of 2, specifically for the exact counter-rotation case.

Bottom-driven flow with free surface

To create a rotating flow with a free surface, we just remove the top disk in figure 1.4 and we get the configuration presented in figure 1.1(a). Before discussing the bottom-driven configuration, we recall here the classical case of Newton's bucket.

Newton's bucket A cylindrical container (a bucket) partially filled with water is suddenly driven into rotation around its axis by a motor (twist chord) to a constant speed. As already known, the surface will finally adopt a parabolic shape that depends on the rotating speed and the initial height -Froude number and height-to-radius ratio when expressed in a nondimensional scale. After a certain time, the spin-up time, the flow reaches a steady state where both the liquid and the container rotate at the same speed, like a solid body. This steady flow is easy to understand. Nevertheless, the spin-up process is more complex. It will be briefly reviewed later in section 1.5.

We now arrive to the configuration (shown in figure 1.1a) where we generate a rotating flow with a deformable free surface by the only rotation of a disk placed at the bottom of the container. We will concentrate on this configuration in the present thesis.

Configuration and parameters

The rotating flow with free surface has drawn much investigation as it gives rise to numerous interesting phenomena and instabilities, a complexity due to the fact that rotation, boundary layers and free surface are involved. The experimental set-up is represented schematically with more details in figure 1.5. As shown in this sketch, a cylindrical tank with inner radius R is filled with a Newtonian fluid of kinematic viscosity ν and density ρ l up to an initial height H. The disk is accelerated from rest to a state with constant speed Ω, which can also be expressed as a frequency f d = Ω/(2π). Under the effect of centrifugal forces, gravity and surface tension, the free surface deforms into a shape represented by equation z = h (r, θ, t). The fluid viscosity, the initial height and the disk angular speed are the three main control experimental parameters that we could vary. In order to characterize the rotating flow created in this configuration, the following non-dimensional parameters are defined:

G = H R , Re = R 2 Ω ν , F r = RΩ 2 g , W e = ρ l Ω 2 R 3 σ , (1.3.1)
which are respectively the aspect ratio, the Reynolds number, the Froude number and the Weber number. For a given tank, the aspect ratio is a direct indicator of the filling level.

The Reynolds number characterizes the flow regime (viscous-dominated/inertia-dominated, which gives an indication on the laminar/turbulent nature of the flow). An alternative definition based on H is sometimes used for shallow-water flows. As the surface deforms mainly under the effect of gravity and centrifugal forces, we propose to use the Froude number to characterize the surface deformation as it measures the ratio between centrifugal and gravity accelerations. Finally, the Weber number evaluates the importance of surface tension. In fact, noticeable effects of surface tension on surface deformation are always restricted to the meniscus region at the triple line along the side wall, a region with small radial extent with respect to that of the container (typically R ∼ 10 cm).

At low Reynolds number, the flow remains axisymmetric. As it is increased, two kinds of symmetry breaking phenomena may be observed depending of the parameters [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF], Vatistas, 1990, Jansson et al., 2007]. They both give rise to polygonal patterns travelling azimuthally at an angular speed smaller than that of the disk. The first one is reported to occur for rather small disk speeds, so that the free surface remains almost flat. These polygonal patterns can thus be detected through visualization of the bulk flow, as shown in figure 1.6. Hereafter this instability is named as bulk instability (BI). On the contrary, as shown in figure 1.7, the second kind of patterns occurs for large disk angular speeds, when the free surface is strongly deformed. Sometimes the speed is so large that the disk center becomes dry (dewetted) † . In these cases, travelling polygonal pattern can be directly observed and measured at the free surface. Such rotating wave is called rotating polygon (RP) in the thesis.

The ultimate objective of the present PhD is to describe and interpret the development of rotating polygons. In order to draw a global picture of symmetry-breaking developing in the bottom-driven free-surface flow, a state of the art on bulk instability, rotating polygons, as well as switching and sloshing, two phenomena also reported in the literature, is given in the following.

Bulk instability

When the disk rotates at low speeds, the centrifugal force and the radial pressure gradient are not strong enough to cause a significant deviation of the free surface from horizontal. The Froude number remains very small F r 1 ‡ . By increasing the speed for a given initial height, the flow axisymmetry may break and rotating waves may grow through instability and saturate. Yet, the surface remains almost flat. This fact calls for the use of flow tracers, such as Kalliroscope mica flakes, ink or fluorescein, so that these three-dimensional structures can be made visible. Figure 1.6 gives an example of the large vortices created by this instability, as observed by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] in experiments with water and Kalliroscope. In their study, patterns with a number of vortices varying from 2 to 8 were observed for shallow layers at G ≤ 0.107. The experiments also showed that this number decreases when the rotating speed is increased until the flow axisymmetry is recovered. [START_REF] Hirsa | Symmetry breaking to a rotating wave in a lid-driven cylinder with a free surface: Experimental observation[END_REF] observed experimentally a similar pattern with wave number 4 for a much higher aspect ratio G = 2 in a small set-up. This wave number 4 was successfully reproduced in [Lopez et al., 2004] by numerical simulations with flat stress-free surface model for the air/water interface. However the same model fails to predict the wavenumber 3 for an aspect ratio G = 0.25 [START_REF] Miraghaie | Flow induced patterning at the air-water interface[END_REF]). This is explained by the fact that the imperfect free surface in experiments raised questions for the flat stress-free model to capture odd modes. Considering the instability for deep systems is related to a near-wall jet which is not sensible to the surface imperfection, the numerical and the experimental results are consistent. For this small aspect ratio, they considered that the shear-layer between a core in †. The dry center region is normally circular. In fact, the polygons shown in figure 1.1(b) and figure 1.7 have all circular dry center regions. There is a thin fluid layer between the border of the dry disk and the rotating polygonal contour.

‡. The Froude number will be always calculated using our definition (1.3.1) throughout the thesis, unless stated differently. solid body rotation and an outer meridional circulation is unstable to azimuthal disturbances and results in the rotating wave. This scenario of shear-layer instability is also shared by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. At the same time, [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] have determined experimentally the critical Reynolds number for each mode observed in their range of aspect ratios.

This instability was also studied numerically by [Kahouadji, 2011] during his thesis. Linear instability analysis for axisymmetric base flow with flat surface have been performed for a large range of aspect ratios: G ∈ [0.0280.121]. The numerically determined critical Reynolds numbers were found consistent with those of [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] for aspect ratios G > 0.07 while the numerical results underestimated the experimental values for smaller aspect ratios. Disk vibration and surface imperfection seem to be possible causes for the discrepancy. This is the ongoing work in the context of Antoine Faugaret's thesis at LIMSI.

Rotating polygons

Characteristics When the disk rotates at larger angular speed, the Froude number may become O(1). On strongly deformed free surfaces, rotating polygonal structures may occur (see figure 1.7). Rotating polygons were first observed by [Vatistas, 1990] in a water experiment using a set-up close to the one of figure 1.5 but with a gap present between the disk and the inner wall. Stable polygons varying from triangle to hexagon were generated by increasing the rotation rate of the disk for aspect ratios around G = 0.5. These patterns rotate in the same direction as the disk but with a smaller angular velocity. [Vatistas, 1990] announced also a periodic "sloshing" between a state of elliptic cross section and a pair of vortices for the surface shape for two large aspect ratios G = 1.2 and 2.3, which seems to resemble the switching phenomenon observed by [START_REF] Suzuki | Surface switching of rotating fluid in a cylinder[END_REF] (see below). The patterns were proven to be exceptionally robust as they reappear after they have been perturbed or even destroyed by inserting an obstruction inside the flow. Besides, [Vatistas, 1990] presented phase diagrams for two specific aspect ratios, illustrating that the number of polygonal corners increases with the rotating speed. This is in sharp contrast with the bulk instability discussed previously. The range of rotating speeds in which a polygon exists decreases with the number of corners and the initial liquid height delays the emergence of a given pattern.

As an exploratory investigation of this fascinating rotating polygons, [Vatistas, 1990] gives a rather complete description of the characteristics. Similar experiments were later conducted by many other studies [Jansson et al., 2006, Vatistas et al., 2008, Bergmann et al., 2011, Bach et al., 2014, Iga et al., 2014]. The experimental set-ups used in these works are generally of the same type, as shown in figure 1.5, but may differ in sizes. A summary of the set-ups is given in [Mougel, 2014]. Despite these differences, the above characteristics of the rotating polygons were verified in a qualitative way.

(a) (b)

Figure 1.8 -Illustration of switching phenomenon extracted from [START_REF] Tasaka | Surface switching statistics of rotating fluid: Disk-rim gap effects[END_REF] observed in a small experimental set-up (R = 4.2cm) for an aspect ratio G = 0.95, F r ∼ 30.

Sloshing phenomenon observed by [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] in an experimental set-up with cylinder radius of R = 168 mm (R d = 150 mm) in a configuration G = 0.476, F r = 3.4.

Switching and sloshing Beside the steady rotating polygons, a non periodical switching between a symmetric state and an ellipsoidal non-axisymmetric state was also observed wit strongly deformed surface ( [START_REF] Suzuki | Surface switching of rotating fluid in a cylinder[END_REF]). The switching is shown in figure 1.8(a). However, the switching is observed in configurations with high aspect ratio (around G = 1) and a narrow range of disk speed, namely f d ∈ [11.5, 13.9] rad/s (Froude number F r ∈ [22,33]). Like the rotating polygons, this rotating wave also has a large amplitude and it seems to be related to the laminar-turbulent transition, as found by analyzing the intensity of flow velocity perturbations and the dynamics of surface deformation ( [START_REF] Suzuki | Surface switching of rotating fluid in a cylinder[END_REF], Tasaka and Iima, 2009, Iima and Tasaka, 2016, Tasaka and Iima, 2017]). When using a set-up with a disk radius smaller than the cylinder inner radius, [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] reported a sloshing phenomenon involving an alternation between an almost axisymmetric state and another sate with a large-amplitude azimuthal wave reaching the side wall, as shown in figure 1.8(b). The range of angular speeds over which the sloshing takes place is also narrow (f d ∈ [2, 2.58] rad/s), leading to a Froude number range F r ∈ [2.7, 4.5].

To conclude with the switching and sloshing phenomena, they were all observed in large Reynolds number regimes (Re ∼ 10 5 ). However, their occurrence is restricted to a narrow range of Ω, thus of Froude number, as is shown in the phase diagram of figure 1.9. Both of them present an unsteady surface evolution, which is different from the case of steadily rotating polygons. It is worth noting that the sloshing occurs at Froude numbers of unit order, while for switching, F r is one order of magnitude larger.

Phase diagram

We mention here the phase diagrams established by [Jansson et al., 2006] and [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] using their respective experimental set-ups in the parameter plane of the initial height and the rotating speed to get a direct look at the occurrence and transition of these rotating waves, reproduced in figure 1.9. It is worth noting that the vertical coordinate in figure 1.9(a) is 'ω/(g/R) 1/2 ', the square root of the Froude number defined in equation (1.3.1). This alternative definition of Froude number is much used by [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF], Mougel, 2014]. As a reminder, the disk radii R d are smaller than that of the inner cylinder wall for the experimental set-ups used in [Vatistas, 1990, Vatistas et al., 2008] (R d = 126 mm, R = 142 mm) and in [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] (R d = 150 mm, R = 168 mm). In fact, these set-ups introduced another non dimensional parameter R d /R which necessarily influence the formation of rotating polygons. In our set-ups, the gap between the disk and the cylinder inner wall has been reduced to a negligible level in order to get rid of this additional parameter.

From these diagrams, we observe that the rotating polygons are present for G > 0.15 and F r > 0.25 (the lowest limits). These ranges, especially for the Froude number are completely different from the range of bulk instability waves observed by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]. However, we still do not know if a bulk instability pattern, like the mode 3 observed at moderate aspect ratio G = 0.25 (and F r ≈ 0, Re = 2000) investigated by [START_REF] Miraghaie | Flow induced patterning at the air-water interface[END_REF], can resist somehow to the surface deformation when increasing Ω. This is a question that we will address in the present thesis.

Except the dashed lines (discussed later) shown in figure 1.9(a), the other observations in figure 1.9(a,b) were all conducted using water as a working fluid. Amongst all the set-ups to develop rotating polygons (as listed table 2.1 in [Mougel, 2014]), the smallest disk rotating speed giving rise to a rotating polygon is around f d = 1 Hz. Taking a typical cylinder radius of 14 mm, we can evaluate the critical Reynolds number for the first bifurcation: Re = 1.23 × 10 5 . Clearly, the rotating polygons with water are mostly observed at large [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] in the plane ( √ F r, G). (b) Summary done by [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] for the previous experimental studies, represented in the plane ( √ F r, G). In this latter graph, we added the points relative to mode 3 and axisymmetric states discussed in [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] for comparison.

Reynolds numbers O(10 5 ). The high Reynolds numbers prevailing in these experiments precludes any realistic numerical simulation of such configurations. Ethylene glycol with its viscosity 16 times larger than that of water was also tested in [Jansson et al., 2006] and experimental results correspond to the dashed lines in figure 1.9(b). Actually, only polygons with wave number m ≤ 3 were found and the transition lines are represented by the dashed lines in figure 1.9(a). This liquid gives rise to Reynolds numbers one order of magnitude smaller than that of water experiments. It seems that the effect of viscosity is small on the transitions from wavenumber m to m + 1. However, viscosity plays a damping effect on higher modes, which has been confirmed in a quite recent work [START_REF] Abderrahmane | The effect of viscosity on the rotating waves and polygonal patterns within a hollow vortex core[END_REF]. The increase of viscosity seems a good idea to numerically investigate such a transition. [Jansson et al., 2006] used two set-ups of different size for their experiments: R 1 = 145 mm (the original reported radius 131mm was corrected in [START_REF] Bach | From Newton's bucket to rotating polygons: experiments on surface instabilities in swirling flows[END_REF] to 145mm) and R 2 = 194 mm. The original phase diagram presented in [Jansson et al., 2006] in the dimensional plane (H, f d ) shows similarities in terms of transition lines between the two set-ups but large differences are present in the non-dimensional phase diagram shown in figure 1.9(b). We have also included the two triangular patterns and two axisymmetric configurations discussed in [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] into figure 1.9(b). Actually these four cases are consistent with the findings of [Jansson et al., 2006] (they used different set-ups but the same disk radius R = 194 mm).

The fact that the emergence of rotating polygons depends strongly on the experimental protocol and the size of the experimental set-up level makes it difficult to obtain an agreement between all studies at a quantitative level, even these studies agree qualitatively. We can also notice that the three-dimensional surface shape is seldom measured precisely and analysed quantitatively. A vertical laser slice was proposed in [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] to reconstruct 3D shapes but the authors reported that this method was not accurate enough. As the present instability induces large amplitude surface waves, time and space resolved surface measurements should allow us to calculate the growth rate and the phase velocity, as was done in [Cobelli et al., 2011a, Cobelli et al., 2011b, Paquier et al., 2015] who investigated other kinds of surface waves. A part of the present work is devoted to the accurate time and space resolved measurents of surface deformations.

After this panorama on rotating polygons, we now focus on physical interpretations and on the different modellings in the literature.

From description to understanding

Kelvin's equilibria Considering the fact that only polygons with wave number smaller than 6 are observed in bottom-driven experiments, [Vatistas, 1990, Vatistas et al., 2008] conjecture that stable polygons are related to Kelvin's equilibrium of a ring of m point vortices (m ≤ 6). The vortices, observed in some cases, are considered by [Jansson et al., 2006] to be possibly caused by the strong shear layer due to the stationary wall. They speculate that there is a phase locking between the rotating pattern and the disk rotation in their experiments. [START_REF] Bach | From Newton's bucket to rotating polygons: experiments on surface instabilities in swirling flows[END_REF] showed later that there is no phase locking. As for the emergence of instability patterns, in a total misunderstanding between two teams, [Jansson et al., 2006] states that the surface shape can spontaneously break the axisymmetry, while [START_REF] Vatistas | Experimental confirmation of Kelvin's equilibria[END_REF], Abderrahmane et al., 2009, Abderrahmane et al., 2011] showed that the transition between two types of polygons occurs progressively. For instance, the intermediate states between a clear triangle and a clear square involve a superposition of the two patterns.

To verify the conjecture of [START_REF] Vatistas | Experimental confirmation of Kelvin's equilibria[END_REF], [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] calculated the velocities of the vortices from PIV surface velocity measurements and expressed some doubts whether the vortices are strong enough to induce the observed angular speed of the entire pattern [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF], Bergmann et al., 2012]. Quickly, the group of Vatistas published their new evidence [START_REF] Abderrahmane | A note on relative equilibria in a rotating shallow water layer[END_REF] to validate the conjecture of point vortices. These authors concluded that their own results and also those of [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF], Bergmann et al., 2012] are indeed consistent with predicted results by the point vortex theory. Nevertheless, the point vortex approach is not predictive in terms of pattern formation. To conclude with, all the discussions on the origin of rotating polygons do not converge to a uniform one.

Base flow: theoretical models

To understand the first bifurcation to a rotating polygon, an investigation of the basic state followed by a linear stability analysis is a classical and necessary approach. In their paper, [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] outlined the importance of the axisymmetric state and its experimental characterization as a starting point for linear instability studies. [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] measured the flow velocities on the free surface for two axisymmetric states: one with a dry region around the disk center (hereafter denoted as AD for convenience, see the upper red circle in figure 1.9b) and another one with the entire disk covered by liquid (hereafter denoted as AW, see the lower red circle in figure 1.9b). As shown in figure 1.9(b), AW is located before (in terms of F r) the first bifurcation to rotating polygon and the AD is above the region where polygons are observed. The surface shapes of the two cases are shown in figure 1.10(a) and (c) respectively. The authors assumed that the flow is purely inviscid except in thin boundary layers. A potential flow and a Rankine vortex were then used respectively to calculate the surface velocities of the two axisymmetric state -AD and AW -based on measured height profiles [Tophøj and Bohr, 2013]. These analytical velocities were compared to measured ones by PIV, which is shown in figure 1.10(b,d). The configuration of the case AW will be much studied in the present work since quantitative surface height and velocities are available for comparisons. Also, these measurements have been used to derive the simplified base flow models used in [Tophøj et al., 2013] and [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF]. [Tophøj et al., 2013] then assumed the whole flow to be potential and independent of z in the pure dry case of figure 1.11(c). We will verify this point by numerical calculations of the flow and by measuring the velocity at different heights. [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF] the model of [Tophøj et al., 2013] to the wet (figure 1.11a) or dry-composite cases (figure 1.11b), using a Rankine vortex in the core. This inviscid model uses an angular momentum balance to find the radius of the Rankine vortex model (the point r = r s in figure 1.11 a,b) in order to close the equation system. Actually, to establish the balance of angular momentum, the authors assume that the angular momentum injected at the bottom disk (the region from r s to R for wet and dry composite cases, and from r c to R for the dry potential case) is all absorbed at the side wall and that the transfer takes place in two thin channels along the disk and the side wall boundary layers which are assumed turbulent. [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF]. Theoretical models of the base flow used for linear stability analysis: (a,b) Rankine vortex proposed by [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF] for wet and dry-composite cases and (c) potential flow by [Tophøj et al., 2013] for dry potential case. The top graphs represent surface shape and the bottom graphs show the modeling of azimuthal velocity in each cases. For these models, the fluid follows the rotation of the disk at angular speed rΩ in the inner solid body rotation (SBR) region and supposes a potential flow (PF) with flow circulation Γ in the outer region bordered by two thin turbulent boundary layers (represented by the red bands in the top graphs).

The linear stability analysis of these models explains the emergence of rotating polygons and sloshing as a result of a resonance between different kinds of waves: gravity waves and centrifugal waves or Kelvin waves on the vortex core. When compared to the experimental results, the phase diagram obtained from this simplified model explains the wave numbers observed in experiments and also their variation with the parameters. Such models are strongly simplified: they neglect any meridional flow, boundary layers... This gives a motivation to reduce the discrepancies by characterizing the basic flow in a more complete and precise way, using numerical simulations. In a more recent work, [START_REF] Mougel | On the instabilities of a potential vortex with a free surface[END_REF] performed a global stability study on the potential vortex by taking into account the whole free-surface shape.

New additional wave resonances were reported, supported by experimental evidence.

In a recent asymptotic analysis of the base flow, [Iga, 2017] gives an in-depth characterization of the internal and boundary layers: velocity profiles, scaling laws based on an assumption that the flow is laminar. The author questions on the suitability of the turbulent model used by [Tophøj and[START_REF] Tophøj | Rotating polygon instability of a swirling free surface flow[END_REF]Mougel, 2014] to calculate the angular momentum as he points out that the flow is not turbulent when the flow is axisymmetric according to measurements of [START_REF] Tasaka | Flow transitions in the surface switching of rotating fluid[END_REF]. The flow regime, especially in the boundary layers, will be checked through numerical simulations and accurate velocity measurements.

All these models have a great physical interest, but they contain drastic simplifications which makes quantitative agreement with experiments difficult. One goal of this thesis is to precisely characterize the actual base flow using numerical simulations and experiments.

Base flow: numerical approach

Numerically, it is easier to obtain the flow solution of a configuration in which the whole fluid body covers the disk, meaning no dry center. So a large part of the present work will be devoted to the configurations without dewetting.

We can first evaluate the regime where the first bifurcation to rotating polygons occurs from the existing phase diagrams. Taking the configuration G = 0.1856 of [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] as an example, we can herein only compare amongst similar experimental set-ups (namely with the same radius R ≈ 140 mm and working fluid : water). For an aspect ratio G ≈ 0.19, [Jansson et al., 2006] found (Re, F r) ≈ (1.3 × 10 5 , 0.58) when the first unstable wavenumber m = 2 is triggered. [START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] observed the first unstable wave number m = 2 at (Re, F r) ≈ (1.6 × 10 5 , 0.56). Meanwhile, [START_REF] Bach | From Newton's bucket to rotating polygons: experiments on surface instabilities in swirling flows[END_REF] found also m = 2 at (Re, F r) ≈ (1.2 × 10 5 , 0.4) for a slight different G ≈ 0.2. However, for G = 0.2, [Mougel, 2014] observed a first unstable wavenumber m = 3 at (Re, F r) ≈ (3.2 × 10 5 , 3.4) while the dewetting takes place at (Re, F r) ≈ (2.13 × 10 5 , 1.51).

We can see clearly the quantitative difference in terms of the critical value of Froude number to develop the first rotating polygons. This fact is consistent with the conclusion that a quantitative agreement hasn't been obtained yet amongst different research groups. However, we can see that the rotating polygons occur always at high Reynolds numbers in the existing water experiments. Despite Reynolds number has not been considered as a control parameter for rotating polygons in water experiments, it represents a huge challenge for direct numerical simulations, especially in the situation where the surface is strongly deformed. [Iga, 2017] points out that the flow regimes in the circular region of the phase diagrams shown in figure 1.9 are subjected with weak disturbances. However, this flow regime can still be considered to be axisymmetric considering the disturbances have very small effect on the mean flow. This is quite reasonable since the bulk instability has the great chance to settle in at this high level of Reynolds number, even before the surface starts to deform perceptibly. However, we are still not clear how much these disturbances, necessarily three dimensional, influence the the mean axisymmetric flow. In our study, we will characterize the influence of these disturbances through numerical and experimental investigations in this "transitional regime" bridging the bulk instability to the rotating polygons. As the rotating speed increases, turbulence develops on the large structures of bulk instability.

We can have a looking back at the accomplished numerical investigations on the bottomdriven rotating flow with deformable free surface. Different from the numerous experimental investigations, the numerical work is more limited.

In fact, most of numerical works were done for flows in closed cavities or with flat free surface. As for free-surface flows, [Bouffanais and Lo Jacono, 2009b] conducted three dimensional simulation of the current configuration with small free surface deformation. [START_REF] Brady | Two-fluid confined flow in a cylinder driven by a rotating endwall[END_REF] simulated immiscible two-fluid in a confined cylinder, allowing for large surface deformation. The viscosity and density ratio of the two fluids were modified to study their influence on the swirling flow, but they were not reduced to the small level of the air-water case. [START_REF] Herrada | Off-axis vortex breakdown in a shallow whirlpool[END_REF] tackled the large density and viscosity ratio along the air-flow interface. However, the configurations the authors studied were restricted to a Reynolds level less than 2000 corresponding to a very small physical scale R = 1 mm. This little scale is comparable to the capillary scale, in which case surface tension will play an important role in the flow. Therefore, the studied regime is not suitable for studying rotating polygons.

[ Kahouadji and Martin Witkowski, 2014] calculated the steady flow solution in configurations corresponding to real experiments. Preliminary comparison with measured surface shape showed good agreement. The steady flow solver developed in his work will be a starting point for the present thesis to ameliorate the numerical tools in order to investigate configurations with large Reynolds numbers and Froude numbers.

Spin-up of the axisymmetric free surface flow

We want to study the flow motion during the transition of the liquid layer from rest to an established state after the disk is impulsively set rotating. We have chosen this unsteady process to get some insight on the interaction of the flow structure and the surface deformation. Due to the large discontinuity along the air-liquid interface, unsteady simulations of the rotating flow with free surface still remain a great challenge.

Let's start with the spin-up of a fluid in Newton's bucket when the whole container rotates as a single unit. Many studies on the spin-up of a closed cylindrical container filled with an homogeneous liquid can be found in the literature. Two reviews [Benton andClark Jr, 1974, Duck andFoster, 2001] outlined important results of the up-to-then investigations on the spin-up of homogeneous filling completely a closed container § . As a matter of fact, most of the early work used generally asymptotic methods for the specific spin-up problem in which a small speed variation is applied to a flow already established in solid-body rotation. The small Rossby number then allows for the linearization of the Navier-Stokes problem. The first order linear solution shows already interesting results, for example, the spin-up §. Spin-up of stratified fluids was also discussed in these reviews but not within the scope of the present study.

time is determined by Ekman pumping at small Ekman number ¶ . Non-linear solutions of the spin-up were also discussed in the reviews and were developed so as to achieve better agreement with experiments [Wedemeyer, 1964, Watkins and Hussey, 1973, Watkins and Hussey, 1977]. Simulations of this problem were also conducted to assess the theoretical and experimental work [START_REF] Warn-Varnas | Numerical solutions and laser-doppler measurements of spin-up[END_REF], Hyun et al., 1983].

The studies on the spin-up of fluid in solid-body rotation only flourished during the second half of last century but did not continue to the present century. Since we will concentrate on free surface flow, we will not go into further details concerning the state of the art for spin-up in closed configuration. Within the scope of the present work, we will use this configuration to validate numerical and experimental tools. In comparison with the spin-up of fluid in closed containers, very limited studies included the effect of the free surface on the spin-up in the Newton's bucket configuration. The latest work of this subject in the review given by [START_REF] Duck | Spin-up of homogeneous and stratified fluids[END_REF] is the theoretical and experimental work published by [O 'Donnell and Linden, 1991]. We will first summarize the published works during the period up to 1991. [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] studied the effect of free surface for small variation of rotation speed. They pointed out that the introduction of free surface induces tow effects compared with the closed situation: the viscous effects at the top surface are removed to the lowest order. Secondly, the free surface induces a radial motion in addition to that produced by the Ekman layer on the bottom. The Froude number was also proposed to measure the importance of surface effect on the spin-up process.

Later, [START_REF] Goller | Unsteady rotating flow in a cylinder with a free surface[END_REF]] solved the simplified Navier-Stokes equations by extending the model of the confined flow given in [Wedemeyer, 1964] to the free surface problem. The solved surface height showed good agreements with experiments. [CederlÖf, 1988] chose a parabolic bottom for the cylinder in order to match the parabolic surface. [START_REF] Homicz | Numerical model for fluid spin-up from rest in a partially filled cylinder[END_REF] investigated the fluid spin-up with free surface intersecting with the bottom or both end walls based on an extension model of [Wedemeyer, 1964] and [START_REF] Goller | Unsteady rotating flow in a cylinder with a free surface[END_REF]. Their results were later validated by experimental work performed by [START_REF] Choi | Transient free surface shape in an abruptly rotating, partially filled cylinder[END_REF], Choi et al., 1991]. We arrive now to the work of [O' Donnell and Linden, 1991] who presented theoretical and experimental results of spin-up for a full range of Froude numbers up to the one at which the surface meets the disk center. They showed the spin-up time is longer than that in a closed cylinder. Since then, interest on the spin-up with free surface in Newton's bucket configuration seems to decrease in the community[Park and Hyun, 2008]. We end the state of the art on the spin-up of a homogeneous fluid partially filled in a closed cylinder by introducing a most recent paper on this subject [START_REF] Yan | Numerical simulations of flows inside a partially filled centrifuge[END_REF], in which comprehensive axisymmetric simulation was used to solve both the liquid and air phase. A review on a more general situation on the spin-up of two (or more) immiscible fluids with comparable density or viscosity in a rotating cylindrical container gives more references. Furthermore, instead of rotating the whole container, rotating only the lid disk of the container can give birth to interesting topology changes of the surface, like the "Mt.

Fuji", as shown in [START_REF] Fujimoto | Topology changes of the interface between two immiscible liquid layers by a rotating lid[END_REF], Tsai et al., 2015, Carrión et al., 2017].

Contrary to the spin-up problem in Newton's bucket which has been investigated previously, very few experimental or numerical works have been devoted to the spin-up of the rotating free surface flow driven by the bottom disk. This stimulates us to tackle this problem in the present work.

Thesis outline

Through this general review of the rotating polygon problem, the following conclusion can be drawn: experimental observations seem to be strongly dependent on experimental setups, and there is a rather large variety of them in the published works on this subject. As a result, it is difficult to draw a universal theory intuitively just from the observed phenomena.

The approach to elaborate a base flow and then to study the linear stability is a good direction to get close to the real origin of the rotating polygons. Considering that the theoretical models of the base flow, simplified as they are, give already stimulating results, linear instabilities from a realistic base flow, computed numerically, is expected to clarify more about the bifurcation mechanism.

So the first objective of the present work is to implement numerical and experimental investigations in order to get precise knowledge of the axisymmetric base flow with deformable surface, especially in the large Reynolds number regime. The influence of the small threedimensional perturbations (bulk instability) on this base flow is also to be studied. The second objective is to characterize the unsteady flow at moderate Reynolds number and high Froude number. This study explores the transient regime for the bottom-driven flows in the 2D axisymmetric framework both experimentally and numerically. The third objective is to explore the rotating polygon regime through accurate measurements of velocity and surface deviation, and to compare rotating polygons to bulk instability patterns. The linear stability analysis of the base flow with surface deformation is outside the scope of the present work.

The thesis is organized as follows.

After the present chapter which gives a general introduction and a state of the art, chapter 2 presents two in-house numerical codes, which have been used to conduct steady and unsteady simulations of the rotating flow with free surface. The numerical method of the steady code named ROSE, is presented in details. A more general code named Sunfluidh allows for powerful direct simulations, the numerical methods of which will be outlined. Numerical tests will also be conducted to evaluate their precision and their performance in terms of CPU time. Finally, the results of the two codes will be compared, and confronted to other numerical and experimental results.

In chapter 3, we first describe the two experimental set-ups used in the present study to generate bottom-driven flows. We then present the Laser Doppler Velocimetry and the adjustments needed because of the refraction due to fluid and curved geometry. Thereafter, the techniques used to measure surface deformation are presented: two ready-to-use laser techniques for measurements at a point or along a line, and, more importantly, the imple-mentation of a third technique called Fourier Transform Profilometry (FTP) allowing for surface measurement in two-dimensional domains.

Chapter 4 is devoted to numerically and experimentally characterizing the base flow of the bottom-driven flow with a deformable free surface. Steady axisymmetric solutions of simulations by ROSE reveal the flow structure at high Reynolds numbers. The influence of several control parameters are characterized. Results of unsteady 2D and 3D simulations are then presented and confronted to the steady numerical solution. We then focus on the transitional regime between the bulk instability and rotating polygons, in which the influence of three dimensional fluctuations on the mean flow could be determined numerically and experimentally. Finally, we compare the mean flow velocities with existing models.

In chapter 5, we concentrate on the unsteady processes in free surface rotating flows. Firstly, we benchmark the transient regime with axisymmetric or almost axisymmetric free surface through numerical simulations and experiments. Actually, two different flow configurations are considered: Newton's bucket and the bottom-driven flow. Results of unsteady axisymmetric simulations and experimental measurements are discussed in details. We then quantitatively study the three dimensional instabilities by characterizing the relation between the bulk instability and the rotating polygons through precise measurements of flow frequencies and surface heights. This chapter ends with the presentation of a novel instability pattern, which had not been reported in the literature to our knowledge.

At the end of this thesis, a general conclusion and some perspectives are given.

Chapter 2

Numerical tools

Numerical simulation is a powerful tool to get precise and complete knowledge of complex flows. This is one of the two approaches in the present work. Another one is experimental investigation, as will be discussed in Chapter 3. Two in-house codes are used to study the incompressible rotating flow with free surface, generated by a rotating disk at the bottom of a cylindrical tank.

In this chapter, a state of art on the numerical methods used in the simulations of twophase flows will be firstly given. Then an in-house code ROSE (see below), developed for steady simulations will be presented in detail. Subsequently, an outline will be given for another in-house code Sunfluidh, developed by Yann Fraigneau at Limsi. After that, some numerical tests are used to evaluate the accuracy and the performance of each code. Finally, validations for the two codes are given. The two codes allow a cross validation. Besides, other numerical results in publications and experimental results are also included for validation. All the parameters involved in numerical simulations are expressed in non-dimensional scale. Lengths are scaled by the disk radius R, velocities by RΩ, time by 1/Ω and pressure by ρ(RΩ) 2 .

In most configurations studied in this Ph. D thesis, physical length scales are sufficiently large compared to capillary length scale so that noticeable surface tension effects are restricted to the vicinity of triple lines, i.e. the contact region with the side wall or the disk. Numerically, we enforced the surface to meet the solid boundaries perpendicularly or with a constant angle.

State of art on modeling two-phase flows

A general discussion about numerical methods for simulation of two-phase flows will be addressed here. The difficulty of two-phase flows consist in the presence of the phase interface, where discontinuities of density and viscosity are located. There exist generally two main approaches to model the phase interface based on the grid type:

• Grid is predefined and fixed for the whole domain containing all phases. The phase interface is tracked later by an additionally defined function.

In this approach, the Marker and Cell (MAC) method was first developed by [START_REF] Harlow | Numerical calculation of timedependent viscous incompressible flow of fluid with free surface[END_REF]. [Daly, 1969a, Daly, 1969b] adapted this method to two fluid problems and added surface tension to examine the Rayleigh-Taylor instability for finite density ratio. Then [START_REF] Chapman | Nonlinear effects in the collapse of a nearly spherical cavity in a liquid[END_REF] and [START_REF] Mitchell | Asymmetric cavitation bubble collapse[END_REF] used their method to study the collapse of a cavitation bubble. The marker particles used to identify the fluid in the MAC method were replaced by a marker function, which gave birth to the Volume-Of-Fluid (VOF) method, as presented in [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF].

A detailed review about VOF was done by [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF].

New methods kept emerging. [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]] introduced a Front-Tracking method where the interface is marked by connected points instead of the maker function in VOF while a fixed grid is still used for each phase.

The Level-Set method were used in the works of [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], Sussman et al., 1999, Chang et al., 1996, Sussman and Fatemi, 1999]. In this method, a continuous level-set function is advected with the flow, which is different from the discontinuous marker function in the VOF method. The interface is then defined by the zero contour of a distance function from the interface.

• Grids are deformed to be conforming to the phase interface. The interface is thus identified by the boundary of two sub-grids corresponding to two different phases.

In the second category of methods, the flow is resolved in a moving grid, which adapts itself to the shape of the interface. As a result, it needs a grid generation procedure at each time step, which allows a precise description of the interface. Besides, an algorithm to couple the surface evolution to the flow solver needs to be established.

[ Ryskin and Leal, 1984a] developed a technique, in which the grid is generated once the interface is modified and then the Navier-Stokes equations are resolved using the finite-difference methods. Normal stress balance is used to update the interface shape to flow solution. Then they used this technique to compute the steady solution of a buoyant, deformable, axisymmetric bubble motion ([Ryskin and[START_REF] Leal ; Ryskin | Numerical solution of freeboundary problems in fluid mechanics. part 2. buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF]Leal, 1984c]). Their steady algorithm was then generalized by [Kang andLeal, 1987, Kang andLeal, 1989] to an unsteady one in their study of the bubble deformation in a uniaxial and biaxial straining flow respectively.

Steady simulation by ROSE

The code ROSE (short for ROtating Surface Evolution) was initially developed by [Kahouadji, 2011] and Martin Witkowski Laurent based on the methods proposed in [Ryskin and Leal, 1984a[START_REF] Leal ; Ryskin | Numerical solution of freeboundary problems in fluid mechanics. part 2. buoyancy-driven motion of a gas bubble through a quiescent liquid[END_REF], Ryskin and Leal, 1984c]. It is developed for fast steady simulation for rotating flows in cylindrical configurations. It can treat flows with or without surface deformation. I have taken its last version used in [Kahouadji and Martin Witkowski, 2014] at the beginning of my Phd and continued to improve the code to calculate flows with strongly-deformed free surface, especially at large Reynolds numbers. I will give a detailed description of the numerical methods and boundary conditions of this code to simulate the bottom-driven flow.

Curvilinear coordinates

A Newtonian fluid is filled partially into a cylindrical tank. In this cylindrical geometry, the variables can be expressed in a cylindrical frame (e r , e θ , e z ) When the fluid is driven into motion by the bottom disk, the upper surface will be deformed under the effect of gravity and centrifugal forces as well as the surface tension. Here we choose a curvilinear coordinates system to adapt the numerical grid to the deformed surface h(r, t). The mapping function transforming the cylindrical coordinates (r, z)

∈ [0, 1] × [0, h(r, t)] to the curvilinear coordinates (ξ(r, z), η(r, z, t)) ∈ [0, 1] × [0, G] is defined as:    ξ = r η = G h(r,t) z.
(2.2.1) Actually, r is independent of the time while z is a function of time. The projection from the irregular interface-conforming grid to the regular rectangular grid is also illustrated in the figure (2.1). Once this system of curvilinear coordinates is chosen, all the concerned variables can be expressed as a function of ξ and η. Their partial derivatives in the cylindrical system can thus be computed by:

           ∂ ∂r = ∂ ∂ξ ∂ξ ∂r + ∂ ∂η ∂η ∂r = ∂ ∂ξ -η h h ∂ ∂η ∂ ∂z = ∂ ∂ξ ∂ξ ∂z + ∂ ∂η ∂η ∂z = G h ∂ ∂η .
(2.2.2)

We can use again the chain rule to compute the second order derivatives,

                         ∂ 2 ∂r 2 = ∂ 2 ∂ξ 2 -2η h h ∂ 2 ∂ξ∂η + η 2 h 2 h 2 ∂ 2 ∂η 2 + η 2h 2 -h h h 2 ∂ ∂η ∂ 2 ∂z 2 = G 2 h 2 ∂ 2 ∂η 2 ∂ 2 ∂r∂z = G h ∂ 2 ∂ξ∂η - Gh h 2 ∂ ∂η -η Gh h 2 ∂ 2 ∂η 2 .
(2.2.3)

Navier-Stokes equations and Boundary conditions

Navier-Stokes Equations

The dimensionless incompressible Navier-Stokes equations in the velocity-pressure form are written as:

       ∇ • V = 0 ∂V ∂t + (V • ∇)V = -∇P + 1 Re ∇ 2 V (2.2.4)
where V is the vector of flow velocity. Re is the Reynolds number defined by

Re = R 2 Ω ν
in which R denotes the cylinder inner radius, Ω is the disk angular speed and ν is the kinematic viscosity. The pressure term P is the total pressure defined by

P = P f + z F r + c 0 ,
where P f refers to the fluid pressure and F r denotes the Froude number defined by

F r = RΩ 2 g
with g as the gravity constant and c 0 a constant related to the definition of zero pressure.

If we project the vectorial equations into the cylindrical coordinates, we can get the scalar formulation of the Navier-Stokes equations:

                               ∂V r ∂r + V r r + 1 r ∂V θ ∂θ + ∂V z ∂z = 0 ∂V r ∂t + V r ∂V r ∂r + V θ r ∂V r ∂θ + V z ∂V r ∂z - V 2 θ r = - ∂P ∂r + 1 Re (∇ 2 V r - 2 r 2 ∂V θ ∂θ - V r r 2 ) ∂V θ ∂t + V r ∂V θ ∂r + V θ r ∂V θ ∂θ + V z ∂V θ ∂z + V r V θ r = - 1 r ∂P ∂θ + 1 Re (∇ 2 V θ + 2 r 2 ∂V θ ∂θ - V θ r 2 ) ∂V z ∂t + V r ∂V z ∂r + V θ r ∂V z ∂θ + V z ∂V z ∂z = - ∂P ∂z + 1 Re ∇ 2 V z (2.2.5)
where V = (V r , V θ , V z ) are the three components of the velocity in the cylindrical frame (e r , e θ , e z ) and the Laplace operator is

∇ 2 = ∂ 2 ∂r 2 + 1 r 2 ∂ 2 ∂θ 2 + ∂ 2 ∂z 2 + 1 r ∂ ∂r .
For an axisymmetric and steady problem, the variables are independent of time t and the azimuthal coordinate θ, i.e., they are only functions of r and z. Then we have:

V(r, z) =     V r (r, z) V θ (r, z) V z (r, z)     et P (r, z) (2.2.6)
The equations (2.2.5) can thus be simplified as:

                               ∂V r ∂r + V r r + ∂V z ∂z = 0 V r ∂V r ∂r + V z ∂V r ∂z - V 2 θ r = - ∂P ∂r + 1 Re (∇ 2 V r - V r r 2 ) V r ∂V θ ∂r + V z ∂V θ ∂z + V r V θ r = 1 Re (∇ 2 V θ - V θ r 2 ) V r ∂V z ∂r + V z ∂V z ∂z = - ∂P ∂z + 1 Re ∇ 2 V z (2.2.7)
while the Laplace operator turns out to be

∇ 2 = ∂ 2 ∂r 2 + ∂ 2 ∂z 2 + 1 r ∂ ∂r .
The steady axisymmetric flow is governed by the equations (2.2.7) expressed in primitive formulation. With axisymmetric property, the Navier-Stokes equations can be more conveniently expressed using streamfunction ψ, azimuthal vorticity ω, and angular momentum Γ:

                           0 =rω - ∂ 2 ψ ∂r 2 - 1 r ∂ψ ∂r + ∂ 2 ψ ∂z 2 0 = ∂ ∂r 1 r ∂ψ ∂z ω - ∂ ∂z 1 r ∂ψ ∂r ω - ∂ ∂z Γ 2 r 3 - 1 Re ∂ 2 ω ∂r 2 + 1 r ∂ω ∂r - ω r 2 + ∂ 2 ω ∂z 2 0 = ∂ ∂r ∂ψ ∂z Γ - ∂ ∂z ∂ψ ∂r Γ - r Re ∂ 2 Γ ∂r 2 - 1 r ∂Γ ∂r + ∂ 2 Γ ∂z 2 .
(2.2.8)

In this formulation, the continuity condition is validated automatically. Also, the treatment of the coupling between pressure terms and velocities is avoided. The relations between the unknowns of the two systems (2.2.7) and (2.2.8) are:

V r = 1 r ∂ψ ∂z , V z = - 1 r ∂ψ ∂r , ω = ∂V r ∂z - ∂V z ∂r , Γ = rV θ .
(2.2.9)

By projecting the cylindrical coordinates to the curvilinear system, the equations (2.2.8) are rewritten as:

                                                       ξω - ∂ 2 ψ ∂ξ 2 + 2η h h ∂ 2 ψ ∂ξ∂η - η 2 h 2 + G 2 h 2 ∂ 2 ψ ∂η 2 + 1 ξ ∂ψ ∂ξ -η 2h 2 -h h h 2 + h ξh ∂ψ ∂η = 0 ∂ ∂ξ G ξh ∂ψ ∂η ω - ∂ ∂η G ξh ∂ψ ∂ξ ω + Gh ξh 2 ∂ψ ∂η ω - G h ∂ ∂η Γ 2 ξ 3 - 1 Re ∂ 2 ω ∂ξ 2 -2η h h ∂ 2 ω ∂ξ∂η + η 2 h 2 + G 2 h 2 ∂ 2 ω ∂η 2 + 1 ξ ∂ψ ∂ξ - ω ξ 2 + η 2h 2 -h h h 2 - h ξh ∂ω ∂η = 0 ∂ ∂ξ G h ∂ψ ∂η Γ + Gh h 2 ∂ψ ∂η Γ - ∂ ∂η G h ∂ψ ∂ξ Γ - ξ Re ∂ 2 Γ ∂ξ 2 -2η h h ∂ 2 Γ ∂ξ∂η + η 2 h 2 + G 2 h 2 ∂ 2 Γ ∂η 2 - 1 ξ ∂Γ ∂ξ + η 2h 2 -h h h 2 + h ξh ∂Γ ∂η = 0
(2.2.10) The new system of equations is discretized in the curvilinear grid and solved by the Newton-Raphson method which will be discussed later.

Boundary conditions

According to the boundary nature in the bottom-driven flow, the fluid partially filled in the cylindrical tank is constrained to the following boundary conditions.

• The impermeability condition on both the rotating disk and the stationary sidewall;

• Regularity at the axis;

• Continuity of the velocity and stress components at the interface. The tangential stress exerted by the air on the interface is neglected due to its negligible dynamic viscosity compared to that of the liquid.

In the steady state, the surface is no longer evolving. Therefore, the velocity component normal to the surface should be zero. If we want to follow the temporal evolution of the surface in spin-up experiments, the normal velocity component is not zero and corresponds exactly to the velocity of the interface (see Appendix A.2).

The boundary conditions at each side can be translated into the following mathematical relations.

• East border: cylinder's side wall

ξ = 1, 0 < η < G :                  ω = ∂ 2 ψ ∂r 2 = ∂ 2 ψ ∂ξ 2 -2η h h ∂ 2 ψ ∂ξ∂η = ∂ 2 ψ ∂ξ 2 as ∂ψ ∂η = 0 ψ = 0 Γ = 0
(2.2.11)

• West border: axis of rotation

ξ = 0, 0 < η < G :              ω = 0 ψ = 0 Γ = 0
(2.2.12)

• South border: rotating disk

0 < ξ < 1, η = 0 :                  ω = 1 r ∂ 2 ψ ∂z 2 = G 2 ξh 2 ∂ 2 ψ ∂η 2 ψ = 0 Γ = ξ 2
(2.2.13)

• North border: free surface This is where the free surface locates. There are two sets of conditions:

(i) The condition to guarantee a steady shape of the surface and impermeability

V • n = 0 ⇒ V z = h (r) • V r (2.2.14)
which can be reduced to:

ψ = constant = 0 (2.2.15)
(ii) Generally, the stress balance along an interface between two fluids can be expressed as:

(T a -T)n = σn(∇ • n) -∇σ (2.2.16)
where T, T a are the stress tensor of the fluid and the air respectively. n = 1 √ 1+h 2 (-h , 0, 1) t is the unit outward normal vector to the surface z = h(r). The surface tension, σ, will be considered constant in the following.

Tangential stress balance When modeled in the code ROSE, the free surface has zero tangential stress.

           t t • T • n = 0 ⇒ -2h ∂V r ∂r + 2h ∂V z ∂z + (1 -h 2 ) ∂V r ∂z + ∂V z ∂r = 0 e t θ • T • n = 0 ⇒ -h ∂V θ ∂r - V θ r + ∂V θ ∂z = 0
(2.2.17)

for the liquid side. Here, t and e θ are the two tangential unit vectors to the surface:

Chapter 2. Numerical tools

e θ = 1 √ 1 + h 2 (0, 1, 0) t (2.2.18) t = 1 √ 1 + h 2 (1, 0, h ) t (2.2.19)
By following the same step in the proof (A.2.6), the equation (2.2.17) can be simplified as:

           ξω = -2Gh h(1 + h 2 ) ∂ψ ∂η * -h ∂Γ ∂ξ + 2h Γ ξ + G(h 2 + 1) h ∂Γ ∂η = 0
(2.2.20)

Normal stress balance Projecting equation (2.2.16) in the direction n yields the normal stress balance at the interface:

n t • (T a -T) • n = σ(∇ • n) (2.2.21)
The jump of normal stress across the interface must balance the curvature force per unit area. The computation of normal vector and curvature are very well explained in [Bush, 2013]. The curvature is also derived in axisymmetric case. There is however a small typo in his equation 5.41. The corrected expression is :

∇ • n = - h r(1 + h 2 ) 1/2 + h (1 + h 2 ) 3/2 .
(2.2.22) Equation (2.2.21) can thus be developed into: (2.2.23) in dimensional scale. P * f , P * a are respectively the static pressures at the liquid and air sides of the interface. Equation (2.2.23) can be non-dimensionalized to:

P * f -P * a - 2µ 1 + h 2 h 2 ∂V r ∂r + ∂V z ∂z -h ( ∂V r ∂z + ∂V z ∂r ) +σ h r(1 + h 2 ) 1/2 + h (1 + h 2 ) 3/2 = 0,
P f -P a - 2 Re(1 + h 2 ) h 2 ∂V r ∂r + ∂V z ∂z -h ( ∂V r ∂z + ∂V z ∂r ) + 1 W e h r(1 + h 2 ) 1/2 + h (1 + h 2 ) 3/2 = 0
(2.2.24) * . Note that there is a typo in the relation ξω = -2Gh h h 2 (h 2 +1) ∂ψ ∂η in the equation ( 16) of [Kahouadji and Martin Witkowski, 2014] where W e denotes the Weber number defined as :

W e = ρR 3 Ω 2 σ
Since the fluid pressure by P f = P -h F r -c 0 , we obtain

P - h F r -P a -c 0 - 2 Re(1 + h 2 ) h 2 ∂V r ∂r + ∂V z ∂z -h ( ∂V r ∂z + ∂V z ∂r ) + 1 W e h r(1 + h 2 ) 1/2 + h (1 + h 2 ) 3/2 = 0, (2.2.25)
By moving the term h/F r to one side of the above equation, we obtain

h F r = P -(P a + c 0 constant C ) - 2 Re(1 + h 2 ) h 2 ∂V r ∂r + ∂V z ∂z -h ( ∂V r ∂z + ∂V z ∂r ) + 1 W e h r(1 + h 2 ) 1/2 + h (1 + h 2 ) 3/2 (2.2.26)
Projected in the curvilinear system, the equation (2.2.26) is written as:

h F r = P + C - 2 Re G h ∂V z ∂η - Gh h ∂V r ∂η - h h h 2 + 1 V r viscous term + 1 W e h r(1 + h 2 ) 1/2 + h (1 + h 2 ) 3/2 interfacial tension term (2.2.27)
This equation will be used to update the surface height, as will be presented in the section 2.2.3. The constant C can be determined by the volume conservation of the fluid phase

2π 0 1 0 h(ξ)ξdξdθ = πG.
(2.2.28)

Pressure derivative at interface

Since pressure is included in equation (2.2.27), we need also to calculate the pressure to renew the surface. As is known, the velocity components can be calculated directly from the solutions (ψ, ω, Γ) of equations (2.2.10) through the formula (2.2.9). We will take advantage of the relation between the pressure and velocity in equations (2.2.7) which can be rewritten Chapter 2. Numerical tools in the curvilinear system as:

                                   ∂V r ∂ξ + V r ξ -η h h ∂V r ∂η + G h ∂V z ∂η = 0 V r ∂V r ∂ξ + G h V z -η h h V r ∂V r ∂η - V 2 θ ξ = - ∂P ∂ξ + η h h ∂P ∂η + 1 Re ∇2 V r - V r ξ 2 V r ∂V θ ∂ξ + G h V z -η h h V r ∂V θ ∂η + V r V θ ξ = 1 Re ∇2 V θ - V θ ξ 2 V r ∂V z ∂ξ + G h V z -η h h V r ∂V z ∂η = - G h ∂P ∂η + 1 Re ∇2 V z (2.2.29)
where,

∇2 = ∂ 2 ∂ξ 2 -2η h h ∂ 2 ∂ξ∂η + η 2 h 2 + G 2 h 2 ∂ 2 ∂η 2 + 1 ξ ∂ ∂ξ + η 2h 2 -h h h 2 - h ξh ∂ ∂η (2.2.30)
With the help of the steady kinematic condition (2.2.14) and the mass conservation condition (the first equation of the system (2.2.29)), a combination of the second and the fourth equation of the system (2.2.29) leads to

∂P ∂ξ = -V r ∂V r ∂ξ -V z ∂V z ∂ξ advection + V 2 θ ξ rotation + 1 Re -h ( ∂ω ∂ξ + ω ξ ) + G(1 + h 2 ) h ∂ω ∂η viscous .
(2.2.31)

An integration from ξ = 0 to ξ = 1 of this equation gives rise to the pressure distribution at the liquid side along the surface, which is inserted later to equation (2.2.27) to update the surface height.

Algorithm of ROSE

The presence of the surface height h(r) in equations (2.2.10) and the boundary conditions [(2.2.11), (2.2.13), (2.2.20),(2.2.27)] make the calculation difficult to accomplish if all the equations are solved together. Here, we use "a projection method" to uncouple the resolution of the equations (2.2.10) and resolution of the surface shape. Concretely, a prescribed surface height is inserted into the equations 2.2.10, which is then resolved by the Newton-Raphson method. The newly computed flow velocities as well as pressure obtained from the direct solutions (ψ, ω, Γ) are then inserted to equation (2.2.27) to update the surface height. Newton method is also used to resolve equation (2.2.27). This process is included into an outer loop until the flow solution and surface height match each other. As shown in figure 2.2, two Newton loops are implemented to resolve separately the Navier-Stokes equations and the equation of normal stress balance inside an outer loop in which a relaxation factor is introduced to converge to the steady solution. Normally, the two inner Newton loops take a number less than 10. The outer loop bounded by the red rectangle is to repeat the two inner Newton loops so that the residual calculated on the liquid volume .

resid V = 1 0 (h n+1 -h n )dr, (2.2.32)
is small enough. The relaxation factor α needs to be set optimally large in order to converge to the steady solution with least iterations for the outer loop. Usually with a value α = 0.5, only 50-100 iterations for the outer loop are needed to get a residual of machine precision. It is worth noting that α can not be too large to avoid divergence of calculation, especially for strongly deformed surface. The solutions for configurations with lower Reynolds or Froude number are always saved to initialize the calculations of configurations with high Reynolds or Froude number, which will greatly reduce the computation time.

A detailed development about the Newton-Raphson method can be found in [Ypma, 1995]. Here, I just introduced the principle of the method. Generally, a problem can be simplified as A(x) -B = 0

where A and B can also depend on x if the problem is nonlinear. Alternatively, it's to find the fixed point of the function F(x) = Ax -B. Then we need to solve the equation

F(x) = 0
which could be difficult through direct resolution. Newton's method consists in approaching iteratively the solution by solving

J f (x n ) • ∆x n+1 = -F(x n ) (2.2.33)
where J f is the Jacobian matrix of the operator F. A simulation test has been conducted for the configuration G = 0.25, F r = 1.23, Re = 2062. A mesh of 401 × 101 and a relaxation factor α = 0.5 are used. Each iteration of the outer loop is composed of 6 N-S iteration and 10 iterations for the normal stress balance loop. When initialized by the solution for a configuration with a lower Froude number F r = 0.88, the calculation of the test configuration took only 80 seconds to decrease the residual resid V from 10 -2 to 10 -12 within 62 iterations in the outer loop. The calculation was calculated on four Intel(R) Xeon(R) CPU E5-2623 processors as the PARADISO solver is parallelized.

(V r , V θ , V z ) and P (ξ, G) (V r , V θ , V z ), P (ξ, G), h C = 1 in (2.2.27)
compute Jacobian matrix A h and second term b h from Eq.(2.2.27)

b h max < ξ 2 solve A h U h = b h (h, C) = (h, C) + U h
Steady ( ĥ, C) yes no Newton loop for normal stress balance

h new = (1 -α)h + α ĥ (h new , C) Steady solution (ψ, ω, Γ) (V r , V θ , V z ), (h, C) Outer loop

Grid dependence analysis

Since the finite-difference method of second-order central scheme is used to approximate all the spatial derivatives, the simulation results should theoretically have a second-order precision.

The second order can be easily confirmed for simulations with flat surface (see [Kahouadji, 2011]). To complete his study on the grid dependency, a bottom-driven configuration with flat surface at a much higher Reynolds number G = 0.1856, Re = 15 000 is used to study the influence of grid size for solutions of velocity components (V r , V θ , V z ), stream function (ψ) and the total kinetic energy E c . Five grids are used, which are respectively 162 × 50, 324 × 100, 648 × 200, 1296 × 400 and 2592 × 800 to compute a domain [0.4, 1] × [0, 1]. Figure 2.3 shows the differences in terms of different parameters between the computed solutions and the ones deduced by Richarson extrapolations. Fittings from the most resolved points show always a slope equal to 2 in the logarithmic scale. The configuration (G, F r, Re) = (0.248, 1.496, 1047) [Kahouadji and Martin Witkowski, 2014] is also undertaken to study grid influence as strongly deformed surface is present in this configuration compared to the above one. Four grids (N r × N z ) are tested: 100 × 25, 200×50, 400×100 and 800×201. Since the surface is deformed, it's difficult to find a common node for all these grids. Therefore, only the total kinetic energy E c is used to study its error with respect to the grid size. Figure 2.4 shows that the error between calculated values of Ec and the Richardson extrapolation value,which is proportional to the grid spacing with a slope of 2 in the logarithmic scale. This confirm the second-order precision for the spatial discretization even on a curvilinear grid. .248, 1.496, 1047) where curvilinear grid is used for spatial discretization as the surface is strongly deformed.

Unsteady simulation by Sunfluidh

We will present in this section the numerical codes to establish unsteady simulations, always for the rotating flow in a bottom-driven experiment. Attempt has been undertaken to extend the steady ROSE to unsteady simulations for axisymmetric configurations without or with deformed free surface (see Appendix A). Since it's not so efficient using Newton method for unsteady simulations, we used mainly another code, Sunfluidh, to implement unsteady simulations.

Sunfluidh is developed by Yann Fraigneau from the year of 2011 for 2D/3D unsteady simulations of incompressible flows or flows under low Mach number hypothesis. It covers a large range of flows. The part for incompressible two-phase flows based on a level-set method is in progress and gives very satisfactory results in 2D axisymmetric configurations. Orthogonal structured grid is required for Sunfluidh. We use mainly Sunfluidh for unsteady simulations of incompressible flow, one-phase or two-phase, at Reynolds numbers less than O(10 4 ) in cylindrical geometry.

In the following, only the numerical methods used in the present context are outlined. Actually, Sunfluidh offers several other options of numerical methods. The reader is referred to [START_REF] Yann | Sunfluidh description[END_REF] for more details. Evaluation of the time cost will also be given for simulations by Sunfluidh. Finally, I will present the numerical set-ups for simulations established in the following chapters about results and analysis.

Numerical methods in Sunfluidh

The Navier-Stokes equations expressed in (V, P ) form (see equations (2.2.4)) are discretized on a structured staggered grid in order to avoid checkerboard instability in the pressure field. Actually, pressure is defined at the cell center while the velocity components are defined at the center of cell interfaces. The viscous term and the non-linear convection term are discretized using a second order central scheme. A conservative formulation is implemented.

A 2nd order Backward Differentiation Formula is used for the time discretization, that is:

( ∂V ∂t ) n+1 = 3V n+1 -4V n + V n-1 2δt + O(δt 2 ), (2.3.1)
where the exponent n represent the state at discrete instant t n . In order to relax the constraint on the time step for numerical stability, the viscous term will be defined at t n+1 which makes the discretized Navier-Stokes equations semi-implicit. The convective flux, NL (denoting hereafter (V • ∇)V for convenience) is estimated by linear extrapolation of the values at the two precedent time steps:

NL n+1 = 2NL n -NL n-1 . (2.3.2)
The pressure gradient is defined explicitly on P n to uncouple the pressure and the velocity. By considering all the previous discretizations, the original equations (2.2.4) are transformed into its discrete incremental version:

1 - 2δt 3Re ∇ 2 δV * = S V , (2.3.3) with δV * = V * -V n and S V = V n -V n-1 3 -∇P n + NL n+1 - 1 Re ∇ 2 V n • 2δt 3 .
The notation V * is a prediction velocity for V n+1 which does not satisfy the divergence free condition as a result of the explicit representation for the pressure gradient in equation (2.3.3). It is solved by the alternating direction implicit method [Hirsch, 2007]. The next step is to project V * onto a space of divergence-free velocity field to update the velocity and pressure.

Considering ∇ • V n+1 = 0, by subtracting equation

3V n+1 -4V n + V n-1 2δt + NL n+1 = -∇P n+1 + 1 Re ∇ 2 V n+1 (2.3.4) from 3V * -4V n + V n-1 2δt + NL n+1 = -∇P n + 1 Re ∇ 2 V * , (2.3.5)
we can obtain the Poisson equation:

∇ 2 Φ = 3∇ • V * 2δt , (2.3.6) with Φ = P n+1 -P n + 1 Re ∇ • V * . (2.3.7)
The Poisson equation is solved later by Successive Over-Relaxed method (SOR) coupled with a nV-cycle multigrid method, directly implemented in the code in order to accelerate the convergence. The solution for the scalar Φ allows to obtain P n+1 and V n+1 using

P n+1 = P n + Φ - 1 Re ∇ • V * , V n+1 = V * - 2δt 3 ∇Φ.
(2.3.8)

The numerical methods are only presented assuming a constant density and dynamic viscosity. The projection method has been adapted in the new version of Sunfluidh to treat two-phase flows. Further details on the level set method can be found in [Couderc, 2007] and the implementation of the surface tension follows the ideas developed in [M. Coquerelle and Glockner, 2016].

Two-dimensional unsteady simulations

Since only a meridian plane is considered in 2D axisymmetric simulations, the computing domain is a rectangle and the grid used is thus rectangular. Grid can be refined at boundaries if necessary (see figure 2.9(a)).

Unforced flat surface :

To validate the Sunfluidh code, we have chosen a well known configuration, G = 1 and Re = 900 with flat surface, studied before by [START_REF] Piva | Steady axisymmetric flow in an open cylindrical container with a partially rotating bottom wall[END_REF] and later on by [Bouffanais and Lo Jacono, 2009a]. The configuration is calculated with Sunfluidh on a regular grid 96 × 96, modified ROSE (see Appendix A.1) on a regular curvilinear grid 151 × 151 and Gerris on a regular grid 128 × 128. For the Sunfluidh simulation, a CF L = 0.5 constraint was used so that the time step was δt = 0.01. A run of this Sunfluidh simulation on four Intel(R) Xeon(R) CPU E5-2623 processors took 104 seconds to achieve the steady solution at t = 185, when a residual defined as

resid(t) = 1 Vol Vol i (V n+1 i -V n i ) 2 δt dv (2.3.9)
decreases below the value 10 -12 . In the definition of this residual, Vol stands the volume of the simulated physical domain, V i stands for the velocity components (V r , V θ , V z ).

The temporal evolution of the velocity components at r = 0.5, z = 0.5, calculated by the three codes are shown in figure 2.5. We show also the comparison of the steady velocity components in figure 2.6. Good agreements are achieved amongst these three codes, which validates Sunfluidh as well as the unsteady ROSE for flat-surface flow.

Deformed surface

Another simulation on a two-phase flow (air and glycerine/water mixture), with G = 0.25, F r = 1.23, Re = 2063 was also established. A computing domain (r, z)

∈ [0, 1] × (a) (b) (c) 0
(a) V r , (b) V θ , (c) V z .
[0, 0.5] was discretized on a regular grid 128 × 64 and the time step was set to 2 × 10 -3 resulting from a constraint CF L = 0.3 for stable simulation. The speed of disk is set directly to the desired speed at t = 0 in the simulation. It took 72 minutes on a single processor to reach the steady solution at t = 192, where resid < 10 -12 is obtained.

The temporal evolution of the velocity components at three probes (r = 0.4, 0.6, 0.8) at the same height z = G/2 is shown in figure 2.7. The oscillations present at the beginning of V r and V z evolution are enhanced by the abrupt spin-up of the disk in the simulation. The implementation of a linear variation of the disk speed in the simulation (see Chapter 5) will damp these oscillations. There is no possible comparison with modified unsteady version of ROSE for this configuration since it is still in development (see Appendix A.2). However, we can compare the steady solution obtained by Sunfluidh simulation with that obtained by ROSE simulation, as shown in figure 2.8.

No significant difference was found between them. Similar comparisons of steady solutions are shown later in figure 2.12 for larger Reynolds and Froude numbers. V θ with 21 equispaced iso-contours between 0 and 1. Probes are also shown here. (c) V z with 11 equispaced iso-contours between [V z,min (-0.08), 0] and 11 equispaced iso-contours [0, V r,max (0.23). The maxima and minima positions of V r and V z are marked by red '+' and blue '×'. Iso-contours from ROSE results are also presented by magenta lines with the same contour levels.

Three-dimensional simulation with flat interface

We describe here the 3D computation of Chapter 4. Three-dimensional simulations cost much more computing time and should be used for well-designed configurations. In cylindrical geometries, hexahedral grids need to be used. I will take the configuration G = 0.1856, F r = 0, Re = 30 000 as an example to show the grid distribution in a 3D simulation. A grid of size 192 × 448 × 128 in (r,θ,z) directions for a computing domain (r,θ,z) ∈ [0.4,1] × [0,2 3 π] × [0, 0.1856] was used. Figure 2.9 shows the grid by plotting one node out of 5 for clear demonstration. The grid spacing in each direction of the cylindrical frame is shown in figure 2.10.

It is worth taking different measures to reduce the calculation cost. A coarse grid with local-refinement is preferred when spatial resolutions are needed, e.g., in the boundary layers. According to the phenomenon we want to capture by the simulation, a complete computing domain can be reduced to one part of it on condition that no impact would be raised. In the present 3D simulation, a hub was introduced at the center (r ∈ [0, 0.4]) of the computing domain, which is considered to have no effect on the solution compared to a complete domain since the central fluid is shown to be in solid rotation through 2D simulations. Also, two thirds in the azimuthal direction was cut to reduce the simulation cost, provided that expected three-dimensional instability mode can be still captured. Since a complete 3D calculation of the whole spin-up from rest to established regime costs a lot, the established solution for the 2D axisymmetric configuration (which is already unsteady) is used to start the 3D simulation to quickly reach established solutions. Even with all these measures, it took still about one month for the three-dimensional simulation on four processors to reach the established state, during which the non-dimensional time is advanced by 140.

Validation of numerical simulations

Comparison amongst different codes

The two in-house codes allow a mutual validation between them. Generally, difficult configurations with large Reynolds and Froude numbers are taken for comparison to confront the their limits.

Comparison for configurations with flat surface

In the firstly place, we show the comparison between these two codes for configurations with flat surface (F r = 0). Take G = 0.1856, Re = 10 000 as an example since it will be studied further in the following chapters. The computing domain is the same: In the following, only results from ROSE and Sunfluidh simulations obtained on the coarsest grids (250 × 51 for ROSE, 128 × 32 for Sunfluid) are shown for more complete comparison. The comparison of iso-contours of flow velocities and streamlines ψ are given in figure 2.12. Even though these grids are coarse especially for a Reynolds number of 10 000, the flow is well resolved as can be seen in figure 2.12. Good agreements are also achieved between ROSE and Sunfluidh.

(r, z) ∈ [0, 1] × [0,
The code ROSE was also used to calculate the configurations discussed in [START_REF] Herrada | Off-axis vortex breakdown in a shallow whirlpool[END_REF], where Reynolds numbers are less than 2000 while Froude numbers varies from O(1) to O(100). As a matter of fact, their configurations are designed based on taking water as the working fluid and an imaginary cylinder radius of 1 mm, which is difficult to realize in real experiments. To simulate their configurations by ROSE, a grid of 150 × 150 was used. The contour levels presented in figure (6,8,10) of [START_REF] Herrada | Off-axis vortex breakdown in a shallow whirlpool[END_REF] are not communicated. (a) Streamlines of the meridional circulation consisting of 21 equispaced solid contours of ψ between 0 and ψ max with ψ max = 0.00347 (0.00341) for ROSE (Sunfluidh) and one dashed contour at value 1 2 ψ min with ψ min = -2.06×10 -4 for ROSE (-2.03×10 -4 for Sunfluidh). (b) Isocontours of the azimuthal velocity V θ , with 21 equispaced contours between 0 and 1 for the two codes. (c) 21 equispaced solid contours of V z between 0 and V z,max with V z,max = 0.127 (0.125) for ROSE (Sunfluidh) and one dashed contour at value 1 2 V z,min with V z,min = -0.088 for ROSE (-0.085 for Sunfluidh). (d) 21 equispaced solid contours of V r between 0 and V r,max with V r,max = 0.137 (0.141) for ROSE (Sunfluidh) and one dashed contour at value 1 2 V r,min with V r,min = -0.167 for ROSE (-0.164 for Sunfluidh). Positions of maximum('+') and minimum ('×') ψ, V z and V r are also shown: magenta markers for ROSE and black markers for Sunfluidh.

So I considered that the contour levels plotted in their figure 2 are always used for all the rest figures in this paper, that is, 15 positive levels between 0 and ψ max and 15 levels between ψ min and 0. The results of ROSE are then plotted into the same graphs of their results, as seen in figure 2.13 .

The results of ROSE are generally coherent with those in [START_REF] Herrada | Off-axis vortex breakdown in a shallow whirlpool[END_REF]. However, the difference of the interface profiles can be remarked for Re = 2000. In fact, the free surface is slightly more deformed in ROSE simulation. To confirm, a two-phase simulation by Sunfluidh for the configuration having Re = 2000 was also carried out using a grid 128 × 128 for the whole domain [0, 1] × [0, 1]. The results are also plotted in the same graph of figure 2.13(c). We can see ROSE agrees well with Sunfluidh for the surface height as well as for the positive streamlines.

Comparison for deformed surface

We have shown the good agreement between ROSE and Sunfluidh calculations for configurations without surface deformation and will continue to compare them in simulation for flow with deformed surface. The configuration G = 0.1856, Re = 10 000 is taken again but with a large Froude number, F r = 1, to create strongly deformed surface. For ROSE, the solution is sought in a monophasic deformed domain on grid N r × N z = 501 × 101 while for Sunfluidh, the steady two-phase flow solution is determined in the domain (r, z) ∈ [0, 1] × [0, 0.5] on a regular mesh N r × N z = 256 × 128 from an initially rest state with a liquid-air interface at z = G. Figure 2.14 shows the comparisons of results of these two simulations by presenting in the same way as in figure 2.12. Here, a strongly-deformed surface is present, which almost touches the bottom disk. Good agreements between ROSE and Sunfluidh are always achieved. The shift of positions of some extrema of velocity components and stream function is just a direct result of taking into consideration the air phase when searching the extrema. If only liquid phase is checked, the extrema are found to differ only by about 1% between the two codes, which is fully satisfactory as these two techniques to deal with surface deformation are completely different.

Comparison with experimental measurements

Numerical results are also compared to measurements obtained from experiments. Some comparisons of surface heights between ROSE simulation results and experimental measurements were given in [Kahouadji and Martin Witkowski, 2014]. In the following, comparison between numerical and experimental velocities is presented. The configuration with G = 0.186, Re = 10 000 with flat surface can be realized by experiments in the large experimental setup (see Chapter 3). The rotating speed of the disk was set to 0.51 rad/s and water was filling to an initial height 26mm. As a result, we have G = 0.1856 and Re = 10 000 and a very small Froude number (F r = 0.0037). No surface deformation was observed during the experiment. LDV (discussed in Chapter 3) was used to measure the azimuthal and axial velocities at the mid height z = G/2.

The measured velocity values are then compared to the previous numerical solutions 12 is used except in this figure: ψ max = 0.0035 and ψ min = -1.83 × 10 -4 are both found in the liquid phase for ROSE while ψ max = 0.0038 and ψ min = -0.0016 are both located in the air phase for Sunfluidh.

V z,max = 0.136 (0.134) for ROSE (Sunfluidh) and V z,min = -0.086 (-0.085). V r,max = 0.147 (close to the disk) and V r,min = -0.168 (under the surface) are found in the liquid phase for ROSE while for Sunfluidh, V r,max = 0.152 is found in air-phase and V r,min = -0.168 on the surface. The markers V z,max and V r,min of ROSE are covered by those of Sunfluidh.

of ROSE and Sunfluidh, as seen in figure 2.15. As seen in figure 2.15(a), the numerical azimuthal velocity components are generally coherent with the measurements and the solid body rotation close to the axis and the side boundary layer are well captured. The difference of V θ is found around the position r = 0.5, the exact position where the overshoot appears in numerical solutions. As for the axial velocity V z , a remarkable difference can be seen from figure 2.15(b). We do not have yet an explanation on this large difference. This will be studied by Antoine Faugaret during his PhD thesis.

Conclusion

In this chapter, the numerical methods used by the two in-house codes are presented: ROSE and Sunfluidh. Both of them have second order precision. Tests have shown that the code ROSE is able to compute efficiently the steady solution of a rotating flow with deformable free surface. Actually, ROSE can calculate configurations with Reynolds number up to 20 000, which is a great improvement compared to the previous version. Simultaneously, Sunfluidh shows great interest to accomplish reliable unsteady simulation of two-phase flows with deformable interface in axisymmetric configurations. Also, the computing cost is quite reasonable considering a full knowledge of the unsteady transition is then acquired. In the following chapters, numerical results of simulations by these two codes will be presented, together with experimental measurements.

Chapter 3

Experimental techniques

This chapter presents the experimental set-ups and methods. Beside flow visualization used to quickly determine the flow state, quantitative measurements of velocities and surface deformation were carried out to get precise flow field information. We first present the two experimental set-ups used in this work, which mainly differ by their sizes. The LDV technique used to measure flow velocities inside cylindrical containers will then be presented, with a focus on the corrections induced by the geometry and the fluid used. Finally, three kinds of techniques for measuring surface deformation will be presented, two ready-to-use laser techniques for measurements at a fixed position and along a line, and a third technique called Fourier Transform Profilometry (FTP) allowing for two-dimensional domains. A detailed presentation of this latter technique developed during the PhD is made.

Experimental set-ups

The rotating flow generated by a rotating disk placed at the bottom of a cylindrical tank is already presented schematically in figure 1 The DC motor model is a RX 320E produced by Parvex. It is mounted directly to the shaft of the disk with a gearbox reduction of 12 to 1. The speed loop control is provided by a servo amplifier RTS also manufactured by Parvex. This servo amplifier that drives the DC motor has an input voltage ±10V. The DC motor can be controlled in two ways, manually with a potentiometer or automatically. In automatic mode, the servo amplifier receives analog voltage signal that controls the rotating speed. Practically, the input of the servo amplifier is provided by a programmed Mbed card, which is a micro-controller powered by a USB port. Since the analog output of the Mbed card is limited to 3.3 V, an amplifier with a gain of 3 is necessary. The Mbed card is also useful to generate TTL * signals with rising or falling edge to synchronize the velocity and surface measurements, as explained later.

In experiments, the Mbed card triggers the measurements in the following way: the measurement is started ∆t w1 seconds earlier than the start of disk rotation and lasts ∆t w2 seconds after the complete stop of disk rotation. The disk is not started abruptly but accelerated from rest to the desired speed within ∆t r seconds with a linear variation of speed , for mechanical reasons. The disk is also stopped using a linear ramp within ∆t e from the constant speed to zero. In practice, ∆t w1 = 5 s, ∆t r is aways set as 0.5 s or 1 s. ∆t p is the duration of the plateau at constant rotating speed Ω, which is set from 30 s to 240 s according to fluid dynamics phenomenon we wish to describe. The descending ramp duration ∆t d is set always as 10 seconds. At last, the duration to continue measurements after the disk has stopped. ∆t w2 is set to 5 s or 10 s. * . TTL: Transistor-transistor logic

Measurements of velocities

Laser Doppler Velocimetry: basics

Laser Doppler Velocimetry has been largely used to establish precise real-time measurement of velocity components in transparent or semi-transparent flows. In our case, LDV is carried out using a Dantec BSA † A signal beam of continuous argon-ion laser with wavelength 660 nm and power 25 mW is emitted and splitted into two coherent ones of equal intensity and with shifted frequencies. These two beams are focused by transmitting optics. In the intersection probe volume, parallel fringes are generated due to interference between the coherent beams. Hollow silver-coated glass spheres with a diameter of 10 micrometers are seeded in the fluid. When seeded particles carried by the fluid pass through the probe volume, the laser light is scattered all around. A photo-detector converts the fluctuating intensity of the scattered light to an electrical signal containing the Doppler frequency proportional to the velocity component perpendicular to the fringe direction.

The frequency shift of the two laser beams makes the fringe pattern in the probe volume move at the shift frequency. As a result, particles moving with positive and negative velocities will generate Doppler frequencies larger and smaller than the shift frequency, respectively. The sign of the flow direction is then determined.

The fringe spacing d f depends on the wavelength λ of the laser light and the angle θ between the beams:

d f = λ 2 sin(θ/2) . (3.2.1)
When a particle goes from one fringe to the other, it travels on a distance d f orthogonally to the fringes, and the Doppler frequency f D provides information about the travel time. Therefore, the normal velocity component V is calculated by

V = d f • f D = λf D 2 sin(θ/2) .
(3.2.2) †. BSA: Buest Spectrum Analyses Flow explorer. The optical head is fastened to a horizontal-vertical displacement platform (see figure 3.3). This optical head is connected to the BSA system thanks to an optical fiber.

The present system gives very precise measurements. For example, the measured value of velocity in a still fluid is 10 -4 m/s, which is considered as the precision of the velocity measurements. For steady flows, the root mean square can also be measured. This parameter characterizes the turbulence level of the flow.

Laser Doppler Velocimetry: corrections

LDV is initially calibrated for measurement in a specific fluid, notably in air, see figure 3.4. The two beams focus and cross directly at the measurement point, located at the focus length f LDV away from the laser head. In the present study, LDV is used to measure flow velocities in liquids (water or glycerine), the refractive index n of which is different from that of air. Consequently, refraction will take place at each interface between two transmission media due to different refractive indices. Table 3.2 lists the indices of the media involved in the present experimental study. The refractive index of glycerine-water solutions depends on the glycerine mass ratio. Details can be found in table 23 of [Glycerine Producer's Association, 1963]. The mixture of 80% glycerine and 20% water (by weight) is denoted hereafter by G8W2.

Laser Head

As the liquid is stored in a tank with a certain wall thickness, the laser beams, after transmission through air, pass through the cylinder wall (PMMA or glass) and finally focus in the liquid. As a result, the lasers will bend twice when passing the cylinder outer and inner borders. The focus point of the bended beams will be moved from its original point in air. Moreover, when doing measurements in the horizontal (r, θ) plane, the cylindrical border introduces a curved interface, which systematically modifies the incident angles of the beams when the laser head is moved in the radial direction. I use the ray-tracer method proposed in [START_REF] Huisman | Applying laser doppler anemometry inside a taylor-couette geometry using a ray-tracer to correct for curvature effects[END_REF] to correct LDV measurements for refraction and curvature effects.

In the following, I will discuss separately the correction used when measuring axial velocity V z in a vertical plane and measuring azimuthal velocity V θ in a horizontal plane. Correction for axial velocity V z Figure 3.5(a) shows a side view of the laser transmission when installed vertically to measure V z . Since the two laser beams meet the cylinder outer wall by two points aligned vertically, the incident angle θ i is always constant(θ i = θ a /2) when the laser head is displaced along the radial direction. Equation (3.2.1) can be rewritten, according to Snell's law n a sin(θ a /2) = n l sin(θ l /2), as:

f D 2V real = sin θ l λ l = sin θ a λ a = f D 2V measure ⇒ C z = V real V measure = 1. (3.2.3)
The subscript a refers to air while l referes to liquid. So n a and n l are respectively the index of air and liquid while θ a /2 and θ l /2 are respectively the angle of refraction in air and in the liquid. So, the measured velocity can be correctly calculated from the Doppler frequency and the correction coefficient C z for the measure of V z is always 1. However, it is clear that the real measurement position is shifted to r f from the theoretical position r air = r LH -f LDV (see figure 3.5(a)). In fact, optical relations prove that r f is always a linear function of r air , for example the blue curve in figure 3.6(a) for the situation where water flow is filled partially inside the large experimental set-up (r f /R = 1.338 r air /R -0.3615).

Correction for V θ When laser head is installed so that the two beams are in the horizontal plane to measure the azimuthal velocity V θ , the incident angle of the beams changes always owing to the curved interface when the laser head is moved horizontally, and the relation in equation (3.2.3) does not hold anymore. A coefficient defined as

C θ = V real V measure = λ l sin(θ l /2) • sin(θ a /2) λ a = n a sin(θ a /2) n l sin(θ l /2) (3.2.4)
needs to be employed to correct measured velocity values. The refractive angle can be obtained using the ray tracer method presented in [START_REF] Huisman | Applying laser doppler anemometry inside a taylor-couette geometry using a ray-tracer to correct for curvature effects[END_REF]. In the situation of measuring water flow in the large set-up, position correction and correction factor are represented by the green curves in figure 3.6(a) and (b), respectively. Two campaigns of measuring true focus positions r f of LDV in water/large set-up situation were conducted to confirm the theoretical position correction and they are respectively represented by the black '×' and red '+' in figure 3.5. An error of 1.3 mm at worst is found between the theoretically-corrected positions and the measured ones. The cylinder walls of the two experimental apparatus have two different refractive indices and thicknesses. Also, two different liquids are used. So specific correction needs to be applied to each combination of liquid and apparatus. Figure 3.6 shows the LDV correction for the combination of water and large tank. Figure 3.7 shows the correction for the combination of G8W2 and large tank. The linear relation for correcting vertical position is r f /R = 1.446 r air /R-0.4704 in this situation. Figure 3.8 shows the correction for LDV measurements in flow of water as well as in flow of G8W2 inside the small set-up. The linear relations for correcting vertical position are respectively r f /R = 1.337 r air /R -0.3707 for water case and r f /R = 1.447 r air /R -0.4841 for G8W2 case in the small set-up. The positions of V θ measurements differ slightly between the two liquids while the values of the correction factor significantly change from water to G8W2. Similar figures can be plotted for other combinations of liquid and apparatus when needed. Overall, these corrections for the position and value of LDV measurements are noticeable and thus absolutely necessary for the set-ups and working fluids we have used. When the azimuthal velocity is measured in the horizontal plan, the position shift is of order 0.1R due to the curvature of the tank wall and the increased liquid refraction inside. And in the vertical plan, the real measurement position is respectively about 1.34 times in water and 1.45 times in G8W2as large as that without correction. In parallel, the correction on the measured value of the azimuthal velocity component can make a difference of 30% in the worst case for radii around the axis.

Measurements of surface deformation

As the present work concerns the rotating flow in an open cylindrical tank, we are also interested in characterizing the free surface. The free surface of liquid will be deformed under the effect of the gravity, centrifugal force and surface tension in the co-rotating frame of the rotating tank. Depending on the situations, the deformed free surface can be either axisymmetric or three-dimensional.

Non-intrusive methods are usually favored to measure the surface deformation as they do not modify the original flow. [START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF] proposed a technique, named Free-Surface Synthetic Schlieren, to measure the free surface deformation by comparing the refracted images of a flat and a deformed interface through which a set of random points printed at the bottom of a channel is observed by a camera from above. The correlation between the two digital images allows for the computation of the height slope, which is integrated at a second step based on a least-square inversion of the gradient operator to obtain the surface deformation height. This technique is not suitable to measure strongly deformed surfaces encountered in our experiments. [START_REF] Cobelli | Global measurement of water waves by fourier transform profilometry[END_REF] proposed an optical profilometry method using Fourier transform, named Fourier transform profilometry (FTP). This technique measures instantaneous 3D surface deformation of a non-transparent liquid with large slope. So we decided to use this technique in our experiments. During my first year of PhD, I implemented a FTP system in the lab.

We will first describe two ready-to-use laser sensors for point and line measurements then three-dimensional measurements of a given surface by FTP in section 3.3.3.

In order to measure the liquid surface height, non-transparency is needed for laser and FTP techniques. Therefore, specific material should be used to seed and color the liquid on one hand but should not change dramatically the liquid properties on the other hand. [START_REF] Przadka | Fourier transform profilometry for water waves: how to achieve clean water attenuation with diffusive reflection at the water surface?[END_REF] found TiO 2 as coloring pigment, which makes the water surface more light diffusive while minimizing surface effects. This material is used for water and glycerinewater mixtures in our experiments for measurement of surface deformation either by laser method or by FTP.

Point measurement

Keyence LK-152 is a laser triangulation sensor to determine the position of a target by sending a laser light and then receiving the reflected light. In conjonction with a LKG3001 controller, it was used for measuring the height evolution of a point on non-transparent surface in the present study. The installation of the point laser head is shown in figure 3.9(a).

Line measurement

LMI3D Gocator 2350 laser profile sensor is used to measure along a line the position of a non-transparent target. Its laser emitter projects a laser line onto the target and the reflected light is captured from an angle by the sensor's camera installed beside the laser emitter. The sensor's emitter, its camera and the target form a triangle, which helps to determine the distance from the sensor to the target using the known characteristics. The measurement range of this model is 400 mm and the field of view is 158 -368 mm beyond a clearance distance of 400 mm. The resolutions in x and z direction are announced as 0.15 mm and 0.019 -0.06 mm, respectively. The installation of this sensor in the laboratory is shown in figure 3.9(b).

3D surface measurement

As mentioned before, FTP has been widely used in free surface deformation measurements: water waves [START_REF] Cobelli | Global measurement of water waves by fourier transform profilometry[END_REF], Cobelli et al., 2011a, Cobelli et al., 2011b], drop impacts on surfaces [START_REF] Lagubeau | Spreading dynamics of drop impacts[END_REF], Le Goff et al., 2013], deformation of particle monolayer at a liquid-fluid interface under compression [START_REF] Jambon-Puillet | Wrinkles, folds, and plasticity in granular rafts[END_REF]]. In the following section, the principle of this technique will be presented in the first place. Then effects of involved parameters on measurement precision are investigated. Finally, experimental tests on solid and liquid surfaces are discussed.

Fourier transform profilometry: principle FTP principle

We use the parallel FTP configuration described in [Maurel et al., 2009, Takeda andMutoh, 1983], which is shown schematically in figure 3.10. The present FTP set-up is composed of two optical devices. A Canon XEED SX 60 digital video projector (resolution: 1400 × 1050), controlled by computer and installed vertically, projects a predefined one-dimensional sinusoidal fringe pattern on a flat surface as reference at one time and on the deformed surface at another time. These two images are registered by a camera (UI CP 3370, resolution: 2048 × 2048), which is also installed vertically so that its optical axis is parallel to that of the video projector and its entrance pupil is at the same height as the exit pupil of the video projector. The distance between the camera entrance pupil and the flat reference plane is denoted by L and the distance between the two optical axis by D. The recorded images display a fringe pattern whose intensity has two sources of noises. One (B) results from illumination inhomogeneity or background variations over the field of view and is an additive term. The other one (A) is due to the local surface reflectivity, which modulates the intensity of the projected fringes. We denote I(x, y) and I 0 (x, y) the signals of images taken on the deformed surface and on the flat one, respectively. Their general forms are then expressed as

I 0 (x, y) = A(x, y) cos(2πy/p 0 + ϕ 0 (x, y)) + B(x, y) I(x, y) = A(x, y) cos(2πy/p 0 + ϕ(x, y)) + B(x, y) (3.3.1)
in a configuration where the fringes are parallel to the x direction, with wavelength p 0 . When I(x, y) is compared to I 0 , we find that there exists a phase shift (∆ϕ = ϕ(x, y) -ϕ 0 (x, y)) in I(x, y) arising from the variation of surface elevation. An extraction of the phase-shift map gives rise to the reconstruction of the surface height through geometrical optics. The mechanism of FTP can be easily understood by the demonstration in figure 3.11.

The phase map ∆ϕ is extracted through the following two steps: and

•
I m = I -B.
Then a band-pass filter centered on the carrier frequency f 0 = 1/p 0 is applied. Inverse Fourier transform is then applied so as to get complex signals in spatial domain, of the form:

I 0f (x, y) = Ã(x, y) exp{i[2πy/p 0 + ϕ 0 (x, y)]}, I f (x, y) = Ã(x, y) exp{i[2πy/p 0 + ϕ(x, y)]}.
(3.3.2)

•
Step 2 is aimed at recovering the phase shift ∆ϕ between I 0f and I f . Quantity ∆ϕ is obtained by taking the imaginary part of

log(I f • I * 0f ) = log(| Ã|) 2 + i∆ϕ. (3.3.3)
However, ∆ϕ can only be determined with values in an interval of extent 2π, thus leading to discontinuities. Phase unwrapping needs to be applied to yield a continuous phase shift map. The unwrapping algorithm unwrap of MATLAB is used in our postprocessing.

The unwrapped phase-shift map is inserted into equation

h(x, y) = L∆ϕ(x, y) ∆ϕ(x, y) -2πD/p 0 (3.3.4)
to yield the height distribution h(x, y) ( [START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-d object shapes[END_REF]).

Influence of parameters on FTP precision

Before using FTP for real measurements, we tested firstly in a numerical approach the whole post-processing algorithm of FTP using synthetic signals. This kind of numerical test for the FTP system were previously conducted in [START_REF] Cobelli | Global measurement of water waves by fourier transform profilometry[END_REF] In their studies, the synthetic signal is computed out of a height variation due to a gaussian shape. Our approach is similar to these studies but we apply it to a right triangular prism of prescribed height h(x,y) which height is 3 cm and the triangle base length is 6 cm. The discontinuity of slope at the crest of the prism is this quite challenging. A right triangular prism [START_REF] Cobelli | Global measurement of water waves by fourier transform profilometry[END_REF] of prescribed height h(x, y) has been taken as a test deformed surface. Within this framework, the influence of each involved parameter on reconstruction precision is studied. From the above section, the following parameters can affect the precision of height reconstruction: fringe wavelength p 0 , pixel size p x , filter width σ c . The pixel size refers to the length in the reference plane corresponding to one pixel in the camera sensor. To simplify the simulation, the background is composed simply of uniform gray level 128 with a random noise (function rand in Matlab) superimposed :

I 0 (x, y) = 128 + w 0 (x, y) + 120 cos(2πy/p 0 ), w 0 ∈ [-5, 5] I(x, y) = 128 + w(x, y) + 120 cos(2πy/p 0 + ∆ϕ(x, y)), w ∈ [-5, 5] (3.3.5) where ∆ϕ(x, y) = - 2πDh(x, y)/p 0 L -h(x, y) ,
and w 0 , w denoting the noises. The noises are imposed to have 10 gray levels, which is typical for the camera (not thermally regulated) we have been using. An amplitude of 120 gray levels for the sinus function holds in the quite favorable case where the entire scale over the 8bits resolution of the camera is exploited.

The 2D Gaussian band-pass filter employed in FTP measurement is of form

F (f x , f y ) = exp - (f x -f 0x ) 2 2σ 2 cx - (f y -f 0y ) 2 2σ 2 cy . (3.3.6)
In the following numerical simulation, we pick the dimensions of the real experiment, the distance between the camera's entrance pupil (or the projector's exit pupil) and the reference plane is set as L = 1.159 m. The distance between the two parallel optical axis is set as D = 17.9 cm. The fringes are along x direction. If two-dimensional FFT are done on the signals I m0 and I m , spectrum peaks are expected at (0, f 0 ) in the spatial frequency domain. As explained in [START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-d object shapes[END_REF], the height information is stored in this main peak whose width is proportional to the height gradient along the direction normal to the fringes, i.e, ∂h ∂y in present situation. A filter is then designed to widely select this main peak: we set f 0x = 0 and f 0y = f 0 as center frequencies on one hand; on the other hand, the standard deviation σ cy needs to be set for optimal extraction. Nevertheless the choice of σ cx is much less constrained and it just needs to be large enough, say 10σ cx , to avoid information loss yet remove undesired small scale perturbation at high spatial frequencies. The same definition of relative error as introduced in [START_REF] Cobelli | Global measurement of water waves by fourier transform profilometry[END_REF] is used in the present study to evaluate reconstruction precision:

ξ = |h r (x i , y j ) -h(x i , y j )|/ |h(x i , y j )| (3.3.7)
where h r is the reconstructed surface height for the reference one h.

Gradient of height

When the triangular prism is placed horizontally, the height of the upper surface with respect to the bottom one varies only in one direction and remains constant in the other direction. We choose here two different cases. In the first one, the crest is parallel to the fringe direction, i.e., along x (see figure 3.12(a1)) while in the second one, it is perpendicular to the fringes (see figure 3.12(b1)). The surface shape is encoded differently into the modulated fringe pattern for these two positions (see figure 3.12(a2, b2)). Figure 3.12 shows the reconstructed profiles for these two orthogonal positions under the same conditions as follows: L = 1.159 m, D = 17.9 cm, p 0 = 3 mm, p x = 1.72 • 10 -4 m. We have thus 17.4 pixels in a fringe. The reconstructed profile h r in the perpendicular position is better and has a relative error of 0.24% against 0.58% in the parallel position, as shown in figure 3.13 (The curve in the perpendicular case coincides with the reference one). This result is not necessarily intuitive but can be explained by the fact that the main peak at f 0 , shown in figure 3.14(b) for the perpendicular case, is much thinner than that for the parallel case shown in figure 3.14(a). This results in an extraction with more height information through filtering.

The nth component of the spectrum of I m (x, y) along the f y direction covers local frequencies f n (see equation ( 23) in [START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-d object shapes[END_REF]):

f n = f 0 + n 2π • ∂ϕ ∂y , (3.3.8)
indicating that the width of the main peak depends on the gradient of the phase-shift and thus the height gradient along the y direction. For the case of figure 3.12(b1), ∂h/∂y = 0, while it can be nonzero in the case of figure 3.12(a1). This is the reason why the main peak for the parallel case shown in figure 3.14(a) is wider than that for the perpendicular case shown in figure 3.14(b). In the following, only the parallel case is taken in the simulations to study the influence of other parameters.

Fringe wavelength By imposing p x = 1.72 • 10 -4 m and σ y /f 0 = 0.3, the fringe wavelength p 0 was varied from 1 mm to 8 mm. The reconstructed curves and the corresponding errors with respect to fringe wavelength are shown in figure 3 We can infer from figure 3.16 that the smaller the wavelength of the fringe pattern, the more precise the reconstruction. The decrease of p 0 moves the main peak even further away from zero, same for the harmonics. This makes filtering easier and gives more precise results. However it should be noted that the number of pixel in a wavelength decreases also when the wavelength is reduced, which can lead to insufficient sampling.

Pixel size

We set this time p 0 = 0.002 m, σ cy /f 0 = 0.3 and vary the pixel size p x from 10 -4 mm to 8 • 10 -4 mm, see figure 3.17. The point with p x = 8 • 10 -4 m corresponds to only 2.5 pixels in a wavelength.

When the wavelength p 0 is fixed, the augmentation of the pixel size p x , equivalent to a diminution of resolution, reduces the number of pixel in a fringe period. The error thus increases but remains below 0.5% until the pixel size is so large that it leaves only 2 pixels in a period.

Filter width

We define the cut-off frequency f c of our band-pass filter as the frequency at which the filter attenuates the input by -3 dB in the power spectrum (or 1/ √ 2 in the amplitude spectrum). The filter width is equal to 2(f c -f 0 ) = 2 √ ln 2σ c in the x and y directions. I will use σ c as an indicator of filter width for presentation convenience. The 2D spectrum map will have its principal peak at the fundamental frequency pair (f 0x , f 0y ) = (0, f 0 ). Consequently, the filter width in the x-direction σ cx just needs to be large enough in x-direction as discussed previously. However, the filter width in y-direction σ cy = σ c needs to be selected carefully to get the principal spectrum at f 0 . We then set σ cx = 10σ c .

As discussed above, the wider a filter is, the more of the main peak will be selected and the more precise will be the measurement. But too wide a filter will possibly capture some part of the zero spectrum component and also some of the harmonics. As a result, there exist an optimal filter width with minimum reconstruction error.

We can estimate the error with respect to filter width by conducting simulations with constant p 0 = 0.003 m, p x = 1.72 • 10 -4 m. From figure 3.18, we find the optimal filter width at σ c /f 0 = 0.4, with error less than 0.5% in the range of 0.33 -0.66.

The above simulations are done in ideal conditions: the background image is set constant, the fringes have a sinusoidal shape (or close to sinusoidal since noise has been added) and the mean values of I m0 and I 0 are almost zero so that no zero-frequency component is introduced. In reality, the signal quality of images taken by the camera can be largely degraded by shadows, reflections, electronic noise from projector and camera. A pair of polarized filters are installed to reduce reflection: one at the exit of the projector and another at the entrance of the camera. Light is found colored and inhomogeneous in space after passing through the filter at the projector's exit, creating a gray-level gradient in background image. All these constraints lead to a smaller filter width than the optimal one previously found numerically. In practice, σ c /f 0 = 0.3 was used to obtain a filter widths equal to f 0 /2.

Range of measurement and uncertainty in profile reconstruction

As discussed in the original paper [START_REF] Takeda | Fourier transform profilometry for the automatic measurement of 3-d object shapes[END_REF] about FTP, this technique has a maximum range of measurement in terms of the height gradient in the direction normal to the fringe in order to make the main peak to be separable from all other spectra. The slope limitation for FTP is

| ∂h(x, y) ∂y | max < 1 3 • L D . (3.3.9)
Based on this point and the results of the previous simulation results, many studies have tried solutions to circumvent such limitation. Generally, they can be classified into three types.

• Larger L/D

• Lower fringe wavelength

• Remove zero spectrum A FTP set-up with larger L/D can increase the measurement range. However, L should not be too large so that enough fringes are presented on the target surface and enough pixels in images. In the present study, the height is set L ∼ 1m and the distance D ∼ 18 cm. Then theoretically, the maximum measurerable slope is about 1.85 for the present installation. The maximum slope for the deformed surface in Newton's bucket with G = 0.49, to be discussed in subsection 3.3.3.4, is 1.4. Whereas rotating polygons with large amplitude discussed in Chapter 5 can challenge the limit of 1.85 locally: for the monogon-1 with G = 0.25 a maximum slope of 1.8 is found; the rotating ellipse with G = 0.5 has a maximum slope of 1; the maximum slope for the rotating triangle with G = 0.1856 is 2 which slightly exceeds the limitation.

The wavelength can be adjusted in the computer by setting the number of fringes projected by the video projector. However when we can increase the fringe number, the contrast and the pixel number inside a wavelength is reduced as well. In our measurements, we chose to project a pattern of 100 fringes.

A π phase-shifting was proposed in [START_REF] Li | Improved fourier transform profilometry for the automatic measurement of three-dimensional object shapes[END_REF] to remove zero component of the signal spectrum. Theoretically, the elimination of zero spectrum component extends the limitation for an extractable width of the main peak to a range of [0, 2f 0 ]. The maximal measurable slope of height is thus increased to three times of that without elimination of the zero component. In this method, two deformed fringe patterns with π phase difference are captured sequentially. Therefore, their technique is limited for static objects.An improved π phase-shifting FTP method is proposed by [START_REF] Hu | Surface profile measurement of moving objects by using an improved π phase-shifting fourier transform profilometry[END_REF]] by using two synchronized line-scan cameras for dynamic measurements.

Besides all these techniques implemented on the measuring stage, the zero spectrum component can also be eliminated at the stage of post-processing using still the standard FTP. Empirical mode decomposition is proposed by [START_REF] Li | Eliminating the zero spectrum in fourier transform profilometry using empirical mode decomposition[END_REF] to eliminate the zero spectrum component, or wavelet transform by [START_REF] Gdeisat | Eliminating the zero spectrum in fourier transform profilometry using a two-dimensional continuous wavelet transform[END_REF]. Since the fringe wavelength used was small enough to leave the main peak rather far from zero in the present study, the method of eliminating the zero spectrum component didn't bring any evident advantages for height reconstruction thus was not pursued.

The relative error of the reconstructed h r can be deduced from equation (3.3.4):

∆h h ≤ ∆L L + ∆D D + ∆(∆ϕ) ∆ϕ + ∆p 0 p 0 (3.3.10)
Therefore, the uncertainty on h r comes from two sources: the error in the determination of phase shift and in the measurement of the geometrical parameters (L, D, p 0 ). On one hand, the contribution from the phase uncertainty results mainly from the signal processing (Fourier transform, filtering) and electronic noise, which generate a zero mean contribution but with a non-zero variation Local shadows and reflection may result in parasite points in height reconstruction, which should be suppressed manually. On the other hand, the errors on the geometric parameters are systematic and can be minimized by calibration on known heights.

Experimental tests

Video projector and camera were installed and aligned with caution to ensure necessary alignments. Then the present FTP was firstly tested on solid surfaces of 3D printed objects shown in figure 3.19. The characteristics of the objects will be given in the following. Liquid surfaces in Newton's bucket are also used to test our FTP system.

Measurement on triangular prism

A triangular prism with the same dimension as that of the prism described in subsection 3.3.3.2 has been printed out using a 3D printer, as shown in figure 3.19 and the FTP system has been tested on it.

Firstly, we select the lengths L = 0.957 mm, D = 17.9 cm and p 0 = 3.0 mm. The filter width was set to be σ cx /f 0 = 0.3 for all measurements. Figure 3.20(a1) shows the height profiles at different positions along the ridge direction, i.e., at different x for this parallel case, and figure 3.20(b1) those at different y positions for the perpendicular case. The profile Figure 3.19 -3D printed objects for testing the FTP system. From left to right: reference plane plate, triangular prism and plate with a Gaussian shape. measured by the line laser displacement sensor (Gocator 2350) is also plotted as a reference. We can see that FTP captures generally well the triangle shape except in the vicinity of the top where a gap of 1.5 mm is present between the FTP profiles and the reference for the parallel case and a gap of 1 mm for the perpendicular case. The error at this sharp point is reasonable as a discontinuity of gradient is present here, which strongly widens the peak which in turn is severely truncated by the band pass filter. The fact that the reconstruction for the perpendicular position is slightly better than that for parallel position is coherent with the results of the simulations above. Zoomed views these profiles on a segment of 0.2 mm in y and x direction respectively, shown in figure 3.20(a2, b2), indicate that the fluctuation of measurements by FTP is less than 0.3 mm for the parallel case and 0.2 mm for the perpendicular, even at the top position.

Apart from the above camera-surface plane distance L = 0.957 m, two other vertical positions are also tested: L = 1.092 m and 1.159 m by moving the prism and the reference plane vertically. For these two L values, the fringe wavelength is respectively 3.4 mm and 3.6 mm. Figure 3.21(a, b) shows the profiles averaged along the ridge direction of the prism measured at the three values of L, respectively for the parallel and perpendicular prism position. By comparing the curves for these three L values, we can infer that within this range of variation of L, the FTP measurements are independent of the distance from the object to the camera.

Measurement on a Gaussian shape

We have tested the FTP system with a smooth 3D printed Gaussian shape (as shown in figure 3.19) of equation:

h(x, y) = 0.03 • exp - x 2 + y 2 2.5 • 10 -4 (m) (3.3.11)
Figure 3.22(a) shows the measurement of the Gaussian shape by FTP under conditions: L = 0.957 mm, D = 17.9 cm and p 0 = 3.0 mm. Profiles across the maximum position are plotted in figure 3.22(b), together with the one measured by laser. FTP generally correctly captures the height except again near the top with an error less than 1.5 mm.

The Gaussian shape could be placed top down in order the FTP could be tested with a surface of negative elevation. Same procedure was used for measuring this surface with negative height. Results are plotted in figure 3

Test on liquid surface

The objective in implementing the FTP system is to measure liquid surfaces. We thus performed some tests using the axisymmetric parabolic surface of solid rotation (Newton's bucket) flows.

The Newton's bucket flow is generated in a set-up with its side wall and bottom rotates as a whole, which is schematically presented in figure 3.24. The working liquid in this test is chosen as G8W2. The initial liquid height is 34.2 mm. The motor is piloted by the Mbed card using the following protocol: a linear ramp of 0.5 s to accelerate, a plateau of 30 s at constant speed Ω = 14.56 rad/s and then another linear ramp of 10 s to stop the cavity.

Once all adjustments of camera and projector settings are done, the liquid surface is filmed by the camera with a frame rate of 30 Hz, from 5 s before the start of the cavity rotation to the re-establishment of the flat surface. The series of images are then post-processed in Matlab. It should be noted that the phase difference method presented in [START_REF] Su | Dynamic 3-d shape measurement method based on ftp[END_REF] was used to compute more precisely the instantaneous phase shift for each registered image along time. If each deformed pattern is compared directly to the initial reference one, it's more difficult to unwrap correctly the phase shift, especially for images associated to strongly deformed surfaces. Figure 3.25(a) gives a snapshot of the liquid surface at steady stage. Since a cylindrical geometry is not suitable for Fourier transform, we will try to obtain a maximum rectangular region. In order to reconstruct mostly the deformed liquid surface, a maximum area covered by fringes is firstly chosen in order to get rid of the shadows and the cylinder wall, as shown the region surrounded by the blue solid circle in figure 3.25(a). A square zone (shown by the red square) containing exactly the blue solid circle is extracted from the original image and we obtain figure 3.25(b). But the four corners of the square are still not completely covered by fringes. Another circle with a radius of 10 pixels smaller that of the blue solid circle is also drawn, as shown by the blue dash circle in figure 3.25(b). The mean value in x direction of the gray levels between these two circles are extended horizontally from the solid blue curves outwards to the left and right borders. In this way, we get an image full of fringes, shown in figure 3.25(c). Finally, the post-processing can be executed directly on the whole image of figure 3.25(c). It is worth noting that the fringes outside the solid blue circles are just created artificially to facilitate the post-processing of FTP. So the heights reconstructed on outside the blue solid circle don't have physical meaning.

As a reminder, L = 0.9543 m and p 0 = 3 mm were measured out based on the liquid surface at rest. Figure 3.26(a) shows temporal height evolution at three pair of symmetric radial positions, measured by FTP as well as by the Gocator (with a sampling frequency of 100 Hz). We can see that FTP values vary slightly (less than 0.3 mm) at steady state, which is in the same order of the laser measurements. The symmetry of surface deformation is also verified for the FTP measurements since the evolution at symmetric radial positions coincide perfectly. The steady profiles across the rotating axis (0, 0) along x and y axes are plotted To conclude for the present FTP system, it can be used for measurements with generally precision of precision 0.4 mm while local errors of 1.5 mm can be present for some local positions with largest deformations. Even though it can not achieve the precision of laser measurement, FTP has the advantage of measuring instantaneously a three-dimensional surface distribution. As a result, we may combine these two methods to measure free surface deformation owing to occurrence of instabilities.

Since the error of the above measurements by FTP is proportional to reference heights, I will try to calibrate the relation between the height measured by FTP and the reference one by the laser in order to ameliorate the precision of the present FTP system.

Amelioration of FTP precision

As mentioned above, the heights measured by FTP seem to be proportional to the reference heights. This conclusion can be clearly shown if h r is plotted as a function of h laser (see figures 3.27-3.30).

First of all, still liquid surface was studied for its perfect horizontalness. Practically, the volume of liquid is adjusted so as to displace the free surface by regular steps. One of these heights is considered as the reference height h = 0 and the others can be reconstructed using FTP techniques. Once this z origin is chosen, we measured out L = 0.9543 m, D = 17.9 cm and fringe wavelength on this surface p 0 = 3 mm. After post-processing the images, surfaces at different heights are reconstructed. The problem that arises is that the measured height of a single horizontal surface at h = 0 is no longer uniform is space since the abrupt height change at the borders will spread to all the frequency domain in after Fourier transform. The variation of the measured height can be 2 mm when scanning along y and is found symmetric about the mean value while it is smaller than 0.3 mm along x direction. So the height h r0 at the image center is plotted as a function of the reference height h laser in figure 3.27. By fitting these points, we can find a linear relationship: where, α, β represent the correlation and the offset coefficients.

Consequently, we may think that the error of FTP measurements may simply come from the wrong measurement of the parameter L because h r is exactly proportional to L according to the reconstruction formula (3.3.4). Then we can choose to simply correct L by a quantity ∆L such that:

L + ∆L L = 1 α =⇒ ∆L = (1/α -1)L (3.3.13)
A value of ∆L = 2.9 cm is obtained in the present case. We may ask whether this correction ∆L is uniform for all above measurements. The same procedure was followed to compute the correction ∆L for the previously measured solid triangle, Gaussian shapes and liquid parabolic surfaces.

The measured heights for the triangular prism presented in figure 3.21 are plotted as a function of the laser results in figure 3.28 (a, b) for the parallel and perpendicular cases, respectively. Linear regression was also done for these curves. The fitted linear lines, as well as the resulting correction ∆L or ∆L ⊥ , were listed in table 3.3. A mean value of the correction values for the parallel and vertical cases is also given for each L. We can see that the slopes of the fitted lines are almost independent of L. Therefore, we cannot obtain a uniform value of ∆L for different measured L according to the relation (3.3.13). But it's worth noting that the correction value of 2.9 cm for L = 0.957 m based on the triangular prism is the same as that found on calibrating still liquid surfaces with an almost identical L. h r (m), parallel h r (m), perpendicular ∆L(cm) L=0.957 m h r = 0.9708 h laser -2.64 • 10 -5 h r = 0.97086 h laser -1.9 • 10 -5 2.9 p=3 mm -→ ∆L = 2.9 cm -→ ∆L ⊥ = 2.9 cm L=1.092 m h r = 0.9669 h laser -3.3 • 10 -5 h r = 0.97141 h laser -1.5 • 10 -4 3.45 p=3.4 mm -→ ∆L = 3.7 cm -→ ∆L ⊥ = 3.2 cm L=1.159 m h r = 0.9664 h laser -1.09 • 10 -4 h r = 0.96622 h laser -1.9 • 10 -4 4.05 p=3.6 mm -→ ∆L = 4.0 cm -→ ∆L ⊥ = 4.1 cm Table 3.3 -Linear regression relations of the curves presented in figure 3.28(a, b) for the triangular prim and resulting corrections ∆L for each L.

For the Gaussian shape with positive elevation, only the line and the column crossing the maximum position of the whole measured surface were studied. Figure 3.29 shows the heights of FTP as a function of laser heights and table 3.4 lists the regression functions and the resulting values of ∆L. ∆L seems to be coherent for the triangular surfaces and this Gaussian surface. It remains almost the same for the same value L but increases with L due to a slightly-decreased α. Since the above corrections ∆L are all calculated based on reconstruction of surfaces with positive elevation, What about the correction for surfaces with negative elevation? So we continue to study the inverse Gaussian shape and the parabolic surface of Newton's bucket.

h r (x, y = 0)(m) h r (x = 0, y)(m) ∆L(cm) L=0.957 m h r = 0.96764 h laser -1.9 • 10 -4 h r = 0.97299 h laser -1.43 • 10 -4 2.95 p=3 mm -→ ∆L 1 = 3.2 cm -→ ∆L 2 = 2.
Like the straight Gaussian shape, the line and the column profiles h(x, 0) and h(0, y) of the Gaussian surface with negative elevation measured by FTP were also plotted as a function of Gocator-measured heights (not presented here). The linear regression relations of these curves, as well as the resulting ∆L, are listed in table 3.5. h r = 0.95283 h laser + 2.65 • 10 -4 6.27 p=3.6 mm -→ ∆L l = 6.8 cm -→ ∆L c = 5.74 cm Table 3.5 -Linear regression relations for the curves presented in figure 3.23(b) for the Gaussian shape with negative elevation and resulting corrections ∆L for each L.

Similarly, the temporally-evolving heights at given radial positions for the Newton's bucket, shown in figure 3.26(a) were also plotted as a function of the corresponding red curves measured by the laser, see figure 3.30. By fitting these curves, we obtain table 3.6, listing all the linear regression relations and resulting ∆L. The values of ∆L calculated on the basis of height evolution at the positions of r = 0 and ±0.25 are coherent with those for the triangular prism and the Gaussian shape with positive elevation for 0.95 m The value of ∆L calculated on the basis of h(r = 0.5, t) is larger than the other two values, which may be explained by the small range [0, 1.3 cm] of the height evolution at r = 0.5. Consequently, h(r = 0.5, t) is not suitable for the fitting.

The calibration on the Gaussian shape with negative elevation seems not coherent with other cases, specially for the value of ∆L at L ≈ 0.95 m. This could be caused by the non perfect fabrication of the back surface of the printed Gaussian model. Consequently, it is not axisymmetric. Then the plot of the FTP measurements is not exactly at the same r/R h r (r, t)(m) ∆L(cm) 0 h r = 0.9695 h laser + 1.04 • 10 -4 3 -0.25 h r = 0.9713 h laser -5.6 • 10 -5 2.8 0.25 h r = 0.9732 h laser + 2.3 • 10 -5 2.6 -0.5 h r = 0.9485 h laser -3.9 • 10 -5 5.2 0.5 h r = 0.9469 h laser -2.3 • 10 -5 5.4

Table 3.6 -Linear regression relations of the height evolution at given probes measured by FTP as a function of laser measured one in the spin-up experiment of solid body rotation at14.56 rad/s and the resulting corrections ∆L for each L.

position where the laser measured. Moreover, the negative surface of the Newton's bucket does not cause an incoherent value of ∆L for L ≈ 0.95 m, when compared to the surfaces with positive elevation. So the correction is not a problem of positive or negative elevation. It seems that the correction ∆L is uniform for the same value of L but it increases slightly with L. However when used in the correction of a theoretical height equal to 30 mm, the variation of ∆L from 3 cm to 4 cm only leads to a difference of 0.3 mm for the reconstructed height when L = 0.95 m or a difference of 0.26 mm when L = 1.16 m. In other words, an error of 1 cm on the correction ∆L only gives rise to a difference of about 0.3 mm on the final reconstructed height. Therefore, the method of calibrating on the measured camera-reference plane distance L works to enhance the precision of the FTP system. Finally, we chose to correct L by adding to it ∆L = 3.5 cm to get a maximum error of 0.5 mm (corresponding to 1.5 mm without correction) at worst for a height of 30 mm.

Introduction

A disk rotating at the bottom of a fixed cylindrical tank partially filled with a liquid induces a flow that potentially shows various instability patterns, such as rotating polygons, switching and sloshing phenomena [Vatistas, 1990, Jansson et al., 2006, Suzuki et al., 2006, Tasaka and Iima, 2009, Iima and Tasaka, 2016, Tasaka and Iima, 2017]. Out of these experimental studies, phase diagrams have been established using two parameters, namely, the initial liquid height and angular speed of the disk [Jansson et al., 2006, Iga et al., 2014, Bach et al., 2014]. The high Reynolds numbers prevailing in these experiments preclude a realistic numerical simulation of such configurations.

An inviscid model for the axisymmetric base state on which these patterns grow has been proposed by [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] and later on improved by [Tophøj et al., 2013] and [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF]. Despite their simplicity, these models were able to capture the occurrence of patterns in terms of wave resonances [START_REF] Mougel | On the instabilities of a potential vortex with a free surface[END_REF]. In a recent asymptotic analysis of the base flow, [Iga, 2017] gives an in-depth characterisation of the internal and boundary layers: velocity profiles, scaling laws. . . This analysis assumes that the flow is laminar which is a strong assumption. In fact, instabilities can appear at Reynolds number values much smaller than that where the rotating polygons settle in. In most experiments, having a typical length scale of the order of 10 cm and using water as a working fluid, these instabilities break the rotational symmetry of the velocity field long before the surface even starts to deform perceptibly. An azimuthal wave propagates in the same direction as that of the disk [START_REF] Young | Period-doubling route to chaos for a swirling flow in an open cylindrical container with a rotating disk[END_REF], Lopez et al., 2004, Poncet and Chauve, 2007, Kahouadji et al., 2010]. As the disk rotation speed is increased, turbulence develops on these large scale structures. At sufficiently high disk-rotation speed, the Froude number is no longer small: the free surface is strongly deformed and may form polygons. There is thus a transitional regime that bridges both types of instabilities.

The aim of the present work is to gain knowledge on this transitional regime by studying the mean axisymmetric flow and the fluctuations using direct numerical simulations and experimental measurements. The flow parameters chosen in the present work mostly match those of the experiment of [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] for which velocity measurements are available. The paper is structured as follows: the experimental setup and the numerical methods are introduced in section 4.2. In section 4.3, the steady axisymmetric flow solution is presented, and the effects of Froude and Reynolds numbers are discussed. Section 4.4 describes how unsteadiness and three-dimensionality affect the axisymmetric mean flow structure. Experimental measurements are carried out to assess the relevance of the numerical simulations. Discussions and comparisons with previous available models are given in section 4.5.

Configuration and methods

The configuration is represented in figure 4.1(a). The cylindrical tank of inner radius R is filled with a layer of liquid with density ρ l and kinematic viscosity ν. At rest, the layer has a uniform height H above the bottom disk. When the disk rotates at angular speed Ω, the surface deforms and may become non axisymmetric and time-dependent.

Our experimental set-up (see figure 4.1b) is a modified version of the one used in [START_REF] Moisy | Experimental and numerical study of the shear layer instability between two counterrotating disks[END_REF]. The tank is made of Plexiglas and the heavy brass alloy rotating disk is driven by a DC brushed motor with a tachometer closed-loop speed control. The inner radius of the cylindrical tank is 140 mm and the gap between disk and cylinder is 0.7 mm. Laser Doppler Velocimetry (LDV) is carried out using a Dantec BSA system fastened to a horizontal displacement platform. It consists in a continuous argon-ion laser with wavelength 660 nm and power 25 mW. Hollow silver-coated glass spheres with a diameter of 10 micrometers are seeded in the fluid. The temperature is monitored during the experiment and its variation did not exceed 2 • C.

In the following, variables R, RΩ and ρ l are respectively chosen as reference length, velocity and density scales to make all quantities dimensionless. The flow inside the tank is characterized by four dimensionless parameters, namely the aspect ratio G, the Reynolds number Re, the Froude number F r indicative of the deformation of the surface and the Weber number W e which quantifies inertia with respect to surface tension:

G = H R , Re = R 2 Ω ν , F r = RΩ 2 g , W e = ρ l Ω 2 R 3 σ , (4.2.1)
where g stands for the gravity acceleration and σ for the surface tension between gas and liquid. All physical properties of the liquid are assumed constant. In all cases investigated here W e is found large, which indicates that surface tension effects hardly affect the surface profile and will not be considered in the following.

Beside experiment, two in-house numerical codes are used to study axisymmetric flow states. The first one, ROSE (ROtating Surface Evolution) computes steady states of the liquid phase with surface tension. Using cylindrical coordinates (r, θ, z) with the origin at the center of the bottom disk, the steady Navier-Stokes equations are written as

                       0 = rω - ∂ 2 ψ ∂r 2 - 1 r ∂ψ ∂r + ∂ 2 ψ ∂z 2 0 = ∂ ∂r ( 1 r ∂ψ ∂z ω) - ∂ ∂z ( 1 r ∂ψ ∂r ω) - ∂ ∂z ( Γ 2 r 3 ) - 1 Re ( ∂ 2 ω ∂r 2 + 1 r ∂ω ∂r - ω r 2 + ∂ 2 ω ∂z 2 ) 0 = ∂ ∂r ( ∂ψ ∂z Γ) - ∂ ∂z ( ∂ψ ∂r Γ) - r Re ( ∂ 2 Γ ∂r 2 - 1 r ∂Γ ∂r + ∂ 2 Γ ∂z 2 ) (4.2.2)
using Stokes' streamfunction ψ, azimuthal vorticity ω and angular momentum Γ. These latter quantities are linked to velocity

V = (V r , V θ , V z ) via V r = 1 r ∂ψ ∂z , V z = - 1 r ∂ψ ∂r , ω = ∂V r ∂z - ∂V z ∂r , Γ = rV θ . (4.2.3)
Mainly for numerical reasons, surface tension is introduced in the normal stress balance at air-liquid interface * through n(T air -T)n = σ∇ • n where T (resp. T air ) is the stress tensor on the liquid (resp. air) side, n = (-h , 0, 1) t /(1 + h 2 ) is the unit outwards normal vector at the surface z = h(r). The deformed meridional domain (r, z) ∈ [0, 1] × [0, h(r)] is adapted to the free surface and mapped to a rectangular one (ξ, η)

∈ [0, 1] × [0, G],
allowing for a Cartesian mesh of N r (resp. N z ) points in the radial (resp. axial) direction in this transformed coordinate setting [Kahouadji and Martin Witkowski, 2014]. The solution is obtained by iterations over the two following steps: (i) system (4.2.2) is mapped to the rectangular domain, discretized, and solved using a Newton-Raphson procedure with a prescribed distribution h(r); (ii) the velocity field obtained at convergence is inserted into the normal stress balance equation which is solved to yield a corrected h(r) distribution.

As a second numerical tool, we use the finite-volume DNS (Direct Numerical Simulation) code Sunfluidh that simulates 2D and 3D unsteady incompressible flows. The Navier-Stokes equations are discretized on a staggered grid with 2nd order accuracy in both time and space and the zero velocity divergence is ensured by an incremental projection method. Further details can be found in [START_REF] Tuerke | Nonlinear dynamics and hydrodynamic feedback in twodimensional double cavity flow[END_REF]. Sunfluidh also implements a level-set method to tackle interfaces.

Steady axisymmetric flow

This section is aimed at describing the structure of the steady axisymmetric solutions (system 4.2.2) and investigating the effects of varying Froude and Reynolds numbers. The results are obtained by ROSE using N r ×N z = 501×101, 401×201 and 401×401 equispaced grid points for aspect ratios G = 0.1856, 0.5 and 1 respectively.

In order to evaluate how spatial resolution affects numerical solutions, we consider the configuration G = 0.1856, Re = 10 000 and first check flat surface computations, by setting F r = 0 in ROSE and using the monophasic version of Sunfluidh. ROSE calculations were conducted on four different grids N r ×N z = 251×51, 501×101, 1001×201 and 2001×401 with equispaced grid points. Maximum or minimum values of the stream function, radial and axial velocity components were shown to vary less than 1.95%, 0.44% and 0.09% respectively from the first three grid systems to the finest one, so that the grids selected in the present paper for ROSE computations ensure sufficient precision. Sunfluidh calculations were conducted on the three uniform grids N r × N z = 128 × 32, 256 × 64 and 512 × 128. The maximum and minimum velocity components were found to differ by less than 2.6%, 0.95% and 0.16% respectively from those of the most resolved ROSE computation.

Validations on configurations with a deformed surface at F r = 1 were also performed. * . This was found to solve the convergence problem discussed in [Kahouadji and Martin Witkowski, 2014]. For ROSE, the solution is sought in a monophasic deformed domain (see section 4.2) on grid N r × N z = 501 × 101 while for Sunfluidh, the steady two-phase flow solution is determined in the domain (r, z) ∈ [0, 1] × [0, 0.5] on a regular mesh N r × N z = 256 × 128 from an initially rest state with a liquid-air interface at z = G. The extrema of V r and V z were found to differ only by about 1% between the two codes, which is fully satisfactory as these two techniques to deal with surface deformation are completely different. for G = 0.1856, Re = 10 000 and W e = 1263.6, for F r = 0.01 (red solid), F r = 0.2464 (blue dash) and F r = 1 (brown dash-dot). Radial distributions of (a) V θ and (b) V r at the free surface; radial distributions of (c) V θ and (d) V z at half the interface height z = h(r)/2; axial distribution at r = 0.8 of (e) V θ and (f) rV r , plotted as functions of z/h(0.8).

Effect of Froude number

is varied from F r = 0.01 (top row in figure 4.2) where surface deformation is negligible, to F r = 1 (bottom row in figure 4.2) where the surface is strongly deformed and almost touches the bottom disk at the center. For the Newton bucket (i.e. solid body rotation), the deviation h(r)-G is proportional to the Froude number. Using numerical simulations with an undeformed free surface, [START_REF] Piva | Steady axisymmetric flow in an open cylindrical container with a partially rotating bottom wall[END_REF] proposed a first order approximation for the surface elevation h(r) based on the normal stress balance, yielding a scaling proportional to F r. [Kahouadji and Martin Witkowski, 2014] observed and interpreted such a scaling for weaker deformations and moderate Reynolds numbers. Figure 4.3 shows the rescaled surface deformation (h(r) -G)/F r for the set of parameters of figure 4.2: this scaling remains here a fair approximation, even though at F r = 1 the surface deformation is of the same order of magnitude as the fluid layer thickness. The main effect of free surface deformation is thus to constrain the flow field: this is quantitatively shown in figure 4.4 where some velocity components at the surface, at half-depth or at fixed radius are extracted for different Froude numbers and almost superpose when z is rescaled by the local depth h(r). This feature is found to be robust for Reynolds numbers greater than Re = 1 000.

Effect of Reynolds number

As the Froude number influences the flow structure only marginally, flat-surface cases (h = G) are conveniently chosen to study the effects of Re at lower numerical cost. Figure 4.5 illustrates how the flow evolves as Re is increased from low values.

One first observes the gradual formation of boundary layers at the rotating disk and at the side wall, but also at the free surface (hereafter called top layer) and at the edge of the solid-rotation core (core layer). These layers then become thinner as Reynolds is increased. The selected grids ensure at least 15 grid points in each layer viewed in the meridional plane even at the highest Reynolds number Re = 19 500. An overshoot of the azimuthal velocity arises at mid-radius in the top layer, similar to the one observed by [Iwatsu, 2004]: the fluid locally spins faster than the disk at the same radius (see the contour of azimuthal velocity in figure 4.5 at Re ≥ 10 000) as a consequence of the inward flow convecting angular momentum with weak dissipation. This overshoot at this location is captured for any of the above grid resolutions.

The existence of two regimes prevailing at low and large Re is best illustrated on figure 4.6. The maximum ψ max of Stokes' streamfunction quantifies the strength of the meridional recirculation: at G = 0.1856, this quantity (figure 4.6a) is first found proportional to Re, and so do the maxima (rV r , V z ) max of rV r , V z (figure 4.6b). Indeed, in the viscous regime up to Re = 300, length scales are of order 1 as well as the azimuthal velocity; balancing the two last terms of the vorticity equation in system (4.2.2) leads to a scaling proportional to Re for the meridional flow.

For Re ≈ 1 500, this trend stops: (rV r , V z ) max saturate at a constant value close to 0.1. When the Reynolds number is further increased, the above structure non longer holds: an Ekman-like layer with typical thickness δ ∝ Re -1/2 forms on the rotating disk at the same time as other layers as will be described in section 4.3.3. In the Ekman layer, V r scales as r [von Kármán, 1921] and thus (rV r ) max is of order 1, independent of Re -same for (V z ) max in the side layer by flow conservation arguments. As δ ∝ Re -1/2 , the meridional recirculation strength ψ max ∝ δ × (rV r ) max is thus expected to decrease as Re -1/2 as well, due to the thinning of the layers inside which most of the circulation takes place. This is observed in figure 4.6a: for G = 0.1856, the interpolation gives an exponent of -0.49. This scaling is fully coherent with the in-depth asymptotic analysis of the bottom and side-wall layers performed by [Iga, 2017]. Figure 4.6 also shows the results for another aspect ratio, G = 0.5: the same scalings are observed. Quantities (rV r , V z ) max saturate around the same value 0.1, which is not surprising as the scalings for the Ekman boundary layer do not involve G. However the recirculation in the low-Re regime is far more intense at G = 0.5, which is due to a larger recirculation zone. This mechanically leads to a shift of the critical Reynolds number for the change of regime towards smaller values, around 400 for G = 0.5.

Flow structure at large Reynolds number

Figure 4.7 gives a more general picture of the steady flow obtained in the axisymmetric framework at large Reynolds number Re = 10 000 (i.e. above the critical Reynolds number introduced previously), obtained for aspect ratio G = 0.5 and G = 1.0. The main features already observed at G = 0.1856 are also present here: a solid body rotation (hereafter denoted as SBR) in the central part and a meridional recirculation (hereafter denoted as MR) at the periphery, with four layers. However the core layer has an intricate sinuous shape and feeds a lower region stretching down to the disk. Further increasing the aspect ratio to G = 1 (see figure 4.7a) shows that the size of the upper sub-cell changes only slightly: its axial extent increases from 0.2 to 0.25. Most of the increase in fluid depth is passed to the extension of the lower mostly z-independent region, while the bottom layer with an Ekman-like structure is hardly modified [Iga, 2017]. For G = 0.1856, this lower region is only slightly visible at the highest Reynolds number (see figure 4.5).

The key point to close existing theoretical models (see discussion in section 4.5) is to determine the radius r s of the boundary between SBR and MR regions [Tophøj et al., 2013, Fabre and Mougel, 2014, Iga, 2017]. This parameter can be tentatively determined from simulation results. To do so, we define r s as the location of the first maximum of V θ (r, z = G/4), as depicted in figure 4.4(c): we chose z = G/4 in order to get rid of the influence of the upper cell and the bottom layer. Figure 4.8 shows the decrease of r s as G is increased, which confirms previous studies [Iga, 2017, Tophøj et al., 2013]. This trend is a consequence of the balance between sidewall dissipation and energy injection at the rotating disk below the MR region.

The axial velocity component V z at r = r s has been plotted as a function of G -z in figure 4.9. Coordinate G -z has been used so that the free surface at G -z = 0 coincide for the three aspect ratios. It is observed that the size of the upper cell indeed weakly varies with G as stated above. Moreover, the graph reveals a flow reversal (V z < 0) similar to that observed in vortex-breakdown bubbles for rotating-lid experiments [Escudier, 1984]. This analogy has already been pointed out by [Iwatsu, 2004, Piva andMeiburg, 2005] or [START_REF] Herrada | Off-axis vortex breakdown in a shallow whirlpool[END_REF] and is sometimes referred to as "off-axis vortex breakdown".

The axisymmetric flow solutions obtained here are physically relevant up to a critical value of the Reynolds number for which the flow becomes unsteady. Using Newton's method implemented in ROSE allows to go beyond this critical value by continuation technique. However, exploring further these steady axisymmetric unstable branches of solutions at higher Reynolds number values is eventually limited by a nontrivial dynamical behavior. This is the reason why the curves in figure 4.8 could not be pursued beyond Re = 19 000-23 000. An alternative way is then to use unsteady computations of the three-dimensional flow.

Three-dimensional flow

For a given aspect ratio, there is a large range of Reynolds numbers for which the velocity field is yet three-dimensional and unsteady while the free surface deforms but remains almost axisymmetric. We focus here on this transitional regime that extends up to the occurrence of rotating polygons. The critical Reynolds number value that determines the lower bound of the transitional regime and the associated critical azimuthal wavenumber strongly depend on the fluid aspect ratio. These critical parameters were determined for G varying from 0.036 to 0.107 experimentally by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] and numerically by [START_REF] Kahouadji | Seuils de stabilité pour un écoulement à surface libre engendré dans une cavité cylindrique tournante à petit rapport de forme[END_REF] using linear stability analysis. For G = 0.25 and G = 2, [Lopez et al., 2004] performed a study using nonlinear numerical simulations and experiments. [START_REF] Cogan | Symmetry breaking and instability mechanisms in medium depth torsionally open cylinder flows[END_REF] extended the computations to the range G = 1.5 to G = 3.5 and [START_REF] Serre | Vortex breakdown in a cylinder with a rotating bottom and a flat stress-free surface[END_REF] studied G = 4. All the numerical simulations have assumed a flat horizontal interface.

In this section, we focus on the specific value G = 0.1856. Assuming a monotonic variation with aspect ratio, we can expect a critical Reynolds number between Re ∼ 1 450 given by [Lopez et al., 2004] for G = 0.25 and Re ∼ 10 000 given by [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] for G = 0.107. Therefore, we choose the values Re = 30 000 and 81 400 in the following as they lie far beyond threshold and are typical of the transitional regime. Furthermore, the value Re = 30 000 allows for comparisons between fully resolved simulations and experiments in which velocities are large enough for accurate LDV measurements.

At such Reynolds number values, the flow is turbulent with coherent structures (see movies 1-3 in the supplementary material of [START_REF] Yang | Axisymmetric rotating flow with free surface in a cylindrical tank[END_REF]). The present section characterises the axisymmetric mean flow and fluctuations in the transitional regime, and gives a description of the unsteadiness (section 4.4.3). -3 (resp. 1.98 • 10 -3 ) and one dashed contour at value 1 2 ψ min with ψ min = -1.12 • 10 -4 (resp. -0.924 • 10 -4 ). The locations of the points with maximum and minimum ψ are respectively indicated by the red and the blue cross.

Simulations: mean flow and fluctuations

A 3D simulation is performed at Re = 30 000 using Sunfluidh. The chosen numerical domain (r, θ, z) ∈ [0.4, 1] × [0, 2 3 π] × [0, 0.1856] covers only a part of the physical one: a solid body rotation is assumed for r < 0.4 which is not simulated, and only a third of the full azimuthal extent is considered using periodic boundary conditions along θ, a choice motivated by the prevalence of the m = 3 mode in our experiment (see section 4.4.2. The 3D mesh consists of 192 × 448 × 128 cells and is refined along the side wall, above the disk and below the free surface. The radial grid spacing decreases from δr = 4.17 • 10 -3 for r ∈ [0.4, 0.8] continuously to δr = 1.4 • 10 -3 at r = 1. Along the axial direction, the grid spacing varies from δz = 3.1 • 10 -3 in the bulk to δz = 2 • 10 -4 at the disk and at the surface. The grid spacing in the azimuthal direction is uniform. We impose a CFL constraint of 0.2, so that the time step is δt ≈ 7.7 • 10 -4 . Sunfluidh also allows for unsteady 2D simulations. A first simulation is led in 2D in the domain (r, z) ∈ [0.4, 1] × [0, 0.1856] at Re = 30 000 using the same grid spacing as above along r and z. Even in the 2D framework, the flow is found unsteady, but it evolves to a statistically permanent regime with mean velocity V 2D , defined by averaging V over a large time interval ∆t = 200.

In order to save some computation time, the 3D simulation is started from an instantaneous 2D field in the permanent regime. After a time period T e , the 3D flow reaches another statistically permanent regime. A mean velocity V 3D is defined by averaging over the same time interval ∆t and over the entire azimuth of the computational domain ∆θ = 2 3 π:

V 3D (r, z) = 1 ∆t∆θ ∆θ 0 Te+∆t Te V(r, θ, z, t) dt dθ . (4.4.1)
Due to Reynolds stresses, the mean flow may differ significantly from that arising from system (4.2.2). Figure 4.10 shows the meridional circulations obtained when considering either V 2D or V 3D at Re = 30 000. Some subtle changes are observed in the sub-cell shapes. Concerning the value of ψ max however, an extrapolation of the axisymmetric steady solution (curve in figure 4.6) would lead to predict a value 1.93 • 10 -3 at Re = 30 000, in

Experiments: mean flow and fluctuations

Velocity measurements have been conducted with a water layer of initial height H = 26 mm so that G ≈ 0.1856 as in [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] (see their figure 12). Two angular speeds of the bottom disk have been chosen in order to match either the Reynolds number used in the numerical simulation of previous section or the Froude number of Bergmann's experiment. In the meridional plane, mean and rms values of the azimuthal and axial velocity components were measured along one or several lines of constant z.

The first case investigated has an angular speed Ω = 1.53 rad/s, leading to Re = 30 000 and F r = 0.0335 with no noticeable surface deformation. Figure 4.12 displays the mean and rms azimuthal velocity component just below the surface as well as the azimuthal and axial components at several heights. The corresponding quantities obtained via the 3D DNS have also been plotted. Concerning the azimuthal component, the agreement is good except in specific regions: near the side boundary at mid-height (figure 4.12b), and near the maximum at the surface, where experimental results do not show any overshoot (figure 4.12a). As for the axial velocity component ad mid-height (figure 4.12c), negative values are found experimentally around r = 0.83 and numerically around r = 0.6, indicating that the shape of the meridional circulation significantly differs from the experiment to the numerics.

We also measured the mean axial velocity profile at other heights (see figure 4.12d). The range for which axial velocity measurements can be performed is limited as both incident beams in the tank must remain within the fluid layer axial extension. For 0.38 ≤ z/G ≤ 0.69, the profile was found independent of z and so do the rms values (not shown). The mismatch between numerics and experiments remains unexplained at this stage. A possible reason could be the fact that surface pollution, especially when water is used as a working fluid, may affect the flow dynamics; its mathematical modeling is still an open issue [START_REF] Peaudecerf | Traces of surfactants can severely limit the drag reduction of superhydrophobic surfaces[END_REF], Moisy et al., 2018]. However, the orders of magnitude for both maximum mean values are in agreement, which is satisfactory since experimental measurements of such a small velocity component are extremely delicate and numerical simulations are demanding.

For the second case shown here, the rotation speed is higher Ω = 4.15 rad/s, which corresponds to Re = 81 400 and the same Froude number F r = 0.2464 as in [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] (their figure 12). Results are displayed in figure 4.13. At such a high Reynolds number, azimuthal velocity profiles are found independent of z in the range 0.192 ≤ z/G ≤ 0.75 investigated (the upper limitation is due to the strong surface deformation). Moreover, the axial velocity profile at mid-height seems to be robust as it is found to be very close to the one measured at Re = 30 000.

Analysis of the unsteady flow

Sections §4.4.1 and §4.4.2 describe the mean axisymmetric flow in the transitional regime and quantify the fluctuation amplitude via rms values. Hereafter we briefly describe the structure and the frequency spectrum of these fluctuations on the case G = 0.1856, Re = 30 000, for which both experiments and simulations are available.

Even though the free surface remains almost flat in this configuration, a flow structure is evidenced when Kalliroscope flakes are added to water. Figure 4.14 reveals a mode with azimuthal wavenumber m = 3, with a large amplitude at the periphery of the solid body rotation zone. This robust structure rotates in the same direction as the disk, however at a lower angular speed (see movie 1 in supplementary material in [START_REF] Yang | Axisymmetric rotating flow with free surface in a cylindrical tank[END_REF]). A quantitative characterization is performed by extracting the gray levels from successive video images along the circle of radius r = 0.8. A spatio-temporal picture is obtained and plotted in figure 4.15(a). The diagram contains inclined stripes from which we can deduce an angular phase velocity of = 1/0.98 rad • s -1 = 1.02 rad • s -1 . This corresponds to a pattern rotating at an angular velocity close to 2/3 of that of the disk Ω = 1.53 rad • s -1 , or equivalently to a frequency f = m /(2π) close to twice the disk frequency f d = Ω/(2π). This is best seen in the spectral domain when applying a 2D Fourier transform to the spatialtemporal signal: in figure 4.15(b), the maximum of the spectrum is located at m = 3 and f /f d = 1.98. Note that a peak at (m, f /f d ) = (1, 1) is also visible, presumably associated to the disk rotation.

Numerical simulations bring some useful information on the minimum ingredient to capture this instability. Figure 4.16(a), shows the evolution of the azimuthal velocity at some probe location for the chained 2D and 3D simulations with Sunfluidh. As mentioned earlier, the temporal mean is barely shifted from 2D to 3D. However the rms value is strongly enhanced in the 3D simulation. The spectrum associated to the 3D established state shows a peak at f /f d = 2.07 in full agreement with the above experimental determination. The above findings on the frequency spectrum were eventually corroborated by experimentally measuring the azimuthal velocity using a LDV probe located at the very same location. Signal and spectrum are displayed in figure 4. 16(b). A marked frequency peak is found at the same frequency -this peak was also found at all other locations investigated for both axial and azimuthal velocity components. The fluctuation amplitude differs between experiments and numerics depending on the probe location. For the case of figure 4.16(b), the discrepancy is relatively strong, but, as illustrated in figure 4.12(a)-(b), better agreement can be found at other locations. On the whole, the 3D simulation restricted to the sector r ∈ [0.4, 1] and θ ∈ [0, 2 3 π] is able to catch most of the instability features of the flow in the transitional regime.

This azimuthal modal structure at such relatively large Reynolds number is likely to be related to the primary instability at threshold. Indeed, the angular speed of the observed pattern normalised by that of the disk is close to 2/3, a ratio that we found to be robust with respect to changes in Re; at threshold, this same ratio (within roughly 10%) has been reported from linear stability analyses, whatever the critical azimuthal wavenumber [Kahouadji, 2011].

Discussion and concluding remarks

Which base flow should be used for linear stability analysis is a major question as it is usually a necessary requisite for accurate predictions of bifurcation thresholds. This study gives a better understanding of the base flow structure. Increasing the rotation speed leads, for Reynolds numbers exceeding a value of typically several hundreds, to a regime where most of the meridional circulation takes place within boundary layers. Increasing the rotation speed also enhances surface deformation, with a deviation h(r) -G almost proportional to the Froude number. Numerical simulations reveal that this deformation however preserves the overall flow structure. Unsteady 3D DNS results show little and localised influence on the mean fields. Along with rms data, these characterisations provide a clear picture of the flow at high Reynolds numbers, far beyond the primary instability threshold. A comparison between the numerical simulations and existing models is of interest. [Iga, 2017] suggested that the rotating polygons grow on a laminar base flow at least in the bottom and side wall boundary layers, while the model presented by [Tophøj et al., 2013] assumes fully developed turbulence. This may deeply impact the angular momentum exchanges between the flow and the walls. As a matter of fact, to close the models, the balance between the angular momentum fed by the bottom disk and dissipated by the side wall is used to find the position of the solid body rotation radius r s . In our experiment and simulation, the flow is clearly turbulent.

Nevertheless, in the transitional regime, the bottom and side-wall layers studied by [Iga, 2017] are confirmed and do not differ significantly from those obtained by the mean flow as the rms is quite small in those layers. In particular the side-wall layer could be both accurately measured and simulated with an excellent agreement (figure 4.12a,b). Overall, our study s . Blue dash-dot: the TMBF model in [Mougel, 2014] with correction of the dissipation, leading to values r TMBF s = 0.57, 0.49 and 0.43. Without corrections, this latter model would predict values 0.65, 0.58 and 0.52 respectively for the three aspect ratios.

shows that the boundary layers are mildly modified by turbulence and most of its structure is captured in the laminar regime. In addition to the boundary layers described by [Iga, 2017], simulations reveal a top layer along the free surface, an overshoot of the azimuthal velocity at the intersection of the top and the core layers, and a sinuous deformation of the core layer. These features of the steady two-dimensional solutions are also present in the mean unsteady 2D and 3D flow simulations. The 3D unsteadiness is found to weaken the sharp overshoot present in 2D (figure 4.11a) without smoothing it out completely. The presence of the top layer in simulations implies that the azimuthal and radial velocity profiles at the surface are not representative of the profiles in the bulk. This could be of importance as most simplified models assume a z-independent azimuthal velocity profile obtained through surface measurements.

We now compare our computations with two existing models namely those presented by [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF] (hereafter the TMBF model) and [Iga, 2017] for three aspect ratios. In these existing models, as stated above, r s is a key parameter. Once r s is known, the azimuthal velocity field can be derived. An additional matching between the solid body rotation region and the peripheral region, assumed potential, was derived by [Iga, 2017] in order to remove the unphysical slope discontinuity at r s † . Azimuthal velocity profiles are compared in figure 4.17. For G = 0.1856, there is a fair match between the 2D steady simulation and both theories. The radius r s is predicted within a few percent when compared to the value r s = 0.58 given in figure 12 of [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] and our experiment where r s is slightly below 0.6. As the aspect ratio increases, the results issued from theory deviate significantly from the numerical results: Iga's model underestimates both the SBR radius and the azimuthal velocity amplitude while the TMBF model overestimates them. Note that †. It's worth noting that there is a typo in the formula (3.64) of [Iga, 2017] where the term

2Ωδs/ 2 √ 2/π + 1 should be 2Ωδs/ 2 √ 2/π + 1 .
in order to improve the theoretical models, their authors suggested to tune the respective dissipations at the side wall and at the disk. In his thesis [see details in chapter 6 of [Mougel, 2014]], Mougel introduced a friction coefficient ratio, thus increasing the dissipation by a factor 3 along the side wall. Similarly, in order to match experimental results, [START_REF] Iga | Axisymmetric flow in a cylindrical tank over a rotating bottom. Part II. Deformation of the water surface and experimental verification of the theory[END_REF] reduced the friction along the disk. In figure 4.17, we only applied such correction to the TMBF model. The numerical simulations and experiments presented in this study are a step toward a more complete understanding of the base flow prevailing in the high Reynolds-number regime that leads to the growth of rotating polygons in water. Specifically, the azimuthal velocity profiles and the structure of the flow in the meridional plane are given in detail. Visualizations and LDV measurements reveal a robust rotating azimuthal modal structure rotating at angular speed 2Ω/3 while the free surface deformation barely deviates from axisymmetry.

This contrasts with rotating polygons where the angular speed is close to Ω/3 [START_REF] Bach | From Newton's bucket to rotating polygons: experiments on surface instabilities in swirling flows[END_REF] and the deviation of the free surface along the azimuth is significant. Evaluating the role of the meridional flow, usually ignored in instability analyses, as well as the potential influence of azimuthal modal structures on the emergence of rotating polygons is left for future studies.

Chapter 5

Unsteadiness: spin-up, bulk instability and rotating polygons

In this chapter, two topics are considered. The first one is a benchmark for the transient regime in rotating flows with axisymmetric or almost axisymmetric free surface. Two different flow types are considered: Newton's bucket and the bottom-driven flow. In both flow types, large enough angular speeds are chosen to create a strongly deformed air-liquid interface. Using experimental measurements and numerical calculations, we accurately characterize the transient spin-up evolution as well as the established state in such two flows. Comparisons of the results obtained from both techniques allow for cross validations.

The second topic is a quantitative study on 3D instability, more specifically on the relation between bulk instabilities and rotating polygons in the bottom driven flow. Experiments are conducted in two experimental set-ups using water and G8W2 as working liquids. Velocities and surface deformation have been measured, which is an original contribution. Finally, the broad variety of patterns observed in this system is illustrated by presenting a regime that, to our knowledge, was not reported yet in the literature: a mode 1 rotating in the opposite direction with respect to the disk rotation.

Newton's bucket: spin-up

Newton's bucket has been discussed in Chapter 3 as it has been used as a benchmark for FTP measurements. In this flow, the velocity field and the surface deformation are both axisymmetric. Using non-dimensional variables, the steady parabolic free surface can be predicted by the formula:

h(r) = G + F r 2 r 2 - 1 2 (5.1.1)
The whole fluid rotates with the cylinder as a solid body once the flow is established leading to v θ = r. This configuration with very basic analytical solutions is chosen in order to assess our ability to measure precisely both a surface deformation and an azimuthal velocity field.

It is a necessary step before considering the bottom-driven flow. By the way, note that 98 Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons the later flow has also a core in solid-body rotation, see Chapter 4. As for the transient regime, its main characteristic is the establishment time, which is the duration for the flow to reach the established regime from rest. The established regime can be either steady or periodically-evolving. Since Newton's bucket has been proven to be stable with respect to several wave families through viscosity damping (see [START_REF] Mougel | Waves in Newtons' bucket[END_REF]), the established state of Newton's bucket is steady. We will take advantage of the well-established properties of Newton's bucket flow to check the limits of achievable comparisons between numerical solutions and experimental measurements as both approaches have their own limitations. In the following Sunfluidh is used for unsteady simulations and ROSE for steady simulations, both in 2D configurations.

Spin-up time: simulations

During spin-up, local variables (velocity components,...) have a transient behavior that strongly depends on the spatial location. In order to characterize the whole spin-up with a single establishment time, a global quantity such as kinetic energy may be used.

Energy is difficult to measure in experiments. A numerical approach is thus undertaken to simulate the transient regime. The temporal evolution of volume-averaged kinetic energy is recorded. The flow establishment is estimated using a residual, defined here as

resid(t) = 1 Vol Vol i (V n+1 i -V n i ) 2 δt dv (5.1.2)
where Vol stands for the volume of the simulated physical domain and V i for the velocity components (V r , V θ , V z ). The flow is considered as established at t = τ s when the residual resid(t) drops below a given threshold . This criterion detects steady states as well as periodic states oscillating at low amplitude, if a suitable value of is chosen (see below).

Simulations with enforced flat surface

In the first place, we will calculate the configurations with a flat surface where symmetry boundary conditions are imposed. We first consider the aspect ratio G = 0.25 and study the influence of Reynolds number on flow establishment. In this situation, only the liquid phase domain

[r, z] ∈ [0, 1] × [0, 0.25] is calculated on a regular grid of N r × N z = 256 × 64.
Figure 5.1 (a) shows the temporal evolution of the volume-averaged kinetic energy e c from t = 0 to τ s for different values of Reynolds number, with a threshold set at a level = 10 -10 . All the variables are presented in dimensionless form. Clearly, the establishment time increases with Reynolds number. The volume-averaged kinetic energy seems to evolve with a law e c (t) = e 0 (1 -exp(-t τ )) where e 0 is the energy of the solid rotation flow. This allows us to define a characteristic time parameter τ which is equal to the time it takes for the flow to reach 63.2% (= 1 -e -1 ) of its steady-state kinetic energy e 0 . The values of τ are drawn versus Reynolds number in figure 5.1 (b) while the values of τ s are drawn in 5.1 (c). The same procedure is applied for the configuration G = 0.5 using a grid 128 × 64. The G=0.25 G=0.5 (a, n) 1st regime 2nd regime 1st regime 2nd regime τ (0.017, 0.984) (0.254, 0.538) (0.0325, 0.974) (0.4, 0.552) τ s (0.196, 0.882) (1.806, 0.475) (0.462, 0.794) (4.179, 0.465) Table 5.1 -Fitting dependency of the establishment time with Reynolds number using a law a • Re n . Data extracted from numerical simulations of Newton's bucket spin-up with flat surface. Two aspect ratios G = 0.25 and G = 0.5 are investigated.

extracted values of τ and τ s are also plotted in figures 5.1(b, c). Results of simulations with an even coarser grid 64 × 32 for G = 0.5 show no significant differences compared to the numerical results using the grid 128 × 64.

From the figures 5.1 (b,c), we can identify two different regimes for the establishment time as function of Reynolds number (see [Wedemeyer, 1964, Watkins andHussey, 1977]). Function of type a • Re n is used to fit the points of the two regimes for the two concerned aspect ratios.

The first regime takes place at low Reynolds number where the flow is dominated by viscous diffusion whose characteristic time is Ω -1 Re (resp. Re) in dimensional (resp. nondimensional) form. The second regime is dominated by Ekman pumping during the transient regime with a dimensional time scale Ω -1 Re 1/2 (resp. Re 1/2 ). These exponents are indeed roughly obtained in the numerics (see table 5.1).

Simulations with deformed surface

The previous study shows that the flow development is governed by the Ekman transport even for moderate Reynolds numbers of a few hundred. Let us focus now on this Ekman pumping regime and investigate if the presence of a free surface affects the establishment time. For a fixed Re = 2 500, the Froude number is varied from 0 to 2 in two-phase flow simulations with the initial interface set to z = 0.25. A grid of N r × N z = 128 × 64 is used to discretize the two-phase domain [r, z] ∈ [0, 1] × [0, 0.5], of which the upper-phase is set to have air viscosity and density. It is inferred from figure 5.2 that the establishment time increases with the Froude number with an almost constant slope, and then saturates. It's worth noting that, for G = 0.25, the surface at center begins to touch the bottom disk at F r = 1. A closer look at the evolution of resid(t), we found that its value begins to oscillate when the liquid surface touches the bottom from F r = 1 and eventually the top boundary from F r = 1.5. The configurations with G = 0.5 are computed on a grid N r×N z = 128×128 for a two-phase domain [r, z] ∈ [0, 1] × [0, 1]. Similar to the behaviour observed at G = 0.25, τ keeps increasing as a function of F r for G = 0.5 before the disk dewetting which takes place at a value between F r = 1.75 and 2.

Overall, at large Reynolds numbers, one can consider that the Froude number does not modify the transient duration significantly. 5.2 -Parameters of three spin-up experiments for generating Newton's bucket, at small (SRS), medium (MRS), and large angular speed (LRS).

Surface deformation, local velocity: experiments and simulations

G8W2 is used as working fluid in Newton's bucket experiments. The main reason to use this liquid is that it has a dynamic viscosity about 60 times that of water, which makes comparisons between numerical and experimental results easier. Increasing viscosity decreases the Reynolds number and consequently shortens the dimensionless establishment time τ , as seen in the previous section. Actually, for the same angular speed (Ω), the dimensional establishment time τ /Ω will also be reduced in dimensional scale. We can give an example here. Assuming the law τ = 4.2 Re 0.465 for aspect ratio G = 0.5 with flat surface thus discarding the influence of the Froude number on the establishment time, using water in our small set-up (R = 7 cm) for Newton's bucket leads to an establishment time of 58 seconds when the angular speed is set to Ω =11.83 rad/s (this angular speed correspond to F r = 1 and Re = 57 583) while using G8W2 leads to an establishment time of 9.4 seconds for the same angular speed (Re = 1153). This reduced establishment time is however still large enough to allow for well resolved sampling of surface height and velocity evolution since our measurement tools have high sampling frequencies. There is another constraint from the mechanical point of view, that comes from the fact that the disk can not be started abruptly. Consequently, the angular speed of the disk is increased linearly from zero to a prescribed value within a certain duration ∆t r in experiments. However, this duration is not negligible compared to the establishment time, as will be shown later. To circumvent this drawback, we have implemented the linear variation of the angular speed used in the experiment into the unsteady simulations by Sunfluidh. This is crucial for subsequent comparisons between numerical and experimental results.

We filled initially 34.2 mm (G = 0.49) of G8W2 in the cylindrical container of the solidbody rotation set-up as shown in Chapter 3 section 3.3.3.4. Three different angular speeds were tested, as shown in table 5.2. The open cavity is accelerated to prescribed speeds within ∆t r = 0.5 s in the spin-up experiments. The liquid temperature was T = 28.6 • C during the whole campaign of measurements. As a reminder, the numerical simulations take into account both air and liquid phases. A regular grid N r × N z = 128 × 128 is used to compute the whole domain [r, z] ∈ [0, 1] × [0, 1] with the air-liquid interface set at z = 0.49. The CF L number is imposed to be 0.3, so that the dimensionless time step is δt = 1.15 × 10 -3 , 1.71×10 -3 and 2.13×10 -3 for the simulations corresponding to respective experiments SRS, MRS and LRS. 102Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons SRS: small rotation speed Figure 5.3 shows the transient evolution of surface height and azimuthal velocity for a case: H = 34.2 mm and Ω = 7.84 rad/s. The Gocator line sensor is set to work at a sampling frequency of 100 Hz for measuring surface elevation and will also be set at this frequency for the following measurements unless specified otherwise. This sampling frequency is much larger than that of any phenomenon we wish to describe in the present work. For a better visualization, the points of experimental measurements in figure 5.3(a) have been re-sampled each 0.25 seconds instead of using the original temporal resolution of 0.01 seconds. In figure 5.3(b), azimuthal velocity measurements by LDV at two representative probe locations are shown. The start of acquisition with the LDV system is synchronized with the start of disk rotation in order to capture precisely the transient regime and to compare with numerical results. As shown in figure 5.3(b), the LDV system captured very low speeds at the beginning of the spin-up. The limited number of measurements in the early transient stage is due to the fast flow dynamic at the measured points. We can see a good agreement between the experimental and numerical results for both the surface and velocity evolutions. Figure 5.4 (a) shows the numerical and measured profiles of surface height every ∆t = 15.68 (2s) from t = 0 to t = 78.4 (10s) as well as the steady profile. The experimental surface profiles are also sub-sampled with one tenth of the original resolution (0.15 mm). The numerical solutions predict reasonably well the surface evolution with an accuracy better than 1%. The profile of V θ at mid-height is presented in figure 5.4 (b). This profile also holds at other z values since the whole fluid rotates at the same speed as the disk at steady state. The rms values are also shown but they are very weak as the flow is laminar. The uncertainty for the radial positions in LDV measurements is also shown, which is related to the error on the identification of the rotation axis in laboratory frame. As the LDV system has a probe volume measuring about 3 mm long, the exact probe position should also have an uncertainty of almost 3 mm. As shown in figure 3.6(a) of Chapter 3, a maximal error of 3 mm for the localisation of measuring positions can be induced by an erroneous identification of the axis position. MRS: medium rotation speed By increasing the angular speed of the open cylindrical cavity to Ω = 11.65 rad/s, we obtain the case MRS with larger surface deformation. Similar figures of transient evolution and steady regimes as the previous case are also drawn for this case, see figure 5.5-5.6. Good agreement between numerical and experimental results is also achieved.

LRS: large rotation speed

The experiment with an angular speed of Ω = 14.56 rad/s, denoted by LRS is now considered. Similar figures of transient evolutions and steady regimes as the previous two cases are also drawn for this case (see figures 5.7-5.8), leading to the same conclusion.

These experimental results give estimations of the establishment time, respectively ≈ 100, 150 and 200 for SRS, MRS and LRS in their respective non-dimensional scales. These values are compatible with the ones obtained by the numerics, respectively τ s = 136.5, 176.3 and 202.8.

Bottom-driven flow: spin-up

The investigation on spin-up dynamics in Newton's bucket serves as a preparatory step to assess numerical and experimental tools now for the free surface bottom-driven flow in a cylindrical tank with stationary side wall. After a transient evolution, the flow obtained in Newton's bucket experiment is always stable and steady [START_REF] Mougel | On the instabilities of a potential vortex with a free surface[END_REF], while the bottom-driven flow is able to develop three-dimensional instabilities, such as rotating polygons and switching/sloshing. Former studies on these instabilities concentrated mainly on the description of phenomena and determination of critical conditions for instability occurrence or transitions by experimental methods. [Vatistas, 1990, Jansson et al., 2006, Vatistas et al., 2008, Abderrahmane et al., 2009, Iga et al., 2014]. Comprehensive phase diagrams of this kind of instabilities have also been drawn from experimental results, see [Jansson et al., 2006, Bach et al., 2014, Iga et al., 2014]. However, a precise knowledge of the rotating free surface flow itself is still missing. In order to obtain a complete and precise characterization of the flow displaying instabilities, we shall conduct direct simulations and quantitative measurements. In the first place, we will investigate the axisymmetric configuration, a regime before 3D instabilities settle in. We have discussed the steady flow structure of axisymmetric configuration in Chapter 4 and will continue in the following section to study the transient regime from rest to established flow, in situations where the surface although strongly deformed, remains nearly axisymmetric.

We have already mentioned the fact that numerical simulations are always restricted by limitation in Reynolds number. The small value of the kinematic viscosity of water leads to large Reynolds numbers when the regime of rotating polygons is reached. The configuration at Re = 30 000 presented in the Chapter 4 corresponds, in the water experiment to a small angular speed which induces a weak surface deformation. Furthermore, two-phase flow simulation, necessary for strongly deformed interfaces, has large discontinuities in density and viscosity across the interface. This is also quite challenging for numerical simulations. 106Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons Instead of implementing turbulence models or employing prohibitive fine grids, we have decided to act on experiments by changing the working fluid from water to a more viscous liquid, for example the glycerine-water mixture as we did for Newton's bucket. A more viscous liquid helps to reduce the Reynolds number (which makes numerical simulation easier) yet maintaining a large enough Froude number so that the free surface is strongly deformed. Fortunately, the increase in viscosity does not prevent the rotating polygons from eventually settling in, as shown in [Jansson et al., 2006, Abderrahmane et al., 2017]. This is discussed in section 5.4.

We show here that Sunfluidh is able to compute two-phase flow in axisymmetric configuration at moderate Reynolds number and large Froude number, even for very large deformations associated to dewetting around the center of the disk. Numerical and experimental results for two different aspect ratios, G = 0.25 and G = 0.5, are considered in the two following subsections.

Aspect ratio G = 0.25

We will firstly present the results for cases with G = 0.25 realized experimentally in the large set-up then present the results obtained for cases with the same aspect ratio but in the small set-up. A potential interest to use both experimental set-ups of different sizes is to check whether the three non-dimensional parameters we have chosen (G, F r, Re) are the only relevant parameters.

Experiments in the large set-up

An initial height H = 35 mm of G8W2 is poured into the large set-up so as to obtain the aspect ratio value G = 0.25. The angular speed of the bottom disk was set to Ω = 7.854 rad/s. The protocol for the experiment is the same as for Newton's bucket spin-up experiments in section 5.1.2: a constant acceleration of the bottom disk from 0 to a prescribed rotation speed during ∆t r = 0.5 s. The non-dimensional parameters are respectively: G = 0.25, F r = 0.88, Re = 3063 (and W e = 3153). Axisymmetric unsteady simulation under Sunfluidh is also carried out in a computational domain (r, z) ∈ [0, 1] × [0, 1] meshed by 128 × 128 cells. Initially, the interface is set to z = 0.25. A CFL number of 0.2 is enforced, corresponding to a time step δt = 1.7×10 -3 . A steady computation under ROSE in the domain [0, 1]×[0, 0.25] mapped by a regular curvilinear grid 401 × 101 is also run.

The flow structure at steady state, obtained from numerical simulations by the two codes, is shown in figure 5.9 which presents isocontours of the stream function and of the azimuthal velocity. The two codes converge almost towards the same steady solutions in the liquid phase. The sub-cell division (SBR and MR) discussed in Chapter 4 is also found for this case, with an overshoot zone just below the free surface around r = 0.5. ). 21 equally spaced contour levels between 0 and ψ max = 0.006 (red '+') and one dashed contour at 1 2 ψ min with ψ min = -0.0024 (located in the air phase, blue '×') are plotted for the Sunfluidh results. Similarly for solutions of ψ computed by ROSE in which ψ max = 0.006 (magenta '+') and ψ min= -1.89 × 10 -4 (magenta '×' in the liquid phase) as only one phase is calculated by ROSE. V θ values are presented using 21 equally spaced iso-contours between 0 and 1 for both numerical codes.

108Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons state obtained using ROSE and Sunfluidh, which are fully consistent. At steady state, the surface deviation at the axis represents 55% from its initial height. The probes for LDV measurements were chosen in such a way that they remain always in the liquid phase during the whole spin-up experiment (see their locations in figure 5.10(b)). Figure 5.11 shows the comparison between numerical and experimental time evolutions of V θ at several chosen probes located at one of the four heights investigated, z = 0.22. A good agreement is achieved between numerical and experimental results which requires the synchronization between the LDV system and the disk rotation. Spectra were also checked for the velocity signals of each probe in established regime. The corresponding power spectra of the established state are plotted in figure 5.12. Only peaks at multiples of the disk rotating frequency with weak amplitudes were found for all the measured probes, as shown in figure 5.13. As a matter of fact, the largest peak of (power) spectra for each probe was found at f = f d . In non-dimensional scale, the largest spectral amplitude is smaller than 4 × 10 -3 and the largest power spectral peak is smaller than 4 × 10 -4 for the height z = 0.22. The same order of magnitude was also found for probes located at other heights. This validates that this experiment is either axisymmetric or that the 3D deviation from axisymmetry does not modify strongly the base flow.

Figure 5.13 shows the azimuthal velocity V θ as function of radius r in steady regime. Apart from the local difference at the overshoot zone between numerical and experimental results, good agreement is generally achieved. .13 -Comparison between numerical and experimental steady azimuthal velocity profiles V θ (r, z) at four different heights z = 0.05 (red), 0.14 (blue), 0.17 (magenta) and 0.22 (black), for the case G = 0.25, F r = 0.88, Re = 3063. Colors are coherent with those used for the probes in figure 5.10(b). Solid curves refer to numerical solutions while symbols are for measurements.

Experiments in the small set-up

Similar cases with G = 0.25 have also been studied using the small experimental set-up, filled with G8W2 to an initial height H = 17.5 mm. The speed of the disk is controlled as shown in figure 3.2 of Chapter 3 and we set here ∆t r = 1 s, ∆t p = 60 s and ∆t d = 10 s. The Gocator line profiler measures the surface deformation h(r, θ 0 , t) along a diameter during each run of experiments. Different angular speeds have been used, as listed in table 5.3. Surface shapes have also been monitored visually to detect the emergence of specific patterns. The second column of table 5.3 lists the visually observed azimuthal wavenumber m for each case investigated. .5 mm (G = 0.25) and different angular speeds in the small experimental set-up. For all cases, the duration of constant acceleration for the disk is ∆t r = 1 s. The corresponding non-dimensional parameters and the main characteristics of surface deformations determined from observation and spectral analysis are also listed. The mode named -1 corresponds to a mode 1 rotating in the direction opposite to that of the disk. Mode column and the last three columns will be exploited in sections 5.3 and 5.4.

N • Ω (rad/s) mode F r Re f max /f d f p /f d |H s (
Here we can characterize the dynamics of the free surface for the different angular speeds in the space-frequency domain. A run of Gocator measurements of the surface deformation during a spin-up experiment yields a two dimensional matrix containing discrete values of h(r, t). We can extract the submatrix corresponding to the established regime (t > τ s ), h s (r, t). The Fourier transform of h s (r, t) along the time dimension gives rise to a complex matrix H s (r, f ). The peaks of its module |H s (r, f )| in the spatial-frequency domain (r, f ) reveal the dominant frequencies and their spatial locations along the free surface. More specifically, we introduce the dominant frequency f m (r) for which |H s (r, f )| is maximum over the whole spectrum at a fixed r. One may also characterize the global maximum (r max , f max ) for which |H s (r, f )| is maximum over all locations r and frequencies f . The fourth column of table 5.3 lists the ratio of f max /f d . The last column shows the values of |H s (r max , f max )| scaled by the cylinder inner radius R.

A priori, f m is a function of r and could differ from one radius to another. In present experiments however, f m is generally found to be uniform with respect to r. We may thus identify this frequency as that of an existing instability wave f i . Together with the wavenumber m obtained from visualization, we can calculate the angular phase velocity 112Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons 2πf p of the wave with f p = f i /m. The ratios of f p and the disk rotation frequency f d for the cases with observable patterns are listed in the 7th column of table 5.3. Besides, f max is generally the most frequent value in the vector f m (r). It should be however noted that exceptionally, a localized region can be identified with a frequency different from that of the instability wave. This will be described later. It is worth noting now that a peak associated to the disk rotation can be present but remains weak with respect to peaks associated to an instability.

We can see that |H s (r max , f max )| of the first 5 cases are all very small. Consequently, the cases from N • 1 to 5 are all considered as almost axisymmetric cases even patterns of mode 3 with weak amplitudes along the azimuth were observed. These weak deviations from axisymmetry will be discussed in detail in section 5.3. Now we begin to discuss the axisymmetric and almost axisymmetric cases, corresponding to the first seven cases in table 5.3.

Axisymmetric flow with wet center region

Since the Reynolds numbers of the (almost) axisymmetric cases in table 5.3 remain moderate, corresponding numerical simulations could be conducted in the axisymmetric framework. The time evolution and the steady radial profiles of measured and numerical surface height are shown in figure 5.14. The maximum difference between experimental and numerical results is found to be 1 mm for case 3 at Ω = 13.11 rad/s. As to be discussed in section 5.3, this may be linked to the fact that this case is the "least steady" case with a maximum value of |H s (r max , f max )| amongst all axisymmetric or almost axisymmetric cases, while axisymmetric numerical solutions remain steady * .

In the following, we concentrate on the axisymmetric case with a dry center, that is case 7 in table 5.3.

Axisymmetric flow with dry center region

When the angular speed of the disk exceeds a critical value depending on the initial liquid height, a dry region occurs at the axis. Such regime is called DC or DP by [START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF]. For case 7 at Ω =19.00 rad/s, the deformed surface is axisymmetric as mentioned above, with a dry region of radius r c = 0.32 around the axis.

This latter experiment has been completely reproduced by numerical simulations with the same non-dimensional parameter values. The computational domain (r, z) ∈ [0, 1] × [0, 0.5] is meshed by N r × N z = 128 × 64 cells with the liquid-air interface initially set to G = 0.25. Figure 5.15 shows the comparison between numerical and measured surface heights as functions of time. The agreement is excellent, with a maximal difference of about 0.5 mm (for r = 0.5), not far above the experimental precision. It is concluded that Sunfluidh is * . This is no longer the case for the highest angular rotation speed Ω = 22.2 rad/s, where the experiment reveals a m = -1 mode. Unsteady axisymmetric simulations give rise to an instability where surface oscillations have the same order of magnitude as the experimental ones but at a very different frequency 1.05f d instead of 0.64f d found in experiments. Of course, we cannot expect to capture the three-dimensional mode -1 by an axisymmetric simulation, but it is interesting to note that this case is already unstable with respect to a mode 0. able to well predict the surface evolution for this axisymmetric case in which a dry region is present. The flow structure for this case is shown in figure 5.16. The core layer with a sinuous shape mentioned in Chapter 4 has disappeared, probably due to the strong surface deformation and to the small Reynolds number. The surface layer is observed and confirmed experimentally (see later figure 5.19c) by the inverted hump formed from r = 0.6 to 0.75.

Azimuthal velocities at different positions in the liquid have also been measured. Since LDV is able to measure velocity at one point at a time, the same spin-up process needed to be repeated as many times as the number of locations. This led to lengthy measurements during which the temperature of the liquid varied from 22.4 • C to 27.9 • C. Temperature is measured by a digital thermometer with a precision of 0.5 • C. Table 5.4 lists physical properties of the liquid and the corresponding non-dimensional parameters for three typical temperatures encountered ( [Cheng, 2008]). The values at T = 20 • C are also given for comparison. Despite the variation of the Reynolds number due to temperature variation, the surface deformation did not show significant deviation among the different runs of the spin-up. Simulations at the four Reynolds numbers were all conducted to evaluate the influence of temperature on the flow velocities.

Transient azimuthal velocity at numerous locations with the liquid are displayed in figure 5.17. Firstly, we observe that the numerical results reproduce better the transient evolution at the two lower heights z = 0.018, 0.125 than at z = 0.24. The simulations run at several Reynolds numbers indicate that the temperature variation cannot be the reason for the mismatch at z = 0.24. Indeed, the time evolution of the volume-averaged energy is shown 114Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons Noting that the difference between numerical and measured curves is maximum near (r, z) = (0.49, 0.125) suggests that the overshoot region at r = 0.49 is more pronounced in the numerical solution than in the experiment. This is clearly visualized on the steady velocity profiles shown in figure 5.19: a good agreement between the numerical and experimental velocity is achieved except in the overshoot region around r = 0.5 at z = 0.125 (figure 5.19b).

Aspect ratio G = 0.5

We consider then an aspect ratio value G = 0.5. A layer of 35 mm of G8W2 is initially filled in the cylinder of the small experimental set-up in the bottom-driven flow configuration (see figure 3.1 of Chapter 3). The control of disk is the same as the one used for the experiments with G = 0.25 in the small set-up (see section 5.2.1.2).

Six disk speeds have been tested, which are listed in table 5.5 together with the corresponding non-dimensional parameters. The first three speed levels Ω = 11.65 rad/s, 13.11 rad/s and 13.84 rad/s lead to an axisymmetric surface deformation, at least up to measurement precision. However, when the speed of the bottom disk is increased further, stable rotating ellipses are observed, for instance in the last three cases of table 5.5. The discussion on the 3D patterns is postponed to section 5.4.2. It should be noted the bottom disk are always wet for all these six cases. 5.5 -Cases of spin-up experiments in the bottom-driven flow with aspect ratio G = 0.5. The fluid is a 35 mm deep layer of G8W2 poured in the small experimental set-up. Six disk angular speeds are tested. The first three cases display an axisymmetric free surface while rotating ellipses are observed for the last three cases . The rotating ellipses are exploited in section 5.4.2.

Unsteady axisymmetric simulations of the first three axisymmetric cases are performed using show the temporal evolution of surface heights at r = 0.25, 0.5 and 0.75 for each case. The difference between the numerically predicted and the measured height is marginal for the two outer positions r = 0.5 and 0.75 in all the three axisymmetric cases, while that at r = 0.25 is larger but still remains relatively small (1.6%). This is more evident when the steady profiles are plotted, (see figure 5.20(d-f)). Actually, the largest difference between measured and numerical steady profiles of the three cases is about 0.5 mm. Steady surface profiles from ROSE simulations are also plotted in figure 5.20(e-f) to confirm the solutions of unsteady simulations by Sunfluidh. Both numerical curves are indistinguishable one from the other.

From the experimental profiles, we observe that the Gocator sensor detects a local cusp at r = 0.33, a radial position below which the difference between numerical and measured profiles become larger. As a matter of fact, the cusp is a measurement artifact created by the presence of a local transparent circle (no light diffusion) at this position, due to a local depletion of pigments (TiO 2 ), as shown in figure 5.21(a). Despite the incorrectly measured height at this fake cusp position, the observed transparent circle is physical and can be related to the velocity field. 120Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons

Bottom-driven flow: bulk instability

As explained in section 5.1.2 for Newton's bucket, changing the working fluid from water to a more viscous liquid G8W2 allows to reduce the Reynolds number value from O(10 4 -10 5 ) to O(10 3 ) in our experimental setups. The latter Reynolds number is a typical value to trigger the primary bifurcation of non axisymmetric instabilities. When using G8W2 as compared to experiment with water, the critical angular velocity for which the flow becomes unstable is thus increased. With negligible surface deformation, i.e. small Froude number, it has been shown thanks to numerical simulations [START_REF] Kahouadji | Seuils de stabilité pour un écoulement à surface libre engendré dans une cavité cylindrique tournante à petit rapport de forme[END_REF] that the dynamics of the free surface does not account for the instability mechanism. We may called them bulk instabilities. To our knowledge, these instabilities have always been observed in experiments with an almost flat interface ( [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF]). The goal of this section is to check the existence of such instabilities when the interface is significantly deformed and if these instabilities may be tracked down by an azimuthal dependence of surface deformation in the azimuthal direction. We will firstly have a look at the signals of the surface heights and their spectra for cases listed in table 5.3. The measured surface heights at r = 0.5 and r = 0.75 are plotted with respect to time in figure 5.22 (a, b) for all the spin-up experiments. Only for cases 1, 6 and 7 at Ω =7.85 rad/s, 15.54 rad/s and 19 rad/s, the surface height seems to become steady, with variations smaller than the measurement precision range of 0.4 mm. By constrast, the heights for cases 2-5 and 8 at angular speeds Ω = 11.65 rad/s to 14.57 rad/s and Ω =22.2 rad/s evolve to a periodic state, with an amplitude captured by the measurement tool. However, the spectrum H s (r, f ) for cases 6-7 at Ω = 15.54 rad/s and 19 rad/s reveals that f m (r) takes the value f d characteristic of the disk rotation for most radial positions, while the overall maximum amplitude H s (r max , f max ) remains very small, as shown in table 5.3. We conclude that cases 6-7 are stable to azimuthal perturbation and are thus axisymmetric.

5.3.1

For case 1 at Ω = 7.85 rad/s, the value of f m (r) is found close to 1.74f d for most radial positions. Figures 5.23 for cases 2-5 at angular speeds between Ω = 11.66 rad/s and 14.57 rad/s (see also column f max /f d in table 5.3).

It has been shown by [Kahouadji, 2011] that the waves growing in the bulk instability have a phase velocity 2πf p close to 2 3 Ω. We may thus infer that this instability dominates the dynamics for the cases 1-5 of table 5.3, as f p = f m /m is also found close to 2 3 f d . In the present study, we have succeeded in detecting the signature of the bulk instability on the free surface, through the precise measurement of the very small deviations from axisymmetry.

Despite this tridimensionalization of surface height for the first five cases in table 5.3, the variation range of the surface deformation in θ (equivalent to variation in time at a given radius) remains less than 1.5 mm, as shown in figure 5.22. The spectral amplitudes of these waves (discarding the harmonics) correspond to the value |H s (r max , f max )| listed in the last column of table 5.3. The variation in θ will be clearly shown later in figure 5.25. Therefore, they can be considered as almost axisymmetric cases, introduced in Chapter 4. However, the oscillation range of surface height for the case at 22.2 rad/s is almost 4 mm (coherent with |H s (r max , f max )| = 0.04R = 2.8 mm as shown in table 5.3 which discards the other harmonics). It is one order of magnitude larger than those of the first five cases. It has its principal peak located at 0.64f d . Since it is mode 1 as observed by visualization, its phase velocity is also 0.64Ω which is similar to that of the first five cases. But the most important fact is that it rotates in the opposite direction of the bottom disk. We will denote this case 122Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons as "monogon-1" in section 5.4.3.

We can try to illustrate the above change of flow regime, determined by the phase velocity, in the diagram of (Re, F r). When only the bottom disk's speed Ω is chosen to modify nondimensional parameters' values, both the Froude and Reynolds number evolve together with Ω. More specifically, the couple (Re, F r) increases along a parabolic line when Ω increases, as indicated by the dashed-red line in figure 5.24. From bottom to top along this parabolic line, we get through three different regimes, as shown in figure 5.24. Firstly, we encounter almost axisymmetric flows where the bulk instability (BI) settles in. After passing through axisymmetric flows with dry center, likely to be a barrier, the flow takes up rotating waves with large amplitude, specially monogon-1 here.

Figure 5.24 -Representation of experimental cases on G8W2 with G = 0.25 from table 5.3 into the diagram of (Re, F r), denoted by the symbols. Filled symbols for wetted cases and empty ones for cases with a dry disk center, except the point of monogon-1, marked specially by a red circle filled with black color. Three different regimes are identified. The filled triangles in the regime of almost axisymmetric with BI refer to bulk instability of mode 3.

The start of the first regime is not determined in this thesis. For cases with flat surface at G = 0.25, the linear stability analysis shows that the first bulk instability occurs at Re c = 3500 with a critical azimuthal wavenumber m = 2 and the mode 3 becomes unstable at Re c = 4600, see [Faugaret, 2014]. However the small-amplitude waves of wave number m = 3 observed in cases 1-5 of table 5.3 have the same characteristics as the bulk instability: same direction of rotation, speed lower than that of the disk. However, the Reynolds numbers of these five cases are all smaller than 3000, thus well below the critical Reynolds number of mode m = 3 in the configuration with flat surface. The conclusions drawn from the stability study with a flat surface must be used with caution: large surface deformations may affect the thresholds of the bulk instability. Furthermore, even for tiny surface deformation, there is a large discrepancy between the numerical prediction and the experimental observations [Faugaret, 2014, Lopez et al., 2004]. The effect of the Froude number on the bulk instability will be studied in future work.

Deviation from axisymmetry surface deformation measurements

In order to construct a three-dimensional map of the surface deformation using the unsteady Gocator line measurements, we have to do some post-processing on the data in established regime h s (r, t) in a manner as explained in the following. The instability frequency f i is first determined from the inspection of the dominant frequencies in f m (r). Given the azimuthal wavenumber m deduced from visual observation, we can deduce the angular phase velocity of the rotating waves, which is 2πf p = 2πf i /m. Since the patterns rotate rigidly at their phase velocity [Jansson et al., 2006, Bach et al., 2014], we can use the transformation θ = 2πf p t to deduce the instantaneous surface elevation map h(r, θ) from the data h s (r, t) over one period 1/f p . We can then introduce the deviation from the time-averaged (or equivalently from the azimuthally averaged) profile:

δh(r, θ) = h(r, θ) - 1 2π 2π 0 h(r, θ)dθ. (5.3.1)
Using this technique, we can reconstruct a snapshot of the surface deformation h(r, θ) in the established regime and deduce the spatial perturbation map δh(r, θ) as defined in equation (5.3.1). The perturbation maps δh(r, θ) for cases 1-5 and case 8 are shown in figure 5.25(a-e). The difference between the weak amplitude of bulk instability waves and the strong amplitude of the monogon-1 of case 8 is clear. Regarding the variation of this amplitude among cases 1-5, as Ω is increased: it reaches its maximum value at disk speed around Ω = 13.11 rad/s and then decreases as the surface moves closer to the disk until an axisymmetric surface is recovered for case 6 at Ω = 15.54 rad/s. This will be also observed later for water experiment in section 5.4.

Aspect ratio G=0.5

For cases 1-3 of table 5.5 at G = 0.5, no surface deformation associated to rotating waves is visible, contrary to the experiments with G = 0.25. In order to check if the bulk instability is present before rotating polygons occur, we focus on surface elevation and azimuthal velocity measurements for case 1 of table 5.5 (G = 0.5, F r = 0.97).

The space-frequency diagram of surface height for this case is displayed in figure 5.26. First of all, the amplitudes of the spectrum is generally small, given that (|H s (r max , f max )| = 0.001 (scaled with respect to R). At most radial positions, the dominant frequency f m (r) equals 2f d . Nevertheless, the whole diagram is dominated by a marked peak at (r max , f max ) = (0.32, 0.81). As shown in section 5.2.2, this is exactly the position of the cusp at the free surface, so the frequency 0.81f d is presumably linked to a local dynamics in the cusp region.

A further investigation on the velocity spectrum is necessary to confirm this assumption. Azimuthal velocities are measured at the probe locations shown in figure 5.27(a) at fixed radial position r = 0.48. The measured steady velocities are compared to the numerical ones in figure 5.27(b). There is only a qualitative agreement for this measurement, as a difference of about 10% is observed close to the disk. The difference reaches 20% near the surface presumably because it is located in the overshoot region. Figure 5.27(c) shows the power spectrum of the measured velocities in the steady regime, with marked peaks at f = 0.84f d . This frequency differs slightly from the frequency observed at the surface at the fake cusp location r = 0.32, which is 0.81f d .

It can be concluded that the bulk instability is present in this case, but the signature on the free surface cannot be measured (at least with our apparatus). The only frequency observable by the laser sensor is at the fake cusp, and may be a consequence of this instability... or not! 126Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons

Bottom-driven flow: rotating polygons

For three different configurations, we present and characterize the rotating polygon instability. We first focus on the aspect ratio G = 0.1856 for which the base flow has been studied in Chapter 4. As water is used as a working fluid, this configuration is quite similar to the one studied by [START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF]. However, the very large values of Reynolds number reached in the experiment when the rotating polygons appear, preclude a DNS as it would be prohibitive. We then fully document two other configurations with aspect ratios G = 0.5 and G = 0.25 for which rotating polygons can be observed with G8W2 used as working fluid. This high viscosity induces thus moderate Reynolds number, which will allow in the future for 3D numerical investigations.

Water experiments: rotating triangles

For the aspect ratio G = 0.1856, a linear stability analysis of a baseflow with flat surface predicts that the first instability occurs at a critical Reynolds number Re c = 6033 and is characterized by an azimuthal wavenumber m = 2 and a frequency of 1.15f d . A second bifurcation involves a mode 3 at Re c = 6622 with a frequency of 1.85f d .

A series of spin-up experiments have been done in the large experimental set-up at aspect ratio G = 0.1856 (26 mm water). The angular speed of the bottom disk has been varied in the range of 1.53 -14.57 rad/s. For each spin-up process, the disk was always accelerated within 0.5 seconds to the desired speed, then it was kept rotating at the prescribed speed during 120 seconds in order to get an established flow state. For this range of angular speeds, the Froude number varies over 2 decades from 0.035 to 3. Meanwhile, the Reynolds number increases roughly by a factor 10, from 30 000 to 290 000, thus always far above the critical Reynolds number of bulk instabilities under flat surface condition. However the Froude number is likely to affect the development of the bulk instabilities and the critical Reynolds number may vary when surface deformation is taken into account.

Figure 5.28 lists the different cases investigated and shows a top view of the flows at steady state visualized using Kalliroscope mica flakes. The bulk instability of wavenumber 3 observed for case (a) at Ω = 1.53 rad/s (Re = 30 035) without no evident surface deformation, is found robust with respect to Froude number changes, at least up to F r = 0.64. Even though the triangular pattern related to the bulk instability is clearly visible and linked to a symmetry breaking of the velocity field, the surface remains always almost axisymmetric in cases (a)-(c). For case d, a significant dewetted central region is present, and the strongly deformed surface still remains axisymmetric. When the angular speed is further increased (case e), a rotating polygon is triggered and settles in, as shown in figure 5.28(e).

The surface deformation is measured for the above cases along a fixed diameter using the Gocator line sensor. LDV has been used to measure the instantaneous azimuthal velocity at (r, z) = (0.8, 0.0928) (i.e. at mid height and at a radius large enough so that it is always immersed in water) during the spin-up and the established regime. c) with bulk instabilities all present a main peak around frequency f i ≈ 2f d (slightly decreasing when Ω increases) while for case (e) with the rotating polygon, it is situated at f i ≈ 1.13f d . Consequently, the angular speed ratio f p /f d is roughly 2/3 for the bulk instability and 1/3 for the rotating polygon. Concerning the intermediate case (d) at F r = 1.94, only a weak peak at f = 0.7f d is slightly visible in the spectrum.

It is interesting to note that the magnitude of the main peak varies as a function of the Froude number. As F r is increased, the amplitude of the saturated bulk instability reaches a maximum and thereafter decreases. An attempt of explanation is that the flow is more and more confined at the periphery and that the region where the instability develops is progressively replaced by air. At larger F r, the flow recovers axisymmetry. The size of the dewetted central region then increases until rotating polygons are triggered.

As done for the experiments on G8W2 with G = 0.25, we can also show all the cases (a)-(e) of figure 5.28 in the diagram of (Re, F r), as shown in figure 5.32. The axisymmetric regime with dry center separates the axisymmetric regime with bulk instability (BI) and the regime of rotating polygons at higher Froude numbers. The transition along the parabolic line here in figure 5.32 is similar to the one found in figure 5.24. For both of these two configurations with small aspect ratios (G = 0.1856, 0.25), the disk dewetting seems to play a role of barrier between the bulk instability and the rotating polygons if the monogon-1 is also considered amongst the family of rotating polygons. Figure 5.32 -Representation of experimental cases on water with G = 0.1856 from figure 5.28 into the diagram of (Re, F r), denoted by the symbols. Filled symbols for wetted cases and empty ones for cases with a dry disk center. Three different regimes are identified. In the almost axisymmetric regime with BI, the bold '×' refers to a BI of mode 2 and the filled triangle refer to BI of mode 3. The empty triangles in RP regime refer to rotating polygons of mode 3 with a dry center.

Aspect ratio G=0.5: rotating ellipses

We move now to the characterization of the rotating ellipse (m = 2) of cases 4-6 of table 5.5. The time-averaged surface height profile is plotted in figure 5.33(a), in order to get an order of magnitude of the surface deformation levels involved. show the perturbation maps δh(r, θ) of surface deformation for these three cases. The spectral analysis of surface deformation reveals marked peaks at f i /f d = 0.825, 0.794 and 0.77 respectively for these three cases. Therefore, they all have a phase velocity such that 2πf p /Ω ≈ 0.4. Different from the campaign of experiments with smaller aspect ratios (G = 0.1856, 0.5), the liquid surface does not touch the bottom disk for all the cases having G = 0.5 listed in table 5.5. When we present these cases in the diagram of (Re, F r) as shown in figure 5.35, there does not exist an axisymmetric regime with dry center. Rotating polygons are observed far before the free surface touches the bottom disk. As talked in section 5.2.2, the bulk instability was not observed or at least locally for the first three cases. We thus consider them always as axisymmetric regime. Different from smaller aspect ratios, the transition from the axisymmetric regime to RP at this large aspect ratio can not be related to dewetting of the bottom disk. Further work needs to be done for cases with small and large aspect ratios to have a complete characterization.

Aspect ratio G=0.25: a new state, the monogon-1

We finally characterize the state called monogon-1, which is observed rotating rigidly in the direction opposite with restpect to that of the disk in the experimental case 8 of table 5.3 with Ω = 22.2 rad/s. The surface deformation has been measured using the Gocator line sensor, and a map of the surface perturbation has been shown in figure 5.25(f). It has 132Chapter 5. Unsteadiness: spin-up, bulk instability and rotating polygons

Conclusion

In this chapter, the transient regime from rest was firstly characterized for Newton's bucket flow and then for the nearly axisymmetric bottom-driven flow, with wetted or unwetted center at moderate Reynolds number. On the experimental side, having large deformations at moderate Reynolds numbers required to use a mixture of glycerine and water as a working fluid. On the numerical side, axisymmetric unsteady simulations with strongly deformed interface were made possible. In all cases investigated, comparisons between numerical and experimental results for axisymmetric cases showed good agreement.

The bulk instability, generally observed in flows with a flat surface, was found robust to surface deformation up to a certain Froude number regardless of the Reynolds number. Precise surface measuring techniques allowed us to visualize the signature of the bulk instability on the surface deformation in some cases, which show actually slight oscillation of surface height in the azimuthal direction. But the interaction of the instability and the surface is still unclear. The study of critical conditions for bulk instability in flows with strongly deformed free surface should be the subject of future work. Eventually, this kind of instability pattern was found more detectable on the spectra of velocity field measurements than on spectra of surface deformation signals.

Rotating polygons are proven to take up even at large viscosity thus moderate Reynolds number, which will facilitate the study of its mechanism in the numerical approach. Precise measurements of surface height and velocity were conducted, which will serve as reference for future numerical simulations.

Finally, a new state has been found experimentally, namely a "mode 1" wave rotating in the direction opposite to that of the bottom disk, which had never been reported for the bottom-driven flow with fixed side wall.

Chapter 6

General conclusion

Conclusion

In this thesis, the rotating flow with deformable free surface was studied by numerical and experimental investigations. The flow is generated by a rotating disk at the bottom of a cylindrical tank which is partially filled with a viscous Newtonian fluid. When the fluid is driven into motion, one may consider intuitively the surface will be deformed axisymmetrically. But experiments using water thirty years ago have shown the axisymmetry can be broken under specific conditions of initial liquid height and rotating speed. In other words, the liquid surface can take up three-dimensional waves with large amplitude, such as rotating polygons (RP for short) and switching/sloshing. However, the mechanisms to develop these instability patterns is still not clear. In this context, my PhD is a step toward a better understanding of the mechanisms of these instabilities in a precise and complete way. To do so, direct steady and unsteady numerical simulations as well as precise experimental measurements were undertaken.

In Chapter 2, I presented the two in-house codes used to simulate the axisymmetric flow. The first code ROSE, using the Newton-Raphson method, allows to calculate at low cost (O(1) minutes) the steady solutions of configurations with a Reynolds number up to 20 000. Another code, Sunfluidh, is able to run unsteady simulation of axisymmetric two-phase flows by using the level-set method and its physical time cost is quite satisfactory (≈ 1 hour). Three-dimensional unsteady simulation is also realizable through the monophasic version of Sunfluidh in order to capture three-dimensional flow structure.

Chapter 3 is devoted to presenting the experimental techniques applied in the lab. Actually, two experimental set-ups that differ in size by a factor two for the radius have been used. The LDV technique is adapted to take into account optical corrections because of the refraction due to fluid and curved geometry in order to obtain accurate measurements of flow velocities. Then three techniques to measure the surface deformation are presented. The first two are ready-to-use laser sensors designed to measure deformations at a point and along a line with the laser triangulation method. In order to acquire the full three-dimensional shape of the free surface, the technique FTP (Fourier Transform Profilometry) was deployed in the lab. Satisfactory precision was obtained for FTP, which may allow real-time measurement of the dynamic surface shape and to track the whole process to develop 3D waves in an upcoming work.

Chapter 4 concentrates on the search of a precise base flow in a high Reynolds-number regime for future linear stability analysis of the rotating polygons or switching/sloshing. This high Reynolds number regime is related to water experiments and we identify it as a transitional regime bridging the BI (short for bulk instability) and RP occurring in bottomdriven flow. Firstly, axisymmetric steady solutions by ROSE reveal the structure of the rotating flow: a central solid body rotation region and a toroidal flow cell located at the periphery. When viewed in a meridional plane, the latter is found bordered by four layers, two at the solid boundaries, one at the free surface, and one located at the edge of the central region, which possesses a sinuous shape. The cell intensity and geometry are determined for several fluid-layer aspect ratios; the flow is shown to depend very weakly on Froude number (associated to surface deformation) or on Reynolds number if sufficiently large. Then this transitional regime demonstrates unsteady and three dimensional properties while the surface is still almost flat. Direct numerical simulations show that the averaged flow shares many similarities with the above steady axisymmetric solutions. At last, experimental measurements of velocities are given, which corroborate most of the numerical results. The numerical simulations and experiments presented in this chapter are a step toward a more complete understanding of the base flow prevailing in the high Reynolds-number regime that leads to the growth of rotating polygons in water. Specifically, the azimuthal velocity profiles and the structure of the flow in the meridional plane are given in detail. Flow visualization and LDV measurements reveal robust rotating azimuthal modal structures (considered as type BI) in this transitional regime. The phase velocity is measured precisely, being 2/3Ω while that of most RPs is close to Ω/3. As a step toward realistic numerical simulations for rotating polygons, we concentrate on the unsteady flow in the last chapter. Since a high Reynolds number leads to difficulties in numerical simulations, a direct three dimensional simulation of a RP experiment using water as working fluid would be very difficult to conduct. We circumvent this problem by replacing the working fluid from water to more viscous liquids, e.g., a mixture of (80%, by weight) glycerine and (20%) water, considering the fact that RPs can always occur even in a more viscous liquid. On the one hand, the decrease of viscosity from O(10 4 , 10 5 ) to O(10 3 ) allows easier simulation, the results of which will be compared directly to experimental results. On the other hand, the Froude number is maintained at a high level so that large surface deformation can always be created. We have tackled the unsteady flow with deformable surface by beginning with axisymmetric configurations as the latter remains unexplored up to now. In this context, numerical simulations and experiments are firstly carried out for Newton bucket flow as a preliminary step. Subsequently, we compare the numerical and experimental results for axisymmetric bottom-driven flows. Good agreements in terms of the surface height and velocities are achieved. This investigation of the unsteady axisymmetric flow with deformed surface is a very good benchmark for the transient regime (thus unsteady flow) and we are now fully confident on our ability to simulate such flow.

The BI observed in water in previous studies [START_REF] Poncet | Shear-layer instability in a rotating system[END_REF] occurs also in experiments with the glycerine-water mixture. For the aspect ratio G = 0.25, the small height deviation (O(0.1 -1) mm)under the effect of these modal structures along the azimuth are measurable by the precise line laser sensor. Spectrum analysis of these height perturbations shows always a phase velocity around 2Ω/3, the same as that in the experiments G = 0.1856 with water. It's worth noting that the velocity measurement by LDV can be relatively more precise to identify the phase velocity in signal's spectrum but it remains local -the velocity signal at well chosen position needs to be used to establish spectrum analysis. LSA (short for linear stability analysis) for flat surface flow also gives rise to close phase velocities for the same aspect ratios. This instability is thus considered to be BI. If the rotating speed keeps increasing, the Froude number will be large enough to dry the disk center. An axisymmetric state is recovered. A further increase of rotating speed will trigger rotating polygons. Here for the configuration of G = 0.25 with G8W2, we observed a stable surface wave of same order amplitude (O(10mm)) as that of most RPs but it rotates in the opposite direction of the disk at an angular speed of 0.64Ω. This new regime is to be investigated further in future work. As a matter of fact, a normal RP with an angular speed around Ω/3 in the same direction of the disk can be observed more commonly after disk-dewetting. That's the case for the experiments G = 0.1856 with water. The question about the relation between the BI and RP can not be answered completely in the present work. However, the characterizing of the two instabilities through precise measurements seems a good start. Preliminary results presented at the end of the last chapter shows their large difference in the phase velocity and oscillating amplitude (expressed in surface deformation or velocities).

According to the experimental phenomenon, the following scenario may be proposed. As stated in the literature (see [Mougel, 2014]), the rotating polygons seem to be controlled mainly by the Froude number while BI are controlled by Reynolds number. For the experiments we analyzed, when the rotating speed is increased, BI occurs first and rotation create large surface deformation. BI persists to a certain level of surface deformation before the surface touches the disk at bottom. During this period, the surface remains almost axisymmetric with slight perturbation along the azimuth under the effect of BI. When the surface gets close to the disk, BI is damped and disappears. After the disk dewetting, the rotating polygons are triggered and larger modes can be developed further if the disk speed keeps increasing. It seems that dewetting plays a role of barrier between BI and RP for small aspect ratios as BI and RP emerge respectively before and after the disk dewetting. However, the transition from BI to RP for configurations with higher aspect ratio like G = 0.5 is not fully understood or described. It should at least be noted that the scenario to describe the transition from BI to RP in configurations with large aspect ratio may be different. The rotating polygons can appear even before the free surface touches the disk. That's the case for the aspect ratio G = 0.5. Above all, the proposed scenario to describe the transition from BI to RP still needs more clarification. Characterization of the two instabilities in a large range of concerned parameters by means of experimental measurements are needed. As a matter of fact, summarizing in figure 6.1, the experiments having G = 0.25 that have been explored experimentally in this Ph. D shows clearly that much of the (Re, F r) space still needs to be covered. Even for this incomplete diagram, various instability patterns have been observed. However, it is worth noting that experimental phenomena could change in different instances of the same experiment. This is the case for the RP2 observed at F r ≥ 2.25 along the blue-dashed parabola in figure 6.1, which were sometimes replaced by axisymmetric states with dry a center. Also, instead of the almost axisymmetric state observed at F r = 0.97 along the same parabola, we could sometimes observe a switching between a circular state and an ellipse. For these delicate cases of (Re, F r) values, influence of the experimental protocol should not be ignored and more experiments should be done to get a statistical evaluation of them.

Perspectives

It took much time to develop the FTP system with a satisfactory precision and this system was not sufficiently used as we end up focusing mainly on axisymmetric comparison. Measurements of RPs or swithing/sloshing by FTP will be done in the future.

We have started investigating the relation between BI and RP in the present work. It will be worth continuing this characterization by means of precise measurements or accurate numerical simulations. For configurations with small aspect ratio, the axisymmetric regime seems to clearly separate the BI from the RP instabilities. However, in the nearly axisymmetric dry center, a small frequency close to the BI could be detected in the experiment. We thus may ask if the BI still remains in the flow with a very slight amplitude even after dewetting. The effect of confinement and dewetting due to this fairly small aspect ratio may also infer in the transition. Therefore, configurations with higher aspect ratios may be easier to investigate as both BI and RP are able to develop before disk-dewetting. Furthermore, configurations without dewetting are more friendly to numerical simulations as only the triple-contact line treatment on the disk is avoided. For the aspect ratio, G = 0.5, the effect of BI on the surface deformation is not so obvious. But velocity signals have shown a clear peak at 0.83Ω. In order to see how the phase velocity shifts from BI to RP, finer rising steps of the disk rotating speed need to be applied in order to track the spectrum evolution in a quasi-static way.

In the numerical aspects, we will implement a three-dimensional simulation of a RP experiment using the glycerine-water mixture since the Reynolds number is successfully decreased to get rid of the turbulence. Meanwhile, a LSA for the rotating flow with deformed free surface seems interesting and should also be done. For configurations with flat surface, LSA predicts correctly the critical Reynolds number for unstable modes. When the surface is deformed, the BI patterns appear. However the conditions to trigger it are modified by the surface deformation. Take the mode 3 for G = 0.25 as an example, linear instability analysis on a baseflow with symmetric flat surface shows Re c = 4600. However, the BI wave appeared already at Re = 1316 in the bottom-driven flow experiments with G8W2. So it is likely that surface deformation changes the condition for BI but how it is changed is still unclear. The linear instabilities for flows with deformed surface can be compared to the characteristics of measured BI and RP to see which one of the two instabilities can be predicted by linear instabilities. Critical values of non-dimensional parameters can also be obtained through LSA. (A.0.1) I will present in the following two sections: the time stepper with flat surface and deformable surface developed on the basis of the original steady code. In fact, the time stepper for simulating configurations with flat surface works and tests have been done to compare with temporal simulations by Gerris [Popinet, 2003] and Sunfluidh. But the one with deformed surface hasn't work yet. Since the code Sunfluidh was mostly used for unsteady simulations later, the attempt to transfer the steady ROSE to unsteady simulation was not pursued further.

A.1 Axisymmetric time stepper with flat surface

When considering a flat surface, we can eliminate all the derivatives of the surface height h as they are all zeros and replace h by the non dimensional initial height G in the system of (A.0.1) and obtain the equations (A.1.1). The spatial discretization is exactly the same as the one in the steady part. First order forward scheme is used for time discretization and an hypothesis V (t = -∆t) = 0 is assumed to start the first time step. And then velocities are solved using Newton-Raphson method at each time step. Algorithm can be clearly understood in the figure A.1.

                           ξω - ∂ 2 ψ

A.2 Axisymmetric time stepper with deformed surface

This section followed the previous one to present the algorithm for extending the original ROSE to unsteady simulation of rotating flows with deformed surface, always in axisymmetric configurations.

+ (1 -h 2 )ω - 2(1 -h 2 ) r ∂ 2 ψ ∂r 2 + 2 r 2
∂ψ ∂r ::::: .2.6) Projected into the curvilinear coordinates, the boundary condition about ω can be written as:

- 2h 2 r 2 ∂ψ ∂r = 0 - 4h r ∂ 2 ψ ∂r∂z - 2 ḣ r + (1 -h 2 )ω - 2(1 -h 2 ) r ∂ 2 ψ ∂r 2 - 2h 2 r 2 ∂ψ ∂r = 0 ⇐= ḣ = V z -h V r = - 1 r ∂ψ ∂r - h r ∂ψ ∂z - 4h r ∂ 2 ψ ∂r∂z - 2 ḣ r + (1 -h 2 )ω - 2 r ∂ 2 ψ ∂r 2 + 2h 2 r ∂ 2 ψ ∂r 2 - 2h 2 r 2 ∂ψ ∂r 2h 2 (ω-1 r ∂ 2 ψ ∂z 2 ) = 0 ⇐= ω = 1 r ( ∂ 2 ψ ∂r 2 - 1 r ∂ψ ∂r + ∂ 2 ψ ∂z 2 ) (1 + h 2 )ω - 2 r ∂ 2 ψ ∂r 2 - 4h r ∂ 2 ψ ∂r∂z - 2 ḣ r - 2h 2 r 2 ∂ 2 ψ ∂z 2 = 0 ω = 2 r(1 + h 2 ) ∂ 2 ψ ∂r 2 + 2h ∂ 2 ψ ∂r∂z + h 2 ∂ 2 ψ ∂z 2 + ḣ (A
ω = 2 ξ(1 + h 2 ) ∂ 2 ψ ∂ξ 2 - Gh h ∂ψ ∂η + ḣ (A.2.7)
The original Navier-Stokes equations in the primitive form (2.2.5) can be developed in the curvilinear system: - .2.11) In summary, we can see the difference in terms of the Navier-Stokes equations and the boundary conditions at south border (the free surface) between the unsteady and steady cases, as shown through the red terms in the equations (A.0.1), (A.2.5), (A.2.7) and the normal stress balance equation (A.2.10), compared to the equations (2.2.10), (2.2.15), (2.2.20) and the equation (2.2.27). Furthermore, the pressure derivative along ξ is also much different

                                   ∂V r ∂ξ + V r ξ -η h h + ∂V r ∂η + G h V z η = 0 ∂V r ∂t + V r ∂V r ∂ξ + G h V z -η h h V r ∂V r ∂η - V 2 θ ξ = - ∂P ∂ξ + η h h ∂P ∂η + 1 Re ∇2 V r - V r ξ 2 ∂V θ ∂t + V r ∂V θ ∂ξ + G h V z -η h h V r ∂V θ ∂η + V r V θ ξ = 1 Re ∇2 V θ - V θ ξ 2 ∂V z ∂t + V r ∂V z ∂ξ + G h V z -η h h V r ∂V z ∂η = - ∂P ∂η + 1 Re ∇2 V z (A.
1 1 + h 2 h 2 ∂V r ∂ξ -h ∂V z ∂ξ = - 1 1 + h 2 (-h ḣ -h h V r ) (A

  Figure 1.1 -Schematic representation (a) of the experimental set-up of the bottom-driven rotating flow in which (b) a rotating triangular pattern is observed (LIMSI experimental set-up).

  Figure 1.2 -The hexagon at the Saturn's north pole and its surroundings. Image in natural color acquired by Cassini wide-angle camera on Nov. 27, 2012.

  Figure 1.3 -(a) Satellite image of Typhoon Mangkhut on September 12, 2018 (UTC). (b) Naruto Whirlpools. Pictures from Wikipedia.

  Figure 1.4 -Configurations of rotating flow in closed cylindrical container. (a): rotor-stator, (b)co-rotation, (c): counter-rotation.

Figure 1

 1 Figure 1.5 -Scheme of configuration for bottom-driven flows.

  Figure 1.6 -Patterns of bulk instability observed experimentally by [Poncet and Chauve, 2007]. (a) G = 0.0714, Re = 49260 and F r = 0.09; (b) G = 0.0714, Re = 24630 and F r = 0.023; (c) G = 0.0429, Re = 61575 and F r = 0.14; (d) G = 0.0429, Re = 34893 and F r = 0.045.

Figure 1 . 7 -

 17 Figure 1.7 -Top view of polygons by [Jansson et al., 2006]: (a) triangle, (b) square, (c) pentagon.

  Figure 1.9 -(a) Phase diagram established in water experiments by[START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] in the plane ( √ F r, G). (b) Summary done by[START_REF] Iga | Various phenomena on a water vortex in a cylindrical tank over a rotating bottom[END_REF] for the previous experimental studies, represented in the plane ( √ F r, G). In this latter graph, we added the points relative to mode 3 and axisymmetric states discussed in[START_REF] Bergmann | Polygon formation and surface flow on a rotating fluid surface[END_REF] for comparison.

  extended

Figure 1

 1 Figure 1.10 -Reproduction of figures 11 and 12 of [Bergmann et al., 2011]. (a) Side view and (b) surface velocity of the axisymmetric case AD: G = 0.2, F r = 9.1, and Re = 8.1 × 10 5 . Quantity r 1 is the radius of the dry center region and corresponds to the notion r c in the present thesis. r 2 = R = 194 mm and Ω = 21.5 rad/s. (c) Side view and (d) surface velocity of the axisymmetric case AW: G = 0.1856, F r = 0.2464, and Re = 132 860. r s is the transition point between the solid-body rotation and potential flow in a Rankine vortex. In the two graph, symbols refer to their PIV measurements while solid curves correspond to model prediction;

Figure 1

 1 Figure 1.11 -Figure reproduced from figure 1 of[START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF]. Theoretical models of the base flow used for linear stability analysis: (a,b) Rankine vortex proposed by[START_REF] Fabre | Generation of three-dimensional patterns through wave interaction in a model of free surface swirling flow[END_REF] for wet and dry-composite cases and (c) potential flow by[Tophøj et al., 2013] for dry potential case. The top graphs represent surface shape and the bottom graphs show the modeling of azimuthal velocity in each cases. For these models, the fluid follows the rotation of the disk at angular speed rΩ in the inner solid body rotation (SBR) region and supposes a potential flow (PF) with flow circulation Γ in the outer region bordered by two thin turbulent boundary layers (represented by the red bands in the top graphs).

Figure 2

 2 Figure 2.1 -Surface-adapted (a) and curvilinear (b) grids for spatial discretization.

  Start h = G; (ψ, ω, Γ) = 0 compute Jacobian matrix A and second term b from N-S equations b max < ξ 1 solve AU = b (ψ, ω, Γ) = (ψ, ω, Γ) + U Steady (ψ, ω,

Figure 2 . 2 -

 22 Figure 2.2 -Diagram of the algorithm for the code ROSE.

Figure 2

 2 Figure2.3 -Grid influence on numerical results of simulations by ROSE for the configuration G = 0.1856, Re = 15 000. The velocity component and the stream function are all taken at the center of the liquid body, (r, z) = (0.5, 0.0928) and the total Energy is calculated on the whole domain.

Figure 2 . 4 -

 24 Figure2.4 -Grid influence on numerical results of simulations by ROSE for the configuration (G, F r, Re) = (0.248, 1.496, 1047) where curvilinear grid is used for spatial discretization as the surface is strongly deformed.

Figure 2 Figure 2 . 6 -Figure 2

 2262 Figure 2.5 -Evolution of velocity components at r = 0.5, z = 0.5 for the configuration G = 1.0 , Re = 900 calculated by three different codes (ROSE, Gerris and Sunfluidh): (a) V r , (b) V θ , (c) V z .

Figure 2 . 8 -

 28 Figure 2.8 -Flow velocities at steady state for the configuration G = 0.25, F r = 1.23, Re = 2063 calculated by Sunfluid: (a) V r with 11 equispaced iso-contours between [V r,min(-0.18), 0] and 11 equispaced iso-contours [0, V r,max (0.155)]. (b) V θ with 21 equispaced iso-contours between 0 and 1. Probes are also shown here. (c) V z with 11 equispaced iso-contours between [V z,min (-0.08), 0] and 11 equispaced iso-contours [0, V r,max (0.23). The maxima and minima positions of V r and V z are marked by red '+' and blue '×'. Iso-contours from ROSE results are also presented by magenta lines with the same contour levels.

  Figure 2.9 -Mesh used for 3D unsteady simulation by Sunfluidh in a domain (r, θ, z) ∈ [0.4, 1] × [0, 2/3π] × [0, 0.1856] for the configuration G = 0.1856, F r = 0, Re = 30 000. (a): meridional plane (r, z) (b): horizontal plane (r, θ). One fifth of cells are shown.

Figure 2

 2 Figure 2.10 -Mesh spacing in radial (a), azimuthal (b) and axial (c) directions of the grid used for 3D unsteady simulation by Sunfluidh of the configuration G = 0.1856, F r = 0, Re = 30 000.

  0.1856]. Grid dependence was studied for these two codes on the basis of this configuration. ROSE calculations were conducted on four different grids N r × N z = 251 × 51, 501 × 101, 1001 × 201 and 2001 × 401 with equispaced grid points. Maximum or minimum values of the stream function, radial and axial velocity components were shown to vary less than 1.95%, 0.44% and 0.09% respectively from the first three grid systems to the finest one. A direct comparison amongst the results on these four grids are shown in figure2.11(a). Sunfluidh (monophasic version) calculations were conducted on three uniform grids N r × N z = 128 × 32, 256 × 64 and 512 × 128. The maximum and minimum velocity components were found to differ by less than 2.6%, 0.95% and 0.16% respectively from those of the most resolved ROSE computation. Some velocity profiles from numerical solutions on the three grids are shown in figure2.11(b) together with those of the ROSE results on the grid 501 × 101.

Figure 2 Figure 2

 22 Figure 2.11 -(a): Results of steady simulations by ROSE on four grid systems N r × N z are plotted. Left: Azimuthal velocity at the surface as function of r; Right: axial distribution of rV r at r = 0.8 plotted as function of z/G. (b) Results of steady state obtained from unsteady axisymmetric simulations by Sunfluidh using three grid systems N r × N z are plotted. Same velocity profiles as in (a). The magenta curves from ROSE calculation on the grid 501 × 101 are also included as cross comparisons.

  Figure2.13 -Reproduction of the streamlines as shown in the figures (6, 8, 13) of[START_REF] Herrada | Off-axis vortex breakdown in a shallow whirlpool[END_REF] for G = 0.5, Re = 300, 1000, 2000, and the corresponding Froude numbers are F r = 9, 100, 400 and Weber numbers are W e = 1.26, 14, 56. The iso-contours of ROSE are plotted at values (i/15) 2 • ψ max (solid lines) and (i/15) 2 • ψ min (dashed lines). ROSE found respectively ψ max = (0.0128, 0.01, 0.078) and ψ min = (-2.63, -1.7, -1.03) × 10 -4 (in the liquid phase) for Re = 300, 1000, 2000. Sunfluidh results for Re = 2 000 are also shown in (c) and it found ψ max = (0.078) in the liquid phase and ψ min = (-0.0059) in the air phase. The interfaces calculated by Sunfluidh and ROSE are superposed.

Figure 2

 2 Figure 2.14 -Comparison of ROSE and Sunfluidh results for the configuration G = 0.1856, F r = 1, Re = 10 000. Same representation as in figure2.12 is used except in this figure: ψ max = 0.0035 and ψ min = -1.83 × 10 -4 are both found in the liquid phase for ROSE while ψ max = 0.0038 and ψ min = -0.0016 are both located in the air phase for Sunfluidh. V z,max = 0.136 (0.134) for ROSE (Sunfluidh) and V z,min = -0.086 (-0.085). V r,max = 0.147 (close to the disk) and V r,min = -0.168 (under the surface) are found in the liquid phase for ROSE while for Sunfluidh, V r,max = 0.152 is found in air-phase and V r,min = -0.168 on the surface. The markers V z,max and V r,min of ROSE are covered by those of Sunfluidh.

Figure 2 .

 2 Figure 2.15 -G = 0.1856, F r = 0 and Re = 10 000. Azimuthal velocity component V θ (r) (a) and axial velocity component V z (r) (b) at mid height z = G/2. The blue o symbols refer to measured steady values and measured rms values are also shown by blue '+' and '×' symbols.

  Figure 3.1 -(a) Large and (b) small experimental set-up.

Figure 3

 3 Figure 3.2 -Temporal line for a typical experiment.

Figure 3

 3 Figure 3.3 -LDV setup.

Figure 3

 3 Figure 3.4 -LDV installed in air.

Figure 3

 3 Figure 3.5 -(a) Measuring V z by LDV (side view ): blue lines show the propagation of laser in the vertical plane; (b) Measuring V θ by LDV (top view): green lines show the propagation of laser in the horizontal plane.

Figure 3

 3 Figure 3.6 -LDV correction when used for measurements in water in the large set-up. Colors of curves are in accordance with figure 3.5: blue lines for V z measurements and green for V θ . (a) Actual focus position versus focus position in air. '×' and '+' markers correspond to two campaigns of in situ measurements of r f in water. (b) Coefficients to correct measured V z and V θ [see equations (3.2.3) and (3.2.4)].

Figure 3 Figure 3

 33 Figure 3.7 -Same representation as in figure 3.6 for LDV correction for a flow of G8W2 in the large set-up.

  Figure 3.9 -Measurement of surface deformation: (a) point measurement set-up by Keyence LK-152 sensor; (b) line measurement set-up by LMI3D Gocator 2350 laser profile sensor

Figure 3 .

 3 Figure 3.10 -Scheme of the FTP set-up. The liquid surface is deformed when the whole cylindrical tank or only its bottom is driven into rotation. The projected fringes by the video projector are perpendicular to the plane of the paper in present configuration.

  Figure 3.11 -Demonstration of FTP principle. Signal 1 and signal 2 are respectively the reference signal I 0 (y) and the deformed signal I(y) along a given x position.

Figure 3

 3 Figure 3.12 -Right triangular prism placed in two orthogonal directions with the crest (a) parallel and (b) perpendicular to the fringes. (a1, b1) reference surface, (a2, b2) corresponding deformed fringes patterns.

Figure 3

 3 Figure 3.13 -Theoretical reconstructed surface for parallel and perpendicular positions as shown in figure 3.12.

Figure 3

 3 Figure 3.14 -Signal transformation in the process of FTP for the parallel (a) and perpendicular (b) cases. Only signals along a column at a fixed x are shown for the reference and deformed fringe images.

Figure 3 Figure 3

 33 Figure 3.15 -Influence of fringe pattern wavelength p 0 on FTP reconstruction.

Figure 3

 3 Figure3.17 -Influence of pixel size p x on FTP reconstruction.

Figure 3

 3 Figure 3.18 -Influence of filter width σ c on reconstructed profile (a) and relative error (b).

Figure 3 Figure 3 Figure 3 Figure 3

 3333 Figure 3.20 -Several profiles measured for the triangular prism placed in parallel (a1) and perpendicular (b1) positions. Red curves indicate the prism profile measured by the laser with a resolution of 0.15 mm. Zoom over a local region to show measurement precision: (a2) for parallel position and (b2) for vertical position.

Figure 3

 3 Figure 3.24 -Scheme of the Newton's bucket set-up.

  Figure 3.25 -(a) shows a snapshot of the deformed fringe pattern on the liquid surface for the steady regime at Ω = 14.56 rad/s and 34 mm. The graphs (a, b, c) show the treatment process to get a maximum rectangular area of the liquid surface in order to facilitate the reconstruction by the FTP. The graph (c) is used for present calculations.

Figure 3

 3 Figure 3.26 -(a): Temporal evolution of the surface height at r/R = 0, ±0.25, ±0.5 ± 0.75 (from bottom to top) for rotating speed Ω = 14.56 rad/s.(b): Steady surface profiles across the rotation axis along the x and y axes. (c) 21 equispaced contours between -0.026 m to 0.014 m of surface height at steady state measured by FTP and the Gocator laser profiler. The measurements of laser were 2D and thus spread along the azimuthal direction to 3D for comparison convenience. The blue circle is the border of the useful region in the images shown in figure 3.25. The red curves in the three graphs are measurements by laser.

Figure 3

 3 Figure 3.27 -Calibration on liquid flat surface at different heights. The heights at the image center of surfaces reconstructed by FTP are plotted as a function of the actual surface heights.

Figure 3

 3 Figure 3.28 -Mean height profiles of the triangular prism measured by FTP, shown in figure 3.21, plotted as a function of the laser-measured results for three values of L: (a) for the parallel case and (b) for the perpendicular case.

Figure 3 .

 3 Figure 3.29 -Height profiles h(x, 0) and h(0, y) of the straight Gaussian shape measured by FTP, shown in figure 3.22(b), plotted as a function of the profile measured by the laser.

  h r (x, y = 0)(m) h r (x = 0, y)(m) ∆L(cm) L=0.957 m h r = 0.95033 h laser + 1.43 • 10 -4 h r = 0.95762 h laser + 9.2 • 10 -5 4.62 p=3 mm -→ ∆L l = 5 cm -→ ∆L c = 4.24 cm L=1.159 m h r = 0.94461 h laser + 3.2 • 10 -4

  Figure3.30h r (r, t) at r/R = 0, ±0.25 and ±0.5 measured by FTP plotted as a function of h laser (r, t) for the Newton' bucket flow at 14.56 rad/s.

  Figure 4.1 -(a) Flow configuration. (b) Experimental set-up.

Figure 4 . 2 -

 42 Figure 4.2 -Effect of the Froude number on the steady axisymmetric flow structure for G = 0.1856, Re = 10 000, W e = 1263.6 and, from top to bottom, F r = 0.01, 0.2464, 1. (a) Streamlines of the meridional circulation consisting of 21 equispaced solid contours of ψ between 0 and ψ max = (0.0034, 0.0035, 0.0035) and one dashed contour at value (1/2)ψ min with ψ min = (-1.993, -1.975, -1.825)×10 -4 . (b) Isocontours of the azimuthal velocity V θ , with 21 equispaced contours between 0 and 1. The thick black line figures the free surface.

Figure 4

 4 Figure 4.3 -Deviation (h(r) -G)/F r of the free surface rescaled by F r. Same parameters as figure 4.2.

Figure 4 .Figure 4 . 4 -

 444 Figure4.2 shows the flow structure obtained at large Reynolds number (Re = 10 000), consisting in a central region with pure solid body rotation at the exact velocity of the disk (see right column) and an outer region with a meridional circulation (see left column). This circulation is relatively weak, as indicated by the low levels of ψ ≤ 0.0035. A striking feature is that the overall flow arrangement is hardly modified when the Froude number

Figure 4 Figure 4 . 6 -

 446 Figure 4.5 -Effect of the Reynolds number on the steady axisymmetric flow structure for G = 0.1856, F r = 0 and, from top to bottom, Re = 300, 1 200, 10 000, 19 500. Same representation as in figure4.2 with here ψ max = (0.0035, 0.0074, 0.0034, 0.0024) and ψ min = (0, 0, -1.993, -1.447) × 10 -4 .

Figure 4

 4 Figure 4.7 -Structure of the steady axisymmetric flow at Re = 10 000 for G = 0.5 (top) and G = 1 (bottom). (a) Streamlines of the meridional circulation with 25 equispaced contours between ψ min = (-1.8181, -2.87) × 10 -4 (dashed line for negative stream function) and ψ max = (0.0038, 0.0043), respectively for G = (0.5, 1.0). (b) Isocontours of V θ , with 25 equispaced contours between 0 and 1.

Figure 4

 4 Figure 4.8 -Location r s of the boundary between SBR and MR regions determined at z = G/4 as a function of Reynolds number Re for 4 different aspect ratios G, deduced from ROSE computations.

Figure 4

 4 Figure 4.9 -Axial velocity component V z (r s , G -z) at Re = 10 000 for G = 0.1856, 0.5 and 1.

Figure 4

 4 Figure 4.10 -Streamlines of the meridional circulation for G = 0.1856, Re = 30 000, computed from (a) V 2D [resp. (b) V 3D ]. The graph consists of 21 equispaced solid contours of ψ between 0 and ψ max = 2.1 • 10 -3 (resp. 1.98 • 10 -3) and one dashed contour at value 1 2 ψ min with ψ min = -1.12 • 10 -4 (resp. -0.924 • 10 -4 ). The locations of the points with maximum and minimum ψ are respectively indicated by the red and the blue cross.

Figure 4

 4 Figure 4.11 -Profiles of mean velocities V 2D (black dotted lines) and V 3D (red continuous lines) obtained by DNS for G = 0.1856, Re = 30 000. Velocity profiles as in figure 4.4.

Figure 4

 4 Figure 4.12 -Experimental and numerical mean and rms velocity distributions for the case G = 0.1856, Re = 30 000 and F r = 0.0335 (0 for numerical simulations). Mean values obtained by LDV (blue circles) and 3D DNS (black solid line) of (a) V θ at z/G = 0.96, (b) V θ at z/G = 0.5 and (c) V z at z/G = 0.5. In each case the rms amplitude is indicated by two curves above and under the mean value (thin dashed and dot dashed lines). Thick dashed line: disk azimuthal velocity. (d) Measurements of V z as function of r at different heights.

Figure 4 .

 4 Figure 4.13 -Experiments at G = 0.1856, Re = 81 400 and F r = 0.2464. Mean distribution of (a) V θ at different heights and (b) V z at z/G = 0.5 (blue circles) with indication of the rms level.

Figure 4 Figure 4

 44 Figure 4.14 -(a) Experiment at G = 0.1856, Re = 30 000 and F r = 0.0335: top view of the water layer seeded with Kalliroscope flakes. The red circle at r = 0.8 indicates the position of the line captured for the spatial-temporal diagram of figure 4.15a. (b) Numerical simulation at G = 0.1856, Re = 30 000 and F r = 0: snapshot of the axial velocity at z = 0.84G. The periodic numerical domain (r, θ) ∈ [0.4, 1] × [0, 2 3 π] has been replicated along the azimuth; for r < 0.4 (inside the black circle) a zero axial velocity has been plotted.

Figure 4

 4 Figure 4.17 -Comparison between azimuthal velocity profiles along the radius for aspect ratios (a) G = 0.1856, (b) G = 0.5 and (c) G = 1. Magenta solid: ROSE calculation at Re = 10 000 and z = G/4. Red dotted: distribution from [Iga, 2017]'s model (the thin black dashed line shows the position of r I s as determined in this model respectively at 0.57, 0.43 and 0.34). Black dashed: Ideal Rankine vortex model with core size r Is . Blue dash-dot: the TMBF model in[Mougel, 2014] with correction of the dissipation, leading to values r TMBF

Figure 5 Figure 5

 55 Figure 5.1 -(a): Volume-averaged kinetic energy evolution for Newton's bucket with symmetric flat surface for G = 0.25 and for various Reynolds numbers. (b): establishment time τ versus Re and (c) characteristic time τ s versus Re. Lower curve in black : G = 0.25. Upper curve in blue for G = 0.5. Symbols : values extracted from simulations. Dashed lines : fitting curves. Solid red lines with slopes of 1 and 0.5 are also plotted as reference.

Figure 5

 5 Figure 5.3 -Comparison of numerical simulation and measurements for the experiment SRS with H = 34.2 mm (G = 0.49), Ω = 7.84 rad/s. (a) Surface height evolutions at r = 0.25, 0.5, 0.75 and (b) azimuthal velocity V θ for two radial locations at z = 0.25. Symbols refer to measurements while solid lines stand for numerical results. Only 1 out of 25 points are represented for Gocator measurements in (a) while in (b) all data points of LDV are shown.

Figure 5

 5 Figure5.4 -Comparison of numerical simulation and experiment for the Newton's bucket case SRS with H = 34.2 mm (G = 0.49), Ω = 7.84 rad/s : (a) surface height profiles h(r, t) every 2 seconds from 0 s to 10 s as well as those in steady regime. Thin colored lines for simulation and symbols for experiments. The red symbols (resp. line) refer to the measured (resp. numerical) surface profile at steady state. (b) Azimuthal velocity V θ at steady state for the previously chosen probes. rms are indicated by + and × symbols. Note that the horizontal error bars have been drawn to show the potential errors of the measurement positions at worst cases due to an erroneous identification of the rotating axis position.

Figure 5 Figure 5 Figure 5 Figure 5

 5555 Figure 5.5 -Similar representation as in figure 5.3 used here for the case MRS.

Figure 5 .

 5 Figure 5.10(a) shows the surface height evolution at three radial positions, taken from the simulation under Sunfluidh, where the whole spin-up process is simulated including the acceleration of the disk.Figure 5.10(b) shows the numerical surface deformations at steady

Figure 5

 5 Figure 5.10(a) shows the surface height evolution at three radial positions, taken from the simulation under Sunfluidh, where the whole spin-up process is simulated including the acceleration of the disk.Figure 5.10(b) shows the numerical surface deformations at steady

Figure 5

 5 Figure 5.10 -Same case as in figure 5.9. (a) Surface height evolution at r = 0.25, 0.5 and 0.75 computed from Sunfluidh 2D unsteady simulation. (b) Comparison of the steady interfaces obtained using ROSE and Sunfluidh-the two curves are superimposed. Symbols are probe locations for LDV measurements.

Figure 5 Figure 5

 55 Figure 5.11 -Same case as in figure 5.9. Comparison between the numerical and experimental transient evolutions of V θ at the probe locations (a) r = 0.55 (b) r = 0.56 (c) r = 0.72 (d) r = 0.90, for z = 0.22. Black lines: simulation, blue lines: LDV.

Figure 5 Figure 5 Figure 5 Figure 5

 5555 Figure 5.15 -(a) Time evolution of h(r, t) at r = 0.25, 0.5, 0.75 for case 7 at Ω =19 rad/s (G = 0.25); (b) radial profiles of the fluid surface at different instants (red for the steady state). Circle points represent the experimental measurements by Gocator and solid lines come from numerical simulations.

Figure 5

 5 Figure 5.19 -Radial profile of V θ at the three measured heights for the same case as in figure 5.17. Circles indicate measurements and solid lines simulation results.

Figure 5 Figure 5

 55 Figure 5.20 -Comparison between numerics and experiments (a-c) Surface height evolution during spin-up at r = 0.25, 0.5 and 0.75 (from bottom to top) and (d-f) steady surface elevation profiles for the bottom-driven flow, corresponding to experiments n • 1,2 and 3 of table 5.5

Figure 5

 5 Figure 5.22 -Time evolution of the surface elevation h(r, t) for all the angular speeds listed in table 5.3 (a) at r = 0.5 and (b) r = 0.75. Note that dimensional variables are used here. Ω increases from top to bottom in (a) and from bottom to top in (b). Color coherence is kept between (a) and (b).

Figure 5

 5 Figure 5.23 -Spectra of the surface height signals presented in figure 5.22 (a) at r = 0.5 and (b) at r = 0.75. The frequencies have been rescaled using the disk frequency f d and the magnitude by the cylinder inner radius R (70 mm). Colors are in coherence with figure 5.22.

  Figure 5.25 -Perturbation maps δh(r, θ) of surface deformation (a)-(e) for cases 1-5 and (f) for case 8 of table 5.3 (G = 0.25), determined from Gocator line measurements. Colors indicate the perturbation height, the scale is in millimeters.

Figure 5 Figure 5

 55 Figure 5.26 -Space-frequency diagram |H s (r, f )| for case 1 of table5.5 (G = 0.5, F r = 0.97), for which no deviation from axisymmetry is visually observed at the surface.

Figure 5

 5 Figure 5.28 -List of experimental cases investigated in the large experimental set-up initially filled with a layer H = 26 mm of water, so that G = 0.1856, at different angular speeds Ω of the bottom disk. Top view of the flow pattern (the liquid is seeded with Kalliroscope) and nondimensional numbers. Only cases (d) and (e) have a dewetted central region indicated by a red circle.

Figure 5

 5 Figure 5.30 -(a) Temporal evolution of V θ at location (r, z) = (0.8, 0.0928) for different Froude numbers, increasing from bottom to top (cases a-e). Note that dimensional variables are used here. (b) Non-dimensional time-averaged velocities (circles) in established regime at the same location, with rms indications.

Figure 5

 5 Figure 5.31 -(a) Spectrum magnitude and (b) power spectrum of the azimuthal velocity V θ at (r, z) = (0.8, 0.0928) in established regime for cases (a)-(e) of figure 5.28 (G = 0.1856).

Figure 5

 5 Figure 5.33 -(a): Time-averaged surface profiles in established regime for cases 4-6 of table 5.5 (G = 0.5) where rotating ellipses m = 2 are present. (b)-(d): Perturbation maps δh(r, θ) of surface deformation for cases 4-6. Colors represent the perturbation height values in millimeters.

Figure 5 .Figure 5

 55 Figure5.35 -Representation of experimental cases on G8W2 with G = 0.5 from table 5.5 into the diagram of (Re, F r). The filled circular symbols represent the first three axisymmetric cases. The filled diamonds denote the wetted rotating ellipses. No dewetting of the bottom disk takes place for all the tested rotating speeds.

Figure 6 . 1 -

 61 Figure6.1 -Experimental investigations of cases having G = 0.25 with two set-ups and two viscous fluids. The four parabolic lines are defined by the dimensions of our two experimental set-ups and the two chosen liquids: water and G8W2, see legends for detailed combinations. Symbols correspond to different states observed in experiments: filled ones for wet cases (except the monogon-1) and empty ones for dry cases.

  Figure A.1 -Diagram of the algorithm of the unsteady code developed based on ROSE for configurations with flat surface

  Figure A.2 -Diagram of the algorithm of the unsteady code developed based on ROSE for configurations with deformed surface

is the vector ∆x n+1 = x n+1 -x n .

	It is a sparse matrix if discretized
	by finite difference. The variable of equation (2.2.33) Repeat the resolution of equation (2.2.33) until a certain precision level is reached, i.e.,
	x n+1 -x n < . The Intel MKL PARDISO solver is used to solve the matrix equation
	related to the equations (2.2.10) while MKL LAPACK routine is used to compute the matrix
	equation related to the equation (2.2.27).

Table 3 .

 3 1 -Geometry and materials for the two experimental set-ups generating the bottomdriven rotating flows.et al., 2004] while the small one is newly designed. Either the large or the small set-up can be mounted on the experiment table so that their respective bottom disks are driven by a DC brushed motor with a tachometer closed-loop speed control. The disk rotation speed can reach 250 rpm in our experiments.

	Parameter	Symbol	large set-up	Small set-up
	Cylinder inner radius	R (mm)	140	70
	Cylinder wall thickness	e (mm)	6.8	5.6
	Rotating disk radius	R d (mm)	139.6	68
	Disk material		brass	brass
	Side wall height from disk	H w (mm)	71.5	114
	Cylinder wall material		PMMA	glass

Table 3

 3

	7 cm

.4 -Linear regression relations of the curves presented in figure

3

.29 for the Gaussian shape with positive elevation and resulting corrections ∆L for each L.

  Figure 5.14 -Time evolution (top) and steady radial profile (bottom) of the surface height in almost axisymmetric cases 1, 3 and 5 (from left to right) of table 5.3. Top: symbols are for measurements with an interval of 0.1 s and solid colored lines are extracted from simulations. Bottom: red curves are for measurements and black curves for numerical results.

				Ω = 7.85 rad/s					Ω = 13.11 rad/s				Ω = 14.57 rad/s
	0.28								0.3							0.3
	0.26								0.25							0.25
	0.24															0.2
									0.2
	0.22															0.15
	0.2								0.15							0.1
	0.18	0		20	40	60	80	100	0.1	0	20	40	60	80	100	0.05	0	20	40	60	80	100
	0.35							0.5							0.5
	0.3							0.4							0.4
									0.3							0.3
	0.25
									0.2							0.2
	0.2							0.1							0.1
	0.15	0	0.2	0.4	0.6	0.8	1	0	0	0.2	0.4	0.6	0.8	1	0	0	0.2	0.4	0.6	0.8	1

Table

	N • Ω (rad/s)	F r	Re	mode
	1	11.65	0.968 1486.2	axi
	2	13.11	1.2265 1701.4	axi
	3	13.84	1.3662 1836.5	axi
	4	14.57	1.5143 1955.2	2
	5	19.00	2.5770 2564.8	2
	6	22.21	3.5210 3031.6	2

  2.8)Combining the second and the last equation of system (A.2.8), we can obtain Consequently, the equation (2.2.26) will be modified in the unsteady case as the kinematic relation V z = h V r to deduce it in the steady case is not true in the unsteady case.

	∂P ∂ξ	=	∂V r ∂t	-h	∂V z ∂t	-	G h	V z -h V r		∂V r ∂η	+ h	∂V z ∂η	-V r	∂V r ∂ξ	-V r h	∂V z ∂ξ
			temporal					ḣ					advection
										+	V 2 θ ξ	+	1 Re	(-h (	∂ω ∂ξ	+	ω ξ	) +	G(1 + h 2 ) h	∂ω ∂η	)
									rotation
	h F r	= P + C -	2 Re	G h	∂V z ∂η	-	Gh h		∂V r ∂η	-	1 1 + h 2 h 2 ∂V r ∂ξ	-h	∂V z ∂ξ
															viscous term
																(A.2.10)
												+	1 W e	h r(1 + h 2 ) 1/2 +	h (1 + h 2 ) 3/2
																interfacial tension term
	We can simplify again the term -1 1+h 2 h 2 ∂Vr ∂ξ -h ∂Vz ∂ξ	by the relation (A.2.3) like:

viscous

(A.2.9)

Remerciements

already mentioned that the oscillation amplitude of this pattern is an order of magnitude higher than that of the almost axisymmetric cases.

Azimuthal velocities have also been measured for this new regime. The chosen probe locations are shown in figure 5.39(a), at two heights z = 0.125 and 0.24 (we recall that G = 0.25 for this case). The figure also displays the time-averaged surface profile as well as two most perturbed instantaneous profiles. The time evolution of the azimuthal velocity at these probe locations during the transient process from rest is shown in figure 5.37 for z = 0.125 and in figure 5.38 for z = 0.24. While the spin-up lasts less than 100 time units, the monogon fully establishes itself rather slowly, from t = 300 to t = 550.

The mean and rms velocities measured by the probes yield the velocity profiles at two different heights presented in figure 5.39(b). On the basis of the surface height spectra, the phase velocity of this pattern as been determined: 2πf p = 0.64Ω. This value is confirmed by the spectra of velocity signals in the established regime, which show peaks at f = (0.65 ± 0.01)f d .

A more detailed study would be needed to understand the nature of this instability in the future, characterize its dependence with respect to the flow parameters. Such a counterrotating monogon has been reported in [START_REF] Bach | From Newton's bucket to rotating polygons: experiments on surface instabilities in swirling flows[END_REF], but in a configuration where the side cylinder wall rotated in the direction opposite to that of the disk, while here this cylinder is at rest. The flow regime is also much different: mostly laminar here (the Reynolds number is below 3300), it was fully turbulent in [START_REF] Bach | From Newton's bucket to rotating polygons: experiments on surface instabilities in swirling flows[END_REF].

ROSE: from steady to unsteady

Efforts were made at the beginning of my Phd to extend the steady code ROSE to unsteady simulations by including the temporal terms in equations A.0.1. The unsteady Navier-Stokes equations in axisymmetric configurations and in the curvilinear system are

A.2.1 Boundary conditions at free surface

The boundary conditions are changed compared to the steady case because the interface height evolves temporarily. We define a function F (r, z, t) = z -h(r, t) to describe the interface when F (r, z, t) = 0.

Taking the definition of F yields

Noting ḣ = ∂h ∂t and h (r, t) = ∂h ∂r , the above equation can be written as The kinematic condition can thus be deduced [START_REF] Kang | Numerical solution of axisymmetric, unsteady free-boundary problems at finite reynolds number. i. finite-difference scheme and its application to the deformation of a bubble in a uniaxial straining flow[END_REF].

The condition of zero tangential stress in the direction t can be developed in the unsteady case:

in unsteady case (A.2.9) from that of the steady case (2.2.31).

A.2.2 Algorithm of the unsteady code

The algorithm diagram for the unsteady ROSE with deformed surface can be described as in the figure A.2.