# Flow topology changes with bubbly flow around a circular cylinder

Eric Thacher, Bruno van Ruymbeke, Céline Gabillet, Nicolas Jacques, Simo A Mäkiharju

► To cite this version:

Eric Thacher, Bruno van Ruymbeke, Céline Gabillet, Nicolas Jacques, Simo A Mäkiharju. Flow topology changes with bubbly flow around a circular cylinder. International Journal of Multiphase Flow, 2024, 179, 10.1016/j.ijmultiphaseflow.2024.104917. hal-04674245

HAL Id: hal-04674245

https://hal.univ-brest.fr/hal-04674245v1

Submitted on 21 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Flow topology changes with bubbly flow around a circular cylinder

Eric Thacher<sup>a,*</sup>, Bruno Van Ruymbeke<sup>b</sup>, Céline Gabillet<sup>b</sup>, Nicolas Jacques<sup>c</sup>, Simo A. Mäkiharju<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Vortex shedding

Bubbly flow

SPOD

POD

Circular cylinder

Cross-flow

## ABSTRACT

Vortex induced vibration (VIV) experienced during flow past a cylinder can reduce equipment performance and in some cases lead to failure. Previous studies have shown that the injection of bubbles in the flow over a cylinder typically leads to a monotonic increase in shedding frequency with void fraction, however, a satisfactory explanation for this phenomenon has not been proposed. Unexplained scatter in the data exists, including that the increase in shedding frequency is not universal. More research is needed to characterize the influence of bubbles on the wake structure, and subsequent shift in shedding frequency. To this aim, the effect of bubbles on the structure of the wake and VIV was examined over two values of Reynolds number, $Re_D = 100,000$ and $160,000$. Time-resolved particle image velocimetry (TR-PIV), proper orthogonal decomposition (POD) and spectral proper orthogonal decomposition (SPOD) of the wake structures, vibration of the cylinder, and bubble image velocimetry (BIV) were used to assess the flow topology changes under the influence of gas injection. Using SPOD/POD analysis in the near wake, it was found that the primary Karman shedding frequency decreased with the injection of gas, from a Strouhal number of $St = 0.2$ to $St = 0.17-0.18$; the width of the spectral peak was found to increase with void fraction. Notably, the vibration of the cylinder at the primary Karman shedding frequency was suppressed following the injection of gas, even at spanwise-averaged volumetric qualities below 0.01%. This suppression occurred regardless of if gas was concentrated locally near the centerline of the channel, or along the span. BIV data suggests that gas accumulation in the near wake, driven by the high velocity vertical motion of gas, serves to uncouple the cylinder motion from the formation of the vortex street downstream while promoting faster wake recovery.

## 1. Introduction

Multiphase flows are encountered in many applications in chemical engineering, heat transfer, as well as ocean and naval engineering. These include devices such as heat exchangers (Kaneko et al., 2014), vortex flow meters (Sun and Zhang, 2009), surface aerated stirring vessels (Patwardhan and Joshi, 1998), froth flotation (Yianatos, 2007), and control surfaces operating at or near the free surface. In each of these applications, both the structural integrity and performance of the device is dependent on the interactions between the component and the surrounding multiphase flow. However, the behavior of the multi-phase flow is challenging to predict *a priori* due to two-way coupling between the dispersed and carrier phase. Numerous flow modifications are possible, including the addition of energy at a variety of scales present in the turbulent flow (Mathai et al., 2020). Nevertheless, the two-way coupling of the two phases is of significant importance for predicting heat transfer and flow mixing (Mathai et al., 2020). This is particularly true in turbulent dispersed flows, which necessitates a

thorough understanding of this complex flow. In the last 30 years, as was summarized by Murai (2014) and Ceccio (2010), there has been significant study on the bubbly drag reduction process in turbulent pipe flows or turbulent boundary layers. However, much less effort has been devoted to the study of bubbly flow coupling with the hydrodynamic loading of immersed bodies.

When considering external flows, one of the main hydrodynamic mechanisms responsible for inducing vibrations is the trailing edge instability that leads to vortex shedding in the wake (Williamson, 1996; Norberg, 2003; Lienhard, 1966; Achenbach and Heinecke, 1981; Roshko, 1955; Bearman, 1984; Gabbai and Benaroya, 2005; Bearman, 2011). This is particularly true for bluff bodies, and hydrofoils truncated at their trailing edge (Ausoni, 2009). The complex interaction between the dispersed bubbly phase and the vortex induced vibration (VIV) of the structure is of importance for the design of structures immersed near the free surface (where natural aeration occurs) or immersed structures subjected to forced aeration of the flow. In part due

<footnote><sup>a</sup>Department of Mechanical Engineering, University of California, Berkeley, Berkeley 94720, USA</footnote>

<footnote><sup>b</sup>Institut de Recherche de l'Ecole Navale (IRENav), EA 3634 - Ecole Navale, 29240, Brest, France</footnote>

<footnote><sup>c</sup>ENSTA Bretagne, CNRS UMR 6027, IRDL, 2 rue François Verny, Brest Cedex 9, 29806, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: ethacher@berkeley.edu (E. Thacher).</footnote>
to the clear design & structural implications of multiphase coupling, there has been a significant body of research on periodic shedding in two-phase cross flow. In particular, due to the relevance of two-phase shedding in heat exchanger applications, there has been a significant body of research on flow over tube bundles (Pettigrew et al., 1991) or flow through a narrow gap (Callison et al., 2023). For brevity, we will not review this work here. Instead, we will focus solely on vortex shedding from cross-flow over a single cylinder. By focusing on this canonical flow, we can draw upon a wealth of past research on the single-phase case as well as find prior work with multiphase flow.

For single-phase flow past a cylinder over a wide range of Reynolds numbers we see the famed von Kármán vortex street. The structure of the vortex street varies with Reynolds number, which is defined as follows:

$$Re_D = \frac{UD}{\nu} \qquad (1)$$

Where U is the freestream velocity, D is the cylinder diameter, and ν is the kinematic viscosity. The structure of the vortex street with Reynolds number has been detailed elsewhere (Williamson, 1996; Norberg, 2003), but for clarity a brief review is presented here. For $Re_D \approx 1000-200,000$, the flow is considered subcritical and remains largely within the shear-layer transition regime (Williamson, 1996), with the only distinction being a transition from ‘high-quality’ to ‘low-quality’ vortex shedding occurring around $Re_D \approx 5000-8000$ (Norberg, 2003). Norberg (Norberg, 2003) proposed that this transition was caused when the shear layer transition to turbulence begins to precede the wake transition to turbulence, causing the wake transition mechanism to be altered. As Reynolds number increases above $\approx 200,000$, the Strouhal number no longer remains constant, the separation point shifts to much later on the cylinder surface, and cylinder roughness also takes on a more significant role (as shown by Lienhard (Lienhard, 1966) and Achenbach and Heinecke (Achenbach and Heinecke, 1981)). In the range of the present set of experiments ($Re = 100,000$ to $160,000$), we remain in the ‘low-quality’ vortex shedding regime, with a nearly constant Strouhal number of 0.2. The definition of Strouhal number is as follows:

$$St = \frac{fD}{U} \qquad (2)$$

Where f is the frequency of vortex shedding. To characterize the Strouhal number variation with Reynolds number across these regimes, many fundamental studies exist (e.g., Lienhard, 1966; Achenbach and Heinecke, 1981; Roshko, 1955). Despite this range of regimes, and while single-phase flow past cylinders is still of interest, it can be considered to be well understood (Bearman, 1984; Gabbai and Benaroya, 2005; Bearman, 2011; Williamson, 1996). However, multiphase flow over a cylinder is much less understood.

A summary of experimental studies on multiphase cross-flow over a cylinder is given in Table 1. In most of these works, the flow is forced and as such the Reynolds number is given by the freestream velocity. However, Uchiyama and Ishiguro (2016), Lee and Park (2020), and Murai et al. (2005) study the motion of bubble driven flows, and as such the Reynolds number is defined using the superficial gas velocity, $u_g = \frac{Q_g}{A}$, where $Q_g$ is the gas flow rate and A is the cross sectional area. The reported Froude number for horizontal forced flows in Table 1 is computed as follows:

$$Fr = \frac{U}{\sqrt{gD}} \qquad (3)$$

Where g is the acceleration due to gravity. Importantly, many of these studies demonstrated that the Strouhal number is not constant at $St = 0.2$ when gas is injected. This change in Strouhal number with void fraction and bubble size is poorly understood both in terms of predicting the change in Strouhal number as well as the mechanisms leading to these changes. The remainder of this literature review will focus on changes to the Strouhal number, vibration of the cylinder,

and alterations to other aspects of the flow topology under the influence of gas injection. While for brevity only experimental works are summarized in Table 1, selected numerical studies are also reviewed.

A feature of primary importance is an understanding of how gas redistributes in the wake of the cylinder. Hulin et al. (1982) demonstrate that bubbles are trapped within vortices, for which the local phase fraction within the vortex increases with Reynolds number. This effect has also been shown numerically, for which Uchiyama and Degawa (Uchiyama and Degawa, 2006) found that bubbles captured in the shed vortices have a local void fraction that is approximately three times higher than the upstream value (which is varied from 0–0.03). Sugiyama et al. (2001) note that there is a local peak in the local void fraction at a distance of one diameter downstream, and posit that this is due to the preferential accumulation of bubbles in the shed vortex at this location. For bubbles to accumulate in the wake vortices, Sugiyama et al. (2001) state that a sufficiently large bubble relaxation time (i.e., sufficient bubble inertia) is required. Importantly, they noted that bubble accumulation in the vortices increased when accounting for buoyancy effects, as the bubble relaxation time increased. The dependence of vortex capture on bubble size has also been shown experimentally (Meng, 1993). The distribution of bubbles rising solely from buoyancy around a circular cylinder was studied by Uchiyama and Ishiguro (2016), for which it was found that in addition to bubbles entrained in the large scale eddies behind the cylinder, bubbles were also trapped in the stagnation point at the front of the cylinder. In particular, Suzuki et al. (2009) show that microbubbles tend to accumulate in the front stagnation point, then are repulsed from the boundary layer near 90° from the front stagnation point. The presence of microbubbles alters the boundary layer dynamics and reduces the turbulence intensity (‘laminarizing’ the flow). Lee and Park (2020) study the gravity induced liquid motion of an oblate ellipsoid bubble swarm past a circular cylinder. A vortex street is not formed, but Lee and Park (Lee and Park, 2020) note that an elongated region of low void fraction in the wake of the cylinder is reduced in size at increasing void fraction. At higher void fractions, shear induced lift is sufficiently large to cause earlier bubble migration to the cylinder wake. These topics are further discussed in the thesis of Lee (Lee, 2021). Kim and Rau (Kim and Rau, 2020) study the liquid only region in the wake of the cylinder, and note that this region decreases in size with increasing Reynolds number.

In addition to changes in the wake topology, numerous studies have assessed the influence of bubbles on the shedding frequency and cylinder vibration. Hulin et al. (1982) found that there is an increase in both the mean frequency and bandwidth of pressure oscillations on the cylinder. Furthermore, they discovered that at high phase fractions ($\approx 10\%$), vortex shedding from the cylinder becomes less regular. Pascal-Ribot and Blanchet (Pascal-Ribot and Blanchet, 2007, 2011) study the lift force on a single rigid cylinder under cross flow in the range of void fraction from 20%–80%. They note that as void fraction increases, the presence of low frequency components in the signal increases.

Hara (1987) demonstrated that air bubbles can reduce the Kármán vortex-induced vibration in both the in-line and cross-flow directions. It is stated that “the reasons for this effect are not quantitatively clear, and we do not yet know all of the parameters influencing this vibration reduction effect”. (Hara, 1987) With regards to cross-flow, Hara (Hara, 1987) found that small bubbles are more effective than large bubbles in reducing vibration. It was proposed that the stochastic motion of air bubbles in the flow was a factor in producing an unsteady lift force-effectively reducing cross-flow vibration. At $Re_D > 130,000$, Watanabe et al. (1990) demonstrated that even at low void fractions of 0.05%, there is a decrease in the drag force, as well as a decrease in the fluctuating lift and drag forces when compared to the single-phase case. Furthermore, the frequency of vortex shedding was increased to values that are typically seen in the supercritical flow range. Voutsinas et al.
<table><caption>**Table 1**
Parameter range of previous experimental studies.</caption><thead><tr><th>Paper</th><th>Phase fraction [%]</th><th>Bubble size [mm]</th><th>Re<sub>D</sub> = u<sub>D</sub>/v</th><th>Cylinder diameter [mm]</th><th>Cross-Section</th><th>Channel size [mm]</th><th>Flow orientation</th></tr></thead><tbody><tr><td>Uchiyama and Ishiguro (2016)</td><td>Not Reported</td><td>0.02–0.06</td><td>702, 981</td><td>30</td><td>Circular</td><td>440 x 440</td><td>Vert. (bubble-driven)</td></tr><tr><td>Lee and Park (2020)</td><td>0.3–2.1</td><td>3.7</td><td>Not Reported</td><td>20, 30</td><td>Circular</td><td>670 x 240</td><td>Vert. (bubble-driven)</td></tr><tr><td>Murai et al. (2005)</td><td>4–20</td><td>1.5–2.5</td><td>Not Reported</td><td>10, 50</td><td>Circular, Square, Triangular, Ellipsoid (2), Star</td><td>10 x 500</td><td>Vert. (bubble-driven)</td></tr><tr><th>Paper</th><th>Phase Fraction [%]</th><th>Bubble Size [mm]</th><th>Re<sub>D</sub> = u<sub>D</sub>/v</th><th>Cylinder Diameter [mm]</th><th>Cross-Section</th><th>Channel Size [mm]</th><th>Flow Orientation</th></tr><tr><td>Present Study</td><td>0.01–0.14<sup>a</sup></td><td>0.3–1.2</td><td>100,000; 160,000</td><td>20</td><td>Circular</td><td>192 x 192</td><td>Horz. Fr = 11.3,18.1</td></tr><tr><td>Kim and Rau (2020)</td><td>Not Reported</td><td>1–2.4–10</td><td>96–2868</td><td>9.5</td><td>Circular</td><td>136.5 x 44.5</td><td>Vert.</td></tr><tr><td>Habeeb and Al-Turaihi (2013)</td><td>22.2–66.7</td><td>Not Reported</td><td>14,928–33,439</td><td>30</td><td>Circular, Square, Triangular</td><td>30 x 100</td><td>Horz. Fr = 0.9–2.0</td></tr><tr><td>Pascal-Ribot and Blanchet (2011)</td><td>10–80</td><td>3–14.9</td><td>Not Reported</td><td>12.15,15.85,20.7 25.6,31.9</td><td>Circular</td><td>70 x 100</td><td>Vert.</td></tr><tr><td>Suzuki et al. (2009)</td><td>0.03</td><td>0.09–0.43</td><td>46,000</td><td>20</td><td>Circular</td><td>80 x 2</td><td>Horz. Fr = 4.7</td></tr><tr><td>Voutsinas et al. (2009)</td><td>0–5</td><td>0.25–2.6</td><td>5000; 10,000</td><td>30</td><td>Square</td><td>90 x 60</td><td>Vert.</td></tr><tr><td>Pascal-Ribot and Blanchet (2007)</td><td>10–80</td><td>3–14.9</td><td>Not Reported</td><td>12.15</td><td>Circular</td><td>70 x 100</td><td>Vert.</td></tr><tr><td>Shakouchi et al. (2002)</td><td>0–30</td><td>2.48</td><td>5000–15,000</td><td>18, 22.5, 30, 45</td><td>Square</td><td>90 x 45</td><td>Vert.</td></tr><tr><td>Iijima et al. (1995)</td><td>0–33</td><td>3–4</td><td>12,200–40,000</td><td>30</td><td>Circular</td><td>200 x 60</td><td>Vert.</td></tr><tr><td>Meng (1993)</td><td>4–20</td><td>0.5–4</td><td>3000–45,700</td><td>2,6,11</td><td>Circular</td><td>Circular D = 40</td><td>Vert.</td></tr><tr><td>Watanabe et al. (1990)</td><td>0–0.48</td><td>0.25–2.75</td><td>130,000–270,000</td><td>30</td><td>Circular</td><td>200 x 610</td><td>Horz. Fr = 8.3–14.7<sup>b</sup></td></tr><tr><td>Hara and Iijima (1989)</td><td>0–20</td><td>5–6</td><td>2000–32,500</td><td>25</td><td>Circular</td><td>200 x 60</td><td>Vert.</td></tr><tr><td>Joo and Dhir (1994)</td><td>0–50</td><td>Not Reported</td><td>430–21,900</td><td>22</td><td>Circular</td><td>125 x 200</td><td>Vert.</td></tr><tr><td>Yokosawa et al. (1986b)</td><td>0–10</td><td>3–5</td><td>5000–80,000</td><td>10, 20, 30, 40</td><td>Cylinder Flat Plate<sup>c</sup></td><td>60 x 120</td><td>Vert.</td></tr><tr><td>Yokosawa et al. (1986a)</td><td>0–10</td><td>3–5</td><td>4000–300,000</td><td>1, 3, 5, 6, 8, 10, 20, 30, 40</td><td>Circular</td><td>60 x 120</td><td>Vert.</td></tr><tr><td>Inoue et al. (1986)</td><td>0–24</td><td>3–5</td><td>5000–80,000</td><td>10, 20, 30, 40</td><td>Circular</td><td>120 x 60</td><td>Vert.</td></tr><tr><td>Hara (1984)</td><td>0–28.5</td><td>Not Reported</td><td>9000–18,000</td><td>30</td><td>Circular</td><td>200 x 60</td><td>Vert.</td></tr><tr><td>Hulin et al. (1982)</td><td>0–25</td><td>2–4</td><td>0–110,000</td><td>35 upstream, 23 downstream</td><td>Trapezoidal</td><td>Circular D = 150</td><td>Vert.</td></tr></tbody></table>

(2009) demonstrate that changing the bubble size affects both the shedding frequency and amplitude. For the three studied bubble sizes (albeit not well characterized and also far from truly monodisperse (Shakouchi et al., 2007)) in their work, there are bi-modal frequency peaks at both average bubble sizes *d* = 2.6 mm and *d* = 0.27 mm. This is an indication that the shedding frequency depends strongly on the bubble size, and that the governing phenomena describing this relationship is much more complicated than the early proposed relationship between void fraction and shedding frequency that is summarized in the Encyclopedia of Fluid Mechanics (Cheremisinoff, 1986). This is not a surprising result, since it makes sense that bubbles will be entrained in the vortices shed from the cylinder at different rates based on their size. Furthermore, the proposed relationship dictates that the shedding frequency should increase with increasing void fraction (Cheremisinoff, 1986). This has also shown to be false, as studies such as that of Uchiyama (Uchiyama, 1999), who conducted a numerical study on flow over rectangular cylinders of varying thickness-to-width ratios at Re = 20000, found that the Strouhal number decreases monotonically with increasing volumetric gas phase fraction from 0.025 to 0.075, regardless of the thickness-to-width ratio.

Several studies have attempted to use the knowledge of gas redistribution in the wake to propose mechanisms by which the shedding frequency may be altered; however, a consensus for this proposed mechanism has not yet been reached. Uchiyama and Degawa (Uchiyama and Degawa, 2006) propose that the increased shedding frequency is due to increased water velocity from the presence of bubbles in the flow. Cook and Harlow (1986) perform numerical studies on the two way coupled motion of bubbles passing around a rectangular obstacle. They note that bubble concentration centers within the shed vortices move more rapidly than the vortices themselves due to buoyancy, thereby pulling the vortices away from the cylinder more rapidly than they would otherwise be shed. In addition to this important finding, they note that the presence of bubbles alters the effective width of the vortex sheet, thereby providing a mechanism for the shedding frequency to change. Suzuki et al. (2009) reported that the shedding frequency increased with the introduction of microbubbles. Their hypothesis for the increase in the shedding frequency is the expansion of the laminar boundary layer coupled with an earlier separation point. Meng (Meng, 1993) note that coupling between bubbles and vortices can suppress the shedding

<footnote><sup>a</sup> Volumetric quality is equal to phase fraction if no gas hold up exists.</footnote>

<footnote><sup>b</sup> Velocities calculated from Re<sub>D</sub>, assuming v = 1x10<sup>-6</sup>.</footnote>

<footnote><sup>c</sup> Plate is 1 mm thick with chamfered edges.</footnote>
frequency due to momentum coupling, but significant attention is not given to this result.

This review of experimental and numerical works studying bubbly flow around a cylinder shows that:

1. At phase fractions above 10%, a gas pocket builds behind the obstacle and the vortex emission in air-water flows becomes less stable.

2. At low void fractions (below 10%), the Strouhal number increases with phase fraction (in most, but not all studies).

3. Bubbles tend to accumulate in vortices downstream.

4. There is a region with locally low void fraction in the wake of the cylinder which varies in size with obstacle shape, Reynolds number, and upstream void fraction.

5. The changes in shedding frequency and cylinder vibration experienced with multiphase flows are bubble size dependent.

6. Bubbles can alter the separation point.

While the above effects are of interest, there is not yet a complete theory for the mechanism by which the shedding frequency changes with bubble size and void fraction. One would expect that such a theory would incorporate the time and spatially-dependent nature of bubble motion within the wake.

While no definitive mechanism has been identified, the explanations reviewed above provide a marked improvement upon the early work presented in the Encyclopedia of Fluid Mechanics (Cheremisinnoff, 1986), as the bubble size effect is accounted for. However, the underlying mechanism may not be universal across the parameter space spanned by Reynolds number, void fraction, and bubble size. It has been noted by Sugiyama et al. (2001) and Uchiyama (Uchiyama, 1999) that there are cases in which the shedding frequency does not increase with increasing phase fraction. Given this lack of consensus on both the nature of the change in shedding frequency as well as the proposed mechanism, more work is needed to identify what effect bubbles have on the shedding frequency.

In this paper, we present experimental results on flow past a circular cylinder at high Reynolds numbers in the subcritical regime. Experiments are conducted for $Re_D = 100,000$ and $Re_D = 160,000$. As demonstrated in Table 1, this subcritical regime at $Re_D > 100,000$ is under-explored at low void fraction. This flow regime is challenging to explore, as pressurization of the test section is required to avoid cavitation at this Reynolds number range. The present study complements the work of Watanabe et al. (1990), by exploring smaller air volume fraction and smaller bubble size. Indeed, in the context of wake flow or vibration control using bubble injection, reducing the air injection rate can lead to more efficient control methods. Six different effective void fractions are studied, by changing the region of the flow upstream in which bubbles are injected. For single-phase flow conditions, the turbulent flow in the wake is characterized by time-resolved particle image velocimetry (TR-PIV), proper orthogonal and spectral proper orthogonal decomposition (POD and SPOD). The vibration of the cylinder is compared between the single- and two-phase bubbly flows. The topology and temporal periodicity of the bubble arrangement is examined in the near and far wakes of the cylinder through high speed visualization, as well as POD and SPOD of the reflected light intensity of the bubbles. Furthermore, the Eulerian gas-phase averaged velocity is computed.

We demonstrate that even at volumetric qualities of less than 0.01% and high Froude number, the bubbles can interact with the near wake and modify the shedding frequency, while also suppressing the vortex induced vibration of the cylinder. The outline of the paper is as follows. Section 2 describes the experimental facility, flow conditions, instrumentation, and processing methods. Section 3 presents the results. Section 4 is a discussion of the main results. Conclusions are presented in Section 5.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Experimental facility. Bubble injector can be oriented horizontally or vertically.</figcaption></figure>

## 2. Methods

### 2.1. Experimental facility

The experiments were conducted in the IRENav facility, located at École Navale, France. A detailed description of the facility is given by Watine et al. (2023), Pernod et al. (2019), and relevant information is re-summarized here. The test section is 1 m long, and has a 192 mm square cross section. The test section is located downstream of a contraction section of contraction factor 1:9.8.

The tunnel was filled with water, and the flow was studied past a cantilevered circular cylinder. The cylinder had a diameter of $D = 20$ mm, and a span of 191 mm, which was clamped to the backside of the test section. Two upstream velocities were investigated: $U = 5$ m/s and $U = 8$ m/s, which were controlled with an accuracy of 2%. To avoid cavitation in the wake of the cylinder, the test section was pressurized. The pressure at the entrance of the test section was adjusted to 1600 mBar and kept constant with an accuracy of 2.5%.

Bubble were injected using an injector array within the tunnel contraction located upstream of the cylinder (as shown in Fig. 1). The injector array could be oriented either vertically (along the center of the channel in the x-y plane) or horizontally (along the centerline of the cylinder in the x-z plane).

The bubble injector, developed by YLEC-Consultants, consists of a hydrofoil with 288 air injection points. The hydrofoil has a span of 600 mm, which is the width of the upstream section of the tunnel. The profile is based on a truncated NACA0024 hydrofoil with an effective chord length of 67.82 mm. Each injection point is made of a needle with an 800 µm outer diameter and 125 µm inner diameter, which is inserted into a concentric hole with 1 mm diameter. Air is injected through the needle. Water is continuously injected through the annular gap between the needle and the hole, so as to have improved control over the bubble size. The global flow rate of water through the injector was kept constant at 38 L/min.

With the injector inserted upstream within the contraction section, the turbulence intensity is homogeneous and largely unchanged when compared to the case with no injector mounted. With the water injection activated (regardless of the orientation of the injector), the turbulence intensity is 2.5% at the entrance of the test section. This is compared to a turbulence intensity of 2% with no injector mounted. With this in mind, the flow past the cylinder with only water injected through the hydrofoil is considered as our single phase flow of reference.

To characterize the quantity of gas present in the flow, we use the volumetric quality. The volumetric quality is computed as $\beta = \frac{Q_g}{Q_g+Q_w}$ where $Q_g$ is the gas volumetric flow rate and $Q_w$ is the water flow
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Volume weighted histogram of bubble sizes, measured upstream of the cylinder. The labeled volumetric qualities correspond to the vertical orientation of the injector, but the histogram is assumed to be unchanged in the horizontal injection direction.</figcaption></figure>

rate. In the case that the slip velocity, $U_{slip} = U_g - U_{wv}$, is zero (no gas hold up exists), the volumetric quality is exactly equal to the void fraction. The values for volumetric quality are adjusted from standard pressure to account for the pressurization of the test section. We report the volumetric quality as a spanwise-average.

It is important to note that when shifting from horizontal to vertical injection the *effective* volumetric quality changes, as in the case of vertical injection not all gas interacts with the cylinder. The vertical width at which gas is pulled into the wake is estimated as 2D. This width is computed as the distance at which the vertical bubble velocity is not significantly different from the single phase case in the near wake region ($x/D \lesssim 2.5$); this is demonstrated in Fig. 17. To ensure consistency, the volumetric qualities that are reported in Table 2 are computed using both water and gas flow rates passing through this width. It is assumed that all gas injected in the horizontal configuration interacts with the cylinder wake. Therefore, the value for the volumetric quality, $\beta$, increases by a factor of 4.8x when shifting from the vertical injector configuration to the horizontal configuration. The histograms of the produced bubble sizes are given in Fig. 2. The 90% interval of bubble sizes *weighted by volume*, are summarized in Table 2, along with the operating points that were investigated in the experiments. As can be seen in Table 2, the range of bubble sizes within which 90% of the total injected gas is contained decreases with Reynolds number. For each of the conditions given in Table 2, experiments were conducted with both the horizontal and vertical bubble injector orientation. Consequently, there are a total of 16 test conditions.

## 2.2. Laser vibrometry

Laser vibrometry provides a non-intrusive method to measure the vibration of the cylinder. The technique uses the Doppler effect, for which the phase shift between the emitted and reflected laser signals relate to the displacement velocity of the cylinder. In the present experiments, the vibration of the cylinder was measured using the Polytec PDV-100; the vibrometer was oriented to measure the displacement in the direction perpendicular to the flow at one point located on the tip of the cylinder which is free to move. The vibrometer signal was sampled at 2048 Hz for a duration of 8s. The power spectral density of the vibration velocity of the cylinder was computed with Welch's method, using 1024 sample segments, windowed with a Hamming window and 50% overlap. This method yielded a frequency resolution of 2 Hz.

<table><caption>**Table 2**
Operating points.</caption><thead><tr><th>D (mm)</th><th>U (m/s)</th><th>Q<sub>z</sub><sup>a</sup> (mL/min)</th><th>d<sub>q0%</sub> Interval (mm)</th><th>β<sub>V<sub>cert.</sub></sub> (%)</th><th>β<sub>H<sub>orz.</sub></sub> (%)</th></tr></thead><tbody><tr><td>20</td><td>5</td><td>0</td><td>N/A</td><td>N/A</td><td>N/A</td></tr><tr><td>20</td><td>5</td><td>1700</td><td>0.2–1.2</td><td>0.010</td><td>0.046</td></tr><tr><td>20</td><td>5</td><td>3400</td><td>0.3–2.2</td><td>0.019</td><td>0.092</td></tr><tr><td>20</td><td>5</td><td>5100</td><td>0.3–2.9</td><td>0.029</td><td>0.14</td></tr><tr><td>20</td><td>8</td><td>0</td><td>N/A</td><td>N/A</td><td>N/A</td></tr><tr><td>20</td><td>8</td><td>1700</td><td>0.2–1.0</td><td>0.006</td><td>0.029</td></tr><tr><td>20</td><td>8</td><td>3400</td><td>0.2–1.5</td><td>0.012</td><td>0.058</td></tr><tr><td>20</td><td>8</td><td>5100</td><td>0.3–2.0</td><td>0.018</td><td>0.086</td></tr></tbody></table>

## 2.3. Flow measurement techniques

The time-resolved flow was characterized using particle image velocimetry (PIV) of the single phase flow and high speed recording of the bubbles within the two-phase flow. Subsequent analysis on the PIV data was performed using proper orthogonal decomposition (POD) and spectral proper orthogonal decomposition (SPOD). For the two-phase flow, bubble image velocimetry (BIV) was performed using the bubble images, as well as SPOD analysis. Each of these analysis methods are described in detail in the subsequent sections.

### 2.3.1. Time-resolved particle image velocimetry (PIV)

Particle image velocimetry measurements were solely taken in the single-phase case. 2D-2C PIV measurements were taken using an Nd:YLF laser with 527 nm wave length. The laser plan was vertical and aligned with the $z=0$ symmetry plane of the channel. The flow was seeded with Polyamide particles with a diameter of $10 \pm 2$ µm and a density of 1.03 g/cm$^3$. The post-processing was performed using Dynamic Studio 7.6-77.0. The flow field was imaged using a 12-bit 2048 x 1952 px$^2$ Phantom v1840 time-resolved PIV camera with 13.5 µm pixel size, which was equipped with a 100 mm Zeiss Makro Planar lens. The scale factor was 61 µm/pixel. The field of view was 6.3D in the streamwise direction, with a height of 5.6D. The time between pulses was adjusted to 64µs and 40µs for $Re_D=100,000$ and 160,000 respectively. Image balancing was applied on the acquired images before PIV processing. The adaptive PIV algorithm was used with a final interrogation area of $16 \times 16$ px$^2$ and no overlap. An example of the raw particle image with the corresponding velocity field is given in Fig. 3.

The analysis of the PIV velocity field using SPOD requires a long measurement time, particularly in flow fields where low frequency structures are of interest. This is true for the case of the recirculating region in the near wake, however, the spectral analysis must also capture the frequency of the primary vortex shedding. Therefore, a compromise was chosen between the sampling frequency and total measurement time. The flow field was imaged at a sampling frequency of 250 Hz for 19.5 s in the $Re_D=100,000$ case, and at a sampling frequency of 400 Hz for 12.2 s in the case of $Re_D=160,000$.

The right hand coordinate system used in the experiments is as follows:

1. The x-direction is aligned with the direction of the flow, and $x=0$ is located at the rear stagnation point on the cylinder.

2. The y-direction is oriented vertically, opposite to the direction of gravity, and $y=0$ is located at the cylinder centerline.

3. The z-direction is oriented along the span of the cylinder, and $z=0$ is located at the centerline of the channel.

For each of the operating points described in Section 2, two independent planar PIV measurements were taken. The laser planes are referred to as the close and far position, and are described in the sketch in Fig. 4. The close position includes the recirculation region immediately downstream of the cylinder. Specifications for the PIV processing are summarized in Table 3.

<footnote>a Standard gas flow rate, for which the volume is given at 1 bar pressure.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Raw image of particles with overlaid mean-subtracted velocity vectors obtained by adaptive PIV. For clarity, only every 3rd vector is shown. The image corresponds to single phase flow at $Re_D=160,000$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. PIV planes used to characterize the single phase flow, along with coordinate axes.</figcaption></figure>

<table><caption>**Table 3**
Specifications for PIV acquisition and processing.</caption><tr><td>Re<sub>D</sub>=100,000/Re<sub>D</sub>=160,000</td></tr><tr><td>PIV acquisition settings</td></tr><tr><td>• Camera: 2048x1952 px<sup>2</sup></td></tr><tr><td>• Lens: 100 mm</td></tr><tr><td>• Field of view: 6.3D x 5.6D</td></tr><tr><td>• Acquisition rate (double frame mode): 250 Hz/400 Hz</td></tr><tr><td>• Minimum number of double-frame pairs for ensemble averaging: 4881</td></tr><tr><td>Pre-processing of acquired images</td></tr><tr><td>• Image balancing with a lightsheet balance map using a cell size of 5x5</td></tr><tr><td>Adaptive PIV specifications:</td></tr><tr><td>• No overlap</td></tr><tr><td>• Interrogation area size: 16x16 px<sup>2</sup><br>(Spatial resolution: 0.049Dx0.049D)</td></tr><tr><td>Default validation criteria:</td></tr><tr><td>• Peak height ratio: 1.5</td></tr><tr><td>• Universal outlier detection method: 5x5 px<sup>2</sup></td></tr><tr><td>• Wall windowing with cylinder mask</td></tr></table>

### 2.3.2. Bubble high speed recording

The bubbles were imaged in a vertical plane, using a 12-bit Phantom v611 camera with 20 µm pixel size. Two spotlights, positioned one above and one below the flow, were used to illuminate the bubbly flow. The minimum frame rate to capture the two-phase images was 1000 Hz, and was as high as 17,000 Hz for selected cases. The images were taken at two positions downstream. These positions vary based on whether

<table><caption>**Table 4**
Specifications for BIV acquisition.</caption><thead><tr><th></th><th>Re<sub>D</sub>=100,000/160,000</th><th>Far position</th></tr><tr><th></th><th>Close position</th><th></th></tr></thead><tbody><tr><td>Field of View</td><td colspan="2">3.5Dx3.4D</td></tr><tr><td>Acquisition Rate (single frame mode)</td><td>4000 Hz/6400 Hz</td><td>17,000 Hz</td></tr><tr><td>Number of Images</td><td>6400 minimum</td><td>12800 minimum</td></tr><tr><td>Interrogation area size [px<sup>2</sup>]<sup>a</sup></td><td>64x64</td><td>32x32</td></tr><tr><td>Resolution</td><td>0.22D</td><td>0.11D</td></tr></tbody></table>

the bubble injector was positioned vertically or horizontally upstream, as the camera lens that was used shifted from one case to the other.

1. Vertical bubble injection images. Close position is from -0.6D to 2.9D, and far position is from 8.6D to 12.1D.

2. Horizontal bubble injection images. Close position is from -0.6D to 4.3D, and far position is from 8.1D to 13.0D.

The scale factor for the vertical and horizontal injector position was 68.7 µm/px and 94.3 µm/px, respectively. A sample image for the case of $\beta = 0.006$ at $Re_D = 160,000$ is shown in Fig. 5. The raw two-phase images provide the basis for the bubble image velocimetry (BIV) and two phase SPOD discussed below.

### 2.3.3. Bubble image velocimetry

Bubble image velocimetry (BIV) employs similar analysis techniques to PIV; however, bubbles are used in place of seeding particles. BIV has been implemented in a variety of applications in literature, from the study of plunging jets (Ryu et al., 2005), to rising bubble plumes (Cheng et al., 2005). In the context of the present experiments, the cross-correlation was applied directly to the bubble intensity images in order to determine the gas-phase velocity. Bubble image velocimetry was only applied on the images with vertical gas injection. Notably, the image represents a spanwise average over the depth of field instead of highlighting a single plane (as is the case for PIV). Depending on the lighting conditions of the image, either 64 × 64 or 32 × 32 interrogation windows with no overlap were used. Since the images were taken from a single camera, the time between frames was fixed at the repetition rate of the camera. Consequently, a higher frame rate (and subsequent reduction of the image resolution) was required, when compared to the PIV measurements. Table 4 summarizes the acquisition settings for the BIV analysis.

Some of the challenges of applying cross-correlation directly to bubble images are highlighted by Cheng et al. (Cheng et al., 2005), such as (1) bubble overlap, (2) irregular motion which is intrinsic to bubbles, (3) time- and spatially-varying scattering of light at the bubble interface, and (4) bubble size-dependence. These challenges are particularly present in the near wake behind the cylinder, as this region experiences the greatest degree of bubble overlap. Consequently, the application of BIV in the present work is to characterize the statistical motion of the gas phase in the wake, by comparing the 25th, 50th (median), and 75th percentile of velocity. These percentiles are computed as the percent of time at which that velocity is exceeded. These measures are expected to be more robust to individual velocity fluctuations, rather than relying on techniques such as POD/SPOD which incorporate time-resolved measurements. Regarding the influence of bubble size on the computed velocity, an analysis of the expected range of Stokes numbers is presented in Section 3.3, along with discussion on the degree to which bubble motion is expected to represent the underlying flow field. Results from the BIV analysis are presented in Section 3.6.

<footnote><sup>a</sup> With no overlap.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Raw two-phase images at $Re_D$ = 160,000, with vertical gas injection. The gray area is in shadow due to the positioning of the spot lights.</figcaption></figure>

### 2.3.4. POD and SPOD analysis

A review of the Proper Orthogonal Decomposition (POD) and Spectral Proper Orthogonal Decomposition (SPOD) techniques, as well as their application within the context of fluid mechanics, is provided by Towne et al. (Towne et al., 2018). Applying their terminology, within the present work we use POD as shorthand for “space-only POD”. POD is applied to decompose a time-resolved dataset of 2D matrices (either two-phase image data or PIV data) into a series of orthogonal 2D modes, which are modulated in time by a time-series of expansion coefficients. Let us denote $\mathbf{x}$ as the vector position in the 2D plane $(x,y)$. We define $\phi_{j,k}(\mathbf{x})$ as the eigenvector of the $k$th mode of the velocity component in direction $j$, such that for a set of $N$ fluctuating velocity fields we can define the following expansion:

$$u'_{j}(\mathbf{x}, t) = u_{j}(\mathbf{x}, t) - U_{j}(\mathbf{x}) = \sum_{k=1}^{N} a_{j,k}(t) \phi_{j,k}(\mathbf{x}) \quad (4)$$

We can further define $\lambda_{j,k}$ as the eigenvalue which is representative of the average energy captured by mode $k$ in the direction $j$. The 2D modes, $\phi_{j,k}(\mathbf{x})$, are selected in order to capture the greatest fraction of the total variance when compared to any other orthonormal basis. Notably these modes are spatially static, and simply represent spatially coherent structures among the flow; these modes do not necessarily have any coherent evolution in time. For the present work, the time-series of expansion coefficients, $a_{j,k}$, is presented as a discrete Fourier transform, since the dataset is time-resolved and ordered. For POD modes in which there is no temporal coherence in the oscillation of the modal energy, no dominant frequencies will appear in this discrete Fourier transform.

In order to capture spatially- and time-resolved modes, we use SPOD. Beginning with the same time-resolved dataset of 2D matrices, we apply Welch's method and split the dataset into a series of $M$ blocks. A pixel-wise discrete Fourier transform is applied to each of the blocks, and in doing so, each entry in the 2D matrix now will have $M$ complex Fourier coefficients associated with each frequency. For a given frequency we now have $M$ 2D matrices of complex coefficients, for which we can apply conventional POD to determine the optimal set of orthogonal 2D modes to capture the variance across the $M$ realizations. These $M$ realizations are not temporally correlated, and in fact may be constructed from partially overlapping data in the original dataset. Consequently, the output of the SPOD analysis is a set of 2D complex orthonormal modes which correspond to each frequency in the discrete Fourier transform. The complex nature of these orthonormal modes allows for both spatial and temporal coherence across the 2D matrix. Additional detail on the theory behind POD and SPOD is provided by Schmidt and Colonius (Schmidt and Colonius, 2020).

Within this work, POD and SPOD are applied to two types of datasets: (1) single-phase velocity data from PIV, and (2) bubble image intensity data from the two-phase experiments. Notably, in both cases the techniques are identical; simply the size of the 2D matrices and quantity of time-resolved 2D matrices will vary. With that said, the physical meaning of the analysis outputs is different. The POD modal shapes computed using PIV data describe spatial coherent fluctuations in velocity along the centerline of the channel; these fluctuations are relative to the mean. By comparison, SPOD modal shapes of the PIV dataset correspond to regions of spatially and temporally correlated velocity, which collectively fluctuate at a particular frequency. Importantly, the 100 Hz peak which corresponds to the frequency of the lights used to illuminate the two-phase flow is removed from the SPOD energy spectra for analysis of the gray-scale two-phase (bubble intensity) images.

In contrast to the PIV dataset, the bubble intensity images highlight bubble edges, where light is reflected to the camera. The intensity of the image increases with bubbles density, which allows us to draw qualitative information on the alternating bubble accumulation and release in the wake. The depth of focus of the images is on the order of the diameter of the cylinder or larger, which means that the bubble accumulation information is depth-averaged. Therefore, the SPOD modes of the two-phase images correspond to regions of the image that are spatially and temporally correlated in terms of the depth-averaged accumulation of gas; the corresponding frequency describes the depth-averaged periodic gathering/ejection of gas in the wake. Given the differences in the meaning of the SPOD modes between the single-phase PIV data and the two-phase image data, care must be taken when drawing conclusions. The processing settings used for SPOD analysis of the two-phase images are given in Table 5.

## 3. Results

### 3.1. Single-phase flow topology downstream of cylinder

The flow downstream of the cylinder was characterized with PIV. Fig. 6 illustrates the time-averaged near wake distribution of quantities such as the streamwise velocity, $V_x$, the vertical velocity, $V_y$, and the Reynolds shear stress, $u'v'$. One can see the symmetrical structure of the mean streamwise velocity along the $y = 0$ axis, along with the existence of the recirculating region characterized by time-averaged reversed flow in the near wake. Furthermore, we observe a classical anti-symmetric distribution of the Reynolds shear stress.

To demonstrate that the single-phase flow aligns with existing literature, wake statistics are computed; these wake statistics require
<figure><img src="image_8.png" /><figcaption>Fig. 6. Single-phase time-averaged flow wake conditions for $Re_D$ = 160,000 with a horizontal upstream injector orientation.</figcaption></figure>

<table><caption>**Table 5**
Settings for bubble intensity SPOD processing.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Re<sub>D</sub>=100,000/160,000</th></tr><tr><th>Close position</th><th>Far position</th></tr></thead><tbody><tr><td rowspan="2">Field of View</td><td>Vert. Inj.</td><td>3.5Dx3.4D</td><td></td></tr><tr><td>Horz. Inj.</td><td>4.8Dx3.6D</td><td></td></tr><tr><td rowspan="2">Spatial Resolution (undersampled)</td><td>Vert. Inj.</td><td>0.007D</td><td></td></tr><tr><td>Horz. Inj.</td><td>0.009D</td><td></td></tr><tr><td colspan="2">Acquisition Rate (undersampled)</td><td colspan="2">500 Hz-3200 Hz</td></tr><tr><td colspan="2">Number of Snapshots per Block</td><td colspan="2">1024</td></tr><tr><td colspan="2">Block Overlap</td><td colspan="2">50%</td></tr><tr><td colspan="2">Minimum Number of Blocks</td><td colspan="2">24</td></tr><tr><td colspan="2">Minimum Frequency Resolution</td><td colspan="2">3.1 Hz</td></tr></tbody></table>

a sufficient number of independent ensemble members for convergence. To determine the number of ensemble members within the measurement time, $t_{meas}$, the correlation time, $t_{corr}$, is computed. The correlation time is defined as the time at which the auto-correlation of the streamwise velocity signal first reaches zero (Kundu et al., 2016). The number of ensemble members is computed as $N_{members} \approx t_{meas}/t_{corr}$. When computed using the maximum correlation time in the near wake region $x > 0$, $-2 < y/D < 2$, the value of $N_{members}$ ranges from 31-125 across the test cases. When computed using the median wake

<table><caption>**Table 6**
Near wake statistics for single-phase flow.</caption><thead><tr><th>Re<sub>D</sub></th><th>Injector Orient.</th><th>Recirculation length</th><th>Eddy length</th><th>Uncertainty</th></tr></thead><tbody><tr><td>100,000</td><td>Vert.</td><td>1.5D</td><td>1.7D</td><td>±0.025D</td></tr><tr><td>100,000</td><td>Horz.</td><td>1.4D</td><td>1.7D</td><td>±0.025D</td></tr><tr><td>160,000</td><td>Vert.</td><td>1.4D</td><td>1.5D</td><td>±0.025D</td></tr><tr><td>160,000</td><td>Horz.</td><td>1.3D</td><td>1.5D</td><td>±0.025D</td></tr></tbody></table>

correlation time, the number of ensembles increases to a minimum of 1940. With this in mind, the measurement time is sufficient to compute converged wake statistics. The wake statistics are shown in Table 6, with the uncertainty reported as 1/2 of the interrogation window size. The mean recirculation length is determined as the streamwise distance along $y = 0$ between the rear stagnation point of the cylinder and the mean stagnation point of the wake (i.e. the point where the mean streamwise velocity, $V_x$, changes sign). The eddy length is computed as the streamwise location of maximum $u'/u'^2$. Indeed, along the wake axis, the point at which the Reynolds stress reaches its highest magnitude is representative of the vortex formation length (Bearman, 1965).

The recirculation length tends to drop slightly when switching from $Re_D$ = 100,000 to 160,000. With the above in mind,
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Mode 1–3 of POD analysis of single-phase streamwise velocity for $Re_D$ = 160,000 with a horizontal upstream injector orientation. a,c,e) Spatial modes $\phi_{x,k}$. b,d,f) Welch's transform of the expansion coefficient $a_{x,k}$.</figcaption></figure>

the effect of the upstream injector position is expected to be limited to a minor shortening of the recirculation region.

## 3.2. SPOD and POD decomposition of single-phase flow

The topology of turbulent energetic flow structures downstream of the cylinder is evaluated using proper orthogonal decomposition (POD)

and spectral proper orthogonal decomposition (SPOD). A relevant background to these techniques is given in Section 2.3.4. These methods were applied on the single-phase streamwise velocity data to determine the modal shapes (and corresponding energies).

The $Re_D$=160,000 single-phase flow case with a horizontal injector position upstream (with only water co-flow activated) is used as a representative case to illustrate key mode shapes from the POD analysis. The POD analysis was conducted on the mean subtracted flow field. For
<table><caption>**Table 7**
POD modal energies for single-phase streamwise flow.</caption><thead><tr><th>Re<sub>D</sub></th><th>Injector Orient.</th><th>Cum. energy (mode 1-3)</th><th>Cum. energy (mode 1-10)</th></tr></thead><tbody><tr><td>100,000</td><td>Vert.</td><td>39.4%</td><td>57.8%</td></tr><tr><td>100,000</td><td>Horz.</td><td>43.3%</td><td>62.1%</td></tr><tr><td>160,000</td><td>Vert.</td><td>40.0%</td><td>58.2%</td></tr><tr><td>160,000</td><td>Horz.</td><td>43.6%</td><td>62.3%</td></tr></tbody></table>

this case, 62.3% of streamwise turbulent energy is contained in the first 10 modes; the cumulative energy contained in these first 10 modes is shown in Fig. 6d. The cumulative energy for all single-phase flow cases is summarized in Table 7.

The modal shape of the first three modes, together with the frequency spectra of the corresponding expansion coefficients, are shown in Fig. 7. Welch's method is applied to smooth the frequency spectra, using 256 sample segments (400 Hz sample rate), windowed with a Hamming window and 50% overlap.

The first two modes correspond to the traveling Karman shedding mode which is characterized by an asymmetrical distribution of the streamwise velocity along the wake centerline. The spectra are characterized by a broadband peak. This is in agreement with observations by Norberg (2003), who show that for $Re_D$ in the range of 100,000-160,000 the Strouhal number reaches a minimum value of around 0.186, with a maximum bandwidth. In our study, the broadband peak contains two fundamental frequencies, of median value $St = 0.185$. The primary shedding peak occurs at $St = 0.2$ while the secondary peak occurs at $St = 0.17$. The third POD mode contains a single frequency at $St = 0.03$, which is notably the difference between the Karman shedding peaks. By visual inspection of its modal shape, this mode corresponds to an extension/contraction of the recirculation region. Notably, the existence of a primary and secondary shedding peak as well as the low frequency flapping of the recirculation region is present for all of the studied cases.

The existence of a low frequency modulation in the shedding frequency has been reported in literature. Lehmkuhl et al. (2013) performed direct numerical simulation (DNS) of flow past a circular at $Re_D = 3900$, and identify low frequency fluctuations at $St = 0.0064$; these fluctuations correspond to a shrinking/expansion of the recirculation bubble. Experimentally, Miau et al. (Miau et al., 1999) studied the flow past a circular cylinder at $Re_D = 27,000$ and again found a low frequency fluctuation in the shedding frequency. They note that these low frequency fluctuations are present in both the vortex formation length and the base pressure. These effects are correlated, such that an extension/contraction of the vortex formation length corresponds to an increase/decrease in the base pressure. They further state that these low frequency fluctuations are one order of magnitude slower than the Karman vortex shedding frequency. This phenomena is not specific to a circular cylinder, as Miau et al. (Miau et al., 1999) demonstrate the same phenomena for a trapezoidal prism, while Watine et al. (Watine et al., 2023) have reported the phenomena in experiments using a cantilevered blunt plate aligned with the flow. According to Watine et al. (Watine et al., 2023), the primary vortex shedding mode is induced by the trailing edge instability, while the secondary shedding mode is induced by the separated boundary layer instability. The low frequency flapping of the recirculating region comes from the combination of these two instabilities. In the Reynolds number range of the present experiments, Palkin et al. (Palkin et al., 2016) applied unsteady LES and RANS with the second-moment (Reynolds stress) closure model to flow over a circular cylinder at $Re_D = 140000$. For both sets of simulations, they found that a low frequency modulation in the shedding frequency occurred at $St \approx 0.02$. This compares favorably to the $St \approx 0.03$ flapping frequency which is seen in the current set of experiments.

To discriminate between the primary and secondary shedding modes (as POD is a multi-frequency content approach), SPOD analysis

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** SPOD Mode 1 energies normalized by the cumulative energy of all modes and frequencies for single-phase flow in the close position. A comparison is given for the vertical and horizontal upstream injector configurations, showing (a) recirculation region oscillation; (b) secondary Karman shedding peak; (c) primary Karman shedding peak; (d) 2 x Karman shedding peak.</figcaption></figure>

was performed on the streamwise velocity data. The energy of mode 1 normalized by the total energy (sum of energies across all modes and frequencies), is shown as a function of Strouhal number in Fig. 8. To improve the clarity of Fig. 8, a 2 Hz wide moving mean filter is applied to smooth the data. Both the primary and secondary shedding peaks, at $St = 0.2$ and $St = 0.17$ respectively, show up clearly in the SPOD mode energies. The secondary shedding peak as well as the difference in the shedding peaks, corresponding to $St = 0.03$, show up for all cases although it is more prominent with the injector positioned horizontally upstream.

In addition to SPOD analysis on the immediate wake of the cylinder, SPOD analysis on the streamwise velocity was also conducted downstream. The modal shapes in the close and far positions for the Karman shedding modes, as well $St = 0.03$, are shown in Fig. 9. As with the POD analysis, the SPOD analysis corresponds to a horizontal upstream injector orientation.

Unlike in the POD analysis, which could not differentiate between the shape of the two Karman shedding modes, we can now compare the spatial structure of the $St = 0.17$ and $St = 0.2$ modes. In the region of $\frac{\chi}{D} < 3$ there are no obvious differences in the shape of the shedding modes, however, it is clear moving downstream that the $St = 0.17$ mode has a narrower wake than the $St = 0.2$ mode. This may be an indication of a change in the separation point in the secondary shedding mode. Furthermore, we can see that the $St = 0.03$ mode in the far position corresponds to an alternating speed up/slow down of the streamwise velocity along the wake centerline. As corroborated with the POD analysis, we can see that this difference in the two primary shedding modes appears to correspond to a shift in the strength and size of the recirculation region and can be considered as a first sign of intermittent reattached behavior of the boundary layer of the cylinder. For both $Re_D=100,000$ and $160,000$ it confirms the assertion of Norberg (Norberg, 2003) that there is some kind of premature transition to the critical regime in the vicinity of $Re_D \approx 160,000$.

To assess the influence of the injection of bubbles on the wake topology, we can now use SPOD and bubble image velocimetry analysis to assess how bubbles are arranged and move in the wake. Prior to doing so, the Stokes number of the injected bubbles is analyzed. We need to determine to what degree we expect that the bubbles will move independently from the underlying liquid phase.
<figure><img src="image_11.png" /><figcaption>(a) St = 0.17, Close Position</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) St = 0.2, Close Position</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(e) St = 0.03, Close Position</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) St = 0.17, Far Position</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) St = 0.2, Far Position</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) St = 0.03, Far Position</figcaption></figure>

<caption>Fig. 9. Select mode 1 shapes at different Strouhal numbers from SPOD analysis of single-phase streamwise velocity for Re_D = 160,000 (with a horizontal upstream injector orientation). The colors correspond to normalized mode intensity.</caption>

### 3.3. Stokes number analysis

To determine the fluid timescales to which the bubbles will respond, it is necessary to estimate the Stokes number. The Stokes number is defined as the ratio of the bubble response time to a characteristic fluid time-scale, $t_{Flow}$, as follows:

$$Stokes_{Flow} = \frac{t_{Bubble}}{t_{Flow}}$$

$$ = \left[ (\rho_{Bubble} + 0.5\rho_{Water}) \frac{d_{Bubble}^2}{18\mu_{Water}} \right] t_{Flow}^{-1} $$

Where $\rho$ is the density of the fluid, $\mu$ is the dynamic viscosity, and $d_{Bubble}$ is the bubble diameter. The bubble response time is computed assuming that the bubble is spherical with an added mass coefficient of 0.5, and the dominant force governing the dynamics is drag. It is challenging to define a single fluid time-scale, in part because the fluid velocity changes substantially when comparing the free-stream velocity
<table><caption>**Table 8**
Stokes number statistics.</caption><thead><tr><th rowspan="2">Re<sub>D</sub></th><th rowspan="2">Injector Orient.</th><th rowspan="2">β[%]</th><th colspan="3">Stokes number, Stokes<sub>Flow</sub></th></tr><tr><th>Median</th><th>50% Interval</th><th>90% interval</th></tr></thead><tbody><tr><td>100,000</td><td>Vert.</td><td>0.010</td><td>1.4</td><td>0.6–3.0</td><td>0.2–6.8</td></tr><tr><td>100,000</td><td>Vert.</td><td>0.019</td><td>2.3</td><td>0.9–4.8</td><td>0.2–8.5</td></tr><tr><td>100,000</td><td>Vert.</td><td>0.029</td><td>2.9</td><td>1.2–5.5</td><td>0.3–8.8</td></tr><tr><td>160,000</td><td>Vert.</td><td>0.006</td><td>2.0</td><td>0.9–3.6</td><td>0.2–7.4</td></tr><tr><td>160,000</td><td>Vert.</td><td>0.012</td><td>2.7</td><td>1.3–5.1</td><td>0.3–8.6</td></tr><tr><td>160,000</td><td>Vert.</td><td>0.018</td><td>3.2</td><td>1.5–5.7</td><td>0.4–9.0</td></tr></tbody></table>

to the wake. Consequently, the spatially varying time scale of the fluid was chosen as $t_{Flow} = \frac{D}{|u_{x,y,t}|}$, with $|u_{x,y,t}|$ being the instantaneous velocity magnitude at time $t$ and a particular $(x, y)$ location. This is the time taken for the local velocity to move one diameter. This length-scale is chosen as the cylinder diameter is the same order of magnitude as several prominent features in the wake, such as: (1) spanwise cell size, (2) vortex diameter, (3) wake width, all of which will strongly influence the motion of the bubble in the wake. The actual bubble motion within the wake is highly nonlinear and the prediction of individual bubble motion in the wake is the subject of ongoing research by the authors.

Since we do not have access to the water velocities in the two-phase flow, the instantaneous single-phase velocity magnitude data, coupled with volume-weighted upstream bubble size histograms (based on data given in Fig. 2), is used to compute a histogram of Stokes numbers. Specifically, a histogram is produced by aggregated the spatial Stokes number data across the region $0 \le \frac{x}{D} \le 2$, $-1 \le \frac{y}{D} \le 1$ as well as in time. This is typically the region where the flapping of the recirculation region takes place in the single-phase flow. Implicit in this analysis is the assumption that bubbles of all sizes respond equally to all fluid velocities, and while this may not be true in practice it is sufficient to produce an estimate of the order of magnitude of the Stokes number.

The Stokes number is dependent on the upstream bubble diameter histogram, which varies with the gas injection rate. The Stokes number data is contained in Table 8. The 50% interval is the Stokes number range which spans the 25th to the 75th percentile of the total volume of gas, while the 90% interval spans the 5th to the 95th percentile. Defining the percentiles by the gas volume rather than by number density accounts for the greater volume contribution of larger (higher Stokes number) bubbles. For brevity, only data computed using the vertical injector orientation is used; since the velocity data is from the single-phase flow, the injector orientation does not have a substantial influence on the estimated Stokes number distribution.

From Table 8, it is clear that most bubbles have Stokes numbers of $O(1)$, when computed using single-phase velocity magnitudes. Consequently, on average we expect most bubbles to take $O(1D)$ to accelerate to the liquid velocity in the wake; however, from the 90% interval we see that bubbles span nearly two orders of magnitude of Stokes number. Therefore, we expect that within the near wake, in the recirculating region, there are going to be bubbles which more closely act as flow tracers in addition to inertia dominated bubbles which may have significant slip velocities.

Other relevant fluid time-scales include the Karman shedding period, which is important in determining whether bubbles can respond to individual vortices being shed. To assess this response, we instead compute the Stokes number, Stokes$_{KS}$, based on the primary Karman shedding frequency of the single phase flow, $St = 0.2$. In this case $t_{Flow} = \frac{D}{0.2U}$. For the $Re_D = 100,000$ case, the median value of Stokes$_{KS}$ varies from 0.7–1.5, depending on the gas flow rate; for the $Re_D = 160,000$ case, the median value of Stokes$_{KS}$ varies from 0.9–1.6. The range of computed Stokes numbers implies the following:

1. While the Stokes number varies with gas flow rate, in all cases there is a significant fraction of bubbles in the wake with Stokes$_{Flow} < 1$ which means that bubbles have sufficiently fast response times to respond to most wake time-scales. In particular, the median value of Stokes$_{KS}$ ≈ 1 means that bubbles can respond to, and potentially be entrapped by, the Karman vortices.

2. If both the primary and secondary Karman shedding frequencies are present in the underlying liquid flow, it is highly unlikely that bubbles will only be trapped by the lower of these two frequencies.

With the above comments in mind, we can now progress to the analysis of the two-phase data. Since bubbles will respond to most time-scales present in the near-wake, including the vortex shedding time-scales, we can study the aggregation and motion of bubbles to provide some insight into the underlying flow.

### 3.4. Time-averaged topology of the gas phase

The distribution in the close plane of the time-averaged bubble intensity is shown in Fig. 10. To compare the topology of the bubble trapping in different conditions, the distribution is normalized for each case. To reduce outlier effects, the minimum and maximum values correspond to the 2nd and 98th percentiles of the raw image intensity, respectively. Since the images are front lit, and the light from individual bubble scattering is nonlinear with bubble size and position, we cannot draw a direct quantitative relationship between the image intensity and the number of bubbles present in a particular region. However, we can make qualitative arguments, specifically, that higher intensity areas of the time-averaged image correspond to time-averaged regions of gas accumulation.

From Fig. 10, we can see that gas tends to collect in the near wake, preferentially in the recirculating region, up to a distance less than $x = 1.5D$ downstream. The circular pattern that extends to roughly $x = 0.5D$ downstream is simply light reflected off of the far wall of the channel, for which there is a circular mount through which the cylinder passes. Furthermore, the bright spot located at approximately $x = -0.3D$, $y = 0.7D$ that is present in 3/4 images is not a function of collected gas in the channel and should be ignored. For the case of the highest gas flow rate ($Q_g = 5100$ mL/min), the vertical width of the region of locally high void fraction is marginally wider but otherwise there is no strong difference in the normalized bubble intensity seen between the different air injection rates. The fluctuation of gas content within the wake is determined using SPOD analysis.

### 3.5. SPOD decomposition of two-phase flow

Prior to presenting the SPOD mode shapes of the bubble intensity, it is necessary to briefly describe the meaning attributed to these modes. Discussion is provided in Section 2.3.4 but key points are re-iterated here for clarity. The intensity of the image increases in regions of higher bubble density, which allows us to draw qualitative information on the bubble accumulation in the wake. One must be careful when comparing SPOD spatial modes from PIV of the streamwise velocity of the single-phase flow and bubble intensity directly, as the single-phase modal shapes correspond to regions of fluctuating velocity, while modal shapes in the two-phase case correspond to regions of fluctuating bubble concentration. Therefore, while coherent regions of bubble concentration are expected to align with coherent flow structures in the underlying liquid velocity, several main differences must be highlighted:

1. As discussed in Section 3.3, the bubbles in the present experiments are not perfect flow tracers. Consequently, the bubbles can slip relative to the underlying flow and will not necessarily move with coherent structures present in the underlying liquid flow.
<figure><img src="image_17.png" /><figcaption>(a) $Re_D = 100,000; \beta = 0.010\%$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) $Re_D = 100,000; \beta = 0.029\%$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c) $Re_D = 160,000; \beta = 0.006\%$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) $Re_D = 160,000; \beta = 0.018\%$</figcaption></figure>

**Fig. 10.** Time-averaged bubble intensity for selected test conditions, with a vertical upstream injector orientation. The color indicates the normalized image intensity. The gray area is in shadow due to the positioning of the spot lights. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2. The depth of focus of the images, while centered at $y = 0$, is on the order of the diameter of the cylinder. This means that the bubble accumulation information is *depth-averaged*. This is in contrast to modal analysis of the PIV data, for which the modal shapes correspond to centerline, $z = 0$, velocity structures.

3. In the case of the PIV data, coherent structures can be directly identified anywhere within the image. In contrast, the presence of coherent structures in the bubble images relies on time- and spatially-dependent bubble migration to the corresponding structures. Consequently, coherent structures downstream are more likely to be highlighted within the two-phase images. The preferential collection of bubbles within the flow is also dependent on the bubble radius, which is an area of study that warrants further attention and is subject to active research by the authors.

Given the time-dependent capture of bubbles in the shed vortices, we expect that coherent structures will be shown more clearly in the far position, as bubbles have time to coalesce within these structures. For the case of $Re_D = 160,000$ with vertical injection, the SPOD mode 1 eigenvalue normalized by the sum of eigenvalues across all modes and frequencies is given in Fig. 11(a); the ratio of the mode 1 to mode 2 eigenvalues is given in Fig. 11(b). According to Schmidt and Colonius (Schmidt and Colonius, 2020), the first SPOD mode is particularly representative of a dominant physical process at a given frequency if the ratio between the first and second mode eigenvalue

is large. The single-phase Karman peaks of $St = 0.174$ and $St = 0.2$ are indicated by vertical lines. To improve the clarity of the figures, a 2 Hz wide moving mean filter is applied to smooth the data. In cases in which the frequency resolution is greater than 2 Hz, no correction is applied. Recall, the 100 Hz peak which corresponds to the frequency of the lights used to illuminate the two-phase flow is removed from the normalized eigenvalue spectra.

From Fig. 11, it is possible to visualize the amplitude of fluctuations in the bubble intensity at different Strouhal numbers. It is immediately apparent that rather than having two Karman shedding peaks, there is a single dominant peak centered at $St = 0.18$. Interestingly, the ratio of mode 1 to mode 2 eigenvalues does not indicate a physical feature at $St = 0.2$, which means that the Karman shedding mode at $St = 0.2$ either (a) does not exist in the two-phase flow or (b) has no imprint on the bubble trapping. However, as discussed in Section 3.3, option (b) is highly unlikely given the Stokes number of the injected bubbles. For vertical injection, the SPOD mode 1 eigenvalue contained in the Karman vortex shedding bandwidth is amplified in both the close and far positions, with increasing air injection rate. Furthermore, the spectral bandwidth of the Karman shedding peak is reduced when shifting from the close to far position. The modal shapes of SPOD mode 1 at $St = 0.18$ for two volumetric qualities, $\beta=0.006\%$ and $0.018\%$, are shown in Fig. 12. The spatial modes are virtually identical for both gas flow rates. These spatial modes highlight coherent regions of alternating bubble accumulation and release along the streamwise direction, as expected for bubbles being trapped by Karman vortices.
<figure><img src="image_21.png" /><figcaption>Fig. 11. Normalized eigenvalue spectra of the primary SPOD modes of the bubble intensity for Re_D=160,000 with a vertical upstream injector orientation. The primary Strouhal numbers of the single-phase flow are indicated by vertical lines (St = 0.174 and St = 0.2). Values of β are reported as percentages.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 12. Spatial distribution of the bubble intensity at St = 0.18, extracted from SPOD mode 1, for Re_D = 160,000 with a vertical upstream injector orientation.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 13. Spatial distribution of the bubble intensity at selected low frequencies, extracted from SPOD mode 1, for case of Re_D = 160,000 with vertical injection.</figcaption></figure>

Returning to the ratio of the SPOD mode 1 and mode 2 eigenvalues shown in Fig. 11b, the low frequency content in the range $0 \le \text{St} \le 0.03$ does not persist to the far position. This is in contrast to the single phase SPOD analysis, for which low frequency fluctuations in the streamwise velocity persist downstream. This implies that the gas that is accumulated/ejected in the near wake at these low frequencies does not move coherently downstream.

In Fig. 13, we show the SPOD mode 1 bubble intensity distribution corresponding to $\text{St} = 0.01$ and $\text{St} = 0.03$. These modes are shown for the $Re_D = 160,000$ vertical injection condition, with volumetric qualities of $\beta=0.006\%$ and $0.018\%$. These Strouhal numbers contain high eigenvalues in the close position (see Fig. 11b for the $Re_D = 160,000$ case). We can see that these modes correspond to a filling/ejection of gas in the recirculation region behind the cylinder.

The influence of switching the bubble injector from the vertical to the horizontal orientation on the mode 1 SPOD eigenvalue of the image intensity is shown in Figs. 14 and 15. For brevity, and given the reduction in spectral bandwidth of the Karman shedding peak at the far position, only the far position is shown. Again, to improve the clarity of the figures, a 2 Hz wide moving mean filter is applied to smooth the data. From these spectra, it is clear that when switching from the vertical to horizontal bubble injector orientation, the spectral bandwidth of the Karman shedding peak increases. Since an increase in gas flux does not appear to substantially change the width of the spectral peak (in a given injector orientation), it is likely that this change is due to the redistribution of gas along the span.

Although the volumetric quality is higher for $Re_D = 100,000$ than for $Re_D = 160,000$, when comparing the normalized mode 1 SPOD eigenvalues within the vortex shedding bandwidth it is clear that the Karman shedding peaks are less prominent in the $Re_D = 100,000$ case. A proposed explanation for this effect is that the strength of the vortices is smaller in the case of the $Re_D = 100,000$ flow; consequently, bubbles are less likely to be pulled into the shed vortices. In this manner, it is clear that the influence of bubbles on the flow past the cylinder is also a function of Reynolds number, since Reynolds number will in part determine the efficiency of bubble trapping in the wake.

From the SPOD analysis of the bubble images, we can assess the frequency at which gas content gathers and is ejected. With this in mind, several effects are apparent:

1. The gas gathers close to the cylinder at broadband low frequencies, but these frequencies do not persist in a spatially coherent manner downstream.

2. The gas is alternatively ejected and gathered at a single frequency peak in the range of $\text{St}=0.17-0.18$, which corresponds to the Karman shedding mode. This differs from the dual Karman shedding modes seen in the single-phase case.

3. There is no preferential gas accumulation/ejection at $\text{St}=0.2$, which when considering the Stokes number of the injected bubbles, indicates that the $\text{St}=0.2$ Karman shedding mode does not exist in the two-phase flow.

With the above comments in mind, it is clear that the flow is two way coupled. Given this understanding, we can now use bubble image velocimetry to estimate the gas phase velocity, and in doing so gain insight into the bubble influence on the flow topology.

### 3.6. Bubble image velocimetry

In an effort to estimate the gas phase velocity, bubble image velocimetry is used to compute the Eulerian motion of bubbles in the wake. As described in Section 3.3, we see that bubbles span nearly two orders of magnitude of Stokes number. Therefore, we expect that within the near wake, in the recirculating region, there are going to be bubbles which more closely act as flow tracers in addition to inertia dominated bubbles which may have significant slip velocities. Therefore, we do not attempt to use bubbles as flow tracers of the liquid velocity, instead, the BIV analysis is solely limited to estimating the gas phase velocity. Interrogation windows of 64 × 64 pixels are used, with the cross-correlation used to determine the Eulerian velocity of the window. In the far field, the images are backlit which allowed for improved cross-correlation accuracy; consequently, 32 × 32 interrogation windows could be used. There is no overlap in the interrogation windows. The bubble image velocimetry analysis is only performed with images from the vertical injector orientation. The settings for BIV acquisition are summarized in Table 4.

To ensure that the Eulerian analysis is computed based on a representative size distribution of bubbles (e.g. not just bubbles towards the larger end of the spectrum), the bubble size distribution in the wake is compared to the size distribution upstream. Importantly, the wake size histogram is constructed using the bubbles from each interrogation window that are most representative of the computed Eulerian velocity vector. To determine the bubbles which are representative of the computed velocity vector within each interrogation window, the preceding image is shifted according to the computed velocity vector, and only bubbles which overlap in the second image are kept. In this manner, bubbles which do not have a match in the second image (and hence do not contribute to the calculation of the Eulerian velocity vector) are not added to the wake bubble size histogram.

The bubble size distribution measured with BIV in the region $0 < \frac{x}{D} < 2$, $-1 < \frac{y}{D} < 1$ is shown in Fig. 16. The BIV volume-weighted bubble size distribution compares favorably with the histogram upstream,
<figure><img src="image_24.png" /><figcaption>**Fig. 14.** Influence of the volumetric quality (reported as percentages) on the normalized eigenvalue spectra of the primary SPOD modes of the bubble intensity in the far position for case of $Re_D$ = 100,000.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 15.** Influence of the volumetric quality (reported as percentages) on the normalized eigenvalue spectra of the primary SPOD modes of the bubble intensity in the far position for case of $Re_D$ = 160,000.</figcaption></figure>

indicating that the BIV analysis is not biased towards a particular range of bubble sizes. It is not clear whether the presence of slightly more large bubbles present in the wake data (when compared to upstream) may be due to a preferential detection of these bubbles in the BIV analysis, due to the tendency of bubbles to combine and form larger bubbles due to interaction in the wake, or due to the preferential trapping of larger bubbles in the near wake.

To capture the range of velocities present in the wake, and compare those to the single-phase velocity data, we present the median, 25th, and 75th percentile velocity data. The percentiles are computed in a time varying sense, such that velocity is slower than the 25th percentile data 25% of the time. Fig. 17 shows the vertical profiles of streamwise velocity plotted as function of the streamwise location, among the different test cases for the vertical injector. Only the upper profile is plotted as the shadow below the cylinder disrupts the BIV measurements. However, the profiles are largely symmetric about the y = 0 axis.

From Fig. 17, there are several important features in the near-wake:

1. There is a wake velocity deficit, for which the gas phase (BIV) measurements consistently show higher streamwise velocities in the 75th percentile, but consistently show similar velocity profiles in the 25th percentile.

<figure><img src="image_26.png" /><figcaption>**Fig. 16.** Volume weighted histogram of bubble sizes used to compute BIV data for case of $Re_D$ = 160,000 and $\beta$=0.006%.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 17. Offset profiles of streamwise velocity, for which the following percentiles of velocity (in time) are presented: (a) 75th percentile; (b) 50th percentile (median); (c) 25th percentile.</figcaption></figure>

2. The median gas phase (BIV) streamwise velocity tends to exceed that of the single-phase flow near the wake centerline.

3. The BIV data does not contain higher streamwise velocities in the separated shear layer, $0.7 \le \frac{y}{D} \le 1$, that are present for the single phase case.

4. There is no obvious influence of the air injection rate except close to the wake centerline, where an increase in the air injection rate leads to a reduction of the gas phase streamwise velocity deficit.

This wake deficit continues in the far wake, as shown in Fig. 18 for the $Q_g = 1700$ mL/min case. From the two-phase BIV measurements, it is clear that the streamwise gas velocity has largely recovered to the freestream value, while the single-phase flow still experiences a wake deficit.

Importantly, the streamwise velocity does not tell the full story. The vertical velocity at different y locations is shown in Fig. 19.

From Fig. 19, it is clear that there is a large negative vertical velocity in the region $\frac{x}{D} \approx 1.5$, $\frac{y}{D} \approx 0.5$, which is not present in the single-phase case. This modified vertical velocity is largely consistent across gas flow rates, with the exception of $\frac{y}{D} = 0.75$, but is consistently larger for the case of $U = 8$ m/s. At a height of $\frac{y}{D} = 1$ this vertical velocity has largely disappeared. Beyond the recirculation region, at $\frac{x}{D} \approx 1.5$, the gas phase velocity is typically more strongly positive than the single-phase flow.

To illustrate in a qualitative manner the momentum injection due to the presence of bubbles, an estimate of the change in velocity is overlaid on the time-averaged bubble intensity image for the case

<figure><img src="image_28.png" /><figcaption>Fig. 18. Offset profiles of streamwise velocity for the $\frac{x}{D} = 10$ location, for which the following percentiles of velocity (in time) are presented: (a) 75th percentile; (b) 50th percentile (median); (c) 25th percentile.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>(a) $\frac{y}{D} = 0.25$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) $\frac{y}{D} = 0.5$</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(c) $\frac{y}{D} = 0.75$</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(d) $\frac{y}{D} = 1$</figcaption></figure>

<caption>**Fig. 19.** Offset profiles of vertical velocity, for which the following percentiles of velocity (in time) are presented: (a) 75th percentile; (b) 50th percentile (median); (c) 25th percentile.</caption>

of $Re_D=160,000$ and $\beta=0.006\%$. The velocity difference is computed using the median gas velocity (calculated with BIV) and median single-phase velocity (calculated with PIV). Note that this estimate of the velocity difference compares the velocity of the bubbles in the two-phase case to the water velocity of the single-phase data. We do not attempt to comment on the slip velocity between the bubbles and underlying water flow in the two phase case, although the degree to which these two velocities match is addressed in Section 3.3.

As can be seen from Fig. 20, it is clear that the region of largest bubble concentration (in the near wake) is produced by the vertical motion of gas relative to the single-phase liquid velocity. The largest vertical velocity difference is just beyond the region of largest bubble concentration (in the near wake), at $\frac{y}{D} \approx 0.5$. This momentum injection as gas moves across the separated shear layer likely serves to disrupt the regular roll in the near wake behind the cylinder, and shifts the roll up of the shear layer further downstream.

To illustrate the influence of the bubbles on the vibration of the cylinder, vibrometer measurements are presented.

### 3.7. Cylinder vibration results

The results from the vibrometer are shown in Fig. 21 for $Re_D = 160,000$. The vibrometer measurement records the pointwise velocity of the cylinder. The acquisition settings are described in Section 2.2.

From Fig. 21, it is apparent that the injection of bubbles into the flow substantially reduces the influence of the shedding frequency peak

<figure><img src="image_33.png" /><figcaption>**Fig. 20.** Velocity difference between the gas phase velocity (estimated using BIV) and the single-phase flow (computed using PIV) for the case of $Re_D=160,000$ and $\beta=0.006\%$. The velocity difference is overlaid on the time-averaged bubble intensity image.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Fig. 21.** Frequency content of the cylinder vibration, for $Re_D$ = 160,000. $V_{cyl}$ corresponds to the pointwise velocity of the cylinder at the free end.</figcaption></figure>

on the cylinder vibration. This effect is irrespective of the effective void fraction increase when switching from vertical to horizontal injection. The shift in the natural frequency peak was likely due to gas collecting in the mounting mechanism for the cylinder, rather than a physical shift in the natural frequency of the cylinder. However, the successful measurement of the natural frequency peak is an indication that the vibrometer was able to track the cylinder motion despite the presence of gas in the channel. The same trends were observed for $Re_D$ = 100,000.

## 4. Discussion

From the previous results, we can conclude that the flow is two-way coupled. The presence of gas in the wake produces a series of effects:

1. Suppression of vortex-induced vibration of the cylinder. This suppression is irrespective of whether gas interacts with the cylinder along the entire span or simply near the centerline.

2. Increase in the centerline streamwise velocity (promoting faster wake recovery).

3. A shift in the shedding frequency to be closer to the lower of two shedding peaks in the single-phase signal (which is the much less energetic secondary peak).

4. With an increase in the gas flux, keeping the injector orientation consistent, there is a larger change in the gas velocity from the single-phase liquid velocity.

5. With an increase in the spanwise distribution of gas, by moving from a vertical to horizontal gas injector orientation, there is a broadening of the spectral shedding peak.

It is difficult to determine the exact cause of each of these effects as we do not have access to the underlying liquid velocity in the two-phase flow. With that said, potential mechanisms can be proposed.

Given that the flow is within the lower critical regime, an increase in the upstream turbulence intensity can lead to a significant shift in the shedding frequency, particularly if there is a premature shift to the critical regime (Norberg, 1987). There is a significant body of work on bubble-induced turbulence, such as that of Ma et al. (2021) and summarized by Mathai et al. (2020). It is worth discussing the influence of the two-phase flow on the carrier phase turbulence in the present work. Bubbles are commonly considered to increase the turbulence level through their random motion and the contribution of their turbulent wakes. This trend is expected to scale with increasing void fraction, $\alpha$. Mathai et al. (2020) present a comprehensive review of the different mechanisms that are present in bubble-induced turbulence. Recent

studies have shown that the turbulence intensity can be decreased by the presence of bubbles at low void fraction, according to the value of the bubblance parameter. The bubblance parameter, $b$, compares the intensity of the bubble induced turbulence with the turbulence intensity of the single-phase flow. The bubblance parameter, $b$, is given by $b = \alpha V_r^2/u'^2$, where $V_r$ is the bubble slip velocity, and $u'$ is the rms value of the single-phase fluctuating velocity component (Alméras et al., 2017). Alméras et al. (2017) demonstrated that for homogeneous isotropic turbulence, in the range of $b < 0.27$, there is an attenuation in the turbulence intensity. However, Alméras et al. did not characterize the bubble-induced turbulence for values of $b$ below 0.13. For the present work, the bubblance parameter does not exceed 0.017 (obtained for the maximum volumetric quality of $\beta$ = 0.14%). This was obtained using shadowgraphy in the measurement section in the absence of the cylinder. This value is much lower than that studied in existing literature (driven by the low volumetric qualities in the present work), implying that the role of bubble-induced turbulence is small. Moreover, previous research by Watanabe et al. (1990), has shown that at void fractions below 0.3%, as is the case in the current work, the bubble induced turbulence is negligible. Therefore, we can conclude that the modification of the turbulence by bubble motion is not expected to play a role in the observed decrease in the vortex shedding frequency.

We can turn our attention to the influence of bubble injection in the near wake. From the significant change in the vertical velocity between the single-phase liquid velocity and the gas velocity, in conjunction with a region of high gas content in the near wake (implying a comparatively higher momentum injection), a mechanism for the restructuring of the wake streamwise velocity may be proposed. Namely, the separated shear layer in the wake is broken up, which eliminates the higher velocity separated flow region and allows more freestream flow to enter the wake. Consequently, while the slow regions of the flow are still present, there is greater mixing with the freestream flow. This promotes faster wake recovery. Gas injection into the wake also may lead to a delayed formation of the vortex street, as regular roll up of the shear layer is pushed further downstream. The stochastic nature of the high velocity vertical motion may insulate the cylinder from vibrating at a regular shedding frequency.

It is also important to note that the primary shedding frequency peak in the two-phase data was shifted to $St = 0.17-0.18$. While it is challenging to determine with the existing data why the lower of the two Karman shedding peaks were selected, a potential explanation may be a shift in separation point, which is not detectable with the present set of experiments. Of note is that the narrow near wake compares well with the lower of the two shedding peaks in the single-phase data, which implies that if a separation point shift was the cause of the two different Karman shedding frequencies present in the single-phase data, then the presence of bubbles in the wake may serve to stabilize this separation point. Indeed, very few bubbles can interact with the boundary layer and modify its nature (Murai, 2014). Murai (2014) also notes that the motion of bubbles relative to the underlying liquid flow can contribute to a reduction in shear stress. A reduction in the turbulent shear stress in the wake is in agreement with a reduction of the streamwise velocity deficit.

Additionally, the spectral bandwidth of the Karman shedding peaks increased with the volumetric quality $\beta$, it is not clear if this increase is simply due to the presence of gas along the span, or if a similar effect could be produced by introducing an equivalent spanwise-averaged volumetric quality with higher gas content locally near the $z = 0$ centerline. However, the broadly similar spectral width with increasing gas flux in a given injector orientation suggests that this effect is likely driven by the redistribution of gas along the span, rather than the increase in volumetric quality.
## 5. Conclusions and future work

The effect of bubbles on the structure of the wake and vortex induced vibration of the cylinder was examined over two values of Reynolds number, $Re_D = 100,000$ and $160,000$, both of which fall within the subcritical regime. The vibration of the cylinder, PIV, BIV and SPOD/POD of the bubble intensity were used to assess the flow topology changes under the influence of gas injection. Using SPOD/POD analysis in the near wake, it was found that the primary shedding frequency decreased with the injection of gas, from a Strouhal number of $St = 0.2$ to $St = 0.17-0.18$; the width of the spectral peak was found to increase with void fraction. The reduction in the shedding frequency, as well as elimination of the secondary shedding peak, may indicate that the gas injection stabilizes the separation point. Furthermore, the vibration of the cylinder at the primary Karman shedding frequency of $St = 0.2$ was significantly attenuated following the injection of gas, even at volumetric qualities below 0.01%. It appears that stochastic gas accumulation/release in the near wake, driven by high velocity vertical injection events of bubbles into the near wake, serve to uncouple the cylinder motion from the formation of the vortex street downstream. While not possible with the present experimental setup, simultaneous measurement of the gas and liquid velocities through techniques such as particle shadow velocimetry (e.g. such as in Ma et al. (2022) may provide additional insight into the two-way coupling between the gas and liquid phases.

The observed reduction in the shedding frequency was in contrast to the increase reported by the majority of previous literature, when studied elsewhere in the parameter space spanned by Reynolds number, bubble size, and void fraction. The highly nonlinear influence of these parameters on the wake topology is chronically under-explored, and warrants further attention.

In addition to changes in the Strouhal number, it was shown that the Karman shedding peak in the SPOD analysis of the bubble intensity images is less prominent for $Re_D = 100,000$ when compared to the $Re_D = 160,000$ case. It is proposed that this is due to the reduction in efficiency of bubble trapping within these structures at the lower Reynolds number, since the vortices will be less energetic. At sufficiently different Reynolds numbers, the efficiency of bubble trapping by the Karman vortices may be substantially altered in a manner that changes the way in which bubbles alter the underlying liquid flow field. The Reynolds number dependent motion of bubbles of different sizes within the wake of a cylinder (in which individual bubbles are tracked in 3D and overlaid on simultaneous PIV measurements) is the subject of ongoing work by the authors. Furthermore, a companion paper is being prepared which explores the influence of bubble injection on super-critical Reynolds number flows.

Further exploration of the influence of bubbles on the separation point must also be conducted, to better characterize the mechanism by which the Strouhal number is reduced. In particular, the influence of higher void fraction on the flow topology in the wake must be explored, in a manner which will connect the current low void fraction experiments with previous literature that has studied void fractions from $O(1\%)$ to $O(10\%)$. These future experiments are necessary to further explore the parameter range spanned by bubble size, Reynolds number, and void fraction beyond what has been studied to date. With these future studies, a coherent detailing of the mechanism by which bubbles alter the flow around a cylinder can finally be constructed.

### CRediT authorship contribution statement

**Eric Thacher:** Writing – original draft, Visualization, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Bruno Van Ruymbeke:** Visualization, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Céline Gabillet:** Writing – review & editing, Supervision, Software, Project administration, Methodology, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Nicolas

**Jacques:** Writing – review & editing, Funding acquisition. **Simo A. Mäkiharju:** Writing – review & editing, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

The authors gratefully acknowledge the support of National Science Foundation award 2024672, program managers Ron Joslin and Shahab Shojaei-Zadeh. This work was also supported by the ISBlue project, Interdisciplinary graduate school for the blue planet (ANR-17-EURE-0015) and co-funded by a grant from the French government under the program "Investissements d'Avenir" embedded in France 2030.

### References

Achenbach, H., Heinecke, E., 1981. On vortex shedding from smooth and rough cylinders in the range of Reynolds numbers $6 \times 10^3$ to $5 \times 10^6$. J. Fluid Mech. 109, 239-251.

Alméras, Elise, et al., 2017. Experimental investigation of the turbulence induced by a bubble swarm rising within incident turbulence. J. Fluid Mech. 825, 1091-1112.

Ausoni, Philippe, 2009. Turbulent Vortex Shedding from a Blunt Trailing Edge Hydrofoil. Tech. Rep., EPFL.

Bearman, P.W., 1965. Investigation of the flow behind a two-dimensional model with a blunt trailing edge and fitted with splitter plates. J. Fluid Mech. 21 (2), 241-255.

Bearman, P.W., 1984. Vortex shedding from oscillating bluff bodies. Annu. Rev. Fluid Mech. 16, 195-222.

Bearman, P.W., 2011. Circular cylinder wakes and vortex-induced vibrations. J. Fluids Struct. 27, 648-658.

Callison, Elizabeth, et al., 2023. Balanced multiphase mixing through a narrow gap. Int. J. Multiph. Flow 165, 104481.

Ceccio, Steven L., 2010. Friction drag reduction of external flows with bubble and gas injection. Annu. Rev. Fluid Mech. 42, 183-203.

Cheng, Wen, et al., 2005. Bubble velocity measurement with a recursive cross correlation PIV technique. Flow Meas. Instrum. 16 (1), 35-46.

Cheremisinoff, Nicholas P. (Ed.), 1986. Encyclopedia of Fluid Mechanics, vol. 3, Gulf Publishing Company, Houston, Texas, pp. 771-810 (Chapter 27).

Cook, Thomas L., Harlow, Francis H., 1986. Vortices in bubbly two-phase flow. Int. J. Multiph. Flow 12 (1), 35-61.

Gabbai, R.D., Benaroya, H., 2005. An overview of modeling and experiments of vortex-induced vibration of circular cylinders. J. Sound Vib. 282, 575-616.

Habeeb, L.J., Al-Turaihi, R.S., 2013. Experimental study and CFD simulation of two-phase flow around multi-shape obstacles in enlarging channel. Am. J. Mech. Eng. 1 (8), 470-486.

Hara, F., 1984. Air-bubble effects on vortex-induced vibrations of a circular cylinder. In: Paidoussis, M.P., Griffin, O.M., Sevik, M. (Eds.), Excitation and Vibration of Bluff Bodies in Cross Flow. In: ASME Winter Annual Meeting, vol. 1, ASME, pp. 103-113.

Hara, F., 1987. Vibrations of circular cylindrical structures subjected to two-phase cross flows. JSME Int. J. 30 (263), 711-722.

Hara, F., Iijima, T., 1989. Vibrations of two circular cylinders in tandem subjected to two-phase bubble cross flows. J. Fluids Struct. 3, 389-404.

Hulin, J.-P., Fierfort, C., Coudol, R., 1982. Experimental study of vortex emission behind bluff obstacles in a gas liquid vertical two-phase flow. Int. J. Multiph. Flow 8 (5), 475-490.

Iijima, T., Hara, F., Nishikubo, K., 1995. Unsteady fluid force acting on two tandem circular cylinders subjected to two-phase cross flow. JSME Int. J. Ser. C 38 (2), 219-226.

Inoue, A., et al., 1986. Studies on two-phase cross flow. Part I: Flow characteristics around a cylinder. Int. J. Multiph. Flow 12 (2), 149-167.

Joo, Y., Dhir, V.K., 1994. An experimental study of drag on a single tube and on a tube in an array under two-phase cross flow. Int. J. Multiph. Flow 20 (6), 1009-1019.

Kaneko, S., et al. (Eds.), 2014. Flow-Induced Vibrations: Classifications and Lessons from Practical Experiences, second ed. Elsevier.

Kim, Dohwan, Rau, Matthew J., 2020. Two-phase wakes in adiabatic liquid-gas flow around a cylinder. In: Fluids Engineering Division Summer Meeting, Vol. 83723. American Society of Mechanical Engineers, p. V002T04A034.
Kundu, Pijush K., Cohen, Ira M., Dowling, David R., 2016. Fluid Mechanics. Academic Press.

Lee, Jubeom, 2021. Modeling and Experimental Validation of the Bubble-Induced Reynolds Stress in the Homogeneous Bubble-Swarm Past a Circular Cylinder (Ph.D. thesis). Seoul National University.

Lee, Jubeom, Park, Hyungmin, 2020. Bubble dynamics and bubble-induced agitation in the homogeneous bubble-swarm past a circular cylinder at small to moderate void fractions. Phys. Rev. Fluids 5 (5), 054304.

Lehmkuhl, O., et al., 2013. Low-frequency unsteadiness in the vortex formation region of a circular cylinder. Phys. Fluids 25 (8).

Lienhard, J.H., 1966. Synopsis of Lift, Drag, and Vortex Frequency Data for Rigid Circular Cylinders. Tech. Rep. Bulletin No. 300, Washington State University, College of Engineering.

Ma, Tian, et al., 2021. Scale-dependent anisotropy, energy transfer and intermittency in bubble-laden turbulent flows. J. Fluid Mech. 927, A16.

Ma, Tian, et al., 2022. An experimental study on the multiscale properties of turbulence in bubble-laden flows. J. Fluid Mech. 936, A42.

Mathai, Varghese, Lohse, Detlef, Sun, Chao, 2020. Bubbly and buoyant particle-laden turbulent flows. Annu. Rev. Condens. Matter Phys. 11, 529-559.

Meng, Haoran, 1993. On Dispersed Two Phase Flows Past Obstacles. Technische Universiteit Eindhoven.

Miau, Jiun-Jih, et al., 1999. Characteristics of low-frequency variations embedded in vortex-shedding process. J. Fluids Struct. 13 (3), 339-359.

Murai, Yuichi, 2014. Frictional drag reduction by bubble injection. Exp. Fluids 55, 1-28.

Murai, Y., et al., 2005. Bubble-driven convection around cylinders confined in a channel. J. Fluids Eng. 127, 117-123.

Norberg, Christoffer, 1987. Effect of Reynolds Number and a Low-Intensity Freestream Turbulence on the Flow Around a Circular Cylinder. Chalmers University of Technology.

Norberg, Christoffer, 2003. Fluctuating lift on a circular cylinder: review and new measurements. J. Fluids Struct. 17 (1), 57-96.

Palkin, Egor, et al., 2016. Scrutinizing URANS in shedding flows: The case of cylinder in cross-flow in the subcritical regime. Flow Turbul. Combust. 97, 1017-1046.

Pascal-Ribot, Sylviane, Blanchet, Yves, 2007. Buffeting lift forces and local air-water flow aspects around a rigid cylinder. Int. J. Multiph. Flow 33 (11), 1237-1254.

Pascal-Ribot, Sylviane, Blanchet, Yves, 2011. An improved scaling model of buffeting lift forces in air-water flows. J. Press. Vessel Technol. 133 (2).

Patwardhan, A.W., Joshi, J.B., 1998. Design of stirred vessels with gas entrained from free liquid surface. Can. J. Chem. Eng. 76, 339-364.

Pernod, Laetitia, et al., 2019. Experimental and numerical investigation of the fluid-structure interaction on a flexible composite hydrofoil under viscous flows. Ocean Eng. 194, 106647.

Pettigrew, M.J., et al., 1991. Flow-induced vibration and related technologies in nuclear components. Nucl. Eng. Des. 131, 81-100.

Roshko, Anatol, 1955. On the wake and drag of bluff bodies. J. Aeronaut. Sci. 22 (2), 124-132.

Ryu, Yonguk, Chang, Kuang-An, Lim, Ho-Joon, 2005. Use of bubble image velocimetry for measurement of plunging wave impinging on structure and associated greenwater. Meas. Sci. Technol. 16 (10), 1945.

Schmidt, Oliver T., Colonius, Tim, 2020. Guide to spectral proper orthogonal decomposition. Aiaa J. 58 (3), 1023-1033.

Shakouchi, T., Tian, D., Ida, T., 2002. Behavior of vertical upward gas-liquid two-phase flow past obstacle in rectangular channel (effect of blockage ratio). JSME Int. J. 45 (3), 686-693.

Shakouchi, T., et al., 2007. Flow and aeration characteristics of micro bubble jet flow. Prog. Multiph. Flow Res. 2, 33-38.

Sugiyama, Kazuyasu, Takagi, Shu, Matsumoto, Yoichiro, 2001. Three-dimensional numerical analysis of bubbly flow around a circular cylinder. JSME Int. J. B 44 (3), 319-327.

Sun, Z., Zhang, H., 2009. Measurement of the flow rate and volume void fraction of gas-liquid bubble flow using a vortex flow meter. Chem. Eng. Commun. 197 (2), 145-157.

Suzuki, T., et al., 2009. Wake structure of circular cylinder in microbubble mixture. J. Phys.: Conf. Ser. 147 (1), 012012.

Towne, Aaron, Schmidt, Oliver T., Colonius, Tim, 2018. Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. J. Fluid Mech. 847, 821-867.

Uchiyama, Tomomi, 1999. Numerical simulation of gas-liquid two-phase flow around a rectangular cylinder by the incompressible two-fluid model. Nucl. Sci. Eng. 133 (1), 92-105.

Uchiyama, Tomomi, Degawa, Tomohiro, 2006. Numerical simulation for gas-liquid two-phase free turbulent flow based on vortex in cell method. JSME Int. J. B 49 (4), 1008-1015.

Uchiyama, T., Ishiguro, Y., 2016. Experimental study of flow around a circular cylinder inside a bubble plume. Adv. Chem. Eng. Sci. 6, 269-280.

Voutsinas, A., et al., 2009. Flow visualization and characteristics of vertical gas-liquid bubbly flow around a rectangular cylinder (bubble size effect). J. Phys. Conf. Ser. 147.

Watanabe, H., Ihara, A., Onuma, S., 1990. Effects of a few small air bubbles on the performance of circular cylinder at critical flow range in water. J. Fluids Eng. 112, 67-73.

Watine, Yann, et al., 2023. Vortex-induced vibrations of a cantilevered blunt plate: POD of TR-PIV measurements and structural modal analysis. J. Fluids Struct. 117, 103832.

Williamson, C.H.K., 1996. Vortex dynamics in the cylinder wake. Annu. Rev. Fluid Mech. 28, 477-539.

Yianatos, J.B., 2007. Fluid flow and kinetic modelling in flotation related processes: Columns and mechanically agitated cells - A review. Chem. Eng. Res. Des. 85, 1591-1603.

Yokosawa, M., et al., 1986a. Studies on two-phase cross flow. Part II: Transition Reynolds number and drag coefficient. Int. J. Multiph. Flow 12 (2), 169-184.

Yokosawa, M., et al., 1986b. Studies on two-phase cross flow. Part III: Characteristics of unsteady flow behavior. Int. J. Multiph. Flow 12 (2), 185-202.