1 Introduction
Frictional losses in the form of drag in turbulent flows is a major drain of energy in applications related to process technology, naval transportation and transport of liquified natural gas in pipelines (Ceccio Reference Ceccio2010). It has been known for a long time that the injection of a small concentration of a dispersed phase into a carrier fluid can result in significant drag reduction, making it of interest for fundamental scientific research in order to understand the mechanism and optimise the effect for engineering applications. Drag reduction has been demonstrated in many physical systems in the past, such as bubbles injected into a turbulent boundary layer over a flat plate (Madavan, Deutsch & Merkle Reference Madavan, Deutsch and Merkle1983, Reference Madavan, Deutsch and Merkle1985; Xu, Maxey & Karniadakis Reference Xu, Maxey and Karniadakis2002), bubbles in channel flows (Lu, Fernández & Tryggvason Reference Lu, Fernández and Tryggvason2005; Murai et al. Reference Murai, Fukuda, Oishi, Kodama and Yamamoto2007; Lu & Tryggvason Reference Lu and Tryggvason2008; Dabiri, Lu & Tryggvason Reference Dabiri, Lu and Tryggvason2013; Pang, Wei & Yu Reference Pang, Wei and Yu2014; Park et al. Reference Park, Tasaka, Oishi and Murai2015) and Taylor–Couette (TC) flow (van den Berg et al. Reference van den Berg, Luther, Lathrop and Lohse2005, Reference van den Berg, van Gils, Lathrop and Lohse2007; Murai, Oiwa & Takeda Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008; van Gils et al. Reference van Gils, Narezo Guzman, Sun and Lohse2013); the reader is referred to Ceccio (Reference Ceccio2010) and Murai (Reference Murai2014) for more detailed reviews. The magnitude of reduction in drag or the driving force for a two-phase system when compared with a single-phase system can be massive; van Gils et al. (Reference van Gils, Narezo Guzman, Sun and Lohse2013) showed drag reduction of up to 40 % with just 4 % of bubbles dispersed into TC flow. Various theories have been suggested to explain the origin of this effect; among them are theories based on effective compressibility (Ferrante & Elghobashi Reference Ferrante and Elghobashi2004; L’vov et al. Reference L’vov, Pomyalov, Procaccia and Tiberkevich2005), for low-Reynolds-number flows disruption of coherent vortical structures present in the single-phase flow (Sugiyama, Calzavarini & Lohse Reference Sugiyama, Calzavarini and Lohse2008) and also the effects of bubble deformability (van den Berg et al. Reference van den Berg, van Gils, Lathrop and Lohse2007; van Gils et al. Reference van Gils, Narezo Guzman, Sun and Lohse2013). However, the exact mechanism is still unknown and it may be expected that different mechanisms such as those mentioned above may be dominant in different flow regimes.
TC flow has been one of the classical models to study and understand various concepts in fluid dynamics for the past several decades due to several reasons: (i) TC flow is geometrically simple; (ii) it is a closed flow system with exact global balances between the driving and dissipation; and (iii) it is mathematically well defined through the Navier–Stokes equations and the respective boundary conditions. It is thus an ideal playground to study and understand the phenomenon of drag reduction in shear flows under the influence of a dispersed phase such as bubbles. In this work we perform two-way coupled numerical simulations of a two-phase TC system to investigate the effects of the dispersed phase on the dynamics of the carrier phase and its consequences on the torque required to drive the system.
A TC setup consists of two coaxial independently rotating cylinders with a fluid confined between them. The geometry of the TC system can be described using the radius ratio
${\it\eta}=r_{i}/r_{o}$
and the aspect ratio
${\it\Gamma}=L/d$
, where
$r_{i}$
,
$r_{o}$
are the radius of the inner and outer cylinders respectively,
$d=r_{o}-r_{i}$
is the gap width and
$L$
is the height of the cylinders. Dimensionless radial and axial distances are defined in the form of
$\tilde{r}=(r-r_{i})/d$
and
$\tilde{z}=z/d$
, respectively. The driving of the TC system can be expressed in terms of the Reynolds number of the inner and outer cylinders defined by
$Re_{i}=r_{i}{\it\omega}_{i}d/{\it\nu}$
and
$Re_{o}=r_{o}{\it\omega}_{o}d/{\it\nu}$
, respectively. In these expressions
${\it\omega}_{i}$
,
${\it\omega}_{o}$
are the angular velocities of the inner and outer cylinders, respectively, and
${\it\nu}$
is the kinematic viscosity of the carrier fluid confined in the gap width. In this study the outer cylinder is kept stationary and only the inner cylinder is rotated (
$Re_{o}=0$
).
For a single-phase TC system, at very low driving the flow is purely azimuthal and in a laminar state. Once the driving of the flow is stronger than a critical value, the purely azimuthal, stable laminar flow is disrupted, leading to the formation of large-scale Taylor rolls. Increasing the driving further leads to the onset of interesting flow regimes such as Taylor vortex flow, wavy vortex flow, modulated wavy vortex flow and finally into turbulent Taylor vortices; a detailed phase diagram of these regimes is given by Andereck, Liu & Swinney (Reference Andereck, Liu and Swinney1986).
The global response of the TC system to the driving can be quantified in terms of the friction factor
$$\begin{eqnarray}C_{f}=\frac{{\it\tau}_{w}}{{\textstyle \frac{1}{2}}{\it\rho}U_{i}^{2}},\end{eqnarray}$$
where
${\it\tau}_{w}$
and
$U_{i}$
are the averaged wall shear stress and velocity of the inner cylinder, respectively. The friction factor can be seen as the dimensionless drag of the system on the inner cylinder. For pipe flow the corresponding friction factor is the most common way to express the wall drag in dimensionless form (Pope Reference Pope2000). For TC flow, alternatively, the response of the system can be quantified in terms of a generalised Nusselt number
$Nu_{{\it\omega}}$
, which is the angular velocity transport from the inner to the outer cylinder, non-dimensionalised by its value for laminar flow
$J_{lam}$
(Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007), i.e.
$$\begin{eqnarray}Nu_{{\it\omega}}=\frac{J}{J_{lam}},\quad \text{with}~J=r^{3}(\langle u_{r}{\it\omega}\rangle _{A,t}-{\it\nu}\partial _{r}\langle {\it\omega}\rangle _{A,t}),\end{eqnarray}$$
where
$\langle \cdots \rangle _{A,t}$
represents averaging in the two homogenous (azimuthal and axial) directions and also in time and
$J_{lam}=2{\it\nu}r_{i}^{2}r_{o}^{2}({\it\omega}_{1}-{\it\omega}_{2})/(r_{o}^{2}-r_{i}^{2})$
. The relation between the Nusselt number and the friction factor reads
$$\begin{eqnarray}C_{f}=\frac{2{\rm\pi}LNu_{{\it\omega}}J_{lam}}{\frac{1}{2}r_{i}U_{i}^{2}}.\end{eqnarray}$$
A recent review on single-phase TC flow with the corresponding phase spaces is given by Grossmann, Lohse & Sun (Reference Grossmann, Lohse and Sun2016).
Figure 1. Schematic of a two-phase TC system with a dispersed phase (bubbles or drops) in the gap width. Set-up consists of two coaxial cylinders of length
$L$
with an inner cylinder of radius
$r_{i}$
rotating with an angular velocity
${\it\omega}_{i}$
and an outer cylinder of radius
$r_{o}$
rotating with an angular velocity
${\it\omega}_{o}$
.
In a two-phase TC system a secondary phase is dispersed into the carrier phase (here bubbles dispersed into water) as shown in the schematic of figure 1. The dispersed phase is transported by virtue of various forces acing on them such as buoyancy (
$F_{B}$
), drag (
$F_{D}$
), lift (
$F_{L}$
), added mass (
$F_{A}$
), history forces (
$F_{H}$
), etc. The dispersed phase in turn has a back-reaction force on the carrier flow, thus affecting the angular velocity transport between the two cylinders. This may lead to a change in the torque required to drive the system. A net percentage drag reduction (DR) for a two-phase TC system can be quantified according to
$$\begin{eqnarray}\text{DR}=\frac{\langle Nu_{{\it\omega}}\rangle _{s}-\langle Nu_{{\it\omega}}\rangle _{t}}{\langle Nu_{{\it\omega}}\rangle _{s}}\times 100=\frac{\langle C_{f}\rangle _{s}-\langle C_{f}\rangle _{t}}{\langle C_{f}\rangle _{s}}\times 100;\end{eqnarray}$$
$\langle \cdots \rangle _{s}$
and
$\langle \cdots \rangle _{t}$
correspond to single-phase and two-phase systems, respectively.
The control parameters for the dispersed phase in the case of bubbles are the bubble diameter (
$d_{b}$
), density ratio (in the case of bubbles dispersed into water
${\it\rho}_{p}/{\it\rho}_{f}\ll 1$
), global volume fraction (
${\it\alpha}_{g}=N_{p}V_{p}/V$
, where
$N_{p}$
is total number of dispersed particles,
$V_{p}$
is the volume of individual particle and
$V={\rm\pi}({r_{o}}^{2}-{r_{i}}^{2})L$
is the total volume of the TC system). In dimensionless form the control parameter can be expressed as Froude number
$Fr={\it\omega}_{i}\sqrt{r_{i}/g}$
which is the square root of the ratio of centripetal and the gravitational accelerations. Since it is a ratio of forces in two different directions, it does not describe the equilibrium position of the bubbles in the TC gap width. However, it can be interpreted as the relative strength of the buoyancy of bubbles as compared with the driving in the TC system, as has been done in previous studies (Murai et al.
Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008; Yoshida et al.
Reference Yoshida, Tasaka, Murai and Takeda2009; Watamura, Tasaka & Murai Reference Watamura, Tasaka and Murai2013).
Shiomi et al. (Reference Shiomi, Kutsuna, Akagawa and Ozawa1993) performed one of the earliest experiments on two-phase TC flow by studying the various flow patterns that develop in a two-phase mixture confined in a concentric annulus where they observed various flow patterns such as dispersed bubbly, single spiral, double spiral and triple spiral flows. Djeridi et al. (Reference Djeridi, Fave, Billard and Fruman1999) studied the bubble capture and migration patterns in a TC cell for two different configurations, namely (i) with a free upper surface and (ii) with a top stationary wall. In a subsequent study Djeridi, Gabillet & Billard (Reference Djeridi, Gabillet and Billard2004) found different spatial structures while using air bubbles and cavitation bubbles separately as the dispersed phase. The focus of these studies was primarily on understanding the different flow patterns that developed in a two-phase TC system, but not on the modification of the torque required to drive the cylinder. Murai et al. (Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008) demonstrated experimentally that a tiny percentage of the dispersed phase (0.1 % volume fraction of bubbles dispersed in silicone oil) can reduce the driving torque on the inner cylinder up to 25 %. The maximum inner cylinder Reynolds number reached in these experiments was
$Re_{i}=4500$
and they found that the overall drag reduction decreased with increasing Reynolds number with almost negligible drag reduction at
$Re_{i}=4000$
. With the help of particle tracking velocimetry (PTV) Yoshida et al. (Reference Yoshida, Tasaka, Murai and Takeda2009) studied the relationship between the observed drag reduction and changes in the vortical structures of the bubbly TC system. More recently Watamura et al. (Reference Watamura, Tasaka and Murai2013) investigated the effect of microbubbles on the properties of azimuthal waves found in TC flow (
$\langle Re_{i}\rangle _{max}=1000$
). Fokoua et al. (Reference Fokoua, Gabillet, Aubert and Colin2015) performed experiments on a small gap width (
${\it\eta}=0.91$
) two-phase TC system to study the correlation between bubble arrangement, wavelength of the Taylor vortices and torque modification. They found that the buoyancy of the bubbles plays a crucial role in the behaviour of the two-phase TC system, thus leading to either drag reduction or drag enhancement.
In the highly turbulent regime (
$Re_{i}=10^{5}{-}10^{6}$
), drag reduction of up to 25 % was achieved by injecting millimetric-sized bubbles into the TC flow by van den Berg et al. (Reference van den Berg, Luther, Lathrop and Lohse2005, Reference van den Berg, van Gils, Lathrop and Lohse2007). In this regime, it has been demonstrated that deformability of the bubbles can be a deciding factor in achieving a strong drag reduction of up to 40 % (van Gils et al.
Reference van Gils, Narezo Guzman, Sun and Lohse2013). The torque required to drive the inner cylinder in a two-phase TC system (or a bubbly TC system) depends on various control parameters as mentioned previously. However, in experiments it is not possible to control all of these parameters independently, thus making it challenging to study the individual effect of each control parameter on the system. For example as shown in experiments by Murai et al. (Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008), the Froude number
$Fr$
(which is directly related to the velocity of the inner cylinder) cannot be controlled independently from the inner cylinder Reynolds number (
$Re_{i}$
). In addition, the diameter of the bubbles dispersed into the gap width depends on the local shear in the flow at the point of injection of the bubbles. In contrast to experiments, numerical simulations allow for more flexible control over the parameters and to explore the underlying physics in detail.
While there have been a number of experimental studies demonstrating the effect of bubbles on the drag in two-phase TC system, numerical studies have been fairly limited. The computational burden of fully resolving the interactions between thousands of bubbles and the carrier phase is massive, thus making it extremely difficult to simulate thousands of fully resolved bubbles in a highly turbulent carrier flow (Lu et al.
Reference Lu, Fernández and Tryggvason2005; Lu & Tryggvason Reference Lu and Tryggvason2008; Dabiri et al.
Reference Dabiri, Lu and Tryggvason2013; Tryggvason et al.
Reference Tryggvason, Dabiri, Aboulhasanzadeh and Lu2013). In previous studies, this problem has been dealt with by assuming the bubbles to be point-wise sources of momentum and advecting them in a Lagrangian manner (Mazzitelli, Lohse & Toschi Reference Mazzitelli, Lohse and Toschi2003; Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008; Chouippe et al.
Reference Chouippe, Climent, Legendre and Gabillet2014). Such an approach constrains the size of the particle to the length of the Kolmogorov scale in the carrier flow. Using one-way coupling simulations Climent, Simonnet & Magnaudet (Reference Climent, Simonnet and Magnaudet2007) analysed the migration of point-like bubbles in the Taylor vortex flow and wavy vortex flow regimes. More recently, Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) numerically studied the dispersion of such point-like bubbles in a TC system in which the maximum inner cylinder Reynolds number was
$Re_{i}=5000$
. However, one-way coupling simulations involve only advecting the dispersed phase based on the local flow conditions in the carrier phase without any back-reaction from the dispersed phase onto the carrier phase. This prevents investigation of any form of drag modification mechanisms that maybe active in a two-phase TC system. This problem is overcome by two-way coupled simulations, namely by implementing the momentum exchange between the bubbles and the carrier phase, while ensuring that the simulations are grid independent and not prone to numerical instabilities.
Figure 2. Phase space of a two-phase TC system; ratio of bubble diameter
$d_{b}$
to the viscous length scale
${\it\delta}_{{\it\nu}}$
is plotted against the inner cylinder Reynolds number
$Re_{i}$
. Filled symbols refer to numerical simulations; hollow symbols refer to experimental studies while the lines between the markers indicate a series of experiments. The shaded grey (darker) area indicates the regime of previous two-way coupled simulations on a two-phase TC system, while the shaded orange (lighter) area refers to the current simulations.
Two-way coupling simulations are crucial in order to investigate the effect of a dispersed phase, such as bubbles on the dynamics of a carrier phase. In the context of homogenous isotropic turbulence it has already been shown through two-way coupled simulations, that at large scales the point-like bubbles reduce the intensity in the turbulence spectrum as compared with one-way coupling simulations (Mazzitelli et al.
Reference Mazzitelli, Lohse and Toschi2003). Using a two-way coupled Euler–Lagrange scheme, Sugiyama et al. (Reference Sugiyama, Calzavarini and Lohse2008) reproduced the drag reduction observed in experiments by Murai et al. (Reference Murai, Oiwa and Takeda2005) and concluded that the microbubbles influence the TC system by disrupting the coherent vortical structures and found that the drag reduction decreased at higher
$Re_{i}$
; in those simulations the maximum
$Re_{i}$
reached was 2500 where they observed a drag reduction of about 5 %.
In this present work we use a similar approach by simulating the carrier phase using DNS and tracking the dispersed phase in a Lagrangian manner along with two-way coupling between the phases to extend the
$Re_{i}$
limit to 8000 and gain more insight into two-phase TC drag reduction. We want to find out the effect of bubbles on the driving torque at these higher
$Re_{i}$
. The experiments by Murai et al. (Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008) found almost negligible torque modification at
$Re_{i}=4000$
, but without independent control over the Froude number
$Fr$
due to its implicit dependence on the angular frequency of the inner cylinder and consequently on the Reynolds number
$Re_{i}$
. Will the drag reduction be sustained if
$Fr$
of the bubbles is controlled or would there be drag enhancement? How does the trajectory of a bubble depend on the level of turbulence in the system and the
$Fr$
number? These are the questions we attempt to answer in this work. Figure 2 shows the phase space of a two-phase TC system where the ratio of the bubble diameter
$d_{b}$
to the viscous length scale
${\it\delta}_{{\it\nu}}$
is plotted against the operating Reynolds number
$Re_{i}$
of the inner cylinder. The markers refer to the studies which explore the effect of bubble dispersion and torque modification in the TC system. There are two distinct groupings of studies seen in figure 2; one where the relative particle size
$d_{b}/{\it\delta}_{{\it\nu}}$
is of the order
$O(1{-}10)$
and the rest where
$d_{b}/{\it\delta}_{{\it\nu}}$
is of the order
$O(10^{2}{-}10^{3})$
, which correspond to two-phase TC experiments in the highly turbulent regime. The solid markers and lines are numerical simulations and experiments on a two-phase TC system, respectively. The shaded grey area refers to previous two-way coupled simulations of a two-phase TC system (Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008), while the shaded orange area refers to the
$Re_{i}$
range of current simulations.
This paper is organised as follows. In § 2, the governing equations for the carrier phase and the dispersed phase along with the numerical details are given. Section 3 contains the results of this work: we analyse the net drag reduction, contour plots of mean velocity and fluctuations, bubble trajectories and mean distribution of bubbles in the radial and axial directions. In § 4 we present a summary of the work along with an outlook for further studies.
2 Governing equations and numerical details
2.1 Carrier phase
The dynamics of the carrier phase is computed by solving the Navier–Stokes equations in cylindrical coordinates. For a stationary outer cylinder, i.e.
${\it\omega}_{o}=0$
, the governing equations read
$$\begin{eqnarray}\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\Delta}\boldsymbol{u}+\boldsymbol{f}_{b}(\boldsymbol{x},t),\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$
In (2.1),
$\boldsymbol{f}_{b}(\boldsymbol{x},t)$
is the back-reaction force from the dispersed phase onto the carrier phase; this term is ignored in a one-way coupling simulation. In a two-way coupling simulation
$\boldsymbol{f}_{b}(\boldsymbol{x},t)$
is calculated according to (Mazzitelli et al.
Reference Mazzitelli, Lohse and Toschi2003; Prosperetti & Tryggvason Reference Prosperetti and Tryggvason2007; Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008; Oresta et al.
Reference Oresta, Verzicco, Lohse and Prosperetti2009)
$$\begin{eqnarray}\boldsymbol{f}_{b}(\boldsymbol{x},t)=\mathop{\sum }_{i=0}^{N_{p}}\left(\frac{\text{D}\boldsymbol{u}}{\text{D}t}-\boldsymbol{g}\right)V_{p}{\it\delta}(\boldsymbol{x}-\boldsymbol{x}_{p}(t)).\end{eqnarray}$$
The back-reaction force
$\boldsymbol{f}_{b}(\boldsymbol{x},t)$
is calculated at the exact position of each particle (
$\boldsymbol{x}_{p}(t)$
) which does not necessarily coincide with the location of grid nodes on a fixed Eulerian mesh. A conservative extrapolation scheme is thus necessary to distribute
$\boldsymbol{f}_{b}(\boldsymbol{x},t)$
from the particle position onto the Eulerian grid which will be discussed later.
For a fully coupled true two-way coupled numerical simulation, we would have to take into account the spatial and temporal evolution of the local bubble concentration
$C(x,t)$
in the continuity equation. The continuity equation would then be expressed as
$\partial _{t}(1-C)+\partial _{j}[(1-C)U_{j}]=0$
(Ferrante & Elghobashi Reference Ferrante and Elghobashi2004). However, in our simulations we consider extremely low volume fractions of bubbles (0.1 %); the maximum instantaneous local volume fraction of bubbles is approximately 10 times the global volume fraction (
$\langle C(x,t)\rangle _{max}\sim 1\,\%$
) and in such a case the error in (2.2) is negligible. In cases of bubbly turbulent systems with higher global volume fractions,
$C(x,t)$
can be much higher than 1 % where the bubbles may interact, coalesce and can form dense groups of bubbles in the flow. In such a case the volumetric effect of the bubbles cannot be neglected anymore.
DNS of the carrier phase governing equations are performed by numerically integrating the Navier–Stokes equations with a second-order accurate finite-difference scheme (Verzicco & Orlandi Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015). The equations are integrated in time using a fractional-step approach. This code has been tested and used extensively in the past to study the dynamics of various turbulent flow states in a single-phase TC system for both rotating and stationary outer cylinder configurations (Ostilla-Mónico et al. Reference Ostilla-Mónico, Stevens, Grossmann, Verzicco and Lohse2013, Reference Ostilla-Mónico, Huisman, Jannink, Van Gils, Verzicco, Grossmann, Sun and Lohse2014a ,Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse b ,Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse c ). Grid nodes are spaced uniformly in the azimuthal and axial directions where the boundary conditions are periodic, while a clipped Chebychev type clustering is employed for spacing the radial nodes.
2.2 Dispersed phase
The bubbles dispersed into the TC system are assumed to be clean, non-deformable and are tracked in a Lagrangian manner with effective forces such as drag, lift, added mass and buoyancy acting on them (Mazzitelli et al.
Reference Mazzitelli, Lohse and Toschi2003; Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008; Oresta et al.
Reference Oresta, Verzicco, Lohse and Prosperetti2009; Chouippe et al.
Reference Chouippe, Climent, Legendre and Gabillet2014; Lakkaraju, Toschi & Lohse Reference Lakkaraju, Toschi and Lohse2014). The generalised momentum equation for particles of density
${\it\rho}_{p}$
dispersed into a carrier fluid of density
${\it\rho}_{f}$
which takes into account the drag, lift, added mass and buoyancy forces reads
$$\begin{eqnarray}\displaystyle {\it\rho}_{p}V_{p}\frac{\text{d}\boldsymbol{v}}{\text{d}t} & = & \displaystyle ({\it\rho}_{p}-{\it\rho}_{f})V_{p}\boldsymbol{g}-C_{D}\frac{{\rm\pi}d_{b}^{2}}{8}{\it\rho}_{f}|\boldsymbol{v}-\boldsymbol{u}|(\boldsymbol{v}-\boldsymbol{u})+{\it\rho}_{f}V_{p}C_{M}\left(\frac{\text{D}\boldsymbol{u}}{\text{D}t}-\frac{\text{d}\boldsymbol{v}}{\text{d}t}\right)\nonumber\\ \displaystyle & & \displaystyle +\,{\it\rho}_{f}V_{p}\frac{\text{D}\boldsymbol{u}}{\text{D}t}-C_{L}{\it\rho}_{f}V_{p}(\boldsymbol{v}-\boldsymbol{u})\times {\bf\omega}.\end{eqnarray}$$
Based on the particle velocity, the position of each particle is updated according to
$$\begin{eqnarray}\frac{\text{d}\boldsymbol{x}_{p}}{\text{d}t}=\boldsymbol{v}.\end{eqnarray}$$
Since we consider clean spherical bubbles the effect of history forces is assumed negligible (Magnaudet & Eames Reference Magnaudet and Eames2000; Mazzitelli et al.
Reference Mazzitelli, Lohse and Toschi2003; Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008). In (2.4) and (2.5),
$\boldsymbol{u}$
is the velocity of the carrier phase at the particle position,
$\boldsymbol{v}$
is the velocity of individual particle,
$\boldsymbol{x}_{p}$
is the position of the particle, while
$C_{D},C_{M},C_{L}$
are the drag, lift and added mass coefficients, respectively. Since we consider only very light particles (bubbles in water, i.e.
${\it\rho}_{p}/{\it\rho}_{f}\ll 1$
), in this case (2.4) can be simplified into
$$\begin{eqnarray}C_{M}\frac{\text{d}\boldsymbol{v}}{\text{d}t}=-\boldsymbol{g}-C_{D}\frac{{\rm\pi}d_{b}^{2}}{8V_{p}}{\it\rho}_{f}|\boldsymbol{v}-\boldsymbol{u}|(\boldsymbol{v}-\boldsymbol{u})+(1+C_{M})\frac{\text{D}\boldsymbol{u}}{\text{D}t}-C_{L}(\boldsymbol{v}-\boldsymbol{u})\times {\bf\omega}.\end{eqnarray}$$
In order to close (2.6), along with the coefficients for drag, lift and added mass (i.e.
$C_{D}$
,
$C_{L}$
and
$C_{M}$
), information on the velocity, total acceleration and vorticity in the carrier phase is required at the exact location of the bubble. Since we assume that the dispersed phase is composed of clean spherical non-deformable bubbles, the value of the added mass coefficient is
$C_{M}=1/2$
(Auton, Hunt & Prud’Homme Reference Auton, Hunt and Prud’Homme1988; Magnaudet & Eames Reference Magnaudet and Eames2000). In bubbly turbulent flows the exact value of the lift coefficient is not very well known and could be a source of discrepancy in numerical studies. Climent et al. (Reference Climent, Simonnet and Magnaudet2007), Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) use a lift coefficient which is dependent on both the bubble Reynolds number and the local shear intensity. By systematically varying the lift coefficient acting on the bubbles from 0 to 0.5 Sugiyama et al. (Reference Sugiyama, Calzavarini and Lohse2008) found that the value of the lift coefficient plays a crucial role in the mean bubble distribution. We use
$C_{L}=1/2$
(Auton Reference Auton1987), a simplified approach which has been used in many previous studies (Mazzitelli et al.
Reference Mazzitelli, Lohse and Toschi2003; Ferrante & Elghobashi Reference Ferrante and Elghobashi2004; Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008; Oresta et al.
Reference Oresta, Verzicco, Lohse and Prosperetti2009; Lakkaraju et al.
Reference Lakkaraju, Toschi and Lohse2014). Small changes in the magnitude of the lift coefficient can result in different trajectories of the bubbles but in our simulations we find that the qualitative behaviour of the bubbles remains almost the same.
The drag coefficient
$C_{D}$
is computed for each bubble individually and is dependent on the bubble Reynolds number defined as
$Re_{b}=d_{b}|\boldsymbol{v}-\boldsymbol{u}|/{\it\nu}$
(Mei & Klausner Reference Mei and Klausner1992; Magnaudet & Eames Reference Magnaudet and Eames2000), namely
$$\begin{eqnarray}C_{D}=\frac{16}{Re_{b}}\left[1+\frac{Re_{b}}{8+\frac{1}{2}(Re_{b}+3.315\sqrt{Re_{b}})}\right].\end{eqnarray}$$
The difference in the location of the bubbles and the grid nodes of the Eulerian mesh restricts us from calculating local flow quantities (i.e. velocity, total acceleration and vorticity) directly from the carrier phase solution. Since the velocities are in three different positions in a staggered grid arrangement, all three velocities belonging to any specific computational cell is first interpolated to the cell centre using the values from the nearest grid points. A tri-linear interpolation scheme, containing information from all of the cell centres surrounding a specific particle is used to calculate the carrier phase velocity at the exact particle position. The same approach is employed for the interpolation of the total acceleration and the vorticity in the flow. The error in the interpolation decreases asymptotically as
$({\rm\Delta}x)^{2}$
as the grid spacing
${\rm\Delta}x$
tends to zero (Yeung & Pope Reference Yeung and Pope1988).
In addition to interpolating local flow quantities, the back-reaction force from the bubbles (2.3) needs to be extrapolated to the surrounding grid nodes. The back-reaction force from the particle
$\boldsymbol{f}_{b}(\boldsymbol{x},t)$
is extrapolated to the surrounding grid nodes residing in fixed computational volume (which is smaller than the bubble volume) using an exponential distribution function, which decreases monotonically with the distance between the grid nodes and the particle position. In the context of simulating particle-laden flows such an approach has already been discussed in detail by Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013). The extrapolation scheme ensures second-order accuracy, conservation of the point back-reaction force while extrapolating the data from a Lagrangian location to a Eulerian mesh and is also stable when a non-uniform grid distribution is used. Time integration of (2.5) and (2.6) is done by a second-order accurate Runge–Kutta scheme. The code has been parallelised using MPI and OpenMP directives and has a strong scaling of up to 1000 cores.
Figure 3. Percentage of drag reduction (DR) as a function of the inner Reynolds number
$Re_{i}$
. The results from the current code are compared against the numerical simulations of Sugiyama et al. (Reference Sugiyama, Calzavarini and Lohse2008) and experimental measurements of Murai et al. (Reference Murai, Oiwa and Takeda2005). The global volume fraction of bubbles
${\it\alpha}_{g}=0.125{-}0.670\,\%$
;
$Fr=0.3{-}1.0$
.
In order to validate the two-phase code, we perform simulations with the exact resolutions and parameters as done in Sugiyama et al. (Reference Sugiyama, Calzavarini and Lohse2008) and compare our results in figure 3. Very good agreement is seen between the results of the current simulations and the simulations of Sugiyama et al. (Reference Sugiyama, Calzavarini and Lohse2008) and also the experiments of Murai et al. (Reference Murai, Oiwa and Takeda2005) where the drag reduction (DR) is plotted against the inner cylinder Reynolds number. For the rest of the cases simulated in this paper the parameters chosen are as follows: the geometry of the TC system is fixed by fixing the radius ratio
${\it\eta}=0.833$
and the aspect ratio to
${\it\Gamma}=4$
. The global volume fraction
${\it\alpha}_{g}$
is fixed to 0.1 %, the relative bubble diameter to
$d_{b}/d=0.01$
, while the Froude number
$Fr$
of the bubbles is varied from 0.16 to 2.56.
To reduce the computational cost, the simulation is initially run without any dispersed phase and once the flow in the single phase reaches a statistically stationary state, bubbles are placed at random positions throughout the domain. The velocity of these bubbles is equated to that of the velocity of the carrier phase interpolated at the bubble location and the two-phase flow is allowed to develop. Once the bubbles are introduced into the flow, the previously developed statistically stationary state of the single-phase flow is disrupted and transients develop in the two-phase flow which eventually die out after approximately 100 full cylinder rotations. In addition, it has been ensured that the mean axial velocity of the carrier phase in the domain is equal to zero. Similar to the single-phase system, azimuthal and axial periodicity is employed for the dispersed phase, i.e. when the bubbles exit an azimuthal or axial boundary they are re-entered at the opposite boundary at the exact same location with the same velocity. This approach is obviously different from experimental studies (Murai et al. Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008; van Gils et al. Reference van Gils, Narezo Guzman, Sun and Lohse2013; Fokoua et al. Reference Fokoua, Gabillet, Aubert and Colin2015), where the bubbles are injected at the bottom and are collected at the top and can be a cause of discrepancy in the comparison between numerical and experimental data.
In the current simulations, through some additional test runs we have ensured that the chosen axial extent of the domain does not influence the results. (Details on the Nusselt number for both single-phase and two-phase case with
$Fr=0.16$
is compared for
${\it\Gamma}=4$
and
${\it\Gamma}=8$
in table 3 of the Appendix.) However, it can be expected that after the two-phase system achieves a statistically stationary state, both the experimental and the computational systems behave in a qualitatively similar manner. An elastic bounce model is used for the interaction of the bubbles with the inner and outer cylinder wall as implemented in previous studies (Climent et al.
Reference Climent, Simonnet and Magnaudet2007; Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008; Chouippe et al.
Reference Chouippe, Climent, Legendre and Gabillet2014). While modifying the elastic bounce interaction may result in a change in the local bubble concentration near the wall, small variations in the coefficient of restitution of bubble impact with the wall does not influence the qualitative behaviour of the motion of bubbles which we study in a later section. In addition, since the bubbles are assumed to be Lagrangian points, there is a possibility of bubbles overlapping among each other. Owing to the very low global volume fraction in the system (
${\it\alpha}_{g}=0.1\,\%$
), the fraction of bubbles overlapping onto each other when taking into consideration their physical size is less than
$10^{-3}$
, which we consider to be negligible.
Figure 4. Typical time series of the instantaneous torques at the inner (solid red line) and outer cylinder (dashed blue line). Here
$Re_{i}=2500$
,
$Fr=0.16$
,
${\it\alpha}_{g}=0.1\,\%$
. The temporal origin is arbitrary and after the system has achieved a statistically stationary state.
The procedure for calculating the Nusselt number
$Nu_{{\it\omega}}$
for both the single-phase and the two-phase systems is as follows. The simulation is run for at least 50 full cylinder rotations after statistical stationarity is achieved and convergence is ensured by comparing the mean azimuthal velocity profiles for 50 and 100 full cylinder rotations. In addition, it is ensured that the net angular momentum current in the radial direction is constant. The torques on both cylinders are then averaged in time and it is ensured that the difference in
$\langle Nu_{{\it\omega}}\rangle$
for both inner and outer cylinders is less than 1 %. Figure 4 shows a typical time series of the non-dimensional angular velocity transport
$Nu_{{\it\omega}}$
on both the inner and outer cylinder after the transients in the two-phase system die out and the system reaches a statistically stationary state. The net percentage drag reduction is then calculated according to (1.4). To maintain a global volume fraction of
${\it\alpha}_{g}=0.1\,\%$
with bubbles of size
$d_{b}/d=0.01$
, approximately 50 000 bubbles are required in the simulation. In order to calculate the effective forces acting on the point-wise particles according to (2.6), and advect them with sufficient accuracy, the size of the bubble should be of the same order as the smallest relevant length scale in the carrier phase. It can be seen from figure 2 that the ratio of bubble diameter (
$d_{b}$
) to the viscous length scale
${\it\delta}_{{\it\nu}}$
is of the order of one for all of the simulations considered in this work. In table 1 we give details on the resolutions for various Reynolds number and the corresponding Stokes number of the bubbles (
$St={\it\tau}_{b}/{\it\tau}_{{\it\eta}}$
, where
${\it\tau}_{b}=d_{b}^{2}/24{\it\nu}$
is the bubble response time and
${\it\tau}_{{\it\eta}}$
is the Kolmogorov time scale). In the next section we discuss the main results of this paper, where we study the effect of the operating Reynolds number
$Re_{i}$
and the Froude number
$Fr$
on the global response of the two-phase TC system and also on the local flow dynamics.
Table 1. Numerical details of the simulations. The first column is the operating Reynolds number of the inner cylinder, the second column is the radial grid spacing near the inner wall, the third column shows the radial grid refinement near the inner cylinder, the fourth column shows the diameter of the bubble normalised by the viscous length scale and the fifth column is the bubble Stokes number for each case. These values are for a geometry with radius ratio
${\it\eta}=0.833$
with an aspect ratio
${\it\Gamma}=4$
and a rotational symmetry
$n_{sym}=6$
in the azimuthal direction while the maximum value of the bubble size relative the radial grid refinement in all simulations is
$d_{b}/{\rm\Delta}r_{min}\sim 3$
.
Figure 5. (a) Friction factor
$C_{f}(Re_{i})$
for single-phase TC flow and for bubbly TC flow with bubbles of five different Froude numbers. Only for
$Fr<1$
does drag reduction occur. This drag reduction is much better expressed in terms of the percentage of drag reduction (DR), which is plotted in (b), again as a function of the inner Reynolds number
$Re_{i}$
. Positive DR indicates that the driving torque in the two-phase case is lower than that of the single-phase case. The dashed line shows the path taken by
$Re_{i}$
–
$Fr$
number in the experiments by Murai et al. (Reference Murai, Oiwa and Takeda2008).
3 Results
3.1 Drag reduction
We first focus on the change in the global response of the TC system with the addition of the dispersed phase. In figure 5(a) we compare the friction coefficient
$C_{f}$
for the single-phase and two-phase cases for five different Froude numbers
$Fr=0.16$
, 0.32, 1.28, 2.56 and 1000. When
$Fr<1$
and in particular for low
$Re_{i}$
, we observe a reduction in the friction factor
$C_{f}$
for the two-phase as compared with the single flow case or as compared with cases with
$Fr\gg 1$
. In order to observe this difference more clearly we compute the net percentage drag reduction DR according to (1.4). It is shown as a function of
$Re_{i}$
for the five different
$Fr$
in figure 5(b). Almost 4 % drag reduction is achieved at
$Re_{i}=2500$
and
$Fr=0.16$
. For fixed and small enough
$Re_{i}$
, increasing the
$Fr$
number of the bubbles leads to a decrease in the drag reduction. Vice versa when the
$Fr$
number is kept fixed, the drag reduction DR decreases for increasing
$Re_{i}$
. This holds in particular for configurations with
$Fr<1$
, i.e. when the buoyancy of the bubbles is dominant over the driving of the system. For
$Fr>1$
, for at least within the examined
$Re_{i}$
range, there is overall no systematic trend with
$Re_{i}$
and negligible drag reduction, thus indicating that bubble buoyancy has a strong role in achieving drag reduction. A unique difference in these results when compared with the experimental study of Murai et al. (Reference Murai, Oiwa and Takeda2008) is that here we have independent control over the
$Fr$
number, i.e. independently of the Reynolds number
$Re_{i}$
. Murai et al. (Reference Murai, Oiwa and Takeda2008) found almost negligible drag reduction beyond
$Re_{i}=3000$
; however, in their setup the
$Fr$
of the bubbles is dependent on
$Re_{i}$
and at
$Re_{i}=2000$
the
$Fr$
number was already above one. The path taken by the
$Re_{i}$
–
$Fr$
number in the experiments is shown in figure 5(b) by a dashed line. In the current simulations, in contrast, we observe that by fixing the
$Fr$
number of the bubbles to less than one, drag reduction can be observed even at
$Re_{i}=8000$
. The details on the Nusselt number for the single-phase and two-phase cases for these simulations in given in table 3 in the Appendix.
It can be expected that at the nearly asymptotically large value of
$Fr=1000$
, the centripetal acceleration can play the role of buoyancy for bubbles immersed in a horizontal channel flow (Xu et al.
Reference Xu, Maxey and Karniadakis2002; Murai et al.
Reference Murai, Fukuda, Oishi, Kodama and Yamamoto2007; Harleman Reference Harleman2012; Pang et al.
Reference Pang, Wei and Yu2014) which may lead to drag reduction. However, we do not find any such drag reduction in the TC system as seen in figure 5 and we think this might be due to different physical mechanisms governing drag reduction in both systems. In comparison with the current simulations, experiments with larger bubbles and thus higher gas volume fractions (
$d_{b}^{+}\gg 1$
,
${\it\alpha}_{g}>1\,\%$
) have been performed in (bubbly) turbulent channel flows to obtain sustainable levels of drag reduction (Xu et al.
Reference Xu, Maxey and Karniadakis2002; Murai et al.
Reference Murai, Fukuda, Oishi, Kodama and Yamamoto2007). In experiments with smaller bubbles and lower gas volume fractions (
$d_{b}^{+}\sim 1$
,
${\it\alpha}_{g}\sim 1\,\%$
), Harleman (Reference Harleman2012) and Pang et al. (Reference Pang, Wei and Yu2014) found extremely small levels of drag reduction which is similar to what we find in our simulations. In addition, the drag reduction effect observed in TC flows in the low-Reynolds-number regime is through disruption of coherent Taylor rolls which are fixed in space and time for single-phase flows and are responsible for the majority of the angular momentum transport. In contrast, for channel flows cospectra analysis shows that these wall-attached large-scale structures are inactive (Jiménez Reference Jiménez2011).
3.2 Carrier phase velocity fields
In figure 6 we compare the contour plots of the azimuthal velocity, averaged in the azimuthal direction and over time. The Reynolds number is fixed at
$Re_{i}=2500$
and averaged contours of two systems with
$Fr=0.16$
and
$Fr=1.28$
is compared with that of the single-phase flow. As mentioned earlier (also see again figure 5), the DR decreases with increasing
$Fr$
for fixed
$Re_{i}$
. In figure 6(c) (
$Fr=1.28$
), a clear signature of the Taylor vortex can be observed with a structure very similar to that of the single-phase case as seen in figure 6(a). When
$Fr<1$
(figure 6
b) a strong footprint of the Taylor vortex is not observed anymore. The strong buoyancy of the bubbles as compared with the driving of the system (i.e.
$Fr<1$
) is responsible for disrupting the Taylor vortices, which have concentrated regions of high strain rates and are thus highly dissipative (Sugiyama et al.
Reference Sugiyama, Calzavarini and Lohse2008). With increasing
$Fr$
the Taylor vortex structure resembles more to that of a single-phase system (Ostilla-Mónico et al.
Reference Ostilla-Mónico, Stevens, Grossmann, Verzicco and Lohse2013), resulting in a drop in DR.
Figure 6. Contour plots of azimuthally and time-averaged azimuthal velocity field
$\langle \bar{u}_{{\it\theta}}\rangle _{{\it\theta},t}$
for a fixed
$Re_{i}=2500$
: (a) single-phase, (b,c) two-phase flow with (b)
$Fr=0.16$
and (c)
$Fr=1.28$
. For the low-
$Fr$
case the Taylor rolls are considerably weakened by the strongly buoyant bubbles.
Figure 7. Contour plots of azimuthal and time-averaged r.m.s. of the velocity fluctuations (
$u_{{\it\theta}}^{\prime }$
) for a fixed at
$Re_{i}=2500$
: (a) single-phase, (b,c) two-phase flow with (b)
$Fr=0.16$
and (c)
$Fr=1.28$
. The same colour bar is used for all three plots. Again, the effect of the strongly buoyant bubbles (
$Fr=0.16$
) on the Taylor rolls can be clearly observed.
This is also observed in figure 7 where we show contour plots of azimuthally and time-averaged root-mean-square (r.m.s.) of the azimuthal velocity fluctuation
$u_{{\it\theta}}^{\prime }$
. These are computed in form of the r.m.s. value as
$u^{\prime }(r,z)=[\langle u^{2}\rangle _{{\it\theta},t}-\langle u\rangle _{{\it\theta},t}^{2}]^{1/2}$
. For the single-phase flow (figure 7
a), there exist two localised regions near the inner cylinder which show a peak in the velocity fluctuations. These regions are associated to the sites of plume ejection (Ostilla-Mónico et al.
Reference Ostilla-Mónico, Stevens, Grossmann, Verzicco and Lohse2013, Reference Ostilla-Mónico, Huisman, Jannink, Van Gils, Verzicco, Grossmann, Sun and Lohse2014a
,Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse
b
). In a two-phase system with strong buoyant bubbles (
$Fr<1$
, i.e. figure 7
b) there is a peak in the fluctuation along the complete axial extent of both the inner and outer cylinders. The bubbles rising near the walls of the cylinders by virtue of their buoyancy are responsible for such a behaviour. However, unlike in the single-phase case where the fluctuations extend into the gap width in the two-phase case with
$Fr<1$
they are localised near the cylinder walls giving and indication that the plumes are much weaker. This is not observed in the case where the bubbles are weakly buoyant (
$Fr>1$
, i.e. figure 7
c). Instead, we see that the axial extent of the plume ejection site is extended when compared with the single-phase system. This is a result of bubbles being captured by the Taylor vortices at the sites of plume ejection after sliding on the inner cylinder for a short period. For a single-phase flow, the plumes ejected on the inner cylinder wall are responsible for majority of the angular velocity transport. In addition, they assist in the transition to the ultimate state of turbulence where the turbulence is fully developed in both bulk and boundary layers (Ostilla-Mónico et al.
Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014b
).
Figure 8. Three-dimensional instantaneous snapshot of the azimuthal velocity field
$u_{{\it\theta}}$
for a fixed
$Re_{i}=2500$
: (a) single-phase, (b,c) two-phase flow with (b)
$Fr=0.16$
and (c)
$Fr=1.28$
.
In figure 8 we show three-dimensional instantaneous snapshots of the azimuthal velocity field for the single-phase case and the two-phase cases for both
$Fr=0.16$
and
$Fr=1.28$
. The important observation made here is that while the ejection of plumes from the inner cylinder (
$\tilde{r}=0$
) is very strong for the single-phase case and
$Fr=1.28$
, it is much weaker for
$Fr=0.16$
. The strong buoyancy of the bubbles (
$Fr=0.16$
) causes the plume ejection region to sweep the inner cylinder surface and thus when averaged in time a clear footprint of the ejection region is no longer observed (cf. figure 6
b). The weakening of the plumes and the Taylor rolls is also clearly observed in figure 9 where we plot contours of the normalised wall-parallel strain rate close to the wall. The local strain rate is normalised using the maximum strain rate for each case and these contour plots are a direct indication of the spatial distribution of the wall shear stress on the inner cylinder. As seen in Ostilla-Mónico et al. (Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014b
), the local inner cylinder boundary layer profiles vary depending on whether the region belongs to a plume ejection site or a plume impacting site. This is also reflected in the contours shown in figure 9, where in the case of single-phase flow high
$e_{r,{\it\theta}}$
regions can be observed at the plume impacting sites. The strain rates in these sites are comparatively weaker in the case of
$Fr=0.16$
due to the loss of strength in the plume ejections and as a consequence also the coherent Taylor rolls. With increase in
$Fr$
numbers which is shown in figure 9(c), the plumes ejections regain their strength and the flow becomes similar to the single-phase case with bubbles having a much weaker effect.
Figure 9. Near-wall strain rate contours (
$\hat{e}_{r,{\it\theta}}$
) at a constant radius cut (
$r^{+}=0.25$
) for
$Re_{i}=2500$
. (a) Corresponds to the single-phase, (b,c) to the two-phase case with (b)
$Fr=0.16$
and (c)
$Fr=1.28$
.
Figure 10. Comparison of azimuthal, axial and time-averaged viscous dissipation profiles near the inner cylinder wall.
$Re_{i}$
in both (a,b) is 2500. For the two-phase system, (a)
$Fr=0.16$
and (b)
$Fr=1.28$
.
In order to understand more on how the bubbles affect these plumes, in figure 10 we plot the viscous dissipation profile near the inner cylinder wall. The total dissipation per unit mass is calculated as a volume and time average of the local dissipation rate, and as shown by Eckhardt et al. (Reference Eckhardt, Grossmann and Lohse2007) is directly related to the driving torque. The viscous dissipation is normalised accordingly
$\hat{{\it\epsilon}}={\it\epsilon}/{\it\nu}(U_{i}/d)^{2}$
and in figure 10, we compare the dissipation profiles of a two-phase system with strongly buoyant bubbles (
$Fr<1$
) and weakly buoyant bubbles (
$Fr>1$
) with a single-phase system. When
$Fr<1$
, the viscous dissipation is lower than that compared with the single phase and the difference can be clearly observed in figure 10(a). This is not the case in figure 10(b), where the difference is almost negligible and thus also resulting in minimal drag reduction. Such a behaviour is related to the pattern of distribution of bubbles near the inner cylinder and is discussed in a later section. In addition to a drop in DR with increasing
$Fr$
at a fixed
$Re_{i}$
, we also note a drop in DR with
$Re_{i}$
even for bubbles with
$Fr<1$
. With increase in
$Re_{i}$
, the coherent flow structures become less responsible for the angular velocity transport (Ostilla-Mónico et al.
Reference Ostilla-Mónico, Stevens, Grossmann, Verzicco and Lohse2013, Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014b
). The strong buoyancy effect of the bubbles (i.e.
$Fr<1$
) thus acts on weaker coherent structures and lose their efficiency on affecting the angular velocity transport in the system. This trend was also observed in the experiments by Murai et al. (Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008) who found that the DR dropped to almost zero at
$Re_{i}=4000$
. We show that although there is a drop in DR with increasing
$Re_{i}$
, little DR can be attained even at
$Re_{i}=8000$
provided
$Fr<1$
.
3.3 Motion of the dispersed phase
Figure 11. Full trajectory of a single bubble in the
$r$
–
$z$
plane of the two-phase TC system. (a,c)
$Re_{i}=2500$
and (b,d)
$Re_{i}=5000$
. The
$Fr$
number for the bubbles in (a,b) is
$Fr=0.16$
and in (c,d) is
$Fr=1.28$
. The bubble is tracked for almost 20 full inner cylinder rotations in all cases after the system has reached statistically stationary state. The same colour scale is used for cases with the same
$Fr$
numbers and indicates the axial velocity
$v_{z}$
of the bubbles. While the less buoyant bubbles (
$Fr=1.28$
, (c,d)) tend to get trapped in the Taylor rolls, the more buoyant ones (
$Fr=0.16$
, (a,b)) rise through and weaken them.
In this section we look at the different trajectories of bubbles for different
$Re_{i}$
and
$Fr$
numbers. In figure 11, we show four different trajectories of an individual bubble for two different
$Re_{i}$
and
$Fr$
numbers, respectively. There is a strong azimuthal motion imposed on the bubbles by the carrier flow while their interaction with the underlying coherent Taylor rolls can be easily visualised in the
$r{-}z$
plane as has been done in figure 11. For strong buoyant bubbles, i.e.
$Fr<1$
(cf. figure 11
a), the bubble primarily has an upward motion moving through different Taylor vortices as expected. For higher
$Re_{i}$
where the turbulence is more developed and intense, the bubbles have a tendency to get trapped in the small vortical structures inside the large-scale Taylor rolls. This can be observed at an axial height of
$\tilde{z}\simeq 0.5$
and
$\tilde{z}\simeq 2.5$
in figure 11(b). So also for larger
$Fr$
number (here
$Fr=1.28$
) at fixed
$Re_{i}=2500$
as shown in figure 11(c), the bubbles have a tendency to be trapped inside the Taylor vortex (in this case at a height of
$\tilde{z}\simeq 3.0$
) while rising along the inner cylinder wall. This behaviour has been observed also in the experiments by Murai et al. (Reference Murai, Oiwa and Takeda2008) and Fokoua et al. (Reference Fokoua, Gabillet, Aubert and Colin2015). The influence of such a motion can be seen in the contour plots of the velocity fluctuations (cf. figure 7
c). When the
$Re_{i}$
is increased to 5000 and
$Fr>1$
, the smooth structured motion of the bubble inside the Taylor vortex as seen in
$Re_{i}=2500$
is lost. Due to increased levels of turbulence in addition to loss of importance of coherent structures the bubble motion is more erratic as compared with the previous case. For both
$Re_{i}$
, a clear difference can be observed in the trajectories of the bubble depending on whether
$Fr>1$
or
$Fr<1$
. For low
$Re_{i}$
a bubble which smoothly passes through each Taylor vortex tends to get trapped with increase in the
$Fr$
number; while for a higher
$Re_{i}$
the bubbles which spend short periods of time in small-scale vortices moves in a more erratic manner with increasing
$Fr$
. In all of the above cases, it is to be noted that there is also a strong azimuthal motion of bubbles.
As observed in figure 11 the bubbles dispersed into the TC system exhibit various kinds of motion such as sliding along the inner cylinder wall, outer cylinder wall, organised or erratic motion in the Taylor vortices. In order to understand and categorise such bubble motion more precisely, it would be useful to look at the mean radial and axial distribution of the bubbles and also their corresponding Reynolds numbers. For this purpose we divide the domain into three different regions (i) inner boundary layer (IBL), (ii) bulk region (BULK) and (iii) outer boundary layer (OBL). The method of calculating the inner and outer boundary layer thickness is described in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Stevens, Grossmann, Verzicco and Lohse2013, Reference Ostilla-Mónico, van der Poel, Verzicco, Grossmann and Lohse2014b ) which we discuss here in brief. A straight line is first fitted through the first three computational grid points near the inner and respective outer cylinder wall. For the bulk a straight line is fitted through the point of inflection and the two nearest grid points. The intersection of the line drawn through the bulk and the lines near the inner and outer cylinder wall give the inner and outer boundary layer thickness, respectively.
Figure 12. Probability distribution functions (p.d.f.s) of the bubble Reynolds number for (a–c)
$Re_{i}=2500$
and (d–f)
$Re_{i}=5000$
. (a,d) Correspond to the inner boundary layer region (IBL), (b,e) to the bulk region (BULK) and (c,f) to the outer boundary layer region (OBL).
In figure 12, we show normalised probability distribution functions (p.d.f.s) of the Reynolds number of the bubbles (
$Re_{b}$
) in the three different regions of the domain as described above. When the Froude number of the bubbles is less than one,
$Re_{b}$
has a comparatively narrow distribution as compared with higher Froude number systems, regardless of the operating Reynolds number
$Re_{i}$
of the inner cylinder and also of the position of the bubbles in the domain. When
$Fr<1$
, the bubbles primarily have an upward drifting motion through the Taylor rolls (cf. figure 11) independent of their position in the domain which is reflected in these histograms. However, for increasing inner cylinder Reynolds number
$Re_{i}$
the p.d.f. becomes more symmetric as compared to lower
$Re_{i}$
which is a result of bubbles getting trapped in the intense vortical structures (also seen in figure 11
b). The distribution of
$Re_{b}$
becomes wider for higher
$Fr$
numbers (less buoyancy); more in the inner boundary layer than the outer boundary layer. This is a result of combination of two events; namely the ejection of plumes from the inner cylinder where the bubbles are trapped before being pulled into the Taylor roll and a relatively stronger component of the azimuthal slip velocity near the inner cylinder as compared to the outer cylinder where it is close to zero.
Now that we have looked at how the motion of bubbles depends on the
$Re_{i}$
and
$Fr$
numbers, in the next part we study the mean distribution of the bubbles in the domain. In table 2 we show the percentage volume fraction of bubbles accumulated in the inner boundary layer (
${\it\alpha}_{i}$
), outer boundary layer (
${\it\alpha}_{o}$
) and the bulk (
${\it\alpha}_{b}$
). It is clear that the percentage of bubbles accumulated in the inner boundary layer is highest and for all cases it increases with
$Fr$
. With increasing
$Fr$
the bubbles initially residing in the bulk and the outer boundary layer now migrate to the inner cylinder wall. Also the percentage distribution of the bubbles between the three different regions (
${\it\alpha}_{i}$
,
${\it\alpha}_{b}$
and
${\it\alpha}_{o}$
) does not change significantly for a fixed
$Fr$
number with increasing the
$Re_{i}$
number.
Figure 13. Radial profiles of the azimuthally, axially and time-averaged local bubble volume fraction (
${\it\alpha}$
) for (a)
$Re_{i}=2500$
and (b)
$Re_{i}=5000$
. Solid blue lines refer to
$Fr=0.16$
, while the dashed red lines refer to
$Fr=1.28$
. The insets show the same profile near the inner cylinder wall. Here
${\it\alpha}$
is normalised with
${\it\alpha}_{g}=0.1\,\%$
, which is the global volume fraction in this case.
Table 2. Percentage volume fraction of bubbles in the inner boundary layer (
${\it\alpha}_{i}$
), bulk (
${\it\alpha}_{b}$
) and outer boundary layer (
${\it\alpha}_{o}$
) for two
$Re_{i}$
and
$Fr$
numbers.
Figure 14. Axial profiles of the azimuthally, radially and time-averaged local bubble volume fraction (
${\it\alpha}$
) for (a)
$Re_{i}=2500$
and (b)
$Re_{i}=5000$
. Solid blue lines refer to
$Fr=0.16$
, while the dashed red lines refer to
$Fr=1.28$
. Here
${\it\alpha}$
is normalised with
${\it\alpha}_{g}=0.1\,\%$
, which is the global volume fraction in this case. Again, the trapping of the less buoyant bubbles in the Taylor rolls becomes evident (larger
$Fr$
).
In figure 13 we show the azimuthally, axially and time-averaged radial profiles of the bubble volume fraction for
$Re_{i}=2500$
and 5000 which gives a clear indication of accumulation of bubbles near the inner cylinder wall independent of the
$Re_{i}$
number. The local volume fraction of the bubbles near the inner cylinder wall (inset of figure 13) is higher for
$Fr=0.16$
as compared with
$Fr=1.28$
. But table 2 shows that there is a higher percentage of bubble accumulation near the inner cylinder wall for
$Fr=1.28$
, which means that there is a higher percentage of bubbles sticking to the inner cylinder wall. When
$Fr=0.16$
, the local volume fraction (
${\it\alpha}$
) in the bulk of the system is close to that of the global volume fraction (
${\it\alpha}_{g}$
). With increase in
$Fr$
number, the local volume fraction decreases in the bulk considerably and more bubbles stick to the inner cylinder. These bubbles play a major role in influencing the dynamics of plume ejection and also near-wall viscous dissipation (
$\hat{{\it\epsilon}}$
) as seen in figures 7 and 10, respectively. While in the case of
$Fr<1$
, the disruption of plumes and its consequence on weakening the Taylor rolls results in low-strain-rate regions in the plume impacting regions, there is no such effect when
$Fr>1$
, i.e. when the buoyancy of the bubbles is weaker.
Similar to figure 13 we also plot the axial distribution of bubbles in the axial direction in figure 14. The local volume fraction is averaged in the azimuthal, radial direction and in time. For strong buoyant bubbles (
$Fr=0.16$
), the axial distribution displays a largely uniform behaviour while for the increased
$Fr$
number there are two distinct peaks in the local volume fraction of the bubbles for both
$Re_{i}$
numbers. The positions of these peaks correspond to the plume ejection sites (i.e. the outflow region) on the inner cylinder wall and are an effect of bubbles accumulating near these sites, before they are advected into the Taylor vortices. This distribution of bubbles when
$Fr>1$
is similar to the spiral pattern of bubbles observed in the experimental results of Murai et al. (Reference Murai, Oiwa and Takeda2005, Reference Murai, Oiwa and Takeda2008) when the Froude number is higher than one and the operating Reynolds number of the inner cylinder is
$Re_{i}=2100{-}4500$
.
4 Summary
The drag reduction in turbulent bubbly TC flow compared with the corresponding single-phase system depends on various control parameters such as the operating Reynolds number, the diameter of the bubbles, the relative strength of buoyancy compared with the driving and the strength of coherent vortical structures. It is extremely challenging to independently control each of these parameters in experiments and study its effect on the overall drag reduction. In this work we perform numerical simulations of a two-phase TC system using an Euler–Lagrange approach which allows us to track almost 50 000 bubbles simultaneously up to an operating Reynolds number of
$Re_{i}=8000$
and also gives us the freedom to control various parameters independently. We find that the relative strength of the buoyancy as compared to the inner cylinder driving (quantified by the Froude number
$Fr$
) is crucial to achieve drag reduction at higher Reynolds number. When the Froude number
$Fr$
is less than one, drag reduction was observed up to
$Re_{i}=8000$
, the largest
$Re_{i}$
we studied here. When the
$Fr$
number is larger than one, almost negligible drag reduction was observed for all
$Re_{i}$
. We also observe from the averaged azimuthal velocity fields that for a fixed
$Re_{i}=2500$
, the coherent Taylor roll is much weaker for
$Fr=0.16$
due to the rising strongly buoyant bubbles, as compared with
$Fr=1.28$
. The weakened Taylor rolls in turn imply drag reduction. The effect of the bubbles on the Taylor rolls is also observed in the contour plots of the azimuthal velocity fluctuations. By analysing the individual trajectory of the bubbles, we showed that while for a low
$Re_{i}=2500$
the bubbles have primarily either a upward motion (lower
$Fr$
) or coherent motion inside the Taylor roll (larger
$Fr$
). Similar behaviour was observed even in the case of a higher
$Re_{i}$
, but with a more erratic movement.
In this work we have studied numerically the effect of small bubbles (
$d_{b}/{\it\delta}_{{\it\nu}}=0.55{-}1.25$
) on the dynamics of a turbulent TC system using a Euler–Lagrange approach. The next line of simulations should focus on finite size bubbles which are much larger than the Kolmogorov scale. It has already been shown by experiments that almost 40 % drag reduction can be achieved with just 4 % volume fraction of bubbles (van Gils et al.
Reference van Gils, Narezo Guzman, Sun and Lohse2013). However, the Euler–Lagrange model is limited to the use of sub-Kolmogorov and Kolmogorov scale bubbles while techniques such as front tracking, immersed boundaries, etc., are limited by the number of bubbles that can be simulated in a highly turbulent field. Simulation techniques which can model the deformation properties of the dispersed phase using subgrid models are required in order to successfully advect millions of particles/droplets/bubbles in a highly turbulent carrier phase.
Acknowledgements
We would like to thank Y. Yang and E. P. van der Poel for various stimulating discussions during the development of the code. This work was supported by the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation programme funded by the Ministry of Education, Culture and Science of the government of the Netherlands and the FOM-CSER program. We acknowledge PRACE for awarding us access to FERMI based in Italy at CINECA under PRACE project number 2015133124 and NWO for granting us computational time on Cartesius cluster from the Dutch Supercomputing Consortium SURFsara.
Appendix
In table 3 we give the data for the Nusselt numbers used in figure 5 and also reference data for the simulations with a larger aspect ratio (
${\it\Gamma}=8$
).
Table 3. Nusselt numbers for the single-phase and two-phase cases for different
$Re_{i}$
and
$Fr$
numbers. Last two columns show the Nusselt numbers with double the aspect ratio
${\it\Gamma}=8$
.
