1. Introduction
Despite its simple set-up, single-phase flow past a circular cylinder exhibits a diverse range of intricate fluid-flow phenomena that have fascinated researchers for over a century (Strouhal Reference Strouhal1878; von Kármán & Rubach Reference von Kármán and Rubach1912; Berger & Wille Reference Berger and Wille1972; Williamson & Govardhan Reference Williamson and Govardhan2004). The behaviour of isothermal single-phase flow around a circular cylinder is determined entirely by the Reynolds number
${\textit{Re}} = U_0 D / \nu$
, where
$U_0$
represents the velocity of the incoming uniform flow,
$D$
is the cylinder diameter, and
$\nu$
is the kinematic viscosity of the fluid. For instance, the onset of two-dimensional vortex shedding occurs at
${\textit{Re}} \approx 47$
, followed by the development of three-dimensional wake structures in the range
${\textit{Re}} \approx 190$
to
$260$
. Eventually, the flow in the wake region becomes turbulent at higher
${\textit{Re}}$
values (Roshko Reference Roshko1954; Bloor Reference Bloor1964; Grove, Shair & Petersen Reference Grove, Shair and Petersen1964; Gerrard Reference Gerrard1966; Dennis & Chang Reference Dennis and Chang1970; Fornberg Reference Fornberg1985; Dušek et al. Reference Dušek, Gal and Fraunié1994; Norberg Reference Norberg1994; Mittal & Balachandar Reference Mittal and Balachandar1995; Zhang et al. Reference Zhang, Fey, Noack, König and Eckelmann1995; Barkley & Henderson Reference Barkley and Henderson1996; Henderson Reference Henderson1996; Thompson, Hourigan & Sheridan Reference Thompson, Hourigan and Sheridan1996; Williamson Reference Williamson1996b
; Zdravkovich Reference Zdravkovich1998; Rajani, Kandasamy & Majumdar Reference Rajani, Kandasamy and Majumdar2009; Chopra & Mittal Reference Chopra and Mittal2019). Williamson (Reference Williamson1996b
) and Zdravkovich (Reference Zdravkovich1998) systematically reviewed the studies on single-phase flow past a cylinder, covering theoretical analyses, numerical simulations and experimental observations.
The present computational work centres around two-phase liquid–gas flow across a stationary cylinder in a three-dimensional setting, wherein a horizontal circular cylinder is inserted beneath the free surface within the liquid cross-flow. The density and viscosity contrasts between the liquid and gas phases in our system resemble those of water and air. Prior experimental (Sheridan et al. Reference Sheridan, Lin and Rockwell1995, Reference Sheridan, Lin and Rockwell1997) and numerical (Reichl, Hourigan & Thompson Reference Reichl, Hourigan and Thompson2005; Colagrossi et al. Reference Colagrossi, Nikolov, Durante, Marrone and Souto-Iglesias2019; Zhao et al. Reference Zhao, Wang, Zhu, Ping, Bao, Zhou, Cao and Cui2021; Patel et al. Reference Patel, Sun, Yang and Zhu2025) studies involving a similar set-up have reported substantial changes in the flow pattern around the cylinder and organisation of the wake structure in comparison with the single-phase counterpart. These modifications stem from the coupling between the free surface and flow perturbations arising from the cylinder. Compared with single-phase flow, the flow dynamics of two-phase flow is influenced by multiple parameters other than the Reynolds number. To account for the elevation of the free surface, buoyancy and surface tension forces in a two-phase set-up, we define three dimensionless parameters in addition to
${\textit{Re}}$
: the gap ratio
$G=h/D$
, the Froude number
${\textit{Fr}}=U_{0}/\sqrt {gD}$
and the Weber number
${\textit{We}}=\rho _{w}U^{2}_{0}D/\sigma$
. Here,
$h$
,
$g$
,
$\rho _{w}$
and
$\sigma$
denote the distance between the top of the cylinder and the unperturbed free surface, gravitational acceleration, the density of water and the surface tension of the free surface, respectively. Note that the kinematic viscosity in
${\textit{Re}}$
corresponds to that of water,
$\nu _{w}$
, in the current two-phase setting. The Froude number, which reflects the balance between inertial and gravitational forces, plays an important role in free-surface flows dominated by large-scale gravity waves, including phenomena such as sloshing and surfing in open vessels (Yu & Tryggvason Reference Yu and Tryggvason1990; Lugt & Ohring Reference Lugt and Ohring1992; Oguz, Prosperetti & Kolaini Reference Oguz, Prosperetti and Kolaini1995; Ye & Chu Reference Ye and Chu1997; Kiger & Duncan Reference Kiger and Duncan2012; Li, Yang & Zhang Reference Li, Yang and Zhang2024). The Weber number characterised by surface tension mainly influences the air entrainment and bubble fragmentation, dealing with small-scale surface effects such as ripples, droplets and jets. In scenarios with strong air entrainment, the surface tension imposes significant impacts on hydrodynamic loads and noises (Kiger & Duncan Reference Kiger and Duncan2012; Veron Reference Veron2015; Drenckhan & Saint-Jalmes Reference Drenckhan and Saint-Jalmes2015).
Over the years, low-Reynolds-number (
${\textit{Re}} \leqslant 250$
) free-surface flow past a circular cylinder has been studied extensively from various perspectives. Reichl et al. (Reference Reichl, Hourigan and Thompson2005) examined scenarios with marginally deformable free surfaces, classifying wake states into full, partial and no vortex shedding based on different combinations of
$G$
and
${\textit{Fr}}$
. These regimes were reinforced by Colagrossi et al. (Reference Colagrossi, Nikolov, Durante, Marrone and Souto-Iglesias2019), who also paid particular attention to the metastable state – identified by a secondary frequency in the lift-force signal that is typically an order of magnitude lower than the primary vortex-shedding frequency. Later, González-Gutierrez et al. (Reference González-Gutierrez, Gimenez and Ferrer2019) conducted the global stability analysis to assess how free-surface deformation influences the onset of vortex shedding. In our recent work (Patel et al. Reference Patel, Sun, Yang and Zhu2025), we described a broad spectrum of free-surface dynamics in liquid–gas flows past a circular cylinder, ranging from standing and travelling interfacial waves in highly deformable free surfaces to multiscale gas-bubble entrainment caused by the breakup of the free surface. Besides previous studies, researchers have also explored other related configurations, which include a cylinder in shallow free-surface flows (Moballa, Chern & Borthwick Reference Moballa, Chern and Borthwick2020), a rotating circular cylinder (Guo et al. Reference Guo, Ji, Han, Liu, Wu and Kuai2023; Lin & Yao Reference Lin and Yao2023), a pair of stationary tandem cylinders (Zhao et al. Reference Zhao, Wang, Zhu, Cao, Bao, Zhou, Cheng and Han2024), and fixed elliptic (Subburaj, Khandelwal & Vengadesan Reference Subburaj, Khandelwal and Vengadesan2018) and square (Karmakar & Saha Reference Karmakar and Saha2020) cylinders close to the free surface.
High-Reynolds-number (
${\textit{Re}} \gt 2000$
) free-surface flow past a circular cylinder has also attracted research interest over the years. Sheridan et al. (Reference Sheridan, Lin and Rockwell1995, Reference Sheridan, Lin and Rockwell1997) employed particle image velocimetry to measure the free-surface flow around cylinders at medium Froude numbers (
${\textit{Fr}} = 0.6{-}0.8$
) and in the high-Reynolds-number regime at
$5990 \leqslant {\textit{Re}} \leqslant 9120$
. They identified the presence of an oblique downward shear layer originating from the water surface, resulting in a jet-like flow with nearly zero vorticity between the cylinder and water surface. They also described a metastable flow state characterised by two distinct wake states. Each state exhibited limited stability, leading to time-dependent transitions between the two states. Miyata, Shikazono & Kanai (Reference Miyata, Shikazono and Kanai1990) performed experiments to study the flow around a horizontal cylinder near the free water surface at
${\textit{Re}}=4.96 \times 10^4$
. Their result suggests that under the shallow-submerging condition, the difference in the flow between the top and bottom sides of the cylinder causes asymmetric vortex distribution about the central horizontal plane, which further leads to a reduction of drag coefficient in comparison with single-phase flow. More recently, Zhao et al. (Reference Zhao, Wang, Zhu, Ping, Bao, Zhou, Cao and Cui2021) performed large-eddy simulations with the same parameters as the experiments of Sheridan, Lin & Rockwell (Reference Sheridan, Lin and Rockwell1997), enabling a detailed analysis of free-surface deformation and the forces acting on the cylinder. They observed significant free-surface deformation at
$G = 0.4$
and
${\textit{Fr}}=0.6$
due to the local Froude number exceeding unity in the gap above the cylinder. Zhao et al. (Reference Zhao, Wang, Zhu, Cao, Bao, Zhou and Han2022) conducted an unsteady Reynolds-averaged numerical simulation at a higher Reynolds number of
${\textit{Re}}=4.96 \times 10^4$
. Their findings indicate that the decreasing gap ratio not only shifts the front stagnation point towards the free surface, creating an asymmetrical pressure distribution on the cylinder, but also introduces a significant free-surface proximity effect. This effect modulates the hydrodynamics, causing the force time histories to deviate from a sinusoidal form and become irregular. Alzabari, Wilson & Ouro (Reference Alzabari, Wilson and Ouro2023) also conducted large-eddy simulations of shallow free-surface flows past a circular cylinder at
${\textit{Re}} = 13\,333$
. They observed that free-surface deformation at Froude numbers starting from
${\textit{Fr}}=0.53$
leads to a loss of coherence in the wake vortices, which eventually merge with the ground vortices. Moreover, using proper orthogonal decomposition, they showed that the energy contained within the first couple of modes drops significantly at elevated Froude numbers. Other high-Reynolds-number studies include free-surface flows around spheres (Sareen et al. Reference Sareen, Zhao, Sheridan, Hourigan and Thompson2018; Chizfahm, Joshi & Jaiman Reference Chizfahm, Joshi and Jaiman2021; Hunt et al. Reference Hunt, Zhao, Silver, Yan, Bazilevs and Harris2023) and flat plates (Díaz-Ojeda et al. Reference Díaz-Ojeda, Huera-Huarte and González-Gutiérrez2019; Hu et al. Reference Hu, Liu, Zhao and Hu2023).
From the research progress reviewed, it is evident that previous studies focused mainly on either laminar flow with low Reynolds numbers (
${\textit{Re}}\leqslant 250$
) or turbulent flow with high Reynolds numbers (
${\textit{Re}} \gt 2000$
). Furthermore, the Reynolds number was fixed at a constant in most previous studies to investigate the effect of other parameters on the flow dynamics. In two-phase flow, the turbulence can be generated either in the wake of the cylinder or in the air–water mixed region due to air entrainment (Kiger & Duncan Reference Kiger and Duncan2012; Veron Reference Veron2015; Drenckhan & Saint-Jalmes Reference Drenckhan and Saint-Jalmes2015; Hendrickson et al. Reference Hendrickson, Weymouth, Yu and Yue2019b
; Li et al. Reference Li, Yang and Zhang2024). Considering that the Reynolds number imposes significant impacts on the status of the onset and intensity of turbulence, it is useful to investigate the flow dynamics at transitional Reynolds numbers.
In the present study, we focus on two-phase flow past a fixed submerged cylinder at transitional Reynolds number. A fixed gap ratio of
$G = 0.5$
, Weber number of
${\textit{We}} = 1000$
and Froude number of
${\textit{Fr}} = 1$
are maintained. The Reynolds number varies from 400 to 2000, at which wake turbulence occurs in single-phase flows. The relatively low Reynolds numbers allow us to perform direct numerical simulations (DNS), and as such numerical error corresponding to any turbulence models can be avoided. The main purpose of the present study is to examine the effect of Reynolds number on the interaction between the vortex shedding and air entrainment. At
${\textit{Re}} = 400$
, suppressed vortex shedding allows quasisymmetric shear layers to meander considerably downstream. Conversely, increasing
${\textit{Re}}$
from 500 to 2000 promotes synchronous vortex shedding of the separated shear layers and the shear layer originating from the free surface, leading to intermittent air entrainment. The investigation extends to the Reynolds number’s influence on bubble and velocity statistics.
The remainder of the paper is organised as follows. Section 2 provides a description of the numerical methods and case set-up. In § 3, the results, including the hydrodynamic forces of the cylinder, vortex shedding, air entrainment and statistics of the velocity, are analysed. The conclusions are summarised in § 4.
2. Numerical methods and case set-up
The simulations were conducted using an in-house developed code Computational Air-Sea Tank (CAS-Tank, Yang, Lu & Wang Reference Yang, Lu and Wang2021). Let
$t$
and
$\boldsymbol{x}= [x,y,z ]$
denote time and three-dimensional Eulerian coordinates, respectively; an incompressible two-fluid flow is governed by the following continuity and momentum equations:
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0 , \end{align}
\begin{align} \frac {\partial (\rho \boldsymbol{u})}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\rho \boldsymbol{u}\boldsymbol{u}\right ) &= -\boldsymbol{\nabla}\!p + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (2\mu \unicode{x1D64E}\right ) + \boldsymbol{f}_{\!s} + \rho \boldsymbol{g} , \end{align}
where
$\boldsymbol{\nabla }= [\partial /\partial x, \partial /\partial y, \partial /\partial z ]$
is the nabla symbol,
$\boldsymbol{u}= [u, v, w ]$
is the velocity,
$p$
is the pressure,
$\rho$
and
$\mu$
are density and dynamic viscosity, respectively,
$\unicode{x1D64E}= (\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^\textrm T )/2$
is the strain rate tensor, and
$\boldsymbol{g} = [0, -g, 0]$
is the gravitational acceleration. In (2.2),
$\boldsymbol{f}_{\!s}$
is the surface tension force, which is defined as
\begin{equation} \boldsymbol{f}_{\!s} = \sigma \kappa \boldsymbol{\nabla }H(\phi ), \end{equation}
where
$\sigma$
denotes the surface tension;
$\kappa = \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\nabla }\phi )$
denotes the interface curvature;
$\phi$
represents the level-set function, with
$\phi = 0$
representing the air–water interface; and
$H(\phi )$
is the Heaviside function
\begin{equation} H(\phi)= \begin{cases}0,& \phi \lt 0, \\ 1, & \phi \geqslant 0\;.\end{cases} \end{equation}
The density and viscosity are determined using the level-set function
$\phi$
as
\begin{align} \rho &= \rho _a + (\rho _{w} - \rho _a)H(\phi ) , \end{align}
\begin{align} \mu &= \mu _a + (\mu _w - \mu _a)H(\phi ), \end{align}
where subscripts
$a$
and
$w$
denote air and water, respectively. The interface is captured using the Coupled Level Set and Volume-of-Fluid method (Sussman & Puckett Reference Sussman and Puckett2000). The following convection equations for the level-set (LS)
$\phi$
and volume-of-fluid (VOF)
$\psi$
functions are solved:
\begin{align} \frac {\partial \phi }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{u})&=0, \end{align}
\begin{align} \frac {\partial \psi }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\psi \boldsymbol{u})&=0 . \end{align}
The LS function
$\phi$
refers to the assigned distance from the interface, and the VOF function
$\psi$
is the volume fraction of water in each grid cell. The LS function is used to calculate the normal vector of the interface, whereas the VOF function maintains the mass conservation. Details of the numerical method are described in Yang et al. (Reference Yang, Lu and Wang2021).
Table 1 summarises the key parameters of the simulation. In the two-phase flow past a cylinder, the flow features are controlled by six non-dimensional parameters, including Reynolds number (
${\textit{Re}}$
), Froude number (
${\textit{Fr}}$
), gap ratio (
$G$
), Weber number (
${\textit{We}}$
), density ratio (
$\rho _a / \rho _w$
) and viscosity ratio (
$\mu _a / \mu _w$
). To investigate the differences between single-phase and two-phase flow, we conducted case S400 for single-phase flow at
${\textit{Re}} = 400$
. For two-phase flow, we focus on the effect of Reynolds number, and cases T400–T2000 with the Reynolds number ranging from 400 to 2000 were conducted. Other parameters remain constants, including
${\textit{Fr}} = 1.0$
,
$G = 0.5$
,
${\textit{We}} = 1000$
,
$\rho _a / \rho _w = 1.20 \times 10^{-3}$
and
$\mu _a / \mu _w = 1.80 \times 10^{-2}$
. These settings yield strong air entrainment and a representative shear-layer interaction. The present investigation is focused on regimes characterised by relatively low surface tension and, consequently, weak capillary effects. By maintaining constant Weber (
${\textit{We}}$
) and Froude (
${\textit{Fr}}$
) numbers for all cases, the capillary length (
$\lambda_c = \sqrt{\sigma / (\Delta\rho g)}$
, with
$\Delta \rho$
being the fluid density difference) is identical and likewise constant at a value of
$0.032D$
. Since the analysis does not extend into the sub-Hinze regime, resolving the capillary length scale is not strictly necessary.
Simulation parameters for all cases. Case S400 represents a single-phase simulation. All other cases involve two-phase flow simulations with the following fixed parameters:
${\textit{Fr}} = 1.0$
,
$G = 0.5$
,
${\textit{We}} = 1000$
,
$\rho _a / \rho _w = 1.20 \times 10^{-3}$
and
$\mu _a / \mu _w = 1.80 \times 10^{-2}$
. Case T2000F employs a finer resolution of
$D/60$
, in contrast to the
$D/40$
resolution used in the other cases.
Figure 1 depicts the computational and geometrical parameters of the simulations for two-phase cases. The domain size is
$L_x \times L_y \times L_z = 45D \times 25D \times 3D$
, where
$D$
is the cylinder diameter. The origin of the coordinates is placed at the centre of the cylinder in an
$x$
-
$y$
plane, which is
$10.0D$
from the inlet and
$1.0D$
below the initial free surface. The initial liquid depth is
$17D$
. Cartesian grids are used for discretisation. The number of grid points is
$N_x \times N_y \times N_z = 1180 \times 540 \times 120$
. In the streamwise and vertical directions, the grid is refined around the cylinder and the air–water interface. The grid resolution in the subdomain with
$-1.5D \leqslant x \leqslant 20.5D$
and
$-3D \leqslant y \leqslant 3D$
is
$\varDelta _x = \varDelta _y = 0.025D$
. The grid is gradually stretched toward the four boundaries in the
$x$
- and
$y$
-directions, with the coarsest grid resolution near the boundaries being
$\max (\varDelta _x) = 0.19D$
and
$\max (\varDelta _y) = 0.19D$
. In the spanwise direction, the grid is evenly spaced, and the resolution is
$\varDelta _z = 0.025D$
.
Computational domain and geometrical parameters for the present simulations of two-phase flow past a circular cylinder.
The computational domain size and grid resolution were chosen carefully based on a convergence test. We first used the CAS-Tank to conduct a series of two-dimensional simulations to determine the computational domain size. The present domain
$45D \times 25D$
is found to result in the same drag and lift coefficients in the
$x$
-
$y$
plane as the larger domain
$80D \times 80D$
does. Subsequently, the adequacy of the spanwise domain length was assessed by examining two-point correlations of fluctuating velocity components in relevant flow regions. Figure 2 represents these correlations at the point
$(x,y)=(3D,0.1D)$
for case T2000. The correlations are observed to decay to reasonably small values over a spanwise distance of
$1.5D$
, confirming that the spanwise domain is sufficiently large to contain large-scale energetic structures. To ensure that the grid resolution is sufficiently fine, we conducted a grid-dependence test using case T2000. We conducted the simulation with three sets of grids: a base mesh (
$D/40$
) with 79 million grid points, a coarser mesh (
$D/27$
) with 21 million grid points, and a geometrically similar finer mesh (
$D/60$
) with 240 million grid points. It was discovered that when the grid resolution is reduced from
$D/27$
to
$D/40$
, the lift coefficient grows by
$35\%$
, indicating that
$D/27$
is insufficient. As we further refine the grid resolution from
$D/40$
to
$D/60$
, the lift coefficient only changes by
$3.5\,\%$
. The changes in other quantities, including the drag coefficient, root-mean-square (r.m.s.) values of the lift and drag coefficients, and the Strouhal number of vortex shedding, are all smaller than
$5\%$
. Other key parameters, including the mean velocity (figures 18, 21, 23, 24), kinetic energy of velocity fluctuations (figures 25, 26), vorticity and free-surface geometry (figures 8
e, 27), vortex structure (figures 10
e, 28), frequency spectra of the air entrainment statistics (figure 12
e) and bubble-size spectrum (figure 16) are also verified. Even though the grid resolution of
$D/40$
cannot fully resolve down to the Hinze scale – as the computational cost for such a resolution is currently prohibitive – it is sufficient to capture the large-scale bubbles accurately (Each simulation case already consumes approximately 1.8 million CPU hours). This ensures that the statistical trends of the two-phase turbulence as a function of the Reynolds number are reliable and trustworthy. Based on this test, the grid resolution of
$D/40$
is chosen to conduct simulations of other cases. The computational domain for the single-phase case (case S400) is
$L_x \times L_y \times L_z = 45D \times 32D \times 3D$
. The cylinder centre is placed
$10D$
from the inlet and
$16D$
from both top and bottom boundaries. The number of grid points is
$N_x \times N_y \times N_z = 1180 \times 560 \times 120$
. The grid is equally refined in the same region as the two-phase flow. After conducting the simulation, the results of the hydrodynamic forces acting on the cylinder of case S400 were compared with those from previous studies. As presented in table 2, it becomes evident that the results obtained from the CAS-Tank again show a reasonable level of agreement.
Variation of spanwise two-point correlations of the streamwise (
$R_{uu}$
), vertical (
$R_{vv}$
), and spanwise (
$R_{ww}$
) velocity fluctuations with spanwise separation
$\Delta_z$
. The reference point for the correlations is
$(x,y)= (3D, 0.1D)$
.
Comparison of hydrodynamic characteristics for the single-phase flow case (case S400). The time-averaged drag coefficient
$\langle C_d \rangle_t$
, r.m.s. value of lift coefficient fluctuation
$C_{l,\textit{rms}}$
, base pressure coefficient
$ C_{p_b}$
and Strouhal number
$ St$
are presented. The Strouhal number is defined as
$ St = f\!D/U_0$
, where
$ f$
is the vortex-shedding frequency,
$ D$
is the cylinder diameter, and
$ U_0$
is the free-stream velocity.
A uniform velocity profile was prescribed at the inlet, and a zero-gradient condition was applied at the outlet. In the vertical direction, a free-slip condition is adopted at both the top and bottom boundaries. In the spanwise direction, a periodic condition is used. A sharp-interface immersed boundary method is applied to enforce the no-slip condition at the cylinder surface (Cheny & Botella Reference Cheny and Botella2010; Shrivastava, Agrawal & Sharma Reference Shrivastava, Agrawal and Sharma2013; Cui et al. Reference Cui, Yang, Jiang, Huang and Shen2017). The solver’s capability for accurately simulating general free-surface flows has been extensively validated (Yang et al. Reference Yang, Lu and Wang2021), with successful applications to diverse scenarios such as plunging jets (Li et al. Reference Li, Yang and Zhang2024), wave breaking (Lu et al. Reference Lu, Yang, He and Shen2024) and bubbly channel flow (Lu, Yang & Deng Reference Lu, Yang and Deng2025). Crucially for the present study, its applicability to similar configurations was established in our prior work (Patel et al. Reference Patel, Sun, Yang and Zhu2025). Once the flow reached a statistically steady stage, 2000 instantaneous flow fields were stored to calculate turbulent statistics. The sampling rate of the data was
$\delta t = 0.05 D/U_0$
.
3. Results and discussions
Before we start quantitative analysis, we provide a description of the key features of the flow. Figure 3 shows an instantaneous snapshot of the air–water interface, together with the streamlines, coloured by the streamwise velocity. Warm and cold colours represent positive and negative values, respectively. In the upstream of the cylinder, the velocity is reduced and the pressure rises because of the obstruction effect of the cylinder. The water surface forms a bumping-up shape ahead of the cylinder due to pressure rising. In the near-wake region, the recirculation zone is downward inclined, due to the downward motion of the jet-like flow over the cylinder. Above the jet-like flow, the shear induces another recirculation zone, while the air entrainment results in a two-phase region above the jet-like flow. Similar features of the flow have also been described in previous studies (Sheridan et al. Reference Sheridan, Lin and Rockwell1995, Reference Sheridan, Lin and Rockwell1997; Bouscasse et al. Reference Bouscasse, Colagrossi, Marrone and Souto-Iglesias2017; Zhao et al. Reference Zhao, Wang, Zhu, Ping, Bao, Zhou, Cao and Cui2021; Hendrickson, Yu & Yue Reference Hendrickson, Yu and Yue2022; Patel et al. Reference Patel, Sun, Yang and Zhu2025).
Instantaneous air–water interface and streamlines for
${\textit{Re}} = 400$
(case T400). The streamlines are coloured by the streamwise velocity.
3.1. Hydrodynamic forces
We start the quantitative analysis with the hydrodynamic forces exerted on the cylinder. Figure 4 shows the time histories of lift coefficient
$C_l$
and drag coefficient
$C_d$
, defined, respectively, as
\begin{align} C_l &= \frac {F_l}{0.5 \rho _w U_0^2 A}, \end{align}
\begin{align} C_d &= \frac {F_d}{0.5 \rho _w U_0^2 A}, \end{align}
Time history of the lift coefficient (
$C_l$
) and drag coefficient (
$C_d$
) for cases detailed in table 1: (a) S400, (b) T400, (c) T500, (d) T800, (e) T1400 and ( f) T2000.
where
$F_l$
and
$F_d$
are the lift and drag forces, respectively;
$A = L_z D$
is the projected area;
$\rho _w$
is the water density; and
$U_0$
is the free-stream velocity. In all cases under investigation, the cylinder is fully submerged, and therefore the buoyancy remains constant. To focus on the hydrodynamic force,
$F_l$
is defined with the buoyancy excluded as
\begin{equation} F_l = F_y - F_b , \end{equation}
where
$F_y$
denotes the total vertical force exerted on the cylinder, and
$F_b$
is the buoyancy force, defined as
\begin{equation} F_b = \frac {\rho _w g \unicode{x03C0} {D}^2 L_z}{4} \;. \end{equation}
Figure 4 compares the time series of
$C_l$
and
$C_d$
for different cases. For single-phase flow at
${\textit{Re}}=400$
, as shown in figure 4(a), both
$C_l$
and
$C_d$
oscillate around a mean value in a quasiperiodic manner. The period of
$C_l$
is twice that of
$C_d$
. These features of single-phase flow have been reported extensively in the literature. A major difference between two-phase and single-phase flow is exhibited in the mean lift coefficient, which is zero in single-phase flow, but is negatively valued in two-phase flow. As explained by Miyata et al. (Reference Miyata, Shikazono and Kanai1990), Zhao et al. (Reference Zhao, Wang, Zhu, Ping, Bao, Zhou, Cao and Cui2021), Patel et al. (Reference Patel, Sun, Yang and Zhu2025), the negative lift coefficient results from the downward jet-like flow. For two-phase flow at
${\textit{Re}}=400$
, the magnitudes of the oscillation in both
$C_l$
and
$C_d$
are smaller than those for single-phase flow. The quasiperiodic oscillation recovers for
${\textit{Re}}\geqslant 500$
.
To further quantify the effect of Reynolds number on the lift and drag coefficients, we adopt the following decomposition:
\begin{align} C_l &= \langle C_l \rangle _{t} + C_l^{\prime \prime }, \end{align}
\begin{align} C_d &= \langle C_d \rangle _{t} + C_d^{\prime \prime }, \end{align}
where a pair of angular brackets with subscript
$t$
represent time averaging, and the double prime denotes fluctuations with respect to time averaging. The r.m.s. values of fluctuations, denoted using a subscript ‘
$\textit{rms}$
’, are then used to evaluate the magnitude of oscillation. For the single-phase case, the computed mean drag coefficient is 1.12 and the r.m.s. value of the lift coefficient is 0.25. These values demonstrate reasonable agreement with previous studies (Wieselsberger Reference Wieselsberger1922; Zdravkovich Reference Zdravkovich1998; Norberg Reference Norberg2003; Jiang & Cheng Reference Jiang and Cheng2018). Figure 5 compares the mean and r.m.s. values of the drag and lift coefficients for different two-phase cases. It is observed that as the Reynolds number increases from 500 to 2000, the mean lift coefficient
$\langle C_l \rangle _t$
increases, while the mean drag coefficient
$\langle C_d \rangle _t$
, r.m.s. value of lift coefficient
$C_{l,\textit{rms}}$
and r.m.s. value of drag coefficient
$C_{d,\textit{rms}}$
all decrease monotonically. However, the results for
${\textit{Re}}=400$
do not follow the monotonic variation. Particularly, the values of both
$C_{l,\textit{rms}}$
and
$C_{d,\textit{rms}}$
for
${\textit{Re}}=400$
are smaller than those for other Reynolds numbers. This observation indicates that the flow status changes between
${\textit{Re}}=400$
and
$500$
.
The difference in the flow status is also exhibited in the frequency of
$C_l$
and
$C_d$
. To obtain the dominant frequency, the frequency spectra of lift coefficient
$\varPhi _l$
and drag coefficient
$\varPhi _d$
are calculated, and the results are shown in figure 6. Figure 6( f) compares two grid resolutions,
$D/40$
(case T2000) and
$D/60$
(case T2000F), showing reasonably good agreement between them. It is seen from figure 6(a) that in single-phase flow, the dominant peaks of
$\varPhi _l$
and
$\varPhi _d$
occur at
$f\!D/U_0 = 0.210$
and
$0.420$
, respectively. Figure 6(b) shows that in two-phase flow at
${\textit{Re}}=400$
the primary peaks of both
$\varPhi _l$
and
$\varPhi _d$
move to a low frequency
$f\!D/U_0 = 0.034$
, and the magnitudes of the high-frequency peaks are reduced to a low level. At higher Reynolds numbers, as shown in figure 6(c–f), the low-frequency peaks of
$\varPhi _l$
and
$\varPhi _d$
are reduced, and the high-frequency peaks dominate again. This non-monotonic behaviour of the high-frequency peak is consistent with the r.m.s. values of
$C_l$
and
$C_d$
shown in figure 5(c,d), indicating again the change in the flow status between
${\textit{Re}}=400$
and 500. The low-frequency and high-frequency peaks of
$\varPhi _l$
and
$\varPhi _d$
are correlated to the quasiperiodic motions of the bulge-shape surface ahead of the cylinder and the vortex shedding, respectively. Detailed analyses of the mechanism underlying the two peaks of
$\varPhi _l$
and
$\varPhi _d$
are given in §§ 3.2 and 3.3.
3.2. Vortex structures
The high-frequency peaks in the lift and drag spectra are closely associated with vortex shedding. In single-phase flow, each shedding event on either side of the cylinder generates a horizontal pressure imbalance, corresponding to one drag cycle. During the shear-layer roll-up, vorticity accumulates and the accelerated flow creates a local low-pressure region, which induces an asymmetric vertical force. The net lift force is directed toward the side of shear-layer roll-up and opposite to the shedding side. As a result, two shedding events are required to complete one lift cycle, and the dominant frequency of the drag coefficient is twice that of the lift coefficient. In two-phase flow, the pressure field is modified by asymmetry. Throughout a lift cycle, the low-pressure region remains on the upper surface of the cylinder, and its variation in magnitude and position induces in-phase changes in both the vertical and horizontal pressure imbalances. Thus, the lift and drag coefficients oscillate synchronously and share the same dominant frequency.
Values of (a) mean lift coefficient
$\langle C_l \rangle _t$
, (b) mean drag coefficient
$\langle C_d \rangle _t$
, (c) r.m.s. lift coefficient
$C_{l,\textit{rms}}$
and (d) r.m.s. drag coefficient
$C_{d,\textit{rms}}$
for different Reynolds numbers. The results for
${\textit{Re}} = 400$
are shown using blue pentagons, while red diamonds correspond to cases with
${\textit{Re}} \geqslant 500$
.
The frequency spectra of lift coefficient
$\varPhi _l$
and drag coefficient
$\varPhi _d$
for cases detailed in table 1: (a) S400, (b) T400, (c) T500, (d) T800, (e) T1400 and ( f) T2000. Panel ( f) also compares results from two different grid resolutions:
$D/40$
(case T2000) and
$D/60$
(case T2000F).
To further demonstrate the deterministic role of vortex shedding on the high-frequency peak of the frequency spectra of lift and drag coefficients, we examine the Strouhal number, defined as
\begin{equation} St = \frac {f\!U_0}{D} , \end{equation}
where
$f$
represents the vortex-shedding frequency,
$U_0$
is the free-stream velocity, and
$D$
is the cylinder diameter. The value of
$f$
is identified using the dominant frequency of the spanwise-averaged spanwise vorticity
$\langle \omega _z \rangle _z$
at a downstream location of the cylinder (
$x / D = 5, y / D = -1$
) (marked using a cross symbol in figure 8). Note that for
${\textit{Re}} = 400$
, the vortex shedding is suppressed in the vicinity of the cylinder. Instead, the characteristic frequency arises from the periodic meandering of the separated shear layers, which constitutes the weak vortex shedding observed further downstream at their tails. This definition of Strouhal number has been used to quantify the frequency of vortex shedding in both single-phase and two-phase flows (Gerrard Reference Gerrard1978; Norberg Reference Norberg1987; Williamson Reference Williamson1996a
; Williamson Reference Williamson1996b
; Norberg Reference Norberg2001, Reference Norberg2003; Reichl et al. Reference Reichl, Hourigan and Thompson2005; Jiang & Cheng Reference Jiang and Cheng2017). Figure 7 shows that the Strouhal number of vortex shedding is identical to those corresponding to the lift and drag coefficients in all cases. As the Reynolds number gradually increases, the Strouhal number increases monotonically, suggesting a higher vortex-shedding frequency. Similar to single-phase flow (Roshko Reference Roshko1954; Corke Reference Corke2024), the Strouhal number increases in a power law of
${\textit{Re}}^{-1}$
as
$St=0.52-72.7{\textit{Re}}^{-1}$
, although the fitted value of the constants in the two-phase flow varies from the single-phase flow.
The Strouhal (
$St$
) number based on the lift and drag coefficients and vortex shedding for two-phase flow with different Reynolds (
$\textit{Re}$
) numbers. The vortex-shedding frequency is identified using the spanwise-averaged spanwise vorticity
$\langle \omega _z \rangle _z$
at
$(x,y)=(5D, -1D)$
. The dashed line represents a fitted power law of
$St=0.52\mbox{--}72.7 {\textit{Re}}^{-1}$
.
Contours of instantaneous spanwise vorticity in an
$x$
–
$y$
plane for two-phase flows at (a)
${\textit{Re}} = 400$
(case T400), (b)
${\textit{Re}} = 500$
(case T500), (c)
${\textit{Re}} = 800$
(case T800), (d)
${\textit{Re}} = 1400$
(case T1400) and (e)
${\textit{Re}} = 2000$
(case T2000). The thick arrow shows the jet-like flow schematically. The black cross marks the probe at
$(x, y) = (5D, -1D)$
, where the frequency of the spanwise vorticity is calculated.
To further investigate flow status, we present the contours of the instantaneous spanwise vorticity (
$\omega _z$
) in an
$x$
–
$y$
plane in figure 8 for different cases. For all cases, two separated shear layers with positive and negative vorticity are attached to the lower and upper surface of the cylinder, respectively. Between the free surface and the top of the cylinder, the magnitude of vorticity is small inside the jet-like flow, shown schematically using a thick arrow. Above the jet-like flow, air entrainment induces a two-phase shear flow. In the downstream region of the cylinder, flow perturbations emerging from the cylinder and the free surface interact with each other, resulting in a more complex distribution of vorticity. The descriptions about the dominant flow structures are also reported in previous studies (Bouscasse et al. Reference Bouscasse, Colagrossi, Marrone and Souto-Iglesias2017; Zhao et al. Reference Zhao, Wang, Zhu, Ping, Bao, Zhou, Cao and Cui2021, Reference Zhao, Wang, Zhu, Cao, Bao, Zhou and Han2022; Patel et al. Reference Patel, Sun, Yang and Zhu2025).
Similar to a single-phase flow, the Reynolds number, as shown in figure 8, influences the shear layer attached to the cylinder. For the case with
${\textit{Re}} = 400$
, the two shear layers are quasisymmetric in the near-wake region about an oblique line aligned to the jet-like flow, denoted by the dashed line in figure 8(a). The symmetry of the shear layer extends for a considerable distance (
$\approx 3D$
). The minimum-pressure point on the upper surface of the cylinder remains almost stationary, leading to a small magnitude of the high-frequency signal in the spectra of lift and drag coefficients. As the jet-like flow becomes weak in a downstream region, the unsteady flow above the jet-like flow starts to influence the attached shear layer, and ultimately induces weak meandering of the shear layer. Large peaks observed in the time histories of
$C_l$
and
$C_d$
are attributed to the unsteady roll-up of detached shear layers. Specifically, the crest (event
$\text{E}_1$
) and trough (event
$\text{E}_2$
) identified in figure 4(b) are shown magnified in figures 9(a) and (b), respectively. Event
$\text{E}_1$
corresponds to the roll-up of separated shear layers, visualised via instantaneous spanwise vorticity contours in figure 9(c). During this event, the low-pressure region on the cylinder’s upper surface migrates towards its base, resulting in increased drag. Conversely, event
$\text{E}_2$
corresponds to the roll-up of the two-phase shear layer over the jet-like flow, as depicted in the instantaneous spanwise vorticity contours in figure 9(d). In this case, the low-pressure region moves away from the cylinder’s base, leading to a decrease in drag. This roll-up process is unsteady and exhibits intermittent behaviour. In another words, vortex shedding is rapidly suppressed when the separated shear layers and the two-phase shear layer over the jet-like flow do not exhibit synchronised vortex shedding.
(a,b) Enlarged view of the
$C_d$
time history close to events
$\text{E}_1$
and
$\text{E}_2$
marked in figure 4(b) for the T400 case. (c,d) Corresponding instantaneous spanwise vorticity (
$\omega_z / U_0D^{-1} \in [-2,2]$
) contours in an
$x$
-
$y$
plane.
For cases with
${\textit{Re}} \geqslant 500$
, figure 8(b–e) shows that the symmetry of the attached shear layers is broken in a closer region to the cylinder, resulting in vortex shedding. Meanwhile, in the two-phase region above the jet-like flow, the shear layer is shortened in comparison with
${\textit{Re}}=400$
. This shear layer is destabilised, and rolls up to form a small recirculation zone (see the flow region marked with letter A in figure 8
b). Some bubbles are trapped in this recirculation zone. Downstream of this recirculation zone, the scale of the spanwise vorticity decreases as the Reynolds number increases from 500 to 2000. This is a typical Reynolds-number effect.
In single-phase flow, the vortex shedding occurs at
${\textit{Re}}=47$
(Strykowski & Sreenivasan Reference Strykowski and Sreenivasan1990; Williamson Reference Williamson1996b
; Zdravkovich Reference Zdravkovich1998; Rao et al. Reference Rao, Radi, Leontini, Thompson, Sheridan and Hourigan2015). In contrast, the vortex shedding is weak at a much higher Reynolds number (
${\textit{Re}}=400$
) in the present two-phase flow. This indicates a stabilisation effect of the jet-like flow. As shown in figure 8, an additional shear layer occurs above the jet-like flow. This scenario is similar to the single-phase flow past a cylinder near a flat wall, where the shear layer attached to the flat wall also imposes a stabilisation effect. Specifically, when the gap between the bottom of cylinder and the wall is
$0.4D$
, the vortex shedding does not take place at
${\textit{Re}}=400$
. As the gap decreases to
$0.2D$
, the critical Reynolds number for vortex shedding further increases to
$\mathcal{O}(1000)$
(Lei et al. Reference Lei, Cheng, Armfield and Kavanagh2000).
Isosurfaces of
$Q / U_0^{2} D^{-2} = 10$
coloured by streamwise velocity (
$u / U_0 \in [-1,1]$
) for two-phase flows at (a)
${\textit{Re}} = 400$
(case T400), (b)
${\textit{Re}} = 500$
(case T500), (c)
${\textit{Re}} = 800$
(case T800), (d)
${\textit{Re}} = 1400$
(case T1400) and (e)
${\textit{Re}} = 2000$
(case T2000). Each case is shown in a perspective view (left column) and a down-top view (right column). Background contours in left and right columns represent spanwise vorticity (
$\omega_z / U_0D^{-1} \in [-2,2]$
) in the
$x$
-
$y$
and
$x$
-
$z$
planes, respectively.
The vortex structures can be identified using a Q-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988). Figure 10 shows the isosurfaces of
$Q / U_0^{2} D^{-2}=10$
, coloured by streamwise velocity. As shown in figure 10(a), no vortex structures are identified by the Q-criterion inside the dashed ellipse at
${\textit{Re}} = 400$
. The absence of vortex structures suggests a non-turbulent status. Inside the jet-like flow, no vortex structures are identified for all cases. This is similar to a planar jet, where vortex structures are developed on the two sides of the jet. Above the jet-like flow, the location where vortex structures are identified moves upstream as the Reynolds number increases, indicating the occurrence of instability. The flow topology exhibits a clear transition with Reynolds number. In the shear layer, two-dimensional planar vortices gradually emerge, while in the wake, three-dimensional vortices dominate and shift closer to the cylinder as
${\textit{Re}}$
increases, consistent with the upstream movement of the recirculation zone and enhanced mixing. Specifically, at
${\textit{Re}}=400$
, planar vortices in the shear layer are nearly absent, though three-dimensional vortices appear in the downstream wake. At
${\textit{Re}}=500$
(figure 10
b), large-scale vortices begin to form in the shear layer, and three-dimensional vortices extend further upstream compared with
${\textit{Re}}=400$
. At
${\textit{Re}}=800$
(figure 10
c), distinct planar vortices are clearly visible in the shear layer, while the wake contains a greater number of three-dimensional structures. At higher Reynolds numbers (
${\textit{Re}}=1400$
and
$2000$
, figure 10
d,e), elongated planar vortices span nearly the entire body length, indicating strong three-dimensionality, although localised interruptions suggest possible amplitude modulation. Overall, increasing Reynolds number enhances vortex generation, reduces vortex scale, and promotes a transition from non-turbulent to turbulent flow in both the near-surface and wake regions.
3.3. Air entrainment and bubbles
Instantaneous air–water interface and volume rendering of the dimensionless streamwise velocity
$u/U_{0}$
at
${\textit{Re}}=1400$
(case T1400). The cuboid is the subdomain where the VOF is monitored.
In this section, we discuss the formation and transport processes of entrained air bubbles. Before we provide quantitative analysis, we summarise the key processes through an instantaneous snapshot of the interface in figure 11. The contours of streamwise velocity are superimposed. The water surface bumps upward ahead of the cylinder in the streamwise direction. Subsequently, the surface breaks up, air entrainment (indicated by arrow 1) and the wake generated behind the cylinder leading to the formation of bubbles. Most bubbles rise to the surface due to buoyancy, swirling and breaking up in the upper recirculation zone (indicated by arrow 2), while some are delivered downstream by the flow (indicated by arrow 3).
It can be observed from a video (see the supplementary movies available at https://doi.org/10.1017/jfm.2025.10839) that both the air-entrainment rate and bump height oscillate quasiperiodically in time. To quantitatively describe this temporal variation, the integration of VOF
$\mathcal{I}_\psi$
in the region of interest and the instantaneous gap height
$h^*$
are calculated. The integration
$\mathcal{I}_\psi$
is defined as
\begin{equation} \mathcal{I}_\psi (t) = \int _{V_m} \left (1 - \psi \right ) \mathop {\rm d}\!V, \end{equation}
where
$V_m$
is the volume of the monitoring box as demarcated in figure 11. The instantaneous gap height is defined as
\begin{equation} h^*(t) = A_c / L_z, \end{equation}
where
$A_c$
represents the cross-sectional area between the top edge of the cylinder and the water surface.
The frequency spectra of the VOF
$\varPhi _{\mathcal{I}_\psi }$
and gap height
$\varPhi _{h^*}$
for (a)
${\textit{Re}} = 400$
, (b)
${\textit{Re}} = 500$
, (c)
${\textit{Re}} = 800$
, (d)
${\textit{Re}} = 1400$
and (e)
${\textit{Re}} = 2000$
. Panel (e) also compares results from two different grid resolutions:
$D/40$
(case T2000) and
$D/60$
(case T2000F).
The Strouhal number corresponding to the integration of VOF
$\mathcal{I}_\psi$
and vortex shedding for different Reynolds numbers. The vortex-shedding frequency is identified using the spanwise-averaged spanwise vorticity
$\langle \omega _z \rangle _z$
at
$(x,y)=(5D, -1D)$
.
Contours of instantaneous spanwise vorticity (
$\omega_z / U_0D^{-1} \in [-2,2]$
) in an
$x$
-
$y$
plane at different stages of one period of vortex shedding. The results for
${\textit{Re}}=500$
(case T500) are shown. The time instant shown in figure 8(b) is defined as
$t=0$
, and different panels show the contours at (a)
$t=0$
, (b)
$t=0.2T_s$
, (c)
$t=0.4T_s$
, (d)
$t=0.6T_s$
, (e)
$t=0.8T_s$
and ( f)
$t=1.0T_s$
, where
$T_s=2.7D/U_0$
corresponding to the vortex shedding Strouhal number
$f = 0.37U_0/D$
identified from figure 7.
The frequency spectra of the gap height
$\varPhi _{h^*}$
and integration of VOF
$\varPhi _{\mathcal{I}_\psi }$
for different cases are compared in figure 12. Figure 12(e) compares two grid resolutions,
$D/40$
(case T2000) and
$D/60$
(case T2000F), showing reasonably good agreement between them. Both high- and low-frequency signals can be identified from
$\varPhi _{\mathcal{I}_\psi }$
. The high-frequency peak dominates for
${\textit{Re}} \geqslant 500$
, while the low-frequency signals are larger in magnitude for
${\textit{Re}} = 400$
. This feature in the dominant frequency is similar to the lift and drag coefficients shown in figure 6. Figure 13 further compares the Strouhal number corresponding to the high-frequency peak of
$\mathcal{I}_\psi$
and spanwise-averaged spanwise vorticity
$\langle \omega _z \rangle _z$
. It is evident that they are identical to each other in all cases, indicating that the vortex shedding is responsible for the high-frequency oscillation in the air-entrainment rate.
To further visualise the effect of vortex shedding on the air entrainment, we plot the geometry of the air–water interface together with the contours of instantaneous spanwise vorticity at different stages of vortex shedding in figure 14. The results for
${\textit{Re}}=500$
are shown here as a representative, while results for higher Reynolds numbers remain similar. The results of a complete vortex shedding period
$T_s$
are shown, where
$T_s=2.7D/U_0$
is the period of vortex shedding corresponding to the Strouhal number
$St=0.37$
identified from figure 13. As shown in figure 14(a–d), after the rolling up of the shear layer above the jet-like flow, the shear layer grows again from
$t=0$
to
$t=0.6T_s$
. The shear layer attached to the upper surface of the cylinder also grows. Meanwhile, bubbles are trapped inside the high-vorticity region, and evolve downstream together with the spanwise vortex. From
$t=0.6T_s$
to
$t=0.8T_s$
, the shear layer attached to the lower surface of the cylinder starts to roll up, and breaks the upper attached shear layer, causing vortex shedding below the jet-like flow. The shear layer above the jet-like flow is broken synchronously, and meanwhile the spanwise vortex formed at
$t=0$
is diffused to a larger region. At
$t=1.0T_s$
, the two shear layers on the two sides of the jet-like flow roll up again and enter the next period of vortex shedding. It is seen from figure 14(e, f) that a rare bubble occurs in the flow region between the two spanwise vortices as pointed out by the black arrow. This leads to the intermittency in the air entrainment, and thus a high-frequency signal is detected from the spectra of the volume of entrained air in figure 12.
Different from other quantities examined, including the lift and drag coefficients, spanwise vorticity and integral volume of entrained air, the gap height
$\varPhi _{h^*}$
between the cylinder top and water surface does not show any high-frequency signal. As shown in the right column of figure 12, only low-frequency signals with
$f\!D/U_0 \lt 0.1$
are detected. The low-frequency components show a multipeak structure, indicating complex nonlinear interactions. Furthermore, the low-frequency signals detected in the temporal variation of
$\varPhi _{h^*}$
,
$\varPhi _{\mathcal{I}_\psi }$
,
$\varPhi _l$
and
$\varPhi _d$
show an overlap in the range of frequency, indicating that the vertical oscillation in water surface ahead of the cylinder is the dominant factor that induces the low-frequency variation of air-entrainment rate and hydrodynamic forces. The quasiperiodic vertical oscillation of the water surface ahead of the cylinder reflects variations in pressure and gravitational potential energy head. These variations, in turn, drive oscillations in the jet-like flow and the overlying two-phase shear layer. Consequently, the lift and drag coefficients on the cylinder, as well as the air-entrainment rate, are influenced by this fundamental low-frequency surface oscillation. We did not observe a strict periodic switching between the spilling-jet and von Kármán shedding regimes mentioned by Colagrossi et al. (Reference Colagrossi, Nikolov, Durante, Marrone and Souto-Iglesias2019); in particular, the von Kármán vortex shedding is not fully suppressed in a periodic manner. Instead, the low-frequency components appear to result from the coupling between the jet-like flow in the gap above the cylinder, the free-surface deformation and three-dimensional nonlinear interactions. The observed coupling likely signifies the presence of a metastable state, consistent with the analyses of Reichl et al. (Reference Reichl, Hourigan and Thompson2005), Colagrossi et al. (Reference Colagrossi, Nikolov, Durante, Marrone and Souto-Iglesias2019), and our prior work (Patel et al. Reference Patel, Sun, Yang and Zhu2025). The ratio between the low- and high-frequency components further supports this conclusion, with a quantitative estimate of
$\mathcal{O}(10^{-1})$
that mirrors the results of these studies. The spanwise vorticity distributions obtained in our simulations are also consistent with those reported in the mentioned studies.
The bubble statistics are useful to identify the status of air entrainment. In the present study, the bubble identification technique introduced by Herrmann (Reference Herrmann2010) is employed to delineate and characterise individual bubbles.
Projection of bubble coordinates onto an
$x$
-
$z$
plane for (a)
${\textit{Re}} = 400$
(case T400) and (b)
${\textit{Re}} = 500$
(case T500). High-density regions with tightly clustered bubble coordinates are marked in red, whereas those with low bubble density are shown in blue.
In our previous study (Patel et al. Reference Patel, Sun, Yang and Zhu2025), a spanwise periodic distribution of bubbles was observed at
${\textit{Re}}=150$
. The wavelength was approximately
$1.0D$
. Increasing the Reynolds numbers can induce turbulence and enhance bubble fragmentation (Garrett, Li & Farmer Reference Garrett, Li and Farmer2000; Deane & Stokes Reference Deane and Stokes2002). Therefore, it can be expected that the quasiperiodic distribution becomes less significant at higher Reynolds numbers. To examine this point, we conducted a time accumulative projection of bubble coordinates onto an
$x$
-
$z$
plane. The projection included a collection of bubble coordinates recorded over the time window of
$\Delta t D/U_0 = 100$
at a sampling rate of
$0.05D/U_0$
. Figure 15 compares the projection for
${\textit{Re}}=400$
and 500. As expected, the distribution of bubbles detected from figure 15 is almost uniform in the spanwise direction for both cases. As the Reynolds number further increases, the mixing becomes stronger, and the spanwise heterogeneity is further weakened. The images for higher Reynolds numbers remain similar to
${\textit{Re}}=400$
and
$500$
, and as such the figures are not shown.
We also calculated the size spectra of the bubble number density
${S}_b (R_{\textit{eff}} )$
as
\begin{equation} {S}_b\big(R_{\textit{eff}}\big) = \left \langle \frac {n(R_{\textit{eff}},t;b)}{V b} \right \rangle _t , \end{equation}
where
$n(R_{\textit{eff}},t;b)$
is the number of bubbles, whose effective radii fall between
$R_{\textit{eff}}$
and
$R_{\textit{eff}} + b$
in a given fluid volume
$V$
at time
$t$
(Li et al. Reference Li, Yang and Zhang2024). In the present study,
$b = 0.001D$
is selected. The fluid volume
$V$
is a cuboid in a subdomain with
$x/D \in [0.5, 20.5]$
,
$y/D \in [-5.0, 1.0]$
,
$z/D \in [0, 3]$
, which contains most air bubbles. The effective spherical radius is defined as
\begin{equation} R_{\textit{eff}} = \left (\frac {3 v_e}{4 \unicode{x03C0} }\right )^{1/3} , \end{equation}
where
$v_e$
represents the volume of an individual bubble. The size spectra of bubble number density are presented in figure 16. The images for Reynolds numbers of 800 and 1400 remain similar with
${\textit{Re}} = 400$
and 500, and as such the figures are not shown. For
${\textit{Re}} = 2000$
, figure 16(c) compares two grid resolutions,
$D/40$
(case T2000) and
$D/60$
(case T2000F), showing reasonably good agreement between them. An
$R_{\textit{eff}}^{-6}$
scaling law is identified. In a previous study of wave breaking, Deane & Stokes (Reference Deane and Stokes2002) identified a different scaling law from the size spectra of bubble number density. During the active phase of wave breaking, the turbulence generated by the plunging of an overturning jet breaks the entrained air pocket into small bubbles. At this stage, the turbulence is active, and an
$R_{\textit{eff}}^{-10/3}$
scaling law was identified. After the wave plunging, the breaking enters a quiescent phase and without energy being injected, turbulence intensity is weakened due to the spatial diffusion and dissipation. At this stage, bubbles are no longer broken into smaller ones. Large bubbles float up and leave the water, while small bubbles stay in the water because of their smaller terminal velocity and the trapping effect of vortex structures. The negative power of the scaling grows to
$-6$
. The growth of the scaling power was also observed in other flows, including plunging jet (Li, Deng & Yang Reference Li, Deng and Yang2025) and wave breaking in the wake of a stern (Hendrickson et al. Reference Hendrickson, Weymouth, Yu and Yue2019a
). In the present simulation cases, without plunging processes, the system lacks the intense turbulent energy input that drives bubble fragmentation in breaking waves. Therefore, the conditions governing air entrainment and bubble formation are significantly different, likely being more sensitive to the interplay of gravitational potential energy (e.g. in determining interface stability or the energy required for deformation) and surface tension, rather than being dominated by strong turbulent stresses. Therefore, the
$R_{\textit{eff}}^{-10/3}$
scaling detected in the active phase of wave breaking is not identified in the present study. The interaction between turbulent vortices and bubbles is similar to the quiescent phase of wave breaking, and therefore the
$R_{\textit{eff}}^{-6}$
scaling is observed.
Size spectra of bubble number density for (a)
${\textit{Re}} = 400$
(case T400), (b)
${\textit{Re}} = 500$
(case T500) and (c)
${\textit{Re}} = 2000$
. The dashed lines correspond to an
$R_{\textit{eff}}^{-6}$
scaling. Panel (c) also compares results from two different grid resolutions:
$D/40$
(case T2000) and
$D/60$
(case T2000F).
3.4. Statistics of unsteady flow
Two types of unsteady flows are generated in two-phase flow past a submerged cylinder, namely the wake flow downstream of the cylinder and mixed-phase flow near the water surface. To further investigate the characteristics of these two types of unsteady flows, we examine the statistics, including the mean velocity and kinetic energy of velocity fluctuations.
Contours of mean streamwise velocity
$\langle u \rangle _{zt}$
for (a)
${\textit{Re}} = 400$
(case T400) and (b)
${\textit{Re}} = 500$
(case T500). The green dashed line represents the time-averaged interface corresponding to the isopleth of
$\langle \psi \rangle _{zt} = 0.5$
. The red dash-dotted line demarcates the recirculation zone behind the cylinder via the isopleth of
$\langle u \rangle _{zt} = 0$
. The purple arrow denotes the measurement of recirculation zone length.
Vertical profiles of mean streamwise velocity
$\langle u \rangle _{zt}$
at different streamwise locations with (a)
${x/D = 2}$
and (b)
$x/D = 5$
. Results from two grid resolutions are compared:
$D/40$
(case T2000) and
$D/60$
(case T2000F).
Figure 17 compares the contours of the mean streamwise velocity
$\langle u \rangle _{zt}$
for
${\textit{Re}}=400$
and 500, where a pair of angular brackets with subscript
$zt$
represents time- and spanwise-averaging. The green dashed line represents the mean water surface, identified using the mean value of the VOF function as
$\langle \psi \rangle _{zt} = 0.5$
. The results for
${\textit{Re}}\geqslant 800$
are not shown here because the flows for
${\textit{Re}}\geqslant 500$
are in a similar status. To provide comparison of the streamwise velocity for all cases, we plot the vertical profile of
$\langle u \rangle _{zt}$
at
$x=2D$
and
$5D$
in figure 18. Figure 18 demonstrates excellent agreement between the results obtained with the
$D/40$
(case T2000) and
$D/60$
(case T2000F) grids, indicating that the computed streamwise velocity statistics are well converged for these resolutions.
For
${\textit{Re}}\geqslant 500$
, the Reynolds number only influences the flow region near the cylinder. As shown in figure 17, the acceleration on the side of the cylinder and deficit in the wake region is similar to single-phase flow. Behind the cylinder, there is a recirculation zone. Following Massey, Langella & Swaminathan (Reference Massey, Langella and Swaminathan2019), we identify the recirculation zone as the flow region with
$\langle u \rangle _{zt} \leqslant 0$
. The edge of the recirculation zone is demarcated using the red dash-dotted line in figure 17. Outside the recirculation zone, the velocity recovers to the inflow magnitude gradually as the water flows downstream. As the Reynolds number increases from 500 to 1400, it is seen from figure 18(a) that the streamwise velocity is positive for
${\textit{Re}}\geqslant 500$
around
$(x, y) = (2D, -0.7D)$
, with its magnitude increases monotonically . This observation indicates that the point
$(x, y) = (2D, -0.7D)$
is located outside the recirculation zone for
${\textit{Re}}\geqslant 500$
, and increasing the Reynolds number from 500 to 1400 leads to the reduction of the length of the recirculation zone. The length of the recirculation zone,
$L_R$
, can be quantified as the distance from the cylinder surface to the farthest point of the recirculation zone, which is shown schematically in figure 17(a). Figure 19 compares the value of
$L_R / D$
for all cases. The decrease in the length of the recirculation zone is confirmed from the figure.
In the flow region close to the water surface with
$y/D\geqslant 0$
or far away from the cylinder, the profiles of
$\langle u \rangle _{zt}$
shown in figure 18 collapse for
${\textit{Re}}\geqslant 500$
. As is reported in other flows, including plunging jet and wave breaking, the effect of Reynolds number is less significant than the Froude number in turbulence flows close to a free surface (Duncan et al. Reference Duncan, Qiao, Philomin and Wenz1999; Caulliez Reference Caulliez2013; Deike, Popinet & Melville Reference Deike, Popinet and Melville2015; Li et al. Reference Li, Yang and Zhang2024). This point is supported by the results of the present cases with
${\textit{Re}}\geqslant 500$
. The results for
${\textit{Re}}=1400$
and 2000 are close to each other in terms of both the profile of
$\langle u \rangle _{zt}$
and the length of recirculation zone, shown in figures 18 and 19, respectively. This observation indicates that there exists a threshold value of the Reynolds number, beyond which the flow dynamics remains similar.
The flow status for
${\textit{Re}}=400$
is drastically different from other cases in both the wake and near-surface regions. As shown in figure 8, the stabilisation effect of the jet-like flow leads to a significant suppression of vortex shedding at
${\textit{Re}}=400$
in comparison with other cases. This results in a larger recirculation zone in the wake region at
${\textit{Re}}=400$
. As is evident from figure 19, the length of the recirculation zone decreases drastically as the Reynolds number increases from 400 to 500. Thus, negative streamwise velocity occurs in a larger region at
${\textit{Re}}=400$
. As shown in figure 18(a), around
$(x, y) = (2D, -0.7D)$
, the mean streamwise velocity
$\langle u \rangle _{zt}$
is negative at
${\textit{Re}}=400$
, while it is positive at higher Reynolds numbers. Near the water surface, the maximum negative velocity occurs around
$x=5D$
for
${\textit{Re}}=400$
. In contrast, it occurs at an upstream location
$x=2D$
for
${\textit{Re}}=500$
. In other words, the low-Reynolds-number effect leads to the elongation of the recirculation zone near the water surface. This further results in a smaller and larger magnitude of the negative mean streamwise velocity at the upstream location
$x=2D$
and downstream location
$x=5D$
for the case with
${\textit{Re}}=400$
, in comparison with other cases with higher Reynolds numbers (figure 18).
The recirculation zone length
$L_R$
with different Reynolds numbers.
Contours of time- and spanwise-averaged vertical velocity (
$\langle v \rangle _{zt}$
) at different Reynolds numbers: (a)
${\textit{Re}} = 400$
(case T400) and (b)
${\textit{Re}} = 500$
(case T500). The green dashed line represents the time-averaged interface corresponding to the isopleth of
$\langle \psi \rangle _{zt} = 0.5$
.
Figure 20 compares the contours of the mean vertical velocity
$\langle v \rangle _{zt}$
for
${\textit{Re}}=400$
and 500. The two-phase flow is similar to a single-phase flow in the upstream region of the cylinder, where the vertical velocity is positive and negative, respectively, near the upper and lower surfaces of the cylinder. In the downstream region of the cylinder, the jet-like flow breaks the symmetry of the flow about the horizontal central plane of the cylinder, and leads to a downward shift of the direction of the main stream. This further results in a negative mean vertical velocity below the cylinder, which is different from the positive vertical velocity in this region of a single-phase flow.
Vertical profiles of mean vertical velocity
$\langle v \rangle _{zt}$
at different streamwise locations with (a)
${x/D = 1}$
and (b)
$x/D = 3$
. Results from two grid resolutions are compared:
$D/40$
(case T2000) and
$D/60$
(case T2000F).
To quantify the effect of Reynolds number, figure 21 compares the vertical profiles of
$\langle v \rangle _{zt}$
at
$x=1D$
and
$3D$
for different cases. The excellent agreement between the
$D/40$
(case T2000) and
$D/60$
(case T2000F) grids is also observed for the vertical velocity statistics, confirming that the computed vertical velocity statistics are well converged for these resolutions. Similar to the mean streamwise velocity, the Reynolds number mainly influences the vertical velocity in the vicinity of the cylinder. It is seen from figure 21(a) that, as the Reynolds number increases from
$400$
to
$1400$
, the magnitude of the negative vertical velocity above
$y=0$
increases. The jet-like flow is formed due to the lift of surface elevation ahead of the cylinder. As the water flows from a higher elevation to a lower elevation, it is accelerated in the vertical direction by gravity. From figure 20, it is seen that the water elevation at
$x=-1D$
remains almost the same for
${\textit{Re}}=400$
and 500, indicating that the vertical momentum gained from the gravitational acceleration remains close in these two cases. As the Reynolds number increases, the momentum loss due to the viscous effect decreases, and, as a result, the negative velocity magnitude inside the jet-like flow increases. At a downstream location of
$x = 3D$
, the Reynolds number shows a reversed effect on the negative vertical velocity of the jet-like flow. As shown in figure 21(b), the negative velocity magnitude decreases with increasing Reynolds number. This can be attributed to the occurrence of turbulence at higher Reynolds numbers, which enhances the momentum exchange between the jet-like flow and surrounding fluids. It is also seen from figure 21 that the profiles of
$\langle v \rangle _{zt}$
for
${\textit{Re}}=1400$
and 2000 are close to each other, an observation that is consistent with
$\langle u \rangle _{zt}$
shown in figure 18. This observation confirms that the Reynolds number mainly influences the mean flow when it is sufficiently low.
Contours of kinetic energy of velocity fluctuations (
$k$
) at different Reynolds numbers: (a)
${\textit{Re}} = 400$
(case T400), (b)
${\textit{Re}} = 500$
(case T500), (c)
${\textit{Re}} = 800$
(case T800), (d)
${\textit{Re}} = 1400$
(case T1400) and (e)
${\textit{Re}} = 2000$
(case T2000). The green dashed line represents the time-averaged interface corresponding to the isopleth of
$\langle \psi \rangle _{zt} = 0.5$
. The blue crosses mark the locations of peak
$k$
values.
To further examine the intensity of the unsteady flow, we calculated the kinetic energy of velocity fluctuations
$k$
, defined as:
\begin{equation}k = \frac{1}{2}\left(\langle{\rho u^\prime u^\prime}\rangle_{zt} + \langle{\rho v^\prime v^\prime}\rangle_{zt} + \langle{\rho w^\prime w^\prime}\rangle_{zt} \right),\end{equation}
where the single prime denotes fluctuations with respect to time- and spanwise-averaging.
Figure 22 compares the contours of
$k$
for two-phase flows at different Reynolds numbers. The cross symbols in the figure represent local peaks of
$k$
. By contrasting contours of
$k$
with those of mean velocity shown in figures 17 and 20, we observed that the unsteadiness is weak inside the jet-like flow. At
${\textit{Re}}=400$
, the dominant peak of
$k$
occurs at a downstream location of the two-phase region above the jet-like flow. This observation confirms that the shear layer above the jet-like flow is stable. As the Reynolds number increases to
$500$
, this peak occurs at an upstream location above the jet-like flow. This is consistent with the rolling up of shear layer observed from figure 8(b). The location of this peak in the two-phase region remains almost unchanged as the Reynolds number increases from 500 to 2000. This observation suggests again that the effect of Reynolds number is insignificant in a flow close to a free surface.
In the far-wake region of the cylinder, a weak peak of
$k$
is observed at
${\textit{Re}}=400$
. This peak is caused by the weak meandering of the shear layer (see figure 8
a and corresponding discussions). As the Reynolds number increases to
$500$
, the peak of
$k$
in the wake region also moves upstream because of the occurrence of vortex shedding there. Two peaks are observed for
${\textit{Re}}=500$
and 800, located at the edge of the two shear layers. As the Reynolds number further increases to
$1400$
, the two peaks merge into one peak, indicating a stronger mixing effect induced by the vortex shedding.
4. Conclusion
We conducted three-dimensional DNS of a two-fluid–structure interaction problem involving air–water free-surface flow around a submerged circular cylinder. We employed a coupled LS and VOF formulation to capture the free-surface dynamics and a sharp-interface immersed boundary method to resolve the flow around a stationary cylinder in our set-up. The present study was motivated by the limited understanding of flow and interface dynamics in the intermediate Reynolds (
${\textit{Re}}$
) number range, where turbulence starts to develop. To bridge this gap, we probed into a relatively broad range of
${\textit{Re}}$
, varying from
$400$
to
$2000$
. The other dimensionless parameters that dictate the competition among inertia, surface tension and buoyancy were held constant by adopting the Weber and Froude number values of
${\textit{We}} = 1000$
and
${\textit{Fr}} = 1$
, respectively. The initial spacing between the cylinder’s top and the free surface was set equal to one cylinder radius in all cases.
At
${\textit{Re}} = 400$
, the overall organisation of flow structures qualitatively resembles that of low-Reynolds-number set-ups with
${\textit{Re}} = 150$
and
$180$
in prior studies (Reichl et al. Reference Reichl, Hourigan and Thompson2005; Patel et al. Reference Patel, Sun, Yang and Zhu2025). The vortex shedding is suppressed due to a jet-like flow emerging from the gap above the cylinder. This jet-like flow also tilts shear layers originating from the cylinder surface in the downward direction. Even at a relatively high
${\textit{Re}}$
of
$400$
, these attached shear layers remain stable and extend roughly two cylinder diameters downstream. Upon slightly increasing the Reynolds number to
${\textit{Re}} = 500$
, we observe substantial changes in the vorticity dynamics. First, the shear layer originating from the free surface, due to flow separation, becomes shorter and loses its stability. Second, the shear layers at the cylinder surface also contract in length and start shedding vortices. The shedding of vortices from unstable shear layers persists at subsequent Reynolds numbers in our range up to
${\textit{Re}}=2000$
.
The free-surface deformation and the flow pattern around the cylinder have direct consequences on the time variation of the hydrodynamic forces acting on the cylinder. Unlike single-phase flow, the lift and drag forces remain nearly stable over time due to the absence of vortex shedding at
${\textit{Re}} = 400$
. Subsequently, the fluctuations in the lift- and drag-force signals are restored for
${\textit{Re}} \geqslant 500$
. Moreover, due to the broken symmetry of the stagnation point, both lift and drag forces have identical dominant frequencies, which increase with
${\textit{Re}}$
following the relation
$0.52\mbox{--}72.7{\textit{Re}}^{-1}$
. The primary frequency associated with the spanwise vorticity component also follows the same variation. We discovered that the amplitude–frequency spectrum of the lift force contains weak low-frequency components for
${\textit{Re}} \geqslant 500$
, which stem from the coupling between the jet-like flow in the gap above the cylinder and the free-surface deformation, as previously observed in so-called metastable states (Colagrossi et al. Reference Colagrossi, Nikolov, Durante, Marrone and Souto-Iglesias2019; Patel et al. Reference Patel, Sun, Yang and Zhu2025).
From the direct visualisation of instantaneous vortex structures using the Q-criterion, we found that turbulence in the wake amplifies with
${\textit{Re}}$
and gradually propagates upstream towards the cylinder. This was further substantiated by inspecting the distribution of the time-averaged kinetic energy carried by the velocity fluctuations in the wake region. At
${\textit{Re}} = 400$
, the hot spots corresponding to the elevated disturbance kinetic energy are situated away from the free surface and the cylinder, which subsequently move closer to the cylinder and free surface due to the onset of vortex shedding and turbulence for
${\textit{Re}} \geqslant 500$
. The emergence of turbulent flow structures also affects the breakup dynamics of the air–water interface. In our earlier study at
${\textit{Re}} = 150$
(Patel et al. Reference Patel, Sun, Yang and Zhu2025), we detected the entrainment of gas rolls through interface breakup, which then fragment into bubbles in the downstream area. For the current range of Reynolds numbers, no roll-like structures are detected. Instead, bubbles are formed directly from the interface breakup. However, we observe entrainment of streak-like structures during interface breakup due to a stable jet-like flow in the gap at
${\textit{Re}} = 400$
.
We observed that the dominant frequency of fluctuations in the entrained volume near the free surface and the cylinder matches the vortex-shedding frequency. Once bubbles are entrained, they undergo subsequent fragmentation in the cylinder wake. We found that turbulence in the wake tends to smear out the inhomogeneity in the bubble concentration. This is in contrast to the spanwise periodic pattern of bubble concentration observed in the laminar regime at
${\textit{Re}} =150$
(Patel et al. Reference Patel, Sun, Yang and Zhu2025). Lastly, the time-averaged bubble-size density exhibits a
$-6$
power-law scaling, as also observed during gas entrainment in various settings (Deane & Stokes Reference Deane and Stokes2002; Hendrickson et al. Reference Hendrickson, Weymouth, Yu and Yue2019a
; Li et al. Reference Li, Deng and Yang2025).
The present work serves as the first step towards understanding the free-surface flow past a submerged cylinder within the transitional Reynolds-number range. This study will be extended in the future to understand the role of buoyancy (characterised by Froude number) and surface tension (characterised by Weber number), which will provide further insights into the coupled nature of the interface dynamics and turbulence in the wake of the cylinder.
Comparison of contours of mean streamwise velocity (
$\langle u \rangle _{zt}$
) obtained using different grid resolutions: (a)
$D/40$
(case T2000) and (b)
$D/60$
(case T2000F).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10839.
Acknowledgements
J.S. and Z.Y. acknowledge Professor Jianjun Tao from Peking University for the fruitful discussion and helpful suggestions on the results analyses.
Funding
This research is supported by the National Natural Science Foundation of China (NSFC) projects (no. 11972038, 12302291, 12402260) and the Chinese Academy of Sciences Project for Young Scientists in Basic Research (YSBR-087). K.P. and X.Z. acknowledge financial support from the Max Planck Society. X.Z. also acknowledges the support from the German Research Foundation (DFG) through projects 521319293, 540422505 and 550262949. Open access funding provided by Max Planck Society.
Declaration of interests
The authors report no conflict of interest.
Appendix. Additional results from grid-dependence testing
This appendix contains supplementary figures 23–28, which show additional contour plots and statistical data that aid the discussion of the grid-independence study in the main text.
Comparison of contours of mean vertical velocity (
$\langle v \rangle _{zt}$
) obtained using different grid resolutions: (a)
$D/40$
(case T2000) and (b)
$D/60$
(case T2000F).
Comparison of contours of kinetic energy of velocity fluctuations (
$k$
) obtained using different grid resolutions: (a)
$D/40$
(case T2000) and (b)
$D/60$
(case T2000F).
Comparison of vertical profiles of kinetic energy of velocity fluctuations (
$k$
) at
$x/D = 1.5$
, obtained using different grid resolutions.
Contours of instantaneous spanwise vorticity (
$\omega_z / U_0D^{-1} \in [-2,2]$
) in an
$x$
-
$y$
plane for case T2000F.
Isosurfaces of
$Q/U_0^2 D^{-2}=10$
coloured by streamwise velocity (
$u / U_0 \in [-1,1]$
) shown in (a) perspective view and (b) down-top view for case T2000F.
Open access
