Experiments and Simulations of Free-Surface Flow behind a Finite Height Rigid Vertical Cylinder

: We present the results of a combined experimental and numerical study of the free-surface ﬂow behind a ﬁnite height rigid vertical cylinder. The experiments measure the drag and the wake angle on cylinders of different diameters for a range of velocities corresponding to 30,000 < Re < 200,000 and 0.2 < Fr < 2 where the Reynolds and Froude numbers are based on the diameter. The three-dimensional large eddy simulations use a conservative level-set method for the air-water interface, thus predicting the pressure, the vorticity, the free-surface elevation and the onset of air entrainment. The deep ﬂow looks like single phase turbulent ﬂow past a cylinder, but close to the free-surface, the interaction between the wall, the free-surface and the ﬂow is taking place, leading to a reduced cylinder drag and the appearance of V-shaped surface wave patterns. For large velocities, vortex shedding is suppressed in a layer region behind the cylinder below the free surface. The wave patterns mostly follow the capillary-gravity theory, which predicts the crest lines cusps. Interestingly, it also indicates the regions of strong elevation ﬂuctuations and the location of air entrainment observed in the experiments. Overall, these new simulation results, drag, wake angle and onset of air entrainment, compare quantitatively with experiments.


Introduction
The turbulent flow around a fully immersed cylinder has been the object of numerous experimental [1] and numerical studies [2,3]. When a free surface is present, the interaction between the cylinder wall and the flow structure leads to surface Kelvin waves [4] and eventually air entrainment. Turbulent free-surface flows are well documented for vertical flat plates [5,6], hull shapes [7] and hydrofoils [8]. However, the simple configuration of a finite height vertical rigid cylinder partially immersed in water flow has been less studied in the recent years. In 1947, Hay [9] reported an extensive investigation on the flow around a semi-submerged cylinder and provided many photographs from above and below the free surface for a wide range of velocities and diameters. The data uncovered the scaling relationship for the dimension of the various waves ahead and downstream of the cylinder; the wave dimensions seemed proportional to the square of the translating velocity [10,11]. Recent experiments [12] have provided the underneath velocity fields showing the presence of two recirculation regions: the first one below the cylinder's free-end and the second one behind the cylinder. These correspond to the two modes of vortex shedding attached to the cylinder. This complex flow structure affects the cylinder drag so that it is lower than the drag for fully immersed flow [9,10,13].

Experimental Setup
The experiments were performed for cylinders with D = 5 and 11 cm within a flume of 34 m in length, W = 90 cm in width and 120 cm in height filled with water. These are hollow cylinders made of polymethyl methacrylate (PMMA), which has a Young modulus around 3 GPa. Each cylinder was clogged at its free end; the total height was 65 cm and the wall thickness was 5 mm. They were mounted vertically and perpendicular to the water free surface on a carriage. Then, the cylinder immersed height, h, equaled 23 cm and moved along the flume with a velocity U. According to Griffith et al. [30], the critical horizontal blockage width, D/W, for the appearance of confinement effects on three-dimensional flow is 0.2; here D/W = 0.122. The top of the flume was open to the air and a camera was placed above the free surface in order to obtain the wake angles. The determination of the angles and the critical velocity for air entrainment around the cylinder was measured on images with a 2000 × 2000 pixel charged-coupled device camera (Basler).
The drag on the cylinder was calculated from the axial force, F x , acting on the whole cylinder. It was measured using a piezoelectric sensor (Kistler) with a sensitivity of −7.8 pC/N, placed at the top of the cylinder. It was assumed the contribution of air drag can be neglected, as the density of water is larger. The error on C D was estimated around 0.8%. F x was measured over time, and the time averageF x was done over 5 s. C D was calculated using: C D = 2F x /ρhDU 2 . In the range of Re studied, vortex shedding behind the cylinder was evidenced through the Strouhal number, St = f D/U = 0.2 ± 0.025, from f , the peak of the spectrum of the lift-force, F y , signal [31]. Moreover, the natural frequencies of the cylinder are unrelated to the frequency associated with vortex-induced vibration due to the elastic deformation that acts as an oscillator and triggers special modes of vortex shedding [32][33][34]. Further details of the experiment can be found elsewhere [13].

Numerical Method
The results were achieved with the YALES2 software [35], which is briefly presented below. Then a focus was given to the level-set method, which was used to treat the interface. Finally, the numerical configuration to reproduce experiments is presented.

Description of the YALES2 Library
The objective of the YALES2 code is to solve unsteady Navier-Stokes equations with massively parallel computers. It allows the use of unstructured meshes with several billions of elements [36,37], thus permitting the direct numerical simulation and LES of laboratory and semi-industrial configurations. The low-Mach number Navier-Stokes equations are solved with a projection method for constant and variable density flows [38,39]. These equations are discretised with a finite-volume 4th order central scheme in space and a 4th order Runge-Kutta-like scheme in time. The efficiency of projection approaches for low-Mach number flows is usually driven by the performances of the Poisson solver. In YALES2, the linear solver is a highly efficient deflated preconditioned conjugate gradient (DPCG) with a coarse-grid preconditioning [40].

Treatement of Free-Surfaces in YALES2
Free-surface tracking was ensured by the accurate conservative level-set (ACLS) method [41] and an improved re-initialization algorithm [42]. It predicts precisely the interface dynamics while maintaining the liquid mass. A scalar ψ is transported while the interface Γ is located at ψ = 0.5, i.e., Γ = {x ∈ R 3 |ψ(x, t) = 0.5}. In such a case, the ψ scalar has a hyperbolic tangent profile in the direction normal to the interface: where 2ε is the thickness of the profile and φ = ±|x − x Γ | is the signed-distance function. Considering the flow velocity u is divergence-free, the scalar ψ is first advected by the fluid according to: then reshaped using the re-initialisation equation from Chiodi and Desjardins [42]: where φ map = ln[ψ/(1 − ψ)], τ is a pseudo-time, and n = ∇φ is the interface normal. This step ensures that the ψ profile persists as a hyperbolic tangent while convergence criteria is achieved, i.e., ∂ψ/∂τ = 0. The two-phase coupling was provided by the ghost fluid method (GFM) [43]. More details on the principle of the GFM method are described elsewhere [44]. This method has been extended to three-dimensional unstructured meshes [45] and implemented in the YALES2 solver.

Dimensionless Coefficient and Setup for the Simulation
The simulations were performed in the frame of reference attached to the cylinder. The length of the computational domain was only 1 m long; the width was W and the height was 12D around the cylinder. Thus, the inlet velocity was uniform and the slipping walls were moving opposite the reference frame velocity, as depicted in Figure 1. Clearly, the boundary conditions at the free surface, the contact angle on the cylinder and the cylinder wall conditions are critical for the onset of capillary-gravity flows and are known to lead to finite-time dewetting singularities [46]. The cylinder wall was considered as a no-slip boundary condition and a 90 • contact angle. Some simulations (not shown) have been carried out with slip boundary conditions and lead to little changes to the free-surface patterns but large changes in the pressure and velocity fields. A mesh convergence study was performed with the drag coefficient, C D , defined previously. The element's size was changed from 1.5 to 5 mm on the cylinder wall and from 5 to 20 mm far away in the numerical domain. The final mesh was composed of 47.9 million cells. The edges of the elements were fixed to 2.5 mm on the cylinder wall and to 5 mm in the domain with a growth rate of 1.1. The numerical C D values are compared to the experiments [13].
The drag coefficient, C D , was determined by using the axial force F x on the surface of the cylinder, resulting in pressure expression extracted from Equation (3) and T x determined through the viscous term in Equation (4): where τ is the strain tensor.
/F re f , with F re f = ρhDU 2 /2. Moreover, the viscous coefficient, C x , can thus be defined using viscous force T x as: The wake of the cylinder was characterized by a fully turbulent flow. The smallest eddies could not be captured by the mesh size so the LES formalism accounted for the sub-grid scale impact on the largest scales. Indeed, the Kolmogorov scale, η DRe −3/4 , compared to the size of the elements on the cylinder wall, varies between 114 for U = 0.6 m/s and 215 for U = 1.4 m/s. The wall-adapting local eddy-viscosity (WALE) [47], which is based on the Smagorinsky model, was used in the present simulations.
The simulation was computed by setting the maximum Courant-Friedrichs-Lewy number (CFL) and a surface tension stability constraint. In the present simulation, CFL = 0.9, which corresponds to time steps of 1 ms for the simulation at the highest Re. The surface tension stability constraint was 0.5 and the corresponding surface tension time steps were about 0.2 ms (identical for all simulations). Regardless, the time step was always limited by the surface tension constraint. All simulations were performed on 840 Intel Broadwell cores. Each simulation was conducted for 60 Ut/D of physical time, which is achieved in approximately 15 h of wall-clock time.
The simulation results were consist over 6 runs at different velocities ranging from Re = 30,000 to 70,000 or Fr = 0.86 to 2. Specifically, 4 runs at Re = 30,000, 45,000, 55,000 and 70,000, corresponding to Fr = 0.86, 1,28, 1.57 and 2, respectively, are presented in the next section with a special focus on the vorticity and free-surface effects. All statistics subsequently presented were collected for 42 < Ut/D < 60, with a sampling rate of 100 Hz, in order to eliminate the transient initial state. For example, in our simulation at Re = 40,000, the mass lost for 42 < Ut/D < 60 is about 0.012%. Furthermore, the numerical scheme has been validated on test cases for the YALES2 code [48].

Results and Discussion
A first set of numerical results provides the instantaneous pressure and vorticity field around the cylinder for D = 5 cm. The second set of results concerns the wake angles, the free-surface height and its fluctuations in order to clarify the interaction between the wake, the free-surface height and the capillary-gravity wave theory. This theory predicts the position of crests that correspond to regions of free-surface height fluctuation and air entrainment. The last set of results compares experimental data of the drag and the wake angles on two cylinders as a function of the velocity with numerical results.

Pressure and Vorticity Fields
In Figure 2, the evolution of the dimensionless pressure field, P/ 1 2 ρU 2 , is presented for (a) Re = 30,000 or Fr = 0.86, (b) Re = 45,000 or Fr = 1.28, (c) Re = 55,000 or Fr = 1.57 and (d) Re = 70,000 or Fr = 2. The elevated height ahead of the cylinder at the bow wave [10] indicates pressure higher than the hydrostatic pressure. Indeed, when increasing Re and Fr, high pressure upstream of the cylinder remains roughly constant at P/ 1 2 ρU 2 ≈ 1 given by Bernoulli's equation for a stagnation point. However, for large Re and Fr, the pressure field differs from the classical vortex street flow. Strong free-surface deformations reduce the pressure drop downstream, as shown in Figure 2c,d. In addition to the strong pressure gradient upstream and downstream of the cylinder, a high-pressure region is located about 10D from the cylinder downstream of the cavity for Re = 55,000 or Fr = 1.57; see Figure 2c. It appears after the onset of air entrainment and corresponds to relatively strong free-surface deformation and the V-shaped pattern taking place mainly in the cavity downstream of the cylinder.
Furthermore, the flow dynamics are dominated by vortex shedding along the cylinder for small free-surface deformations visible on pressure (Figure 2a,b) and vorticity fields (Figure 3a,b). From Re = 50,000 or Fr = 1.43, stronger deformations of the free-surface occur and vortex shedding is inhibited in a layer region below the interface as shown in Figures 2c,d and 3c,d. As a result, two steady recirculation regions are formed behind the cylinder. Interestingly, a separated region is observed where a vertical vortex shedding is inhibited, in agreement with previous numerical works [16][17][18]. The small vortices formed deeper are transported toward the free surface further in the wake at a distance which increases with Re and Fr. The recirculation region below the cylinder is then pushed downward and deformed by the finite length of the cylinder. Below the free surface around the cylinder, the vorticity is at its maximum value, and it has been reported that secondary flows could be the origin of bubbly wakes [49]. The observed vorticity patterns in Figure 3 at the free end of the cylinder are similar to the ones observed in LES simulations of a finite height cylinder without free surface [50,51]. These detailed simulations shown the presence of: (i) a pair of vortices or tip vortices on the sides at the free-end of the cylinder and (ii) an arch vortex attached to the stagnation point at the bottom of the cylinder, which makes an angle with the horizontal that varies with Re. In fact, this flow structure is responsible for the concentration of vorticity in a layer region behind the free end of the cylinder, in contrast to the re-laminarized region below the free surface.

Deformation of the Free Surface
The flow behind the cylinder concentrates vorticity. In Figure 4, the instantaneous vertical vorticity component at Ut/D = 60 on the free surface is represented for Re = 30,000 and 60,000, corresponding to Fr = 0.86 and 1.57, respectively. The transition from sub-critical (Fr < 1) to super-critical (Fr > 1) flow clearly modifies the free-surface vorticity distribution. The sub-critical vorticity distribution (see Figure 4a) resembles the classical vortex street distribution of the flow past a cylinder, whereas the super-critical vorticity (see Figure 4b) concentrates vorticity in the neighbourhood of the cylinder. For Re = 60,000, Fr = 1.72, the vortices are concentrated in the cavity, as the flow far downstream the cylinder and below the free surface is almost laminar. Consequently, the shape of the wake changes from the Kelvin wake to a narrower V-shaped free-surface pattern [26,52].  ( (c c c c c) ) ) ) ) ( ( ( ( (d d d d d d) ) ) ) ) Figure 4c,d display experimental snapshots of the view from above the free surface at corresponding values of Re and Fr for comparison. The lighting is from above, so the white areas are surfaces directed normal to the horizontal (maximum or minimum elevation height areas). The analysis of the pictures from the top view allows the wake angles, α and β, to be estimated as a function of the velocity U. The wake half-angle α is defined as the radiated wake angle corresponding to the maximum wave amplitude [26]. The half-angle β is defined as the angle made by the cusp attached to the cylinder. An example for the determination of the half-angles α and β are sketched on the photographs in Figure 4c,d.
The free-surface elevations, presented in Figure 5, describe the bow wave, the cavity region downstream of the cylinder and the Kelvin waves [11,18,25]. The oblique views presented in the first column of Figure 5 show the emergence and decrease of the cusp wake half-angle β. The instantaneous free-surface and the time-averaged elevation height highlight the Kelvin wake's stationary patterns generated by the flow around the cylinder. For instance, at Re = 30,000 and Fr = 0.86, the radiation angle is above that which is predicted by Kelvin [4]. The time-averaged elevation heights (second column of Figure 5) were computed for 42 < Ut/D < 60. As Fr increases, the bow wave height increases, and the cavity behind the cylinder gets deeper. For the highest Fr (Figure 5b-d), the radiation angle seems to evolve as 1/Fr, which is consistent with the experiments of Moisy and Rabaud [26], although the present computation domain only contains several wavelengths.  The crest lines are defined from the capillary-gravity waves theory [24,26,53]. They are compared with the time-averaged free surfaces and the associated relative standard deviations (STD/D) obtained from the numerical simulations. Three crest lines, represented by lines of different colours on Figure 5, have been plotted following the equations of Binnie [53] and considering radiated wavelengths of order 2πD. This far-field theory provides the wake pattern as determined by the parameter U/c min , where c min = (4gσ/ρ) 1/4 = 23 cm/s is the minimum phase velocity of the capillary-gravity waves. For the finite ratio U/c min > 1.938, the evolution of the radiation angle with the wavenumber presents two local extrema, leading to the formation of two cusps in the crest lines where energy concentrates [26]. When increasing Re and Fr, the lower extrema of the cusp angle greatly decreases, leading to a narrower V-shaped pattern towards the cylinder wall (see Figure 5d).
The cusps are correlated with strong unsteady fluctuations, as seen in the last column of Figure Figure 5. Although those locations correspond to newly formed bubbles and advected bubbles which where previously formed, they coincide reasonably well with strong unsteady fluctuations in the wake (last column on Figure 5b), then in the cavity behind the cylinder as Fr increases (Figure 5c,d).
Indeed, air entrainment observed in the wake is due to wave breaking located downstream of the stationary V-shaped pattern from Re = 45,000 or Fr = 1.28, whereas air entrainment in the unsteady cavity occurs upstream of the V-shaped pattern close to the cylinder wall. There is a relation between the modification of the flow dynamics below the free surface and the characteristic V-shaped pattern observed for increasing Re and Fr. The V-shaped pattern coincides with vortex shedding inhibition observed on the free surface ( Figure 4) and in the sub-surface flow (plane z/D = −1 on Figures 2 and 3c,d).

Drag and Wake Angles
The vortex inhibition described above should be translated into global drag. Now, Figure 6a,b, presents the values of C D for D = 5 cm and 11 cm as a function of Fr and Re. C D is first found to slightly increase; then, it passes through a maximum corresponding roughly to the onset of air entrainment. The error bars in Figure 6a represent the standard deviation in the simulation run for 42 < Ut/D < 60. Note that the numerical simulations reproduce the experimental air entrainment found in experiments [13]. At the onset of air entrainment, represented by the transition from filled symbols to empty symbols, the appearance of small bubbles is insufficient to modify the pressure and the drag around the cylinder. Then, after a plateau, the evolution of C D seems to decrease with Re or Fr at the same rate, suggesting a linear dependence, which is in agreement with two-phase flow simulations [54]. The numerical results slightly underestimate C D , but reproduce quantitatively the onset of air entrainment and C D at high Fr or Re. A similar evolution of the drag coefficient is observed for D = 11 cm. However, the maximum of C D occurs at a much larger Reynolds number than the onset of air entrainment, suggesting that the entrapped air is not the only cause of C D 's decrease. As the Fr increases, the free-surface evolution modifies the pressure distribution around the cylinder and, thus, modifies the drag coefficient.
Concerning the angles' evolution, for U/c min > 1, where c min = (4gσ/ρ) 1/4 = 23 cm/s is the minimum phase velocity of the capillary-gravity waves, α remains constant and decreases for large U/c min for D = 5 cm (see Figure 6c. For D = 11 cm, the wake half-angle α remains constant at small U/c min ; then, it oscillates at higher values (see Figure 6d). Such high values may be explained by blockage effects, the multiple reflection from the walls and the limited visualization area away from the cylinder. Both simulations and experiments indicate that the wake half-angle α reduces with an increase in velocity. Note that the simulation results predict quantitatively the half-angles α and β for D = 5 cm. The role of free-surface deformation on global forces is highlighted by the appearance of the β angle at the maximum value of the drag coefficient for both diameters.

Conclusions
In this article, we presented experimental and numerical results of three-dimensional simulations on the drag, the wake angles, the pressure fields, the vorticity fields and the free-surface elevation of a vertical finite height rigid cylinder in turbulent free-surface flow. The experiments were conducted on cylinders of D = 5 and 11 cm. The simulations were performed for D = 5 cm in the range of parameters 30,000 < Re < 70,000 or 0.86 < Fr < 2, which coincides with the occurrence of the V-shaped pattern on the free surface and the onset of air entrainment. The numerical and experimental results are in good agreement. They quantify and illustrate the intricate relation between the viscous, inertial, gravity and capillary effects.
Concerning the drag, C D first increases, then decreases, with Re and Fr. Interestingly, this decline of C D coincides with the occurrence of cusps attached to the cylinder associated with vortex shedding inhibition in a layer region below the free surface. The simulations predict quantitatively the onset of air entrainment [13]. The free-surface elevation simulations are compared to the V-shaped patterns predicted by the capillary-gravity waves theory [26]. The crest lines from the theory indicate the presence of cusps, where energy is concentrated, which seems to correspond to unsteady free-surface elevation fluctuations and air entrainment location in the simulations. In the wake, air entrainment is due to wave breaking, while behind the cylinder, air entrainment is caused by the unsteady cavity. From Re = 50,000 or Fr = 1.43, the emergence of a stationary free-surface V-shaped pattern pushes unsteady wave-breaking areas far away from the cylinder. These results are the basis of more realistic situations taking into account wind and surface waves.

Funding:
The authors acknowledge the financial support of the Normandy Region and the European Union ERDF funding's through the NEPTUNE project. Our work has also benefited from the financial support of the Agence Nationale de la Recherche (ANR) through the programme 'Investissement d'Avenir' from the laboratoire d'excellence Energy Materials and Clean Combustion Center (LabEx EMC3).