2.1. Governing equations

We consider two immiscible and incompressible liquids confined between two coaxial cylinders whose radii are $r_i$ (inner) and $r_o$ (outer), and the axial length is $L_z$. In this work, we fix the outer cylinder and let the inner one rotate with a constant angular velocity $\omega _i$. The two liquid phases are governed by the incompressibility constraint

(2.1)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 \end{equation}

and the balance of momentum

(2.2)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\frac{1}{\rho}\,\boldsymbol{\nabla} p + \nu\,\nabla^2 \boldsymbol{u} + \frac{1}{\rho}\,\boldsymbol{f}. \end{equation}

Here, $\boldsymbol {u} = \boldsymbol {u} ( \boldsymbol {x}, t )$ and $p = p ( \boldsymbol {x}, t )$ denote the fluid velocity and the reduced pressure, respectively. Note that $\rho$ and $\nu$, which are the density and kinematic viscosity, respectively, have the same values for both liquids in this work, and thus take constant values. Also, $\boldsymbol {f} = \boldsymbol {f} ( \boldsymbol {x}, t )$ describes the interfacial contribution to the momentum balance. It is written as

(2.3)\begin{equation} \boldsymbol{f} = \sigma\,\kappa ( \boldsymbol{x}, t )\,\delta \boldsymbol{n} ( \boldsymbol{x}, t ), \end{equation}

where $\sigma$, $\kappa$, $\delta$ and $\boldsymbol {n}$ are the surface tension coefficient, local curvature, Dirac delta function and local normal vector, respectively. In addition to those two equations, we consider an advection equation with respect to the indicator function $H$,

(2.4)\begin{equation} \frac{\partial H}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} H = 0, \end{equation}

to distinguish the two liquids. The indicator function $H$ equals $0$ where the primary phase is, and $1$ if the region is occupied by the secondary phase. Equations (2.1) and (2.2) are solved by a second-order-accurate finite-difference scheme in cylindrical coordinates (Verzicco & Orlandi Reference Verzicco and Orlandi1996; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015).

2.2. Volume-of-fluid method

In this subsection, we describe briefly the method employed to integrate (2.4) and to solve the interfacial contribution $\boldsymbol{f}$ in (2.2). Here, $\boldsymbol{f}$, which is described mathematically as in (2.3), is modelled with the continuum-surface-force approach (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992; Ii et al. Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012)

(2.5)\begin{equation} \boldsymbol{f} \approx \sigma \kappa\,\boldsymbol{\nabla} \phi, \end{equation}

where $\phi$ is the so-called ‘volume-of-fluid’ defined as the volume average of $H$ in a control volume $\mathcal {V}$. By averaging (2.4) in $\mathcal {V}$, we obtain the advection equation with respect to $\phi$ as

(2.6)\begin{equation} \frac{\partial \phi}{\partial t} ={-}\int_{\partial \mathcal{V}} \boldsymbol{u} H \boldsymbol{\cdot} \boldsymbol{n} \,\mathrm{d}\mathcal{S} + \phi\,\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}, \end{equation}

where $\int _{\partial \mathcal {V}} (\cdot )\,\mathrm {d}\mathcal {S}$ denotes the surface integral on the control volume. To compute this surface integral, we extended the THINC/QQ scheme (Xie & Xiao Reference Xie and Xiao2017) to cylindrical coordinates (see Appendix B for details).

2.3. Flow parameters

For the TC set-up, we set the curvature as $\eta \equiv r_i / r_o = 0.714$ and the aspect ratio as $\varGamma \equiv L_z / d = 2{\rm \pi}$, where $d \equiv r_o - r_i$ is the gap width. No-slip and impermeable boundary conditions are imposed in the radial direction, while periodicity is imposed in the axial direction. As shown in figure 1, we employ the full azimuthal cylinder domain because large fluid structures composed of the secondary phase can exist. The cylinder is discretised by $( 128, 1920, 720 )$ grid points in the radial, azimuthal and axial directions, respectively. The clipped-Chebyshev clustering is adopted in the radial direction to resolve boundary layers near the walls, while grid points are uniformly spaced in the other homogeneous directions. Details of the spatial resolution effects and code validations and verifications are described in §§ C.2 and C.5.

Figure 1.

Instantaneous interface snapshots for different total volume fractions $\varphi$ and $We$. The iso-surface $\phi = 0.5$ is shown, where the regions surrounded by the bluish and reddish surfaces are filled by the primary ($\phi < 0.5$) and secondary ($\phi > 0.5$) phases, respectively. Two co-axial transparent cylinders denote the boundaries of the geometry. Weber numbers (a–e) $We = 400$ (high $We$), and (f–j) $We = 40$ (low $We$). The volume fraction $\varphi$ is increased from left to right, from $10\,\%$ to $50\,\%$.

We fix the Reynolds number $Re \equiv r_i \omega _i d / \nu = 960$, which is in the steady Taylor vortex regime (Grossmann et al. Reference Grossmann, Lohse and Sun2016). The two control parameters characterising the secondary phase are its total volume fraction $\varphi$ and the (system) Weber number $We \equiv \rho r_i^2 \omega _i^2 d / \sigma$. We consider ten different total volume fractions $\varphi = 5\,\%, 10\,\%, 15\,\%, \ldots, 50\,\%$, and two system Weber numbers $We = 40, 400$. In addition, we impose the Neumann boundary condition with respect to $\phi$, resulting in a $90^\circ$ contact angle for the secondary phase.

It should be noted that as discussed in the Introduction, there are several parameters that are expected to play crucial roles in the fluid dynamics but are missing in this work. We consider that they are the density ratio between two fluids, and the Reynolds number, among others. When the density ratio between two liquids is not unity, one would expect that the phase having smaller density tends to migrate to the inner cylinder because of the centrifugal force, affecting the phase distributions. When the Reynolds number is varied, interactions between interfacial structures and different turbulent intensities could pose many interesting phenomena. However, even for a fixed and relatively small Reynolds number with unitary density ratio, interplay between flow fields (in particular Taylor rolls) and the surface is rarely studied and is worth being examined. Although this is the reason why we do not consider these parameters in this work, they are to be analysed in the near future.

First, we simulate a single-phase case to initialise the velocity field. Once we obtain a well-developed flow having two pairs of Taylor rolls, the simulation is restarted after some spherical droplets with diameter $0.08 d$ are positioned randomly. The droplets coalesce (because of the contact with different droplets) and break up (because of shear) continuously, and eventually adapt to the flow field. All the statistics presented are collected for at least $100$ time units after the statistically steady state is reached, which is determined by monitoring the time evolution of the Nusselt number (whose variation is within $\approx 4\,\%$) and the interface area (see § C.5).

4.1. Modulated Taylor roll structures

Given that the basic flow is in the steady Taylor vortex regime, one natural question is whether the Taylor rolls persist when a secondary phase is introduced. To answer this question, in figure 4, we show the temporally and spatially averaged velocity field in the $r$–$z$ plane. Again, the reader is referred to Appendix A for cases omitted here. Contours indicate the magnitude of the azimuthal velocity $\langle u_{\theta } \rangle _{\theta,t}$, while vectors represent the meridional velocity components $( \langle u_r \rangle _{\theta,t}, \langle u_z \rangle _{\theta,t} )$. As a baseline for comparison, the leftmost plot shows the result for the single-phase case, where two pairs of Taylor rolls are clearly visible. Four ‘plumes’ (two ejecting and two impacting ones from each cylinder wall) are visible where the azimuthal momentum is actively transported in the radial direction.

Figure 4.

Contours of the averaged azimuthal velocity $\langle u_{\theta } \rangle _{\theta,t}$ and arrows of the averaged radial–axial velocity vectors $\langle u_{r,z} \rangle _{\theta,t}$ in the $r$–$z$ plane. An arrow having half the length of one image width corresponds to ${|\langle u_{r,z} \rangle | = 1}$. The leftmost plot (labelled ‘Single-phase’) denotes the reference single-phase result, while the other plots are for two-phase results. (a–e) Results for high $We = 400$ from $\varphi = 10\,\%$ to $\varphi = 50\,\%$ in steps of $10\,\%$. (f–j) Results for low $We = 40$.

For high $We = 400$ (figure 4a–e), the Taylor rolls are less affected by the presence of the secondary phase. This is mainly because the secondary phase is easily fragmented into droplets that are easily strained by the Taylor rolls.

On the other hand, the flow fields are most affected for low $We = 40$ since the secondary phase is more compact and does not break up easily. Consequently, the resulting velocity modifications are noticeably different for different $\varphi$. At $\varphi = 10\,\%$, although a few large droplets are observed, two pairs of Taylor rolls persist in the flow field. This is to be expected since the volume fraction of the secondary phase is insufficient to create large-scale structures as observed in higher $\varphi$ cases described later. However, this secondary phase structure suppresses the fluid motion, which reduces the circulation of the upper Taylor roll pair (see the bottom half of figure 4(f) with smaller velocity vectors).

For intermediate volume fractions $\varphi = 20\unicode{x2013}40\,\%$, the creation of toroidal structures of the secondary phase attaching to the inner cylinder is observed. By comparing the phase distributions (figure 3g–i) with the corresponding flow fields (figure 4g–i), we notice that smaller sub-rolls (circulatory regions) are formed inside these confined spaces. At the highest volume fraction $\varphi = 50\,\%$, a pair of Taylor rolls can exist separately in both phases because of layering (figure 3j), resulting in a velocity pattern similar to the reference single-phase case. Similar layers are observed also at $\varphi = 40\,\%$, where a pair of narrow Taylor rolls is evolved in between.

Overall, we observe that the flow fields are largely modulated when the secondary phase attaches to the inner cylinder. This is to be expected since at the inner wall, the interface forces the boundary layer to separate, and as a consequence, the base Taylor rolls are modulated by the plume ejecting at this point. Thus it can be concluded that the Taylor vortices are significantly modified at low $We$ due to the formation of large coherent secondary phase structures that adhere to the inner cylinder. We note that this scenario can be changed drastically when the boundary condition is changed (e.g. when the walls are super-oleophobic).

4.2. Taylor rolls and velocity spectra

In the previous subsections, we explained how the interface modulates the velocity field, especially the Taylor vortices, through flow visualisations. In order to get deeper insights into the effect of the interface, we now add a quantitative discussion on the velocity field.

The Taylor vortices are well characterised by axially coherent radial jets, which correspond to plume ejecting (or impacting) regions. To characterise these repeating structures, we consider the one-dimensional spectra of the radial velocity in the axial direction at the middle of the channel, $r = ( r_i + r_o ) / 2$, namely $E_{rr}$. They are defined as $\langle u_r^2 \rangle _{\mathcal {A}} \equiv \int _0^\infty E_{rr}\,\mathrm {d}k_z$, where $k_z$ is the axial wavenumber, and $\langle {\cdot } \rangle _{\mathcal {A}}$ denotes averaging in the $\theta \unicode{x2013}z$ plane. We show these spectra in figure 5 for both Weber numbers and various concentrations $\varphi$.

Figure 5.

Temporally averaged one-dimensional spectra of the radial velocity in the axial direction $E_{rr}$ at $r = ( r_i + r_o ) / 2$, for (a) $We = 400$, and (b) $We = 40$. Grey lines show the reference single-phase result with two pairs of Taylor rolls. Different colours are used to distinguish different $\varphi$ cases: $10\,\%$ (red), $20\,\%$ (orange), $30\,\%$ (green), $40\,\%$ (light blue) and $50\,\%$ (dark blue).

When focusing on the reference single-phase case (grey lines in figure 5), we observe that the dominant mode is $k_z = 2$, indicating that the single-phase flow has two pairs of Taylor rolls extending $L_z / 2$ axially. Even wavenumbers, which are the harmonics of the dominant wavenumber $k_z = 2$, have finite and monotonically decreasing energy with increasing $k_z$. On the other hand, odd wavenumbers have zero energetic contributions, which when combined with the even wavenumber contribution, result in typical sawtooth patterns. We note that although this sawtooth trend was already observed in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Lohse and Verzicco2016) for $Re_i = 3.4 \times 10^4$, it is more pronounced here because of the much lower Reynolds number ($Re_i = 960$) considered here. Since the Taylor rolls are quite stable and little velocity fluctuations exist in our case, almost no energy is transferred through the nonlinearity, which is different from the turbulent regime considered in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Lohse and Verzicco2016).

In the two-phase cases, these sawtooth features are still present, but the spectra are distributed among not only even $k_z$ but also odd $k_z$. In addition, the sharp drop-off at higher wavenumbers for the single-phase case is replaced by broader-tailed decays. The energetic contributions at larger wavenumbers imply that the secondary phase induces many more small-scale fluctuations and injects energy to higher wavenumbers.

When $We$ is high ($We = 400$, figure 5a), the dominant wavenumber remains $k_z = 2$, regardless of the total volume fraction. The secondary phase is prevented from coalescing into large-scale structures (which can affect the Taylor rolls) because of the weak surface tension. The single-phase sawtooth pattern is weakened as $\varphi$ is increased, which again shows the energy redistribution effect of the secondary phase.

When considering the lower $We$ case $We = 40$ (figure 5b), on the other hand, we observe a different scenario. For low volume fraction $\varphi = 10\,\%$, although a long-tailed decay is present, the dominant wavenumber is still $k_z = 2$, which reflects the fact that the volume fraction is not large enough to affect the Taylor roll structures. When $\varphi$ is in the intermediate regime ($\varphi = 20\,\%$ and $30\,\%$), toroidal interfacial structures and sub-roll velocity circulations are observed. These structures have shorter length scales in the axial direction than the original Taylor rolls, which is reflected as a shift in the dominant wavenumber from $k_z = 2$ to $k_z = 4$. A similar explanation can be used for $\varphi = 40\,\%$, where a thin secondary phase layer and circulatory velocity inside can be seen. These confined Taylor rolls also have smaller length scale than the original rolls, resulting in the dominance of $k_z = 4,5,6$. For $\varphi = 50\,\%$, the dominant wavenumber again goes back to $k_z = 2$. This is related to the layering observed in figure 4(j), where a pair of Taylor rolls exists in each phase and thus the large-scale velocity field becomes similar to the single-phase one.

5.1. Global response

In the previous sections, our focus was mainly on the flow organisation (i.e. velocity fields and interface structures). Here, we quantify and explain how the global response is changed by the flow with the secondary phase.

The primary global response of the TC flow is the torque $T$, which is given as a product of the cylinder radius and the surface integral of the azimuthal shear stress. In a non-dimensional form, we can define the Nusselt number of the angular velocity transport $Nu_{\omega } \equiv T/T_{lam}$ in analogy with the heat flux enhancement in Rayleigh–Bénard convection (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007), where $T_{lam}$ is the torque for the purely azimuthal laminar flow (Sugiyama et al. Reference Sugiyama, Calzavarini and Lohse2008),

(5.1)\begin{equation} \frac{4 {\rm \pi}\rho \nu L_z}{1 - \eta^2}\,r_i^2 \omega_i, \end{equation}

which is derived by integrating the $r-\theta$ component of the shear stress tensor on the inner (or outer) cylinder wall.

In figure 6(a), we show how $Nu_{\omega }$ varies with the volume fraction $\varphi$. For the lower $We = 40$, which is shown in red, $Nu_{\omega } ( \varphi )$ shows a strongly non-monotonic dependence. Starting from $\varphi = 0\,\%$, $Nu_{\omega }$ remains almost constant up to $\varphi = 10\,\%$, which is followed by an increasing region up to $\varphi = 40\,\%$. After taking the maximum value there, it suddenly drops at $45\,\%$ and recovers the single-phase value at $50\,\%$. Overall, $Nu_{\omega }$ takes the same or larger values compared to the single-phase flow. For $We = 400$, which is shown in blue, we observe a reduction in $Nu_{\omega }$ when the secondary phase is present; $Nu_{\omega } ( \varphi )$ maintains the constant value for $\varphi = 5 - 40\,\%$, and increases monotonically for $\varphi = 40\unicode{x2013}50\,\%$.

Figure 6.

(a) Angular velocity transport Nusselt number $Nu_{\omega }$, and (b) interfacial surface area $\mathcal {S}_{int} / \mathcal {S}_{cyl}$ as a function of the total volume fraction $\varphi$. Red and blue colours represent the results for $We = 40$ and $400$, respectively. The error bars indicate the magnitude of the temporal fluctuation of the signals (i.e. $1 \sigma$).

One might be tempted to explain these non-monotonic trends by inspecting the interfacial surface area (Rosti et al. Reference Rosti, De Vita and Brandt2019a). In figure 6(b), we plot the interfacial surface area $\mathcal {S}_{int}$ normalised by the inner cylinder surface area $\mathcal {S}_{cyl}$ versus $\varphi$. Here, $\mathcal {S}_{int} \equiv \sum _{0.4 < \phi < 0.6} r\,\Delta \theta \,\Delta z$, where $\sum _{0.4 < \phi < 0.6}$ denotes the summation over cells satisfying the condition $0.4 < \phi < 0.6$. We observe that the high $We$ case takes larger interfacial surface areas than the low $We$ case, which is because of the larger interfacial deformability and more fragmented droplets. Overall, the interfacial surface area increases with increasing $\varphi$ for both $We$ values, which clearly indicates the physical picture that a larger interfacial area contributing to larger $Nu_{\omega }$ does not apply to our situation here.

In figure 6(a), we confirm several issues to be explained regarding $Nu_{\omega } ( \varphi )$: (i) an increasing trend for low $We$ at $\varphi = 10\unicode{x2013} 40\,\%$; (ii) a sudden reduction for low $We$ at $\varphi = 45\,\%$; and (iii) an overall reduction for high $We$. The physical explanations for these numerical findings will be elaborated in the next subsection by inspecting the flow visualisations and quantifying the different contributions to $Nu_{\omega }$.

5.2. Nusselt number decomposition

In the previous subsection, we showed that the $Nu_{\omega } ( \varphi )$ dependence is strongly non-monotonic for both $We$ values, and trends in the interfacial surface area $\mathcal {S}_{int}$ cannot explain these results. In this subsection, we reveal the cause of the complex dependencies of $Nu_{\omega }$ on $\varphi$ and $We$ by decomposing $Nu_{\omega }$ into three contributions.

This decomposition is inspired by an analogous situation in turbulent channel flows, where it is conventional to decompose the total shear stress into the advective and diffusive contributions (Pope Reference Pope2000, Chapter 7), by which one can evaluate the share of each term in the total shear stress. This idea has also been extended to various wall-bounded multi-phase flows by incorporating the additional terms participating in the momentum exchange, e.g. particulate flows (Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015) and two-liquid flows (De Vita et al. Reference De Vita, Rosti, Caserta and Brandt2019).

The corresponding conserved quantity in TC flows is the angular velocity flux $J_{\omega }$ (Eckhardt et al. Reference Eckhardt, Grossmann and Lohse2007) defined as

(5.2)\begin{equation} J_{\omega} \equiv J_{\omega,{adv}} (r) + J_{\omega,{dif}} (r) + J_{\omega,{int}} (r) = \text{const.}, \end{equation}

where the three terms represent

1. (i) the advective contribution (red line in figure 16) $J_{\omega, {adv}} ( r ) \equiv r^3 \langle u_r \omega \rangle _{\theta,z,t}$,

2. (ii) the diffusive contribution (blue line in figure 16) $J_{\omega, {dif}} ( r ) \equiv -\nu r^3\,\partial \langle \omega \rangle _{\theta,z,t} / \partial r$, and

3. (iii) the interfacial contribution (green line in figure 16) $J_{\omega, {int}} ( r ) \equiv -(1/\rho ) \int _{r_i}^r {r^{\prime }}^2$ $\langle \,f_{\theta } \rangle _{\theta,z,t} \,\mathrm {d} r^{\prime }$.

Examples showing that the sum of $J_{\omega, {adv}}$, $J_{\omega, {dif}}$ and $J_{\omega, {int}}$ takes a constant value across the channel are described in Appendix D. We note that $\omega$ is the angular velocity, which can be written as $\omega = u_{\theta } / r$. Averaging (5.2) in the radial direction and normalising by the single-phase laminar value $J_{\omega, {lam}}$ yields

(5.3)\begin{equation} \mbox{{Nu}}_{\omega} = \mbox{{Nu}}_{\omega,{adv}} + \mbox{{Nu}}_{\omega,{dif}} + \mbox{{Nu}}_{\omega,{int}}, \end{equation}

through which we can separate the three contributions to the global response $Nu_{\omega }$. We note that each contribution is a function of the radial position, which is averaged here.

In figure 7(a,b), we show the contributions $Nu_{\omega, {adv}}$, $Nu_{\omega, {dif}}$ and $Nu_{\omega, {int}}$ as functions of $\varphi$ for the higher and lower $We$ cases, respectively. Note that the black lines, which are the summations of three contributions, are identical to $Nu_{\omega }$ shown in figure 6(a).

Figure 7.

Decomposed $Nu_{\omega }$ for (a) $We = 400$, and (b) $We = 40$, as functions of the total volume fraction $\varphi$. Different colours are used to distinguish different contributions: red for advective $Nu_{\omega, {adv}}$, blue for diffusive $Nu_{\omega, {dif}}$, and green for interfacial $Nu_{\omega, {int}}$. The black lines show the sums of the three contributions, which are identical to the lines in figure 6(a). (c) Averaged angular velocity $\langle \omega \rangle _{\theta,z,t}$ normalised by the angular velocity of the inner cylinder $\omega _i$ as a function of the normalised radial position $( r - r_i ) / d$. Different colours are used to distinguish different volume fractions, whose usage is same as in figure 5, whereas different line styles are adopted to identify $We$, namely high $We = 400$ (dashed) and low $We = 40$ (dotted). (d) Ratios of $Nu_{\omega, {int}}$ to $Nu_{\omega, {adv}}$ as functions of the total volume fraction $\varphi$. Red and blue colours represent the results for $We = 40$ and $400$, respectively.

For all cases, including the reference single-phase case, the diffusive contributions $Nu_{\omega, {dif}}$ (blue lines in figure 7(a,b)) give similar values, i.e. $Nu_{\omega, {dif}}$ is less affected by the variations in Weber number and volume fraction of the secondary phase. In order to understand this finding, in figure 7(c) we show the averaged angular velocity profile in the radial direction, $\langle \omega \rangle _{\theta,z,t} ( r )$, which appears in the formulation of the diffusive contribution. As elaborated in § 4.1 and shown in figure 4 as contours, when the Weber number and the volume fraction are varied, azimuthal velocity fields $\langle \omega \rangle _{\theta,t}$ in the radial–axial plane show substantial differences to each other because of the modulations caused by the interface. When we further average in the axial direction, however, $\langle \omega \rangle _{\theta,z,t}$ shows similar monotonic reduction from the inner to the outer cylinders overall, and the deviations are less pronounced, giving an almost identical value for $Nu_{\omega, {dif}}$. This finding can be linked to the physical nature of the diffusive term, which is closely connected to the fluid viscosity that is fixed in this work rather than the flow structures. We note that this is not the case if the viscosities of the two phases do not match, where the phase distribution changes the local diffusivity, thus the diffusive contribution is a response of the system.

The interfacial contribution $Nu_{\omega, {int}}$ (green lines in figure 7(a,b)) is of course zero for the single-phase case, and takes positive values when the secondary phase exists. The amount of increase is different for the two $We$ cases that we consider; The contribution to $Nu_{\omega }$ is less than $10\,\%$ for $We = 400$, while it is from approximately $10\,\%$ to $30\,\%$ for $We = 40$, which reflects the different interface stiffness.

The advective contribution $Nu_{\omega, {adv}}$ (red lines in figure 7(a,b)) is the most dominant in the parameter space under investigation. For the single-phase case, $75\,\%$ of the total originates from this term, and for the two-phase case, it is more than $50\,\%$. A similar analysis was performed by Chiara et al. (Reference Chiara, Rosti, Picano and Brandt2020) for particles in plane channel flow at comparable Reynolds numbers. However, these authors found that the advective contribution is much smaller than the diffusive contribution, which is opposite to our findings in TC flow. This shows that the enhancement of the momentum transport in the radial direction by the Taylor rolls is significant, highlighting the important role of curvature in TC flows.

In the presence of the secondary phase, we observe that the advective contributions $Nu_{\omega, {adv}}$ take lower values than the reference single-phase value, which can be understood as follows. Since $Nu_{\omega, {adv}}$ is a function of $u_r \omega$ (see (5.2)), it is reduced if the magnitudes of the velocities get smaller. This reduction is indeed implied in the mean velocity fields (figure 4), where velocity vectors become shorter than in the single-phase case. This is to be expected since the interface requires energy to be deformed, and as a consequence attenuates the convective momentum transport by suppressing the fluid motions.

Based on the observations for these three contributions, now we work out the reason why the two $We$ cases show opposite trends, i.e. $Nu_{\omega } ( \varphi )$ is increased for $We = 40$ compared to the single-phase value, while for $We = 400$ it shows a reduction. Main roles are played by the advective contribution $Nu_{\omega, {adv}} ( \varphi )$ and the interfacial contribution $Nu_{\omega, {int}} ( \varphi )$. In the presence of the secondary phase, compared to the single-phase flow, $Nu_{\omega, {adv}} ( \varphi )$ is reduced, which is observed for both $We$ values and for all volume fractions. On the other hand, because of the presence of deformable interfaces, $Nu_{\omega, {int}} ( \varphi )$ is added and contributes positively to the total $Nu_{\omega } ( \varphi )$. For $We = 400$, the amounts of increase in $Nu_{\omega, {int}} ( \varphi )$ are marginal, resulting in the net reduction in $Nu_{\omega } ( \varphi )$. For lower $We = 40$, on the other hand, substantially large increases in $Nu_{\omega, {int}} ( \varphi )$ compensate for the reductions in $Nu_{\omega, {adv}} ( \varphi )$, leading to the overall increases in $Nu_{\omega } ( \varphi )$.

Based on the decomposition results as well as on the flow visualisations (figures 3 and 4), we can now give a physical explanation of the non-monotonic $Nu_{\omega }$ trends for the lower $We$ cases.

For $\varphi = 20\unicode{x2013}30\,\%$, although the whole volume fraction is varied, $Nu_{\omega } ( \varphi )$ and its three contributions $Nu_{\omega, {adv}}$, $Nu_{\omega, {dif}}$, $Nu_{\omega, {int}}$ are quite similar to each other. This is expected when we look into their phase distributions (figure 3(g,h); only $20\,\%$ and $30\,\%$ are shown), where similar toroidal structures attaching to the inner cylinder are present. Although their volumes are larger at $\varphi = 30\,\%$, the radii of the torus are hardly changed, resulting in similar interfacial patterns and velocity fields, and thus giving similar $Nu_{\omega } ( \varphi )$.

At $\varphi = 40\,\%$, the volume fraction is large enough to create a region covering the whole gap (a narrow secondary phase region in figure 3i), which connects the inner and outer boundary layers, and another pair of Taylor rolls is formed inside (see figure 4i). This structure increases the number of plumes (see figure 5a) and enhances the radial momentum transport, resulting in the increase in $Nu_{\omega, {adv}} ( \varphi )$ (see figure 7b).

For $\varphi$ larger than $45\,\%$, $Nu_{\omega, {int}}$ shows a drastic reduction despite the small change in the interfacial surface area $\mathcal {S}_{int}$ (see figure 6b). The reason for this is that the interfacial term in the momentum equation $\,f_{\theta }$ is defined as $\sigma \kappa \delta n_{\theta }$, i.e. only the surface area normal to the azimuthal direction characterised by $\delta n_{\theta }$ contributes to $J_{\omega, {int}}$. As the phase distribution indicates (figure 3(j); only $50\,\%$ is shown), the two phases are clearly layered in the axial direction, and ring-like surfaces are formed in between, which have a small surface area normal to the azimuthal direction, leading to a much lower $Nu_{\omega, {int}}$ even though $\mathcal {S}_{int}$ is sufficiently large.

In summary, when $We$ is high, the Taylor rolls, which are present when the secondary phase is absent, clearly remain since interfaces are easily deformed and they do not affect the flow field much. The system response $Nu_{\omega }$ is governed mainly by the advection, which can be seen as the smaller values of $Nu_{\omega, {int}}$ compared to $Nu_{\omega, {adv}}$ in figure 7(d), i.e. the advection-dominated regime. For lower $We$, on the other hand, interfaces tend to create large-scale structures since they are less deformable and fragmented. They change the Taylor roll structures by interfering with the boundary layers, and sometimes affect the radial momentum transport by changing the number of Taylor rolls. In this case, the system response $Nu_{\omega }$ is governed by the interfacial structures, which is characterised as the comparable contribution of $Nu_{\omega, {int}}$ to $Nu_{\omega, {adv}}$ in figure 7(d); i.e. this is the interface-dominated regime.