Free-surface flow past a submerged horizontal cylinder: Reynolds-number dependence

Abstract

Direct numerical simulations of two-phase, free-surface flow past a fully submerged, fixed circular cylinder are conducted for transitional Reynolds numbers $400 \leqslant {\textit{Re}} \leqslant 2000$, with Weber number ${\textit{We}} = 1000$, Froude number ${\textit{Fr}} = 1$ and a fixed gap ratio $G = 0.5$. This parameter combination corresponds to the gas entrainment regime characterised by the production of multiscale gas bubbles through interface breakup in the wake, which is of particular interest for its implications in enhancing gas transfer and mixing in environmental and engineering flows, such as air–water gas exchange processes in rivers and oceans, and the design and performance of naval and offshore structures. For ${\textit{Re}}= 400$, the jet forced through the $0.5D$ gap where $D$ is the diameter of the cylinder, efficiently convects opposite-signed vorticity downstream, suppressing the classical von Kármán instability and yielding a quasisteady recirculation bubble. The jet’s stabilising influence, however, breaks down once ${\textit{Re}} \approx 500$: periodic vortex shedding re-emerges and the wake becomes unsteady in spite of the continuing jet. The corresponding dimensionless shedding frequency Strouhal number $St$ grows with ${\textit{Re}}$ as $0.52-72.7{\textit{Re}}^{-1}$. The onset of unsteadiness first shortens the mean separation length but then drives it towards a saturation plateau for higher ${\textit{Re}}$ values. Surface rupture in the turbulent wake fragments entrained air into a multiscale bubble population whose number density follows $S_b(R_{\textit{eff}}) \propto R_{\textit{eff}}^{-6}$, consistent with gravity–capillary breakup in breaking waves, where $R_{\textit{eff}}$ represents the effective radii of the bubbles. Intermittency in entrainment corresponding to vortex shedding contrasts sharply with the finger-like structures observed under laminar conditions, underscoring the role of turbulent mixing. The coupled analysis of vorticity transport, shear-layer instability and bubble statistics elucidates how momentum exchange and air entrainment over a submerged body are governed under non-turbulent and turbulent conditions.

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 1878; von Kármán & Rubach 1912; Berger & Wille 1972; Williamson & Govardhan 2004). 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 1954; Bloor 1964; Grove, Shair & Petersen 1964; Gerrard 1966; Dennis & Chang 1970; Fornberg 1985; Dušek et al. 1994; Norberg 1994; Mittal & Balachandar 1995; Zhang et al. 1995; Barkley & Henderson 1996; Henderson 1996; Thompson, Hourigan & Sheridan 1996; Williamson 1996b ; Zdravkovich 1998; Rajani, Kandasamy & Majumdar 2009; Chopra & Mittal 2019). Williamson (1996b ) and Zdravkovich (1998) 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. 1995, 1997) and numerical (Reichl, Hourigan & Thompson 2005; Colagrossi et al. 2019; Zhao et al. 2021; Patel et al. 2025) 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 1990; Lugt & Ohring 1992; Oguz, Prosperetti & Kolaini 1995; Ye & Chu 1997; Kiger & Duncan 2012; Li, Yang & Zhang 2024). 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 2012; Veron 2015; Drenckhan & Saint-Jalmes 2015).

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. (2005) 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. (2019), 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. (2019) conducted the global stability analysis to assess how free-surface deformation influences the onset of vortex shedding. In our recent work (Patel et al. 2025), 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 2020), a rotating circular cylinder (Guo et al. 2023; Lin & Yao 2023), a pair of stationary tandem cylinders (Zhao et al. 2024), and fixed elliptic (Subburaj, Khandelwal & Vengadesan 2018) and square (Karmakar & Saha 2020) 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. (1995, 1997) 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 (1990) 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. (2021) performed large-eddy simulations with the same parameters as the experiments of Sheridan, Lin & Rockwell (1997), 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. (2022) 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 (2023) 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. 2018; Chizfahm, Joshi & Jaiman 2021; Hunt et al. 2023) and flat plates (Díaz-Ojeda et al. 2019; Hu et al. 2023).

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 2012; Veron 2015; Drenckhan & Saint-Jalmes 2015; Hendrickson et al. 2019b ; Li et al. 2024). 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 2021). 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:

(2.1) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0 , \end{align}

(2.2) \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

(2.3) \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

(2.4) \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

(2.5) \begin{align} \rho &= \rho _a + (\rho _{w} - \rho _a)H(\phi ) , \end{align}

(2.6) \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 2000). The following convection equations for the level-set (LS) $\phi$ and volume-of-fluid (VOF) $\psi$ functions are solved:

(2.7) \begin{align} \frac {\partial \phi }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{u})&=0, \end{align}

(2.8) \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. (2021).

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.

Table 1. 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$ .

Figure 1. 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.

Figure 2. 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)$ .

Table 2. 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 2010; Shrivastava, Agrawal & Sharma 2013; Cui et al. 2017). The solver’s capability for accurately simulating general free-surface flows has been extensively validated (Yang et al. 2021), with successful applications to diverse scenarios such as plunging jets (Li et al. 2024), wave breaking (Lu et al. 2024) and bubbly channel flow (Lu, Yang & Deng 2025). Crucially for the present study, its applicability to similar configurations was established in our prior work (Patel et al. 2025). 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. 1995, 1997; Bouscasse et al. 2017; Zhao et al. 2021; Hendrickson, Yu & Yue 2022; Patel et al. 2025).

Figure 3. 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

(3.1) \begin{align} C_l &= \frac {F_l}{0.5 \rho _w U_0^2 A}, \end{align}

(3.2) \begin{align} C_d &= \frac {F_d}{0.5 \rho _w U_0^2 A}, \end{align}

Figure 4. 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

(3.3) \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

(3.4) \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. (1990), Zhao et al. (2021), Patel et al. (2025), 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:

(3.5) \begin{align} C_l &= \langle C_l \rangle _{t} + C_l^{\prime \prime }, \end{align}

(3.6) \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 1922; Zdravkovich 1998; Norberg 2003; Jiang & Cheng 2018). 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.

Figure 5. 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$ .

Figure 6. 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

(3.7) \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 1978; Norberg 1987; Williamson 1996a ; Williamson 1996b ; Norberg 2001, 2003; Reichl et al. 2005; Jiang & Cheng 2017). 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 1954; Corke 2024), 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.

Figure 7. 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}$ .

Figure 8. 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. 2017; Zhao et al. 2021, 2022; Patel et al. 2025).

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.

Figure 9. (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 1990; Williamson 1996b ; Zdravkovich 1998; Rao et al. 2015). 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. 2000).

Figure 10. 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 1988). 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

Figure 11. 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

(3.8) \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

(3.9) \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.

Figure 12. 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).

Figure 13. 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)$ .

Figure 14. 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. (2019); 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. (2005), Colagrossi et al. (2019), and our prior work (Patel et al. 2025). 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 (2010) is employed to delineate and characterise individual bubbles.

Figure 15. 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. 2025), 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 2000; Deane & Stokes 2002). 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

(3.10) \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. 2024). 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

(3.11) \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 (2002) 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 2025) and wave breaking in the wake of a stern (Hendrickson et al. 2019a ). 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.

Figure 16. 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.

Figure 17. 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.

Figure 18. 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 (2019), 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. 1999; Caulliez 2013; Deike, Popinet & Melville 2015; Li et al. 2024). 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).

Figure 19. The recirculation zone length $L_R$ with different Reynolds numbers.

Figure 20. 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.

Figure 21. 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.

Figure 22. 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:

(3.12) \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. 2005; Patel et al. 2025). 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. 2019; Patel et al. 2025).

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. 2025), 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. 2025). Lastly, the time-averaged bubble-size density exhibits a $-6$ power-law scaling, as also observed during gas entrainment in various settings (Deane & Stokes 2002; Hendrickson et al. 2019a ; Li et al. 2025).

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.

Figure 23. 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).