Numerical study of instability mechanism in the air-core vortex formation process

Abstract

Air-core vortices are a ubiquitous phenomenon in the intakes of hydropower stations. Due to the transient and instability of two-phase vorticial flow, the prediction of air-core vortex formation is challenging, and understanding of the instability mechanism remains elusive. In this study, the large eddy simulation (LES) method and a coupled level-set and volume-of-fluid (CLSVOF) method are performed to study air-core vortex formation in a benchmark reservoir with a horizontal intake pipe. The process of air-core vortex formation can be classified into an inception stage, an instability stage, and a stability stage. In the instability stage, the surface vortex repeatedly goes through the process of inception, enhancement, attenuation, and extinction. The movement of the counterrotating secondary vortices and water surface level fluctuation plays a negative role in air-core vortex formation. The additional motion generated by the counterrotating pair drives the pair to the back wall. The main vortex undergoes attenuation due to the stretching/tilting effects induced by the secondary vortex. Water level fluctuations briefly increase the submergence depth, which in turn reduces the vertical velocity gradient and vertical vorticity, destabilizing the vortex. The perturbation of the air-core vortex by water level fluctuations is present only at the beginning of the instability stage.

KEYWORDS:

• Air-core vortex
• water level fluctuation
• wave
• counterrotating pair
• large eddy simulation

1. Introduction

Free surface vortices can usually be seen at the intakes of inlet structure in the hydropower stations and pumping stations. According to their shape and swirling intensity, Alden Research Laboratory (ARL) qualitatively divided free surface vortices into six stages, ranging from surface swirling to stable air-core vortex (Hecker, 2017). When a vortex is generated at an intake, it reduces the efficiency of hydraulic machinery and the generating capacity of the power station to different degrees (Knauss, 2017). Among the six stages, the air-core vortex, as the sixth vortex type, has the greatest strength and capacity to draw air and floating trash into the intake pipe, which can be most harmful to the operation of the hydropower station or pumping station. Nonuniform flow patterns and two-phase flow caused by an intake air-core vortex seriously reduce the discharge capacity of the inlet structure of the hydropower station, resulting in vibration and cavitation of hydraulic machinery and affecting the safety of pressure systems and surrounding buildings (Li et al., 2021; Möller et al., 2015). Therefore, it is critically important and meaningful to explore the mechanism the formation of intake air-core vortices, to reveal their evolution, and to seek effective measures to control their evolution and movement. This is a difficult and challenging topic in the field of fluid mechanics because of the complex vortex motion (Andersen et al., 2003; Saleem et al., 2020).

The formation mechanism and hydraulic characteristics of an air-core vortex at an intake have been extensively studied theoretically, experimentally, and numerically. The theoretical analysis mainly focues on the analytical solution of the velocity field of the vortex flow by solving the equation directly based on certain assumptions and different degrees of simplification of the Navier-Stocks (N-S) equation. Rankine (1858) proposed a two-dimensional planar model, assuming that the flow field is irrotational without vorticity outside the core radius, and rotational inside the core radius with constant vorticity. Owing to inconsistency in the flow characteristics between the inside and outside of the core radius, Rankine’s model shows a significant discontinuity in the distributions of tangential velocity and vorticity at the radius of the vortex core. Burgers (1948) and Odgaard (1986) assumed that the free surface vortex is constant and axisymmetric, and solved the N-S equation in cylindrical coordinates to derive the tangential velocity that is continuous and more consistent with the physical phenomenon. More studies using theoretical analysis are needed owing to the different assumptions made for various types of vortex models. In addition to the theoretical study, several experimental results contributed to the development of the theoretical model. Based on the experimental results, the tangential velocity and vorticity of the vortex theoretical models were modified. The water surface profile and the distribution law of the axial and radial velocities of the vortex were also analyzed. In addition, the reasonableness and influence of the axial and radial velocities assumed in the vortex theoretical models were evaluated and discussed by Hite and Mih (1994). In hydraulic engineering, establishing the relationship between the hydraulic characteristics of the air-core vortex and operating and boundary conditions has drawn much attention, where the operating and boundary conditions are the submergence depth (Guo et al., 2020), the inlet diameter and profile (Yang et al., 2014), the discharge flow rate (water, air) (Möller et al., 2015), the arrangement type of the inlet pipe (Huang et al., 2022), and the surface tension and viscosity (Suerich-Gulick et al., 2014).

Numerical simulations have also been applied to reproduce the air-core vortex phenomenon and explain its formation mechanism. The numerical simulations used to reproduce the air-core vortex are challenging because of two main difficulties: capturing the air–water interface and selecting a turbulence model (Huang et al., 2022). To overcome these difficulties, different numerical methods were proposed to simulate the air-core vortex at the intake and their accuracies were evaluated (Škerlavaj et al., 2014). The flow velocity around the free surface vortex and the air-core length of the vortex obtained from the single-phase flow simulations had larger errors compared to the results of two-phase flow models, such as the coupled level-set and volume-of-fluid (CLSVOF) method. Therefore, two-phase flow models are seeing increasing use for the simulation of the complex motion of air-core vortices (Guo et al., 2020; Zi et al., 2021). Various types of turbulence models, from simple to complex, have been applied to simulate air-core vortex phenomena. The Reynolds averaged N-S (RANS) has been widely applied to the study of air-core vortices due to its fast computational speed and relative small requirement for computational resources (Chang & Wei, 2022; Sun et al., 2020). Detached eddy simulation (DES) is a hybrid method of large eddy simulation (LES) and RANS, which combines the advantages of the high efficiency of RANS and the high accuracy of LES. Based on the application DES, the effects of decreasing submergence and increasing discharge on the flow field properties were discussed by Wu et al. (2019). Compared with DES, LES can resolve more details of turbulence and provide more comprehensive information about flow fields. Therefore, in this work, LES is applied to simulate the free surface vorticial flow field at an inlet and analyze the vorticity changes during the formation of an air-core vortex.

Many researchers have found that the free surface vortex undergoes a certain time instability behavior before forming a stable air-core vortex, and the vortex periodically disappears and regenerates (Leweke et al., 2016; Zi et al., 2020). Mizushima et al. (2014) studied the vortex instability phenomenon and swirl direction of bathtub flow by numerical linear analysis and confirmed that the vortex instability was induced when the flow rate exceeding a threshold. Turning to other fields of vortex dynamics, more studies on vortex instability can be found, such as helical vortices for wind turbines and helicopters (Brynjell-Rahkola & Henningson, 2020) and the wing-tip vortex of an airfoil (Del Pino et al., 2011; White & Saric, 2005). Although these studies provide useful references to analyze the phenomenon and mechanism of vortex instability during the formation of air-core vortices, helical and wing-tip vortices are different from free surface vortices involving two-phase flows. In previous studies, the instability phenomena of periodic generation and vanishing of the free surface vortices at the intake were mentioned but not sufficiently described and investigated. The fundamental mechanism of the vortex instability phenomenon remains elusive.

Therefore, in this paper, numerical simulations using a combination of LES and CLSVOF methods are performed to reproduce the vortex instability phenomenon occurring during the formation of an air-core vortex, of which the result is verified by experimental data. The vortex instability phenomenon at the intake of a hydropower station is investigated in detail, which is rarely reported before. Based on the instability characteristics of the vortex, the formation process was divided into different stages. Furthermore, the mechanism of vortex instability induced by vortex pair and water surface fluctuation is revealed. The remainder of this paper is organized as follows: in Section 2, the numerical method, simulation setup and validation are introduced. In Section 3, the phenomenon of vortex instability is described, and the mechanism of vortex instability is revealed. Finally, the conclusions are summarized in Section 4.

2. Numerical model

2.1. Governing equation

The basic idea of LES is that large-scale turbulent pulsation is obtained by numerically solving the filtered N-S equations, and using a subgrid-scale model to represent the effects of small-scale turbulent pulsation on large-scale motion (Tyliszczak et al., 2016). In this paper, the calculation addresses an intake air-core vortex that occurs far from the surrounding sump walls and is not dominanted by turbulent boundary layer flow near the walls. Therefore, the refinement of the mesh near the wall is not included, and the wall-modeled LES (WMLES) approach of Celik et al. (2005) is adopted to reduce the number of grids and optimize the allocation of computational resources. This method overcomes the problem that the wall-resolved LES method requirs a very fine mesh scale to simulate near-wall turbulence at high Reynolds numbers. At present, this method has been widely used in turbulent flow studies (Kawai & Larsson, 2012; Mockett et al., 2012; Shur et al., 2008). The momentum equations and continuity equation are written as follows (Espinoza-Jara et al., 2022): (1) $\begin{array}{rl}& \frac{\mathrm{\partial }{\tilde{u}}_i}{\mathrm{\partial }{x}_i}=0\end{array}$(1) (2) $\begin{array}{rl}& \frac{\mathrm{\partial }{\tilde{u}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\tilde{u}}_i{\tilde{u}}_j}{\mathrm{\partial }{x}_j}=-\frac{1}{\rho }\frac{\mathrm{\partial }\tilde{p}}{\mathrm{\partial }{x}_i}+\nu \frac{{\mathrm{\partial }}^2{\tilde{u}}_i}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_j}-\frac{\mathrm{\partial }{\tau }_{ij}}{\mathrm{\partial }{x}_j}+g_i,\end{array}$(2) where $\left(\stackrel{\sim }{\cdot }\right)$ denotes the filtered quantities in LES; u denotes the fluid velocity; ρ and v are the density and kinematic viscosity of the fluids, respectively; p is the pressure; g is the gravitational acceleration; and τ_ij is the subgrid-scale stress tensor, which is expressed as: (3) ${\tau }_{ij}=2{\nu }_{SGS}{\tilde{S}}_{ij}-\frac{1}{3}{\tau }_{kk}{\delta }_{ij},$(3) where τ_kk is the isotropic part of SGS stress; δ_ij is the Kronecker symbol; ${\tilde{S}}_{ij}$ is the large-scale strain rate tensor and ν_SGS is the subgrid eddy viscosity that must be modeled. They are calculated as: (4) ${\tilde{S}}_{ij}=\frac{1}{2}\left(\frac{\mathrm{\partial }{\tilde{u}}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\tilde{u}}_j}{\mathrm{\partial }{x}_i}\right)$(4) and (5) $\begin{array}{rl}{\nu }_{SGS}& =min[{(\kappa d_w)}^2,{(C_{Smag}\mathrm{\Delta })}^2]\cdot \sqrt{2{\tilde{S}}_{ij}{\tilde{S}}_{ij}}\\ & \cdot \{1-\exp [-(y^{+}/25{)}^3]\},\end{array}$(5) where d_w is the wall distance; κ = 0.41 and C_Smg = 0.2 are constants, and y + is normal to the wall inner scaling. The LES model is based on a modified grid scale to account for the grid anisotropies in wall-modeled flows: (6) $\mathrm{\Delta }=min(max(C_w\cdot d_w,C_w\cdot h_{max},h_{wn}),h_{max}).$(6) Here, h_max is the maximum edge length for a rectilinear hexahedral cell (for other cell types or conditions, an extension of this concept is used). h_wn is the wall-normal grid spacing, and C_w= 0.15 is a constant.

To directly resolve the free surface vortex at the intake, the CLSVOF method was applied to track the air–water interface. The CLSVOF method has the advantages of both the level-set method and the volume-of-fluid (VOF) method. The former can accurately discribe the water–air interface, and the latter can maintain mass conservation (Sussman et al., 1994). The level-set method uses a signed distance function ϕ to describe the free surface (Equation 7), where the value of ϕ is zero at the free surface and positive and negative in water and air, respectively. The VOF method adopts the air–water volume distribution function F (Equation 8), which is defined as the volume fraction of water in the simulation grid cell. F = 0 means that the grid cell is full of air, and F = 1 means that there is no air and only water in the grid cell. When F is between 0 and 1, there is an air–water interface within the grid cell. The CLSVOF tracing the interface of the two phases can be obtained by solving equations for ϕ and F: (7) $\begin{array}{rl}& \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j\varphi )}{\mathrm{\partial }{x}_j}=0\end{array}$(7) (8) $\begin{array}{rl}& \frac{\mathrm{\partial }F}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_{j}F)}{\mathrm{\partial }{x}_j}=0.\end{array}$(8)

2.2. Simulation setup

The calculated geometric model is the same as Möller’s test (2015), which is a hydropower station intake with a horizontal inlet pipe (intake pipe direction, φ = 0°). The structure is sketched in Figure 1. The origin of the coordinates is located at the bottom center of the back wall of the reservoir as shown in Figure 1. The inlet pipe diameter D is 0.4 m, which is selected as the characteristic length. The overall size of the calculation domain is 17D × 4D × 3D (length × width × height), the water depth of the reservoir is 2.5D, and the submerged depth of the intake pipe is 1.5D. The length of the inlet pipe in the reservoir is 2D, and the distance between the inlet of the calculation domain and the intake is 10D. The whole computational domain is plane symmetry with respect to the xz plane across the central axis of the intake pipe.

Figure 1. Sketch of the computational domain (D denotes the diameter of the intake pipe).

The upper part of the reservoir, an open rectangular tank, is set as the atmospheric pressure inlet boundary. The boundaries on both sides, the bottom, and the back wall of the reservoir, and the intake pipe wall are all set as solid walls with the no-slip boundary condition. The area-averaged pressure distribution on the outlet of the pipe, corresponding to a certain flow rate, is assumed to be −25,500 Pa (Ahn et al., 2017). The pressure distribution corresponds to the distribution law of hydrostatic pressure on the outlet. Under this pressure condition, although the outlet water flow rate changed constantly during the whole calculation process, it always fluctuated around the average value of 0.65 m³/s in general. The maximum error between the real-time flow rate and the average flow rate was 7.1 × 10⁻⁴, and the standard deviation was 2.62 × 10⁻⁴. The time-averaged velocity v_D = 5.13 m/s of the outlet was selected as the characteristic velocity. The inlet of the reservoir adopted the pressure inlet boundary of complex air–water two-phase flow defined by the user defined function (UDF) (Kan et al., 2021). Below z = 2.5D, the pressure inlet boundary allows only pure water to pass, where the pressure distribution corresponds to the distribution law of hydrostatic pressure. Above z = 2.5D, only air can pass through the inlet boundary at atmospheric pressure. Before starting the LES for the free surface vortex, a simple initialization calculation is carried out for the flow field, and its solution is used as the initial field of the LES, which can greatly improve the simulation efficiency and reduce the calculation time of the LES. The RNG k-ε turbulence model can be used to handle flows with swirls flow, high strain rates, and large curvatures (Li et al., 2022). The initial field of the LES calculation was obtained by RNG k-ε turbulence model after approximately 8000 iterations to reach a statistically steady state. After the corresponding RSM was reduced to 10⁻⁴ and the monitored outlet flow rate was stable, the LES module was activated. The time step of the LES calculation was 6.41 × 10⁻³ (D/v_D), and the maximum number of iterations at each time step was 20. Similarly, RSM reached convergence when it dropped to 10⁻⁴. The governing equation was discretized based on the finite volume method. In the solution process, the PISO algorithm was used for discretizing the convection field equation. The least squares cell-based scheme was employed for gradient discretization. The transient temporal discretization was realized by a second-order implicit scheme. The bounded central differencing scheme was adopted for the momentum equation. The bounded second-order implicit scheme was conducted on the temporal discretization. The air-core vortex simulation was completed based on FLUENT 18.0 and run on a 52-thread dual-core Intel Xeon Platinum 8171M workstation.

ICEM was used to construct a structured grid for the whole computing domain with the meshes near the solid walls refined. The free surface vortex wonders during the formation process of the air-core vortex. Therefore, it is necessary to encrypt the areas where free surface vortices might occur. Referring to the results of Möller’s paper (2015) on vortex wondering, we refine the plane area with 3D × 3D centered on the inlet pipe. Two sets of grids with different grid sizes were generated. The total numbers of coarse and coarse grids were 4,132,976, and 8,516,679, respectively. The result of fine mesh encryption is shown in Figure 2, grid cell size Δx = Δy = Δz = 0.015D in the vortex area. More grid information can be found in Table 1. The shape of the free surface deformed greatly, where the upper part had a sunken funnel shape and the tip had a small hook-like structure. To ensure that the grid resolution was sufficient for vortex formation, a grid sensitivity study was performed for these two sets of grids.

Figure 2. Visualization of the mesh (fine grid with 8.5 million hexahedron cells).

Table 1. Mesh information for the air-core vortex simulations.

Mesh	Total node number	Vortex area node number	Δx,y,z /D in the vortex area	Average Δx,y,z /D
Coarse	4,132,976	1,728,000	0.020	0.033
Fine	8,516,679	3,188,200	0.015	0.026

2.3. Validation of the numerical model

The development of the air-core vortex is a typical air entrainment phenomenon with a strong correlation with instantaneous air–water motion. To obtain more detailed information on vortex development, it is necessary to evaluate the mesh resolution quality of LES. Celik et al. (2005) proposed the LES index of quality (LES_IQ) to evaluate the analytical accuracy of large eddy simulations of spatial-scale turbulent pulsations, and suggested that LES_IQ of the LES calculation with a high Reynolds number should satisfy 75%∼85%. LES_IQ is usually expressed as the ratio of the analytic part of the turbulent kinetic energy (k_res) to the total turbulent kinetic energy (k_tot), as shown in Equation (9). When the LES_IQ is close to 1, this indicates that LES has solved more details of the vortex flow field under the grid. In this paper, we used LES_IQ to estimate the mesh resolution quality of the calculated results on coarse and fine grids. (9) $LES\mathrm{\_}IQ=\frac{k^{res}}{k^{tot}}$(9) Figure 3 shows the distribution of turbulent kinetic energy k_res at the tangent plane z = 2.38D when the stable air-core is formed (tv_D/D = 1282.5). The position of maximum k_res was selected to study the analytical scale variation in turbulent kinetic energy in two directions, streamwise R_x and spanwise R_y, as shown in Figure 4. In this study, the LSE_IQs in both the fine and coarse grid cases significantly decreased near the maximum value of k_res (R_x/D = R_y/D = 0). However, the fine grids could still resolve more than 90% of the turbulent kinetic energy in the core region where vortices occur and the average LES_IQ in the whole computational domain can reach 97.8%. This indicated that the grid resolution for investigating the vortex and turbulence characteristics of the flow field was sufficient, and it was applied to the subsequent statistical analysis of pulsation and vorticity transport research.

Figure 3. Distribution of k_res on the plane z = 2.38D (tv_D/D = 1282.5).

Figure 4. LES_IQ in the directions of streamwise R_x and spanwise R_y.

Figure 5 shows the calculation results for the air-core vortex based on FLUENT. For the air-core vortex experiment results of Möller et al. (2015) in Figure 5a, their experimental result and our numerical result show reasonable agreement with a funnel-shaped air-core vortex generated at the top of the intake, whose structure included an upper funnel section, a middle straight section, and a small hook-like section. Because the flow rate at the outlet was set to 3.5 times that of the experiment, the diameter of the funnel in the numerical simulation is larger than that in the experiment. The air–water interface in Figure 5b is illustrated by the isosurface of the volume fraction function with F = 0.5. In the following analysis, if not state specifically, F = 0.5 was adopted.

Figure 5. Shape comparison of the air-core vortex: (a) experimental result; (b) numerical calculation (F = 0.5, tv_D/D = 1282.5).

Figure 6 was the marginal histogram of the vortex wandering when the air-core vortex was formed. The sampling frequency of the coordinates of the wandering air-core vortex was 5 Hz. The coordinate position of the stable air-core vortex was collected, whose sampling time was 40 s. The wandering of the air-core vortex was roughly centered at (2.0, −0.1). Compared with the distribution of vortex wandering in the x and y directions, the vortex locations were more concentrated in the y direction. The blue line was the experimental result from Möller (2013). Comparing the vortex motion histogram predicted in this paper, it shows that the vortex center distribution in the y direction was highly consistent with the experimental result in a statistical view. Although there was a certain difference in the distribution of the vortex locations in the x direction near x/D = 2, the overall distributions of the two were still similar.

Figure 6. Marginal histogram of the vortex wandering.

Furthermore, the tangential velocity distribution of the air-core vortex on the plane Z = 2.38D near the water surface was selected for quantitative comparison with the results yielded by previous models, as shown in Figure 7. The tangential velocities of the simulation results, Burgers model, and Hite and Mih model all increase first and then decrease along the radial direction, and the calculation results were more similar to that of the Burgers modelthan the Hite and Mih model. When r/D > 0.2, the tangential velocities from the three approaches were similar, and the error was within 0.55%. The difference was in the vortex core radius region. The calculated vortex core radius and maximum tangential velocity were close to the Burgers model and deviate noticeably from Hite and Mih model, and the maximum tangential velocity obtained by numerical calculation was less than that of the two theoretical models (Burgers model and Hite and Mih model). This may have been due to the theoretical models, which omitted higher-order terms and adopted several assumptions (vorticity viscosity coefficient, axial velocity, etc.), resulting in certain differences between the theoretical model and the real field (Hite & Mih, 1994; Odgaard, 1986). In addition, there were some numerical errors in the LES numerical calculation results, so it was expected that the three results for tangential velocities would not be completely consistent. In general, the numerical simulation results adopted in this paper were in good agreement with the theoretical models, which provided some confidence in the following in-depth analysis.

Figure 7. Comparison of LES and theoretical models for tangential velocity distribution (tv_D/D = 1092.7).

3. Results and discussion

3.1. Evolution of the free surface vortex

The vorticity is generated above the intake of the hydropower station. Over time, the surface dimple caused by the vorticity gradually becomes a fully developed air-core vortex. Finally, air is continuously sucked into the inlet pipe through the air-core vortex. Three stages are identified according to the shape of the free surface and the associated vortical structures: surface swirl, surface dimple, and fully-developed air-core (Zi et al., 2020). Four development stages of the vortex have been identified: surface swirl with very small dimple, well-developed dimple, bubble entraining core, and full air core (Cristofano et al., 2014). In both cases, the stages are classified in an aspect of the flow pattern of the free surface vortex. In this study, we focused on vortex instability evolution during the air-core vortex formation process. Thus, based on the visual air–water interface, the evolution was classified into three stages: the vortex inception stage, instability stage, and stability stage. Figure 8 shows the free surface at the three stages during the formation process of the air-core vortex. Figure 9 shows streamlines on the plane z = 2.38D corresponding to the results in Figure 8. The characteristics of each stage are described below:

1. Vortex inception stage. At this stage, a small swirl appeared on the free surface without significant surface depressions (Figure 8a). The symmetry reservoir boundary and pipe intake contributed to the symmetric flow. As shown by the streamlines on the plane z = 2.38D (Figure 9a), although the flow field developed for 134.7 D/v_D, under the influence of axisymmetric inlet boundary conditions, upstream flows (x > 3D) were still uniform inlet flows. A swirl began to appear above the pipe intake (x < 3D), which resulted from the excitation of random pulsation of the high Re number flow (Re = 2.04 × 10⁶) and gradually converges toward the intake with the movement of water flow.

2. Vortex instability stage. When tv_D/D = 307.8, a surface dimple appeared above the reservoir intake for the first time indicating the start of the instability stage (Figure 8b). During this period, the free surface vortex underwent a periodic process with the repeated inception, enhancement, attenuation, and extinction of vortex instability. Although the free surface showed a periodic pattern with vortex formation and extinction, some physical quantities changed gradually over the periods. In particular, the vertical vorticity, the z component of vorticity, gradually increased in a wave-like manner. In this stage, the free surface vortex was extremely sensitive to various disturbances and changes of field quantities in the flow field. The phenomena of sensitivity and instability are described in detail in Section 3.2. Compared with the phenomenon of symmetrical inflow in Figure 9a, many small vortices appeared near the left wall of the wall in Figure 9b during this period. Meanwhile, the intensity of rotation of the inflow streamlines above the intake pipe increased and the inflow leaned to the left side.

3. Vortex stability stage. The stage consisted of a well-developed dimple and a fully developed air-core vortex. Unlike the previous stage, where dimples periodically generated and vanished, the dimple was stably in this stage. Over time, the tip of the dimple gradually extended to the intake pipe, and air was continuously entrained by the air-core vortex. A fully developed air-core vortex formed, as shown in Figure 8c. Compared with Figure 9a and b, the degree of mainstream deflection in the upstream reservoir during this period increased, and the curvature of the vortex streamlines on the plane above the intake pipe was greater, resulting in a greater intensity of the vortex, as shown in Figure 9c.

Figure 8. The free surface at the three stages during the evolution of the air-core vortex: (a) vortex inception stage, tv_D/D = 134.7; (b) vortex instability stage, tv_D/D = 307.8; (c) vortex stability stage, tv_D/D = 1092.7.

Figure 9. Streamlines on the plane z = 2.38D: (a) vortex inception stage, tv_D/D = 134.7; (b) vortex instability stage, tv_D/D = 307.8; (c) vortex stability stage, tv_D/D = 1092.7. The intake pipe is illustrated with the gray shadow.

The Q criterion is the second matrix invariant of the characteristic equation of the velocity gradient tensor, which can display the vortex structure and its intensity (Hunt et al., 1988). Figure 10 compares the vertical vorticity and Q criterion when the fully developed air-core vortex occurred. The distributions of the vertical vorticity and Q criterion were consistent in the region where the air-core vortex occurs. When the vertical vorticity was larger, the Q value in the corresponding region was larger. From Figure 10a, it could be found that the high Q value was concentrated in the air region, thus showing the swirling strength of the air, not that of the water. In Figure 10b, the vertical vorticity did not jump obviously at the air–water interface, and the water region also had high vertical vorticity, indicating that the vertical vorticity could effectively represent the structure of the free surface vortex. In addition, for a bathtub flow with bottom pipe outflow, based on some simplified assumptions, Burgers (1948) stated that the formation and development of air-core vortex were controlled by vertical vorticity ω_z, the component of vorticity in the direction of gravity. Because the tip of the air-core vortex in the reservoir with horizontal intake air was curved to a certain extent, the vortex was affected by other vorticity components. However, previous investigations showed that vertical vorticity was the most critical physical quantity (Burgers, 1948; Zi et al., 2020).

Figure 10. Comparison of the Q criterion and vertical vorticity (tv_D/D = 1092.7): (a) Q criterion; (b) vertical vorticity (ω_z).

3.2. Description of the unsteady vortex

Figure 11 shows the time record of the maximum and minimum vertical vorticity on the plane Z = 2.37D during the whole air-core vortex formation process. The minimum vertical vorticity fluctuated around the mean value of −0.67 v_D/D during the formation of the air-core vortex. When the flow time was less than 1050 D/v_D, the maximum vertical vorticity increased in a fluctuating pattern and then decreased and stabilized after reaching a peak of 5.02 v_D/D. The vortices at the water–air interface of the reservoir showed extremely transient and unstable characteristics from the generation of the surface swirl to the formation of an air-core vortex. In the vortex instability stage, the free surface vortex underwent inception, enhancement, attenuation, and extinction processes, which appeared for three times in total. After that, the free surface vortex enters the stability stage. The start and end times of these three times are marked in the time record of the vertical vorticity and named phases C1, C2, and C3, respectively, as shown in Figure 11. Another dimple did not immediately form above the intake pipe when the dimple disappeared in this phase, but it regenerated after a short period. Therefore, the boundary between phases C1 and C2 means that the vortex generated in the C1 phase disappeared, which was comparable to the boundary between phases C2 and C3. In phase C1 (282.2 < tv_D/D < 404.6), the values and fluctuations of the maximum vertical vorticity were relatively small. In phase C2 (404.6 < tv_D/D < 646.2), the maximum vertical vorticity fluctuated greatly, the value of which increases. In phase C3 (646.2 < tv_D/D < 704.1), the duration decreased significantly compared with the previous two phases, where the maximum vertical vorticity showed the largest values and fluctuations. The variation pattern of the maximum vertical vorticity was very similar in the three phases. In each phase, the maximum vertical vorticity increased first and then decreased. The boundary position of each phase representing vortex formation or disappearance was consistent with the low position of the maximum vertical vorticity. The appearance and disappearance of the free surface vortex were consistent with the synchronous relationship between the vertical vorticity, which is conformed to our previous discussion that the vertical vorticity could accurately reflect the swirling strength of the vortex. In the following discussions, we indiscriminately use surface vortex to represent the free surface vortex in the instability stage, which may include surface swirls and dimple.

Figure 11. The time record of the maximum and minimum vertical vorticities on the plane z = 2.37D.

3.3. Mechanism of the unsteady vortex

Figure 12 shows the free surface flow at different moments in phase C2, reflecting the process of surface vortex generation and elimination in this phase. In Figure 12a and b, the free surface deformation became high, indicating that the swirling strength of the surface vortex was increasing, and it reached the peak at the time of 498.9 D/v_D, as shown in Figure 12c. When tv_D/D = 581.0, there was a small dimple near the surface vortex, which was caused by swirling flow, but the rotation direction was opposite to the main surface vortex. The reverse vortex could also be observed in the experiment. Figure 12f visualizes a secondary vortex from the experiment result of Möller (2013). He stated that the pattern was unstable and remained only a short duration when a secondary vortex was burned. The experimental observations from Möller (2013) strongly support our discussions on the instability of vortex formation process. Figure 13 shows the water surface and streamlines of the counterrotating pair (tv_D/D = 606.0, 611.1, 619.4), where the pair of vortices rotated in opposite directions. The interaction between the opposite rotations was intense, causing large deformation and undulation of the free surface. Downstream, the counterrotating pair generated a crest, and the crest disappeared after reaching its highest value (tv_D/D = 619.4). During this time, both the main vortex and the secondary vortex were moving simultaneously towards the back wall of the reservoir, slowly moving away from the water inlet, and the distance between them decreased. The streamline chart can not only show the motion path of the counterrotating pair in this process but also reflect the size of the vortex region. As the main and secondary vortices approached each other, the swirling region decreased.

Figure 12. Surface vortex in phase C2: (a) inception, tv_D/D = 439.9; (b) enhancement, tv_D/D = 498.9; (c) peak, tv_D/D = 498.9; (d) attenuation, tv_D/D = 581.0; (e) extinction, tv_D/D = 633.6; (f) sencondary vortex in experiment (Möller, 2013).

Figure 13. Water surface and streamlines of the counterrotating pair (tv_D/D = 606.0, 611.1, 619.4).

The trajectories of the main and secondary vortices in phase C2 are shown in Figure 14a, and the time interval of each node in the figure was 6.41 D/v_D. The red and blue curves are the trajectories of the main vortex, and the secondary vortex, respectively. Consistent with the description in Figure 13 above, over time, the counterrotating pair moved gradually toward the back wall, and the vortices approached. In addition, Figure 14a shows that the trajectory of the secondary vortex was longer than that of the main vortex. Figure 14b is a schematic diagram of the interaction between the counterrotating pair with different swirling strengths. The swirl strength of the vortex can be expressed by the circulation in the circle with a radius of the vortex core; the calculation formula is shown in Equation (10). (10) $\mathrm{\Gamma }=2\pi R_c{v}_{\theta },$(10) where Γ donates circulation; v_θ donates maximum tangential velocity; R_c donates characteristic radius of vortex core when tangential velocity has its maximum value. The two counterrotating vortices with different circulations are shown in Figure 14b. The red circle represents the main vortex, whose circulation was larger than that of the blue circle, representing the secondary vortex. The two vortices with opposite rotation directions and different circulations induced motions on each other, making the pair rotate clockwise around the center Xi at rotational angular velocity Ω. The formulas for the calculation of rotation angular velocity and rotation center are given in Equation (11) and (12), respectively. Compared with the main vortex, the secondary vortex with weak swirling strength moved faster and farther. However, the distribution of nodes in Figure 14a is not uniform. When the counter-rotating pair was moving, the energy of each other was constantly exchanged and dissipated, which resulted in the circulations of the pair changing greatly, leading to the variations in the rotational angular velocity Ω and the pair’s center of rotation Xi. (11) $\begin{array}{rl}& \mathrm{\Omega }=({\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2)/(2\pi L^2)\end{array}$(11) (12) $\begin{array}{rl}& Xi=({\mathrm{\Gamma }}_1{X}_1^c+{\mathrm{\Gamma }}_2{X}_2^c)/({\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2)\end{array}$(12) where L donates distance between the pair; the subscript 1 and 2 donates the main and secondary vortices, respectively.

Figure 14. Analysis of the counterrotating pair: (a) the trajectory of the pair; (b) schematic diagram of the interaction between the pair.

In phase C2, the instability of the main vortex disturbed by the secondary vortex was revealed. As stated in the previous paragraph, the morphology was that of a vertical downward air-core vortex, the stability of which depended on the magnitude of the vertical vorticity. After the secondary vortex was generated, the main vortex was affected by the secondary vortex, the vortex center position of the main vortex shifted and the swirling strength decreased significantly. Finally, the dimple above the reservoir intake disappeared as time passed. Therefore, the presence of the secondary vortex inevitably affected the vertical vorticity of the main vortex, and it was necessary to adopt a quantitative method to analyze the evolution of the vertical vorticity of the counterrotating pair in phase C2. The vertical vorticity transport equation was adopted to examine the mechanism of the interaction of the pair. The vertical vorticity transport equation is the z component of the vorticity transport equation, written as Equation (13): (13) $\frac{\text{d}\mathit{\omega }}{\text{d}t}=(\mathit{\omega }\cdot \mathrm{\nabla })\mathit{v}-\mathit{\omega }(\mathrm{\nabla }\cdot \mathit{v})+\frac{\mathrm{\nabla }\rho \times \mathrm{\nabla }p}{{\rho }^2}+\nu {\mathrm{\nabla }}^2\mathit{\omega }$(13) The term on the left-hand side is the material derivative, comprising the time variation and the convection of vorticity. The first term on the right hand represents both stretching and tilting of the incompressible fluid vortex, which was also the main subject of the analysis in this paper. The second term represents the vorticity variation caused by the compressibility of the fluids. Since the two-phase flow simulated in this paper is incompressible, the volumetric expansion/contraction influence is ignored. The third term indicates both the gravitational term and the effect of the pressure relative to the hydrostatic pressure, which was also negligible in this work. The last term is vorticity variation generated by the viscous diffusion of the vorticity. In this paper as well as many other previous studies, in high Re number flows, the last term plays only a minor role in vorticity transport compared to the other three terms (Ji et al., 2015; Martins et al., 2020; Wojcik & Buchholz, 2014). Therefore, the z component of the simplified vorticity transport equation can be simplified as Equation (14): (14) $\frac{\text{d}{\omega }_z}{\text{dt}}={\omega }_x\frac{\mathrm{\partial }w}{\mathrm{\partial }x}+{\omega }_y\frac{\mathrm{\partial }w}{\mathrm{\partial }y}+{\omega }_z\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$(14) Figure 15 shows the distributions of ${\omega }_z$, ${\omega }_z\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$, and ${\omega }_x\frac{\mathrm{\partial }w}{\mathrm{\partial }x}+{\omega }_y\frac{\mathrm{\partial }w}{\mathrm{\partial }y}$ on the free surface at four moments during the interaction of the pair. The deformation of the free surface affected by the vortex pair motion is discussed in detail in Figure 14 and will not be repeated here. The vertical vorticity of the main vortex gradually decreased along with the proximity of the two vortices, which almost disappeared at 655.5 D/v_D. Contrary to the main vortex, the vertical vorticity of the secondary vortex changed from weak to strong during this process. The last term of the right-hand side of the vertical vorticity transport equation ${\omega }_z\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$ represents the vortex stretching motion. A comparison between the distributions of the vertical vorticity and the stretching term at each moment showed that they reached a state of synchronization, indicating that the vorticity stretching term was dominant in the development of the vertical vorticity. The main and secondary vortices with different vorticity magnitudes were in proximity to each other as they moved toward the back wall. The tilting term of the pair changed significantly when the primary and secondary vortices approached each other. At the initial time of vortex pair generation (tv_D/D = 439.9), the tilting term at the center of the vortex core was relatively small, and increased rapidly as the two approached, indicating that the tilt term was the source of power that pulled the vortex pair in opposite directions.

Figure 15. Distributions of ${\omega }_z$, ${\omega }_z\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$, and ${\omega }_x\frac{\mathrm{\partial }w}{\mathrm{\partial }x}+{\omega }_y\frac{\mathrm{\partial }w}{\mathrm{\partial }y}$ on the free surface.

To show the variation in the circulations of the counterrotating pair in more detail, Figure 16 is plotted to show the maximum tangential velocity and circulations of the counterrotating pair in the same period as Figure 14. During the movement of the pair, the maximum tangential velocity v_θmax and circulations Γ of the main and secondary vortices varied significantly. Figure 16a illustrates that the maximum tangential velocity of the main vortex decreased rapidly and rose at the end of the period, while that of the secondary vortex showed the opposite trend. Since the maximum tangential velocity v_θmax is a part for the calculation of the circulation Γ, the variations of the main and secondary vortex circulations Γ are similar to that of the maximum tangential velocity v_θmax, except for the difference in the variation amplitude, as shown in Figure 16b. Decreases of the circulation value of the main vortex over time were larger, and the revaluation at the end of the period was smaller. This was consistent with the evolution of the air-core vortex in the instability stage: the swirling strength of the main vortex began to decrease at a certain moment, until the surface vortex disappeared, and then energy accumulated again (the circulation rises), preparing for the next vortex generation. The circulation of the secondary vortex increased first and then decreased, and the overall variation was not large. In the whole process of pair motion, the two curves crossed, reflecting the motion and energy exchange under the interaction between the counter-rotating pair. When the vortex pair were close to each other, kinetic energy was transferred from the main vortex with the high swirling strength to the secondary vortex with the low swirling strength. At approximately 622.0 D/v_D, the two circulations were similar, reaching equilibrium, and the main vortex returned to its original position at 627.8 D/v_D as shown in Figure 14a.

Figure 16. The maximum tangential velocity and circulations of the counterrotating pair in phase C2: (a) maximum tangential velocity; (b) circulation.

Six points in the reservoir centrosymmetric line are plotted against the direction of water flow to monitor the level of the water surface, as shown in Figure 17. The distance between point 1 and the back wall was 1D. The distance between point 2 and the back wall was 3D. Points 2 through 4 were equally distributed at intervals of 1D, and points 4 through 6 were uniformly allocated at intervals of 2D.

Figure 17. Diagram of monitoring points and development of surface separation wave.

The time records of the water level at the monitor points during the air-core vortex formation process is described in Figure 18. The water level changes at each measurement point were irregular and chaotic, and thus signal analysis was needed to filter out the white noise to obtain realistic and reliable water level variation. All monitoring data were processed using the Savitzky–Golay method, which featured low-pass filtering to effectively maintain wave shape and wave height (Schafer, 2011). Based on the Savitzky–Golay method, the SG-L curve was obtained. It is obvious that during the formation of the air-core vortex, there was always a fluctuation of the water level in the entire reservoir around H = 2.5D. Frequency analysis was conducted on the SG-L curve, and the result showed that the dominant frequency was 0.285 Hz. Guo et al. (2020) studied the formation of an air-core vortex in a sump with an asymmetric vertical outlet pipe and found that waves were generated by reflection from the back wall and propagated from downstream to upstream. This type of reflection wave propagation from downstream to upstream was also present in horizontal inlet reservoirs. Figure 19 shows the water level variation at different points (from point 2 through point 4) in the instability stage, with the water levels at the three points staggered by 0.04D in the vertical direction for ease of observation and discussion. The arrow pointing from point 2 to point 4 indicated that although all points experienced the same water level fluctuation, there was a certain delay in the wave peaks and troughs at different locations. The delay at different points indicated that the wave motion moved in the opposite direction to the inlet flow, from downstream to upstream. By calculating the ratio of the distance between the different monitor points and the time at which the wave crest occurs, the moving speed of the reflection wave was found to be 0.35 v_D, which far exceeded the reservoir’s inlet flow speed (0.08 v_D). The retrograde waves in this paper were also called surface separated waves, as shown in Figure 17. The surface separation wave originated at the current dividing line, where the upstream side of the dividing line was the flow from the reservoir and the downstream side was the flow from the back wall. The two sides of the flow met and consisted of a stationary point with zero velocity. The upstream reservoir flow reached the flow divider and was reflected, forming a surface separation wave. The height of the crest and trough of the wave decreased as it traveled upstream in the reservoir, further indicating that the surface separation wave was a reflection wave propagating from the downstream to the upstream in a continuous decay. Guo et al. (2020) suggested that if steady air entrainment is finally achieved, then the pressure fluctuation (waves) at different locations reach a state of synchronization after some time. However, during the evolution of the air-core vortex, unlike Guo’s conclusion, there was always a delay in the water level fluctuation in our study. The possible reason was that the different approaches of the intake/outlet pipe led to the presence of surface separation waves at different locations, which in turn caused a difference in the wave motion.

Figure 18. Time records of the water level at the monitor points.

Figure 19. Temporal water level variations at different points in the instability stage.

The circles in Figure 19 show the time when the surface vortex disappeared in phase C1–3. When the water level was high, the surface vortex disappeared in phase C1. The increased depth of submergence was not conducive to the development of the air-core vortex. Figure 20 shows the water surface profile variation at different times in phase C1, from the rise of the water level to the disappearance of the surface vortex. As the level rose, the surface vortex gradually disappeared. Figure 21 shows the time evolution of the vertical vorticity and velocity gradient as the depth of submergence increased in phase C1, with the moment consistent with Figure 20. With the increase in submergence, the vertical vorticity near the surface vortex decreased significantly, and the surface vortex formed above the intake pipe gradually changed from a dimple and disappeared. From the analysis of the vertical vorticity transport equation in Equation (14), it can be concluded that the stretching term ${\omega }_z\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$ dominated the development of the vertical vortex, which consisted of the product of two terms ${\omega }_z$ and $\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$, where $\frac{\mathrm{\partial }w}{\mathrm{\partial }z}$ is the vertical gradient of the vertical velocity component. From 336.7 D/v_D to 404.6 D/v_D, as the depth of submergence increased, the vorticity near the surface vortex and the flow velocity gradient above the intake decreased significantly. As the product of the two terms, the stretching term decreased, followed by the decrement of the vertical vorticity, the intensity of the surface vortex decreased, and the dimples disappeared. The surface vortices in phases C2 and C3 also disappeared at smaller submergence depths, indicating that the effect of water level fluctuation on the air-core vortex diminished with time. Compared to phases C2 and C3, the vortex formed in phase C1 was not as strong, and the vortex was more susceptible to external operating conditions.

Figure 20. Water surface profile variation at different times in phase C1.

Figure 21. Temporal evolution of the vertical vorticity and velocity gradient.

4. Conclusions

In this paper, the LES and CLSVOF methods for complex two-phase flows were combined and applied in numerical simulations to study the air-core vortex phenomenon at reservoir intakes. The air-core vortex formation process was divided into three stages, and the vortex instability during the instability stage was described and analyzed. The main results were summarized as follows:

1. In comparison with the model test photographs, the numerical calculation of the water–air interface morphology was in good agreement with the experiment. The air-core vortex formation process was divided into three stages, and the vortex instability during the instability stage was described and analyzed. During the instability stage, the surface vortex was repeatedly subjected to the process of inception-enhancement-attenuation-extinction. Through detailed analysis of periodic vortex generation and extinction during the instability stage, we found that the movement of the counterrotating vortex pair and water surface level fluctuation could lead to the reduction of the vertical vorticity.

2. The flow induced by the main and secondary vortices of the counterrotating vortex pair with unequal amounts of circulation caused the secondary vortex to move toward the back wall faster than the main vortex. At the same time, energy transfer occurred between the main and secondary vortices, where the swirling strength of the main vortex decreased with the increment of the swirling strength of the secondary vortex. The relative motion of the vortex pair was governed by the tilting term of the vertical vorticity, while the reduction in the stretching term eventually led to a reduction in the vorticity and the disappearance of the dimple.

3. Water level fluctuations briefly caused an increase in the submergence depth at the beginning of the instability stage. This disturbance from the rising water level led to a decrease in vortex intensity. As the water level rose, the flow velocity gradient near the free surface vortex above the intake decreased significantly, which in turn reduced the vertical vorticity stretching term and weakened the vertical vorticity.

4. We described the phenomenon of vortex instability during the formation of the intake air-core vortex, and revealed the formation mechanism of the unstable vortex, which enriched the research of free surface vortex in the field of instability. These results could provide new insight into the optimal design and operation of hydraulic engineering and anti-vortex analyses. For example, water surface fluctuations could be artificially created in areas where vortices were prone to occur to destroy the formation of vortices. An injected jet device installed at the intake could generate additional momentum that worked the same as a reverse vortex and suppress the vortex strength. This anti-vortex measure has been verified by Monshizadeh et al. (2017) to be more flexible and more hydraulically efficient than the structural based anti-vortex methods. In addition, the numerical model combining WMLES and CLSVOF free-surface tracking method used in this paper has been proven to have a good calculation effect, which provided a reliable way for researchers to analyze related problems such as free-surface vortices.

5. The instability of the free surface vortex discussed in this paper was examined under a single hydraulic condition with dimensionless numbers of Fr = 2.89, Re = 2.04 × 10⁶, We = 1.46 × 10⁵, S/D = 1.5, and φ = 0°. These dimensionless numbers controlled the generation and development of the air-core intake vortex, so there must be a certain correlation between the vortex instability phenomenon and these dimensionless numbers. Therefore, in the follow-up study, more hydraulic conditions will be designed to construct the relationship between dimensionless numbers and the intake vortex instability characteristic.

Nomenclature

${\tilde{S}}_{ij}$	=	Large-scale strain rate tensor
D	=	Inlet pipe diameter (m)
dw	=	Wall distance
F	=	Volume fraction function
g	=	Gravitational acceleration (m·s−2)
hmax	=	Maximum edge length for a rectilinear hexahedral cell
hwn	=	Wall-normal grid spacing
p	=	Pressure (Pa)
Rc	=	Characteristic radius of vortex core (m)
S	=	Submergence depth (m)
t	=	Time (s)
ui, uj	=	Instantaneous velocity component in the Cartesian coordinate system (m·s−1)
v	=	Kinematic viscosity of the fluids (Pa's)
vD	=	Inlet time-averaged velocity (m/s)
vθ	=	Maximum tangential velocity (m/s)
y+	=	Distance from the first layer grid to the wall surface in wall units
νt	=	Sub-grid eddy viscosity (Pa's)
ρ	=	Density (kg·m−3)
σ	=	Surface tension (N/m)
τij	=	Subgrid-scale stress tensor
τkk	=	Isotropic part of SGS stress
ϕ	=	Signed distance function
φ	=	Intake pipe direction (°)
ω	=	Vorticity (m/s2)
Ω	=	Rotational angular velocity (m/s)
Γ	=	Circulation (m2/s)

Abbreviations

CLSVOF	=	Coupled level-set and volume-of-fluid
Fr	=	Froude number
LES	=	Large eddy simulation
LES-IQ	=	LES index of quality
NΓ	=	Circulation number
Re	=	Reynolds number
UDF	=	User defined functionlfi
VOF	=	Volume-of-fluid
We	=	Weber number
WMLES	=	Wall-modeled LES

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China [grant number 52006053,52009033]; Natural Science Foundation of Jiangsu Province [grant number BK20200508,BK20200509]; Postdoctoral Research Foundation of China [grant number 2022T150185,2021M690876,2022M711021].