1. Introduction
Two-phase liquid–liquid flows, often referred to as emulsions, involve the presence of two immiscible liquids and play a significant role in various natural and industrial processes. Turbulent emulsions, in particular, introduce further complexity due to the presence of dispersed phases. The interaction between the continuous phase and the dispersed phase leads to additional complexities, resulting in a wide range of phenomena and dynamics influenced by the properties of the phases and the explored parameter regimes (Lemenand et al. Reference Lemenand, Della Valle, Dupont and Peerhossaini2017; Yi et al. Reference Yi, Wang, Huisman and Sun2023). In recent decades, turbulent emulsions have found extensive applications in fields such as petroleum, food, pharmaceuticals and cosmetics (Spernath & Aserin Reference Spernath and Aserin2006; Mcclements Reference Mcclements2007; Wang et al. Reference Wang, Li, Zhang, Dong and Eastoe2007; Kilpatrick Reference Kilpatrick2012), garnering significant interest (Rosti, Brandt & Mitra Reference Rosti, Brandt and Mitra2018; Yi et al. Reference Yi, Wang, Huisman and Sun2023; Ni Reference Ni2024). However, our understanding of how turbulence is influenced by the dispersed phase in turbulent emulsions remains limited.
Due to the presence of the two-phase interface and disparities in liquid properties between the two phases, experimental observation of behaviour in turbulent emulsions becomes particularly challenging. Current experimental research primarily focuses on microscopic droplet formation, size distribution and the macroscopic response of global transport (Bakhuis et al. Reference Bakhuis, Ezeta, Bullee, Marin, Lohse, Sun and Huisman2021; Yi, Toschi & Sun Reference Yi, Toschi and Sun2021; Wang et al. Reference Wang, Yi, Jiang and Sun2022b; Yi et al. Reference Yi, Wang, van Vuren, Lohse, Risso, Toschi and Sun2022). Meanwhile, related simulations have been employed to investigate droplet size distribution, primarily in homogeneous and isotropic turbulence, with a specific emphasis on how turbulence affects droplet breakup behaviour (Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Vela-Martín & Avila Reference Vela-Martín and Avila2022). The feedback of droplets on turbulence has recently been studied in homogeneous shear turbulence. Dodd & Ferrante (Reference Dodd and Ferrante2016) investigated how droplet deformation, breakup and coalescence affect in the temporal evolution of turbulent kinetic energy (TKE). They showed that droplet coalescence reduces the total interfacial surface area, causing a decrease in surface energy and an increase in local kinetic energy. The presence of droplets acts as a sink in the TKE of the bulk fluid, as the dispersed phase was found to slow down the dissipation of TKE compared with the continuous phase (Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019). The effect of droplets and the role of their viscosity on turbulence in homogeneous and isotropic turbulent flows have also received recent attention. The energy is transported consistently from large to small scales by the two-phase interface, and the total interface area is directly proportional to the amount of energy transported. Increasing the dispersed phase viscosity would reduce the amount of energy being transported (Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022). Correspondingly, large velocity gradients are found across the two-phase interface and will gradually disappear as the viscosity of the dispersed phase increases (Farsoiya et al. Reference Farsoiya, Liu, Daiss, Fox and Deike2023). These findings provide a deeper understanding of the effect of droplets on homogeneous turbulence. However, since the studied systems are unbounded, the results cannot be applied directly to wall-bounded turbulence, where strong inhomogeneity and anisotropy could develop due to the boundary layer. This introduces the potential for different observations in both phases.
In contrast to the relatively uniform turbulence dissipation observed throughout different regions in homogeneous turbulence, wall-bounded turbulence exhibits a distinct characteristic where approximately half of the dissipation occurs predominantly in the immediate vicinity of walls (Jiménez Reference Jiménez2012). Recent studies in wall-bounded laminar flows have investigated the modulation of the two-phase interface on the system's drag by allowing droplets to coalesce numerically in the Taylor–Couette (TC) device (Hori et al. Reference Hori, Ng, Lohse and Verzicco2023) or prohibiting droplet coalescence in planar Couette flow (De Vita et al. Reference De Vita, Rosti, Caserta and Brandt2019). The drag enhancement caused by the dispersed phase is attributed to the interfacial contribution, and coalescence could effectively decrease the interfacial area, thus weakening the drag enhancement effect. To investigate drag modulation by different dispersed phases, our previous work (Su et al. Reference Su, Yi, Zhao, Wang, Xu, Wang and Sun2024) examined turbulence modulation induced by dispersed phases with varying density and viscosity through three-dimensional direct numerical simulations in turbulent TC flow. We derived a momentum transport formula, revealing that the two-phase interface consistently enhances drag. Reducing the density or viscosity of the dispersed phase decreases the contributions of the advection and diffusion terms, leading to reduced drag. However, that work only encompassed the effects of the dispersed phase on global properties, leaving an unexplored area in understanding the effect on the local statistical properties of turbulent flow. The effect of the dispersed phase on the production and dissipation of turbulence remains unknown. This motivates further exploration in this new work to gain a comprehensive understanding of the modulation of the transport of TKE by examining the effects of the dispersed phase on the local statistical properties.
In this work, we conducted a comprehensive study on the modulation of statistical properties of turbulence induced by dispersed phases in the TC system. To achieve this, we employed an interface-resolved volume-of-fluid (VOF) method, which allows us to resolve the interface between the two phases and solve the governing equations in a single-equation formulation. This approach enables us to perform operations similar to those in single-phase flow, facilitating effective comparisons and determining the specific effects of the dispersed phase.
We aim to explore the behaviour of turbulent emulsions in a semi-dilute regime, with a specific focus on the potential impact of droplets possessing lower density and viscosity compared with the continuous phase. Through numerical simulations, we can disentangle the effects of the two-phase interface, the density and the viscosity of the dispersed phase on turbulence modulation, uncovering their intertwined coupling. The article is organised as follows. In § 2, the numerical method and settings are described. In § 3, turbulence modulation is discussed based on momentum budget analysis, turbulence fluctuation analysis and TKE budget analysis. Finally, conclusions are drawn in § 4.
2. Numerical method and setting
We conducted interface-resolved three-dimensional direct numerical simulations to investigate the two-phase fluid–fluid turbulence in a TC system. These simulations were performed using a VOF method with a piecewise-linear interface calculation (PLIC) algorithm implemented in the interFoam solver of the open-source OpenFOAM v8 (Rusche Reference Rusche2003; Chen, Zhao & Wan Reference Chen, Zhao and Wan2022). The robustness of OpenFOAM in simulating single-phase TC turbulence and two-phase TC turbulence has been demonstrated in our previous works (Xu et al. Reference Xu, Zhao, Sun, He and Wang2022, Reference Xu, Su, Lan, Zhao, He, Sun and Wang2023; Su et al. Reference Su, Yi, Zhao, Wang, Xu, Wang and Sun2024).
We consider two immiscible and incompressible fluids confined between two coaxial cylinders whose radii are 
$r_i$ (inner) and 
$r_o$ (outer). In this work, we fix the outer cylinder while allowing the inner cylinder to rotate with a constant angular velocity 
$\omega _i$. The two-phase incompressible flow is governed by the Navier–Stokes equations
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$
$$\begin{gather}{\partial _t(\rho \boldsymbol{u})} +\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{uu}) ={-}{\boldsymbol{\nabla}{p}} +\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\tau}} +\boldsymbol{f}, \end{gather}$$
where 
$\boldsymbol {u}$ is the velocity field, 
$p$ is the pressure and 
$\boldsymbol {\tau } =\mu (\boldsymbol {\nabla }\boldsymbol {u}+(\boldsymbol {\nabla }{\boldsymbol {u}})^{\rm T})$ is the viscous stress. Here 
$\rho$ and 
$\mu$ are the density and viscosity of the combined phase. The phase fraction 
$\alpha$ is introduced to characterise the variable density and viscosity, i.e. 
$\rho =\alpha \rho _d+(1-\alpha )\rho _f$ and 
$\mu =\alpha \mu _d+(1-\alpha )\mu _f$, where 
$\rho$ and 
$\mu$ with subscripts 
$f$ and 
$d$ denote the density and viscosity of continuous phase and dispersed phase. The continuum surface force method, as proposed by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992), is adopted in this study to describe the interfacial tension, i.e. 
$\boldsymbol {f}=\sigma \kappa \boldsymbol {\nabla }\alpha$, where 
$\sigma$ denotes the surface tension coefficient and 
$\kappa =-\boldsymbol {\nabla }\boldsymbol {{\cdot }}(\boldsymbol {\nabla }\alpha /|\boldsymbol {\nabla }\alpha |)$ represents the interface curvature.
In the VOF method, the phase fraction 
$\alpha$ is utilised in each cell to characterise the distribution of the two phases. The range of 
$\alpha$ is from zero to one, where 
$\alpha =0$ represents the continuous phase, 
$\alpha =1$ represents the dispersed phase, and 
$0<\alpha <1$ represents the interface region. The evolution of 
$\alpha$ is governed by the transport equation
\begin{equation} {\partial _t \alpha}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha \boldsymbol{u})=0. \end{equation}
Because of the continuity of the phase fraction, the interface between the two phases tends to become smeared. To mitigate this issue, a PLIC-based algorithm has been recently implemented to capture the interface accurately. This algorithm involves representing the interface between the two phases by employing surface cuts, which split each cell to match the phase fraction in that cell. The surface cuts are oriented according to the phase fraction gradient. The phase fraction on each cell face is then calculated from the amount submerged below the surface cut. Based on this algorithm, the resolved interface region (
$0<\alpha <1$) could be confined within a single layer of grid cells between the two phases to ensure the sharpness of the interface. It is important to note that this algorithm may encounter difficulty handling certain cells when the cut position is unclear or when multiple interfaces exist. In such cases, the interface compression approach proposed by Weller (Reference Weller2008) is applied to those cells. In this approach, an artificial compression term, which is only active in the vicinity of the interface, is added to the transport equation to prevent interface smearing based on counter-gradient transport, i.e.
\begin{equation} {\partial _t \alpha}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha \boldsymbol{u}) +\boldsymbol{\nabla}\boldsymbol{\cdot}[\alpha(1-\alpha) \boldsymbol{u}_c]=0, \end{equation}
where 
$\boldsymbol {u}_c = c\boldsymbol {u}\boldsymbol {\nabla }{\alpha }/|\boldsymbol {\nabla }{\alpha }|$ with 
$c$ being the compression factor. In addition, the multidimensional universal limiter with explicit solution (MULES) algorithm is implemented to ensure that the phase fraction 
$\alpha$ remains within the strict bounds of 0 and 1. The combination of the PLIC-based algorithm with the interface compression approach allows the present approach to be easier to implement even with an unstructured mesh, thereby increasing the robustness of the solutions. Therefore, this PLIC-based VOF method has been applied in our study to deal with the two-phase turbulence in the TC system.
To minimise computational costs without compromising the accuracy of our results, we selected a rotational symmetry of order 6 (i.e. the azimuthal angle of the simulated domain is 
${\rm \pi} /3$) and an aspect ratio of 
$\varGamma =L/d=2{\rm \pi} /3$ in the simulated TC system, where 
$d$ corresponds to the gap width between the cylinders and 
$L$ represents the axial length. This choice has been validated for both single-phase and multiphase TC turbulence (Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2013; Spandan, Verzicco & Lohse Reference Spandan, Verzicco and Lohse2018; Assen et al. Reference Assen, Ng, Will, Stevens, Lohse and Verzicco2022). The curvature of the TC system is characterised by the ratio 
$\eta =r_i/r_o=0.714$. The Taylor number, denoted as 
$Ta = \chi (r_o+r_i)^2(r_o-r_i)^2\omega _i^2/(4\nu _f^2)$, is fixed as 
$5.49\times 10^7$, where 
$\chi = [(r_i+r_o)/(2\sqrt {r_i r_o})]^4$ and 
$\nu _f=\mu _f/\rho _f$. The system Weber number, denoted as 
$We$, is given by 
${\rho _f}u_{\tau }^2 d/\sigma$, with a fixed value of 10. The 
$u_{\tau }$ is the friction velocity defined as 
$\sqrt {\tau _w/\rho _f}$, where 
$\tau _{w}$ represents the shear stress at the inner wall for single-phase flow. The frictional Reynolds number at the inner cylinder, denoted as 
$Re_\tau ={\rho _f} u_{\tau } d/\mu _f$, is fixed as 295.39. The system Reynolds number 
$Re={\rho _f} u_i d/\mu _f$ is fixed as 6000, where 
$u_i = r_i \omega _i$ is the velocity of the inner cylinder.
Physically, the length scales of fluid–fluid two-phase turbulence involved can range from the largest scale of the problem down to the Kolmogorov scale of turbulence, and even further to the molecular scale of coalescence and breakup events at the interface. Specifically, during the coalescence event, the length scales associated with film drainage can be as small as a few hundred nanometres or less (Kamp, Villwock & Kraume Reference Kamp, Villwock and Kraume2017). Ideally, it would be advantageous to conduct simulations that fully resolve all scales, similar to the study of single-phase turbulence. However, this approach is not feasible for fluid–fluid two-phase turbulence due to the significant separation between the largest flow scale and the smallest interfacial scale, which can span up to eight to nine orders of magnitude. Such a wide range of scales would require tremendous computational resources. As a result, the conventional approach is to avoid resolving the molecular scales at the interface and instead focus on resolving all turbulence scales, from the larger macroscopic scale down to the Kolmogorov length scale (Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019). In the VOF method, the coalescence and breakup of droplets is handled implicitly and two separate interfaces automatically merge when they occupy the same computational cell. This process is commonly referred to as numerical coalescence. Given that fully resolving film drainage and turbulence is prohibitively expensive from a computational point of view. An alternative approach is to use a subgrid-scale model to determine whether the droplets will coalesce. However, such approaches are highly dependent on the underlying film drainage model used and therefore their predictive capabilities are uncertain. As noted in the review of Soligo, Roccon & Soldati (Reference Soligo, Roccon and Soldati2021), there has been no fully validated method to accurately model film drainage in two-phase turbulence. In our simulations we do not use a film drainage model and allow the droplets to numerically coalesce, which is generally acceptable in two-phase flow simulations in dilute and semi-dilute regimes (Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022).
No-slip and impermeable boundary conditions are imposed in the radial direction, while periodicity is imposed in the axial and azimuthal directions. The inner and outer cylinders are subjected to a Neumann boundary condition for the phase fraction, resulting in a default contact angle of 
$90^\circ$. The maximum Courant–Friedrichs–Lewy number is set to be 0.2. First, a single-phase case is simulated to initialise the velocity field. Once a well-developed flow with a pair of Taylor rolls is obtained, the simulation is restarted, in which the spheres of diameter 
$0.2d$ containing the dispersed phase are uniformly positioned in the domain. Two different volume fractions of the dispersed phase, 
$\varphi = 5\,\%$ and 
$\varphi = 10\,\%$, are considered. The dispersed phase undergoes continuous coalescence and breakup, gradually adapting to the flow field. All the presented statistics are collected for at least 
$3 \times 10^2$ large eddy turnover times (
$(r_o-r_i)/(\omega _i r_i)$) after the two-phase flow reaches a statistically steady state.
The TC system is discretised using a collocated grid system consisting of 
$N_\theta \times N_r \times N_z = 336\times 256\times 192$ grids in the azimuthal, radial, and axial directions, respectively. The grids are distributed uniformly in the azimuthal and axial directions but are unevenly spaced and concentrated near the two cylinders in the wall-normal direction. The grid spacing is measured in units of the viscous length scale 
$\delta _\nu = \nu _f/u_\tau$ for single-phase flow. In the radial direction, the grid spacing varies from 
$0.34\delta _\nu$ near the wall to 
$2.73\delta _\nu$ at the centre of the gap. In the azimuthal direction, it ranges from 
$2.3\delta _\nu$ near the inner wall to 
$3.22\delta _\nu$ near the outer cylinder. The grid spacing remains uniform in the axial direction, with a value of 
$3.22\delta _\nu$. The Kolmogorov scale for single-phase flow, denoted as 
$\eta _k$, is determined to be 
$2.11\delta _\nu$ by employing the exact dissipation relationships given by 
$\eta _k/d=({\chi }^{-2}Ta(Nu_{\omega }-1))^{-1/4}$, where 
$Nu_\omega =T/T_{lam}$ (Eckhardt, Grossmann & Lohse Reference Eckhardt, Grossmann and Lohse2007) with 
$T$ representing the torque required to drive the cylinders and 
$T_{lam}$ corresponding to the torque when the flow is purely laminar. The maximum grid spacing is constrained to be nearly 
$1.5\eta _k$ to ensure that the grid length remains comparable to the scale of the local Kolmogorov length.
We utilise a blended scheme with the blending factor being 0.9 for the temporal term discretisation, which lies between the first-order Euler scheme and the second-order Crank–Nicolson scheme. The introduction of the dispersed phase could lead to numerical instability. The blended scheme ensures a balance between numerical stability and numerical accuracy (Greenshields Reference Greenshields2020). For spatial discretisation, we employ a second-order linear-upwind scheme to discretise the advection term in the momentum equation. The PIMPLE algorithm (Holzmann Reference Holzmann2016), which is a hybrid version of the PISO algorithm and the SIMPLE algorithm, is used to handle the pressure–velocity coupling. The PIMPLE algorithm guarantees better stability for problems that involve very large timesteps and pseudo-transient simulation. The pressure equation is solved using the geometric algebraic multigrid (GAMG) solver coupled with the simplified diagonal-based incomplete Cholesky (DIC) method. The GAMG leverages the multigrid approach, which utilises a hierarchy of grids with different resolutions to accelerate the convergence process. This allows GAMG to quickly converge to a solution, reducing computational time compared with other standard methods. In OpenFOAM, the GAMG is commonly coupled with DIC to speed up the computational efficiency in simulating two-phase flow (Scheufler & Roenby Reference Scheufler and Roenby2019; Chen et al. Reference Chen, Zhao and Wan2022). For solving velocity and phase fraction, we use an iterative solver with a symmetric Gauss–Seidel smoother. The Gauss–Seidel method is known for several advantages over other techniques. Its convergence speed and memory efficiency are particularly noteworthy. In the simulation, we maintain a tolerance of 
$10^{-6}$ for all variables to control the residuals, except for the phase fraction, which has a tolerance of 
$10^{-8}$. The computational accuracy of these settings is examined and validated in Appendices A and B, demonstrating that our methods can effectively simulate single-phase and two-phase cases.
3. Results
To examine the effect of various parameters on drag modulation, we conducted a comprehensive investigation by sequentially altering the volume fraction 
$\varphi$, viscosity ratio 
$\xi _{\mu }$ and density ratio 
$\xi _{\rho }$ to study their individual and coupled effects, as outlined in table 1. In this work, we fix the outer cylinder while sustaining the constant rotational velocity of the inner cylinder. The torque 
$T$ required to drive the inner cylinder is examined to study the drag modulation caused by different types of droplets. This set-up is commonly used to study drag modulation caused by droplets, bubbles and particles (Spandan et al. Reference Spandan, Ostilla-Mónico, Verzicco and Lohse2016; Bakhuis et al. Reference Bakhuis, Verschoof, Mathai, Huisman, Lohse and Sun2018; Yi et al. Reference Yi, Toschi and Sun2021, Reference Yi, Wang, van Vuren, Lohse, Risso, Toschi and Sun2022). Our work focuses on four typical cases, including two-phase flows with the dispersed phase being neutral droplets (
$\xi _{\rho } = 1$ and 
$\xi _{\mu } = 1$), low-viscosity droplets (
$\xi _{\rho } = 1$ and 
$\xi _{\mu } = 1/4$), light droplets (
$\xi _{\rho } = 1/4$ and 
$\xi _{\mu } = 1$) and low-viscosity light droplets (
$\xi _{\rho } = 1/4$ and 
$\xi _{\mu } = 1/4$). The torque is collected from about 
$10^5$ time steps after reaching a statistically steady state. An increase in drag is observed for two-phase flow with neutral droplets, consistent with the experimental observation (Yi et al. Reference Yi, Toschi and Sun2021). However, reducing the viscosity of the dispersed phase to 
$\xi _{\mu } = 1/4$ results in almost no change in the drag enhancement effect compared with the case with neutral droplets. On the other hand, reducing the density of the dispersed phase to 
$\xi _{\rho } = 1/4$ leads to significant drag reduction. Moreover, simultaneously reducing the density and the viscosity can result in stronger drag reduction. For low-viscosity light droplets, the kinematic viscosity is the same as that of the continuous phase. Therefore, the corresponding characteristic Reynolds number 
$Re_d = {\rho _d} u_i d/\mu _d=6000$ is the same as the system Reynolds number for the single-phase case. Considering the same Reynolds number and a viscosity ratio of 1/4, a drag reduction of 
$75\,\%$ would be obtained when the total volume fraction of the low-viscosity light droplets is 
$100\,\%$. In our work, only 
$10\,\%$ low-viscosity light droplets could lead to up to 
$30\,\%$ drag reduction, demonstrating the efficiency of the chosen volume fraction on drag modulation. In addition, we observe similar trends in drag modulation for two different volume fractions of the dispersed phase, 
$\varphi = 5\,\%$ and 
$\varphi = 10\,\%$. For given density ratio 
$\xi _{\rho }$ and viscosity ratio 
$\xi _{\mu }$, an increase in the volume fraction 
$\varphi$ leads to the pronounced amplitude of drag enhancement or reduction, which depends on the droplet properties (
$[\xi _{\rho }$, 
$\xi _{\mu }]$). We will conduct a detailed analysis of these results to elucidate the individual and coupled effects of the viscosity and density of the dispersed phase on the system's drag and turbulence properties, aiming to reveal the underlying mechanisms.
Table 1. Drag modulation in two-phase flow. Here 
$T$ represents the torque exerted on the inner cylinder and 
$T_{\varphi =0}$ denotes the torque specifically associated with the single-phase flow condition.
Given the similarity in drag modulation for the two volume fractions, 
$\varphi = 5\,\%$ and 
$\varphi = 10\,\%$, we focus our subsequent investigations on the 
$5\,\%$ volume fraction cases. In a TC system, the spatial distribution of droplets is affected by the background flow field and potentially by the inertial effect due to centrifugal force if there is a density mismatch between the phases. Figure 1(a–h) displays the instantaneous interface snapshots as well as the azimuthally and time-averaged phase fraction. The averaged radial–axial velocity vectors denote the magnitude and structure of the Taylor vortex, which provides information on the dispersed phase spatial distribution.
Figure 1. (a–d) Instantaneous interface snapshots and (e–h) corresponding azimuthally and time-averaged phase fraction 
$\langle \alpha \rangle _{\theta,t}$. The grey arrows in (e–h) denote the direction and magnitude of the averaged radial–axial velocity vectors 
$\langle u_{rz} \rangle _{\theta,t}/u_i$. Here 
$\langle \alpha \rangle _{\theta,t}$ indicates the proportion of dispersed phase (phase distribution). To facilitate the observation of the phase distribution, a logarithmic scale is employed for the legend labels. (i) The normalised azimuthal velocity 
$\langle u_{\theta } \rangle /u_i$ and ( j) the phase distribution 
$\langle \alpha \rangle$ are obtained as a function of radial position 
$r$. Inset: enlarged view of 
$\langle u_{\theta } \rangle /u_i$ near the inner cylinder. The operator 
$\langle {\cdot } \rangle$ denotes the average in the axial and azimuthal directions and over time.
3.1. Drag modulation
For neutral droplets, the effect of centrifugal force is eliminated as the densities of the two phases are perfectly matched. As illustrated in figure 1(e), the neutral droplets predominantly exhibit voids in the region where the plumes are ejected from the inner (outer) boundary layer to the bulk. This pattern resembles that of neutrally buoyant finite-size particles at a system Reynolds number of 6500 in the same definition as here, which is primarily attributed to the flow structures and the finite-size effect of the particles (Wang et al. Reference Wang, Yi, Jiang and Sun2022a). Although the neutral droplets have a weaker finite-size effect due to their capacity for deformation and breakage, the spatial distribution qualitatively aligns with that of finite-size particles. In addition, compared with the single-phase case, the neutral droplets induce a decrease in the azimuthal velocity in the vicinity of the inner cylinder, as illustrated in the inset of figure 1(i). Considering the fixed azimuthal velocity of the inner cylinder, this decrease causes an increase in the azimuthal velocity gradient at the inner cylinder, resulting in drag enhancement.
Low-viscosity droplets exhibit a distribution similar to that of the neutral droplets. This suggests that decreasing the viscosity of the dispersed phase does not introduce new mechanisms affecting the phase distribution, at least in the parameter regime studied. We note that the low-viscosity droplets are distributed mainly in the bulk region where the effect of viscosity on momentum transport is weak, which leads to a negligible effect of viscosity reduction on the flow field. Therefore, when the droplets are neutrally buoyant, the decrease in viscosity ratio does not lead to a significant change in the system's drag.
For light droplets, the dispersed phase experiences two main forces within the 
$r$–
$z$ plane. First, the centrifugal force generated by the density difference between the two phases causes the dispersed phase to migrate toward the inner cylinder. Second, the drag force, resulting from any velocity difference between the dispersed phase and the continuous phase, causes the dispersed phase to move along with the motion of the Taylor vortex. Ultimately, the dispersed phase tends to aggregate in the region close to the inner cylinder, where the drag and centrifugal forces are roughly balanced. Remarkably, a prominent interfacial structure forms in the ejection region of the plumes near the inner cylinder, and the average phase fraction 
$\langle \alpha \rangle _{\theta, t}$ in this region is approximately 1 (see figure 1c,g). Consequently, there is a notable disparity in the distribution of the phase fraction near the inner cylinder compared with other regions. In addition, an increase in the azimuthal velocity compared with the single-phase case is observed in the vicinity of the inner cylinder. Considering the fixed azimuthal velocity of the inner cylinder, this leads to a decrease in the azimuthal velocity gradient at the inner cylinder, which may be due to the density reduction lowering the effective Reynolds number. The decreased azimuthal velocity gradient at the inner cylinder results in drag reduction.
For low-viscosity light droplets, the alteration from 
$\xi _{\mu } = 1$ to 
$\xi _{\mu } = 1/4$ disrupts the balance between the drag force and the centrifugal force, leading to a modification in the phase distribution. As depicted in figure 1(h), the interfacial structure undergoes a slight stretching in the vertical direction, indicating a greater concentration of the dispersed phase near the inner cylinder, as depicted in figure 1( j). This change can be attributed to the reduced viscosity, which weakens the Taylor vortex and, in turn, makes the centrifugal force more significant. It is important to note that, unlike light droplets, low-viscosity light droplets reduce the azimuthal velocity in the vicinity of the inner cylinder compared with the single-phase case, which results in a greater azimuthal velocity gradient. Therefore, the drag reduction induced by low-viscosity light droplets is mainly attributed to their low viscosity and little related to their modulation on the azimuthal velocity.
Based on the results discussed above, we found that neutral droplets distribute mainly in the bulk region. The presence of neutral droplets causes a decrease in the azimuthal velocity near the inner cylinder. This leads to an increase in the azimuthal velocity gradient at the inner cylinder, increasing the system's drag. Low-viscosity droplets exhibit similar distribution behaviour to neutral droplets. As they are mainly distributed in the bulk region where the influence of viscosity is weak, the effect of viscosity reduction on azimuthal velocity is almost negligible. Therefore, low-viscosity droplets exhibit an indistinguishable decrease in drag compared with neutral droplets. In the case of light droplets, the centrifugal forces caused by density mismatch cause the light droplets to aggregate near the inner cylinder. The presence of light droplets leads to a decrease in the azimuthal velocity gradient at the inner cylinder probably because the light droplets lower the local effective Reynolds number, leading to drag reduction. For low-viscosity light droplets, a greater azimuthal velocity gradient is induced at the inner cylinder, which is in contrast to the drag reduction caused by them. The drag reduction is mainly attributed to the low viscosity of the droplets which causes a reduction in the viscous shear stress at the inner cylinder. In addition, it is observed in figure 1(i) that the azimuthal velocity shows a reduction in the bulk region for all two-phase cases compared with the single-phase case. This phenomenon is explained in the next section.
3.2. Momentum budget analysis
To further investigate the turbulence modulation caused by different types of droplets, a momentum budget analysis is conducted based on the conserved quantity 
$J^\omega$ that characterises the radial transport of azimuthal momentum in the two-phase TC turbulence (Su et al. Reference Su, Yi, Zhao, Wang, Xu, Wang and Sun2024)
\begin{equation} J^\omega=J_{adv}^\omega(r) + J_{dif}^\omega(r)+J_{int}^\omega(r), \end{equation}where the three terms on the right-hand side represent the advection contribution, the diffusion contribution and the interfacial contribution, respectively,
$$\begin{gather} J_{adv}^\omega(r)=\langle r^2 \rho u_r u_\theta \rangle, \end{gather}$$
$$\begin{gather}J_{dif}^\omega(r)={-}\langle \mu (r^3 \partial _r \omega + r\partial _\theta u_r)\rangle, \end{gather}$$
$$\begin{gather}J_{int}^\omega(r)={-}\int^r_{r_i} \langle {r^2 f_\theta}\rangle \,{\rm d} r. \end{gather}$$
Here 
$u_r$ is the radial velocity, 
$\omega$ is the angular velocity and 
$f_\theta$ is the azimuthal component of interfacial tension. It is well-established that the conserved quantity 
$J^\omega$ and the torque at the inner cylinder 
$T$ in the TC system are related through the equation 
$T=2{\rm \pi} L J^\omega$. The advection contribution, the diffusion contribution and the interfacial contribution are explicitly related to density, viscosity and interfacial tension, respectively, offering a convenient way to effectively decouple the effects of density, viscosity and two-phase interface on drag modulation.
For the two-phase flow with neutral droplets (
$\xi _{\rho } = 1$ and 
$\xi _{\mu } = 1$), the advection contribution in the bulk region is modulated to account for the influence of the two-phase interface. Specifically, the advection contribution increases in the half of the bulk region closer to the inner cylinder, while decreasing in the other half of the bulk region (see figure 2a). In addition, the presence of the two-phase interface introduces the interfacial contribution, which consistently exhibits a positive value along the radial position, indicating its contribution to drag enhancement (see figure 2c). Due to the lack of significant modulation in the radial-averaged advection and diffusion contributions, the increase in drag is primarily attributed to the interfacial contribution (see figure 2d).
Figure 2. Momentum budget analysis: (a) the advection contribution, (b) the diffusion contribution and (c) the interfacial contribution as a function of the radial position. (d) The three contributions are averaged in the radial direction to characterise the corresponding terms within the whole system and are shown in the stack column. The dashed line in (d) represents the averaged advection contribution for single-phase flow. All the quantities are normalised by the conserved quantity for single-phase flow. The insets in (b) show enlarged views of the diffusion contribution near the inner and outer cylinders.
In contrast to other scenarios involving drag modulation (e.g. the drag enhancement caused by finite-size particles in turbulent channel flow and turbulent plane Couette flow, where the particle-induced stress has zero contributions at the wall; Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015; Wang, Abbas & Climent Reference Wang, Abbas and Climent2017), the interfacial contribution in our cases exhibits a zero value at the inner cylinder and a non-zero value at the outer cylinder (see figure 2c). Note that the interfacial contribution is obtained by integrating the contribution from the azimuthal component of interfacial tension within each cylindrical plane from the inner cylinder to the specific position. Namely, the interfacial contribution at a specific radial position relies on the interfacial tension in the region between the cylindrical plane and the inner cylinder, not just the interfacial tension within the cylindrical plane.
It is found that the interfacial contribution shows an obvious increase with the radial position near the inner cylinder, while there is an obvious decrease with the radial position near the outer cylinder for the two-phase flow with neutral droplets. Based on the distribution characteristics of the interfacial contribution and 
$J_{int}^\omega (r)= -\int ^r_{r_i} \langle {r^2 f_\theta } \rangle \, {\rm d} r$, we propose a model to explain the effect of interfacial tension in shear turbulence. In this model, when the azimuthal velocity of the droplet phase is less than that of the surrounding continuous phase within a cylindrical plane, the droplet phase is subject to a drag force exerted by the surrounding continuous phase. The drag would be balanced by the interfacial tension 
$\langle\, f_\theta \rangle$, which takes a negative value and opposes the direction of the flow. Conversely, when the azimuthal velocity of the droplet phase is larger than that of the surrounding continuous phase, 
$\langle\, f_\theta \rangle$ takes a positive value and shares the same direction of the flow. In the vicinity of the inner cylinder, 
$\langle\, f_\theta \rangle$ takes a negative value and results in an obvious increase of interfacial contribution with the radial position. The interfacial tension opposes the direction of the flow and tends to reduce the azimuthal velocity as shown in the inset of figure 1(i), thus leading to drag enhancement. Near the outer cylinder, 
$\langle\, f_\theta \rangle$ takes a positive value and results in an obvious decrease of interfacial contribution with the radial position, i.e. the interfacial tension tends to increase the azimuthal velocity there. The interfacial contribution can be regarded as the overall effect of interfacial tension in the region between a specific cylindrical plane and the inner cylinder. The positive interfacial contribution at the outer cylinder indicates that the overall effect of interfacial tension within the gap tends to reduce the azimuthal velocity. As there is a positive interfacial contribution at the outer cylinder for all two-phase cases, azimuthal velocity shows a reduction in the bulk region for all two-phase cases compared with the single-phase case as shown in figure 1(i).
For the two-phase flow with low-viscosity droplets (
$\xi _{\rho } = 1$ and 
$\xi _{\mu } = 1/4$), the three contributions are very close to those of neutral droplets since droplets mainly exist in the bulk region where the viscosity has a very weak influence on the momentum transport. Upon closer examination, it is still possible to glean some insights into the effects of viscosity. In comparison with neutral droplets, a minor increase in the radial-averaged advection contribution is observed, as depicted in figure 2(d). This could be attributed to the reduction in viscosity within the bulk region, leading to a diminished inhibitory effect of viscosity on turbulence within this region. In addition, we observed that the radial-averaged diffusion contribution does not decrease despite the reduction in viscosity in the bulk region. This is because the decrease in droplet viscosity leads to an enhanced velocity gradient across the two-phase interface (Farsoiya et al. Reference Farsoiya, Liu, Daiss, Fox and Deike2023), which potentially contributes to the diffusion term (through 
$\partial _r \omega$ and 
$\partial _\theta u_r$). Moreover, reducing viscosity also decreases the interfacial contribution in the bulk region. We speculate that this is due to the viscosity reduction weakening the interface's ability to impede the surrounding flow field.
For the two-phase flow with light droplets (
$\xi _{\rho } = 1/4$ and 
$\xi _{\mu } = 1$), the advection contribution undergoes a significant reduction. The density ratio 
$\xi _{\rho } = 1/4$ promotes the aggregation of the dispersed phase near the inner cylinder, resulting in a substantial decrease in the upstream advection contribution due to its density-related nature. Moreover, in regions away from the inner cylinder where the dispersed phase's volume fraction is low, the significant reduction in the advection contribution can be attributed to a decrease in 
$u_r u_{\theta }$ since the density remains relatively constant. We here emphasise that the decrease in the advection contribution caused by 
$\xi _{\rho } = 1/4$ surpasses the interfacial contribution, leading to drag reduction.
For the two-phase flow with low-viscosity light droplets (
$\xi _{\rho } = 1/4$ and 
$\xi _{\mu } = 1/4$), the diffusion contribution is diminished near the inner cylinder due to its viscosity-related nature. As discussed earlier regarding the azimuthal velocity, we have shown that the drag reduction is mainly due to the low viscosity of the droplets. The advection contribution will be implicitly reduced through the velocity field (i.e. 
$u_ru_\theta$) to maintain the conserved quantity constant along the radial position.
These results demonstrate that decreased viscosity alone does not significantly affect momentum transport, whereas the decreased viscosity in combination with decreased density can reduce momentum transport significantly and, hence, lead to drag reduction.
3.3. Turbulent fluctuation analysis
The momentum budget analysis mainly focuses on the effects of interfacial tension and the fluid properties of the dispersed phase on global transport. However, while it sheds light on these aspects, the details of how the dispersed phase influences turbulence remain elusive. Therefore, it becomes crucial to delve deeper into statistics of turbulence properties. Due to the density difference between the dispersed and continuous phases, the Reynolds average is no longer the best choice for obtaining the Reynolds stress. Consequently, we adopt the Favre-average method (Favre Reference Favre1969), which applies a density-weighted average to the velocity. The Favre-averaged Navier–Stokes equations are formally similar to the Reynolds-averaged equations for the single-phase case and are easily compared across cases. This advantage makes the Favre average widely used in the study of compressible and multiphase flows. The Favre average applies a density-weighted average to the velocity field, resulting in fluctuations in the velocity vector,
\begin{equation} \boldsymbol{u}^{\prime\prime}=\boldsymbol{u}-\tilde{\boldsymbol{u}},\end{equation}
where 
$\tilde {\boldsymbol {u}} = \langle {\rho \boldsymbol {u}} \rangle / \bar {\rho }$ is the Favre-averaged velocity vector and 
$\bar {\rho }= \langle {\rho } \rangle$. On the other hand, the fluctuation in the velocity vector using the Reynolds average is given as
\begin{equation} \boldsymbol{u}^{\prime}=\boldsymbol{u}-\langle {\boldsymbol{u}}\rangle. \end{equation}The velocity fluctuation due to the Favre average and the Reynolds average are related with
\begin{equation} \boldsymbol{u}^{\prime\prime}=\boldsymbol{u}^{\prime}-\boldsymbol{a}, \end{equation}
where 
$\boldsymbol {a} = \langle {\rho }^{\prime }\boldsymbol {u}^{\prime } \rangle /\langle {\rho } \rangle$ is a function of radial position, characterising the overall effect of density fluctuations within each cylindrical plane. When the dispersed and continuous phases possess identical densities, the Favre average simplifies to the Reynolds average.
To explore the effect of the dispersed phase on the background turbulence, we show the Favre-average Reynolds stress 
${\langle {\rho u_r^{\prime \prime } u_\theta ^{\prime \prime }}\rangle }$ as well as the radial component 
${\langle {\rho u_r^{\prime \prime } u_r^{\prime \prime }} \rangle }$, the azimuthal component 
${\langle {\rho u_\theta ^{\prime \prime } u_\theta ^{\prime \prime }}\rangle }$ and the axial component 
${\langle {\rho u_z^{\prime \prime } u_z^{\prime \prime }}\rangle }$ of the TKE 
${k=0.5 \langle {\rho u_r^{\prime \prime } u_r^{\prime \prime }} + {\rho u_\theta ^{\prime \prime } u_\theta ^{\prime \prime }}+{\rho u_z^{\prime \prime } u_z^{\prime \prime }}\rangle }$ as depicted in figure 3. The Favre-average Reynolds stress (hereafter referred to as Reynolds stress) 
${\langle {\rho u_r^{\prime \prime } u_\theta ^{\prime \prime }}\rangle }$ shows a similar profile to the advection contribution 
$\langle r^2 \rho u_r u_\theta \rangle$ in the momentum transport due to
\begin{equation} \langle r^2 \rho u_r u_\theta\rangle = r^2 \bar{\rho} \widetilde{u_r} \widetilde{u_\theta} +{r^2 \langle {\rho u_r^{\prime\prime} u_\theta^{\prime\prime}}\rangle}, \end{equation}
where 
$r^2 \bar {\rho } \widetilde {u_r} \widetilde {u_\theta }$ is the contribution from mean flow and 
${r^2 \langle {\rho u_r^{\prime \prime } u_\theta ^{\prime \prime }}\rangle }$ is the contribution from turbulent flow. The contribution from the mean flow is nearly negligible, with its maximum value being less than 
$1\,\%$ of the turbulence contribution (see the inset in figure 3a). This suggests that the mean flow has minimal influence on the advection term in momentum transport. Since 
$r$ serves solely as a positional parameter, the advection term in momentum transport is determined completely by the Reynolds stress. Therefore, we can reinterpret momentum transport from the perspective of turbulence.
Figure 3. The modulation on the turbulence properties. (a) The Favre-averaged Reynolds stress 
${\langle {\rho u_r^{\prime \prime } u_\theta ^{\prime \prime }} \rangle }/{\tau }_{i,\varphi =0}$, (b) the radial (wall-normal) component of the TKE 
${ \langle {\rho u_r^{\prime \prime } u_r^{\prime \prime }}\rangle }/{\tau }_{i,\varphi =0}$, (c) the azimuthal (streamwise) component of the TKE 
${ \langle {\rho u_\theta ^{\prime \prime } u_\theta ^{\prime \prime }} \rangle }/{\tau }_{i,\varphi =0}$ and (d) the axial (spanwise) component of the TKE 
${\langle {\rho u_z^{\prime \prime } u_z^{\prime \prime }} \rangle }/{\tau }_{i,\varphi =0}$ are shown as a function of the radial position, where 
${\tau }_{i,\varphi =0}$ is the total shear stress on the inner cylinder in single-phase flow. The inset in (a) shows the Favre-averaged term of the mean flow 
$\bar {\rho } \widetilde {u_r} \widetilde {u_\theta }/{\tau }_{i,\varphi =0}$. The maximum value of the ordinate of the inset is set as 0.01, while the maximum value of the ordinate of (a) is set as 1.
Neutral droplets lead to an increase in Reynolds stress in the majority of the bulk region closer to the inner cylinder, whereas it decreases in the other half of the bulk region, thus accommodating the influence of the two-phase interface. However, the total Reynolds stress changes very little, and the drag enhancement is mainly due to the interfacial contribution. For low-viscosity droplets, their lower viscosity weakens the inhibitory effect of the bulk region on turbulence compared with the continuous phase. However, the overall modulation of Reynolds stress is minimal in comparison with neutral droplets. Thus, the interfacial contribution remains the key factor in determining the drag enhancement. In the case of light droplets, the aggregation of the light droplets near the inner cylinder results in a notable decrease in the upstream Reynolds stress (through 
$\rho$), primarily due to their density-related characteristics. Furthermore, the downstream Reynolds stress is correspondingly reduced (through 
$u_r^{\prime \prime } u_\theta ^{\prime \prime }$) to ensure momentum conservation along the radial position, ultimately leading to drag reduction. For low-viscosity light droplets, the alteration from 
$\xi _{\mu } = 1$ to 
$\xi _{\mu } = 1/4$ further reduces the momentum transport near the wall by decreasing the diffusion term (which is explicitly related to viscosity), thereby once again implicitly reducing the Reynolds stress in the bulk region (through 
$u_r^{\prime \prime } u_\theta ^{\prime \prime }$) and causing a stronger drag reduction. In the studied parameter regime, neutral droplets and low-viscosity droplets increase the system's drag due to the interfacial contribution, while light droplets and low-viscosity light droplets induce drag reduction by decreasing the Reynolds stress.
For the three components of TKE, we observe an overall decrease due to the presence of neutral droplets, indicating that the presence of the two-phase interface weakens the total turbulence intensity. Considering that Reynolds stress did not decrease significantly, it is likely that the presence of the two-phase interface leads to a transformation of turbulence or changes in the energy dissipation process within the system. Recent studies (Perlekar et al. Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022) have demonstrated that the two-phase interface results in a decrease in energy for large-scale vortices and an increase in energy for small-scale vortices, suggesting an alteration in the energy cascade process. This alteration in the energy cascade process is likely responsible for the reduction observed in the three components of TKE, as the significant decrease in the total energy of large-scale vortices outweighs the increase in the total energy of small-scale vortices within the system. Furthermore, the slight increase in the three components of TKE for low-viscosity droplets compared with neutral droplets provides additional support for this conjecture. The slight increase indicates that reducing the viscosity of droplets weakens the reduction scale of TKE (or the alteration in the energy cascade process), aligning with the findings from Crialesi-Esposito et al. (Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022). From an alternative standpoint, the slight increase in the three components of TKE compared with neutral droplets can be attributed to the weakened inhibitory effect on turbulence caused by the decreased viscosity in the bulk region.
Light droplets cause a significant reduction in the three components of the TKE due to their explicit density dependence. Near the inner cylinder, this reduction primarily arises from the decreased density, whereas in the remaining region, it stems from the reduction in velocity fluctuations. In addition, it is observed that low-viscosity light droplets result in a further reduction in the TKE components. This can be attributed to the fact that the low-viscosity droplets weaken momentum transport near the inner cylinder, subsequently causing a decrease in the energy available for transfer from the boundary layer to the bulk region.
Based on the analysis of turbulent fluctuations, we establish a connection between turbulence and the system's drag through Reynolds stress. This approach provides more detailed insights into the effect of the two-phase interface, dispersed phase viscosity and dispersed phase density on turbulence and the system's drag. In the case of two-phase flow with neutral droplets, the drag enhancement is primarily governed by the interfacial contribution. To maintain momentum transport conservation (a constant conserved quantity), the Reynolds stress is adjusted to match the presence of the two-phase interface, but the overall change is not significant. Conversely, low-viscosity droplets will cause a slight increase in Reynolds stress within the bulk region since their low viscosity weakens the suppression of turbulence. However, the overall effect on the drag is not significant compared with the neutral droplets, and the interfacial contribution still dominates the drag enhancement. On the other hand, in the case of light droplets, their low density causes them to aggregate near the inner cylinder, resulting in a significant reduction of Reynolds stress in this region due to the explicit density dependence of Reynolds stress. The Reynolds stress in other regions decreases accordingly to ensure momentum transport conservation, ultimately leading to drag reduction. When the low viscosity is introduced into light droplets, a more pronounced reduction in momentum transport occurs as the droplets predominantly aggregate near the inner cylinder, where viscosity plays a dominant role in momentum transport. Consequently, this leads to greater drag reduction.
To summarise, neutral droplets and low-viscosity droplets primarily contribute to drag enhancement through the interfacial contribution, while light droplets reduce the system's drag by reducing Reynolds stress. Furthermore, low-viscosity light droplets contribute to drag reduction mainly by reducing the diffusion contribution and Reynolds stress. It is important to note that viscosity plays different roles in different regions. Low viscosity near the wall results in a reduction in downstream Reynolds stress by decreasing upstream momentum transport, whereas low viscosity in the bulk region may lead to an increase in Reynolds stress.
3.4. TKE budget analysis
Considering the potential alterations in turbulence characteristics resulting from the introduction of a two-phase interface or changes in droplet fluid properties, it is not advisable to directly correlate the TKE with system's drag. System's drag imparts kinetic energy to the system, but a substantial fraction of this energy dissipates through turbulent mechanisms. Hence, examining the dissipation rate of TKE provides a valuable means to establish a significant linkage between system's drag and turbulence properties. The transport equation of the TKE can be written in the form (Besnard et al. Reference Besnard, Harlow, Rauenzahn and Zemach1992; Wong et al. Reference Wong, Baltzer, Livescu and Lele2022)
\begin{align} 0 &= \underbrace{\langle{\rho\boldsymbol{u}^{\prime\prime} \boldsymbol{u}^{\prime\prime}}\rangle:(\boldsymbol{\nabla}\tilde{\boldsymbol{u}})}_P \underbrace{-\boldsymbol{\nabla}\boldsymbol{\cdot}\langle{p^{\prime}\boldsymbol{u}^{\prime}}\rangle}_{PD} \underbrace{-\boldsymbol{\nabla}\boldsymbol{\cdot}\langle{\rho\boldsymbol{u}^{\prime\prime} \boldsymbol{u}^{\prime\prime}\boldsymbol{\cdot}\boldsymbol{u}^{\prime\prime}}\rangle/2}_T +\underbrace{\boldsymbol{\nabla}\boldsymbol{\cdot}\langle \boldsymbol{\tau}^{\prime} \boldsymbol{\cdot}\boldsymbol{u}^{\prime}\rangle}_D \underbrace{-\langle\boldsymbol{\tau}^{\prime}:\boldsymbol{\nabla} \boldsymbol{u}^{\prime}\rangle}_\epsilon \nonumber\\ &\quad +\underbrace{\langle\boldsymbol{u}^{\prime} \boldsymbol{\cdot}\boldsymbol{f}^{\prime}\rangle}_I +\underbrace{\boldsymbol{a}\boldsymbol{\cdot}(\boldsymbol{\nabla}{\langle{p}\rangle -\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\boldsymbol{\tau}\rangle})}_M, \end{align}
where 
$p^{\prime }=p-\langle {p}\rangle$, 
$\boldsymbol {\tau }^{\prime }=\boldsymbol {\tau }-\langle {\boldsymbol {\tau }}\rangle$ and 
$\boldsymbol {f}^{\prime }=\boldsymbol {f}-\langle {\boldsymbol {f}}\rangle$. Here 
$P$ is the production term, 
$PD$ the pressure diffusion, 
$T$ the turbulence transport, 
$D$ the viscous diffusion, 
$\epsilon$ the dissipation term, 
$I$ interface term and 
$M$ the additional production term due to the density difference between the dispersed and continuous phases. Figure 4 shows the distribution of each term in the TKE transport equation as a function of radial position. The energy transfer rate between the mean kinetic energy and the TKE, 
$P$, is consistently positive in our cases. It is mainly determined by Reynolds stress and azimuthal velocity gradient, i.e. 
$\langle \rho {u_r}^{\prime \prime }u_\theta ^{\prime \prime }\rangle \partial _r \widetilde {u_\theta }$. We use 
$\epsilon$ to represent the dissipation rate of TKE and is always a negative value. This value is directly related to the torque required to sustain the constant rotational velocity of the inner cylinder. Note that the positive and negative values represent the gain and the loss in TKE, respectively. The pressure diffusion 
$PD$, turbulence transport 
$T$ and viscous diffusion 
$D$ often indicate the interchange or redistribution of the TKE. The balance term reflects any imbalance due to numerical effects and it is negligibly small.
Figure 4. TKE budget. Each term on the right-hand side of (3.9) is shown for (a) single-phase flow, (b) two-phase flow with neutral droplets and (c) two-phase flow with low-viscosity light droplets. (d) The quantities are averaged in the radial direction to characterise the corresponding terms within the whole system. All the quantities are normalised by the 
$|\epsilon _{\varphi =0,i}|$, where 
$\epsilon _{\varphi =0,i}$ denotes the dissipation rate at the inner cylinder for single-phase flow and the operator 
$|{\cdot }|$ gets the absolute value of the quantity.
For the two-phase flow with droplets, the presence of the two-phase interface introduces an interface term 
$I$ (see figure 4b). By decreasing the dispersed phase density, an additional production term 
$M$ is observed mainly near the inner cylinder (see figure 4c). By comparing with the single-phase case, we find that the pressure diffusion 
$PD$, turbulence transport 
$T$ and viscous diffusion 
$D$ are all reduced in two-phase cases, indicating that the presence of dispersed droplets weakens the interchange or redistribution of the TKE.
To assess the overall effect of each term in (3.9) on the TKE transport, an analysis was conducted by globally averaging these terms, as illustrated in figure 4(d). The findings indicate that TKE is mainly input into the system through the production term (
$P$) and dissipated through the dissipation term (
$\epsilon$). On the other hand, the effects of the pressure diffusion term (
$PD$), turbulence transport term (
$T$) and viscous diffusion term (
$D$) are nearly negligible. In addition, the interface term (
$I$) and the additional production term (
$M$), which are specific to two-phase flow, exhibit a small positive value, indicating that the two-phase interface and the density difference between the phases enhance the TKE transport. Consequently, the focus is on examining the modulation of the production term (
$P$), dissipation term (
$\epsilon$), interface term (
$I$) and additional production term (
$M$) with respect to their counterparts in single-phase flow, as depicted in figure 5.
Figure 5. (a) The production difference, (b) the dissipation difference, (c) the interface term and (d) the additional production term due to density difference between phases are shown as a function of radial position. All the quantities are normalised by the radial-averaged 
$|\epsilon _{\varphi =0}|$ to illustrate their magnitudes. Here 
$P_{\varphi =0}$ and 
$\epsilon _{\varphi =0}$ denote the production term and dissipation term for the single-phase flow, respectively.
As shown in figure 4, the dissipation term is negative and indicates a loss of TKE. To facilitate observation and analysis, the absolute value of the dissipation term 
$|\epsilon |$ is adopted to characterise the magnitude of dissipation. The production difference (
$P-P_{\varphi =0}$), dissipation difference (
$|\epsilon |-|\epsilon _{\varphi =0}|$), interface term (
$I$) and additional production term (
$M$) are used to examine the effect of droplets on TKE transport, as depicted in figure 5. A positive value indicates an enhancement effect on TKE transport, while a negative value suggests a weakening effect. The interface term (
$I$), specific to two-phase flow, is essentially positive, signifying its role in enhancing TKE transport. Moreover, for light droplets and low-viscosity light droplets, the additional production term (
$M$) is predominantly positive near the inner cylinder, whereas it is close to zero in other regions. This term primarily arises from the interfacial structure adsorbed on the inner cylinder and contributes to the enhanced transport of TKE across the interface.
For the two-phase flow with neutral droplets (
$\xi _{\rho } = 1$ and 
$\xi _{\mu } = 1$), we observe net positive value in turbulence production near the inner cylinder, whereas net negative values are observed in the rest of the domain. This indicates that the production term is enhanced near the inner cylinder, while the production term is weakened outside this region. The dissipation difference shows a negative value in the region closest to the inner and outer cylinders and a positive value outside the closest region. Correspondingly, the interface term primarily leads to an increase in TKE transport near the inner and outer cylinders. The locations where the dissipation difference is positive are approximately close to the regions where the interface term is positive, indicating that the interfacial tension enhances the dissipation near the cylinders. However, as depicted in figure 6(a), when considering the global perspective, although the two-phase interface contributes to the TKE transport through the interfacial tension, it also results in a decrease in the production term. This ultimately leads to a slight reduction in the dissipation term, indicating a decrease in the amount of power required to sustain turbulence. Therefore, the drag enhancement induced by the neutral droplets weakly depends on the TKE transport. Furthermore, the moderating effect of low-viscosity droplets on TKE transport is similar to that of neutral droplets. This similarity suggests that the drag enhancement caused by low-viscosity droplets is also weakly dependent on the TKE transport.
Figure 6. (a) TKE budget analysis. The quantities are averaged in the radial direction to characterise the corresponding terms within the whole system. ‘Other’ denotes 
$(PD+T+D)-[PD+T+D]_{\varphi =0}$, where 
$[PD+T+D]_{\varphi =0}$ denotes the sum of 
$PD$, 
$T$ and 
$D$ for single-phase flow. (b) The dissipation difference 
$|\epsilon |-|\epsilon _{\varphi =0}|$ is integrated within each cylindrical plane from the inner cylinder to a specific radial position. All the quantities are normalised by the radial-averaged 
$|\epsilon _{\varphi =0}|$ to illustrate their magnitudes.
In the case of two-phase flow with light droplets (
$\xi _{\rho } = 1/4$ and 
$\xi _{\mu } = 1$), the presence of these droplets weakens turbulence production throughout the entire domain. This weakening is attributed to the reduced Reynolds stress, which leads to a significant reduction in the total turbulence production, as depicted in figures 5(a) and 6(a). From a global perspective, there is a negligible modulation effect on the sum of pressure diffusion, turbulence transport and viscous diffusion. By isolating the effects of the interface term and the additional production term resulting from density differences, it is observed that a total of 
$10.4\,\%$ of the turbulence dissipation is reduced. This reduction indicates that the power required to sustain turbulence in the system has decreased by this amount, ultimately resulting in a drag reduction. It is worth noting that the reduction in the total dissipation rate is primarily due to the decrease in the dissipation rate near the outer cylinder, as shown in figures 5(b) and 6(b). This is because the density contrast between the phases boosts TKE transport near the inner cylinder, thus weakening the reduction in the total dissipation rate there.
In the case of two-phase flow with low-viscosity light droplets (
$\xi _{\rho } = 1/4$ and 
$\xi _{\mu } = 1/4$), a more pronounced reduction in the dissipation rate is observed. Specifically, a total of 
$17.7\,\%$ of the power required to sustain turbulence in the system has been decreased. The transition from light droplets to low-viscosity light droplets leads to a further decrease in the Reynolds stress, resulting in a more significant weakening of turbulence production. Moreover, the reduction in droplet viscosity diminishes the effect of the additional production term while having minimal effects on the other terms. As a result, a larger reduction in the turbulent dissipation rate is achieved, leading to a more substantial drag reduction.
Based on the aforementioned analyses, it is evident that two-phase flow introduces an interface term and an additional production term arising from the density difference between the two phases. While the interfacial tension essentially enhances the TKE transport, the drag enhancement is not strongly correlated with TKE transport. This mechanism holds for both neutral droplets and low-viscosity droplets, as observed in our examination. Light droplets, on the other hand, cause a reduction in the production term mainly by reducing the Reynolds stress. However, both the density difference between the two phases and the interface between the two phases lead to an increase in TKE transport. Consequently, the reduction in the turbulent dissipation rate is mainly attributed to the decrease in the turbulence production term, which in turn causes drag reduction. The production term is further reduced in the case of low-viscosity light droplets due to their greater reduction in Reynolds stress. In addition, the reduction in viscosity weakens the additional production term. As a result, a greater drag reduction is achieved.
4. Conclusions
In this study, we have examined how the dispersed phase influences system's drag and turbulence properties in a two-phase fluid–fluid TC system operating at a system Reynolds number of 
$6\times 10^3$ and a system Weber number of 10. To achieve this, we employ the interface-resolved VOF method, which enables the resolution of two-phase flows through a single-equation formulation. This approach allows us to conduct operations similar to single-phase flow and facilitates the exploration of how the dispersed phase influences the statistical properties of turbulence.
Through our analysis of momentum transport and turbulent fluctuations, we establish a connection between turbulence statistics and the system's drag using Reynolds stress. This provides detailed insights into the effect of the two-phase interface, dispersed phase viscosity and dispersed phase density on turbulence and the system's drag. In the case of two-phase flow with neutral droplets, drag enhancement is primarily governed by the interfacial contribution. The Reynolds stress is modulated to account for the influence of the two-phase interface, but the overall change is not significant. Lowering the droplets’ viscosity will cause a slight increase in Reynolds stress, as their low viscosity weakens the suppression of turbulence within the bulk region. However, the overall effect on drag is not significant compared with neutral droplets, and the interfacial contribution still dominates drag enhancement. On the other hand, in the case of light droplets, their low density causes them to aggregate near the inner cylinder, resulting in a significant reduction in Reynolds stress in this region due to the explicit density dependence of Reynolds stress. The Reynolds stress in other regions decreases accordingly to ensure momentum transport conservation, ultimately leading to drag reduction. When low viscosity is introduced into light droplets, a more pronounced reduction in Reynolds stress occurs, as the droplets predominantly aggregate near the inner cylinder, where viscosity plays a dominant role in momentum transport. Consequently, this leads to greater drag reduction. In summary, neutral droplets and low-viscosity droplets primarily contribute to drag enhancement through the interfacial contribution, while light droplets reduce the system's drag by explicitly reducing Reynolds stress. Furthermore, low-viscosity light droplets contribute to drag reduction by explicitly and implicitly reducing Reynolds stress.
Furthermore, we have investigated the effect of the dispersed phase on the transport of TKE. In two-phase flow, an interface term arises due to interfacial tension, potentially accompanied by an additional production term resulting from the density difference between the two phases. While interfacial tension primarily enhances TKE transport, the correlation between drag enhancement and TKE transport is not strong, as the two-phase interface also weakens turbulence production. This holds true for both neutral droplets and low-viscosity droplets. On the other hand, light droplets reduce the production term primarily by lowering Reynolds stress. However, the density contrast between the two phases enhances TKE transport near the inner wall. Consequently, the reduction in the turbulence dissipation mainly stems from the decrease in the production term, leading to drag reduction. In the case of low-viscosity light droplets, the production term is further reduced due to a greater decrease in Reynolds stress. In addition, the decrease in viscosity weakens the additional production term. This results in a more significant reduction in the turbulence dissipation, leading to stronger drag reduction.
From the perspective of turbulence, we find that the Reynolds stress plays a key role in drag modulation. It participates in momentum transport in the form of an advection contribution and participates in TKE transport as an important component of the turbulence production term. From a drag reduction perspective, we emphasise the critical role that density ratio plays in this process. By reducing the density of the dispersed phase, the aggregation of the dispersed phase near the inner cylinder is promoted, which amplifies the influence of lower density on Reynolds stress, thereby causing significant drag reduction at a very small volume fraction of the dispersed phase. Furthermore, it enables greater drag reduction by promoting the transportation of low-viscosity droplets towards the near-wall region where viscosity plays a dominant role in momentum transport.
Funding
This work is financially supported by the National Natural Science Foundation of China under grant numbers 11988102, 22478421, 12402299 and 12402298, the New Cornerstone Science Foundation through the New Cornerstone Investigator Program and the XPLORER PRIZE and the Science Foundation of China University of Petroleum, Beijing (grant number 2462024YJRC008).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Conservation of momentum transport
In TC flow, 
$J^\omega$ should be constant in the radial direction, but numerically it does show some dependence on the radial direction. Because of numerical errors, 
$J^\omega$ will deviate slightly from being constant. To quantify this difference, Zhu et al. (Reference Zhu, Ostilla-Mónico, Verzicco and Lohse2016) defined
\begin{equation} \varDelta_J = \frac{\max(J^{\omega}(r))-\min(J^{\omega}(r))}{\langle J^{\omega}(r) \rangle_r}, \end{equation}
where the operator 
$\langle \cdot \rangle _r$ denotes average in the radial direction. As illustrated by Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013) and Zhu et al. (Reference Zhu, Ostilla-Mónico, Verzicco and Lohse2016), 
$\varDelta _J \leq 1\,\%$ is the very strict requirement for the single-phase case. We make sure that the single-phase case meets this criterion. For the two-phase cases, due to the additional errors introduced by the two-phase interface and the disparities in liquid properties between the two phases, 
$\varDelta _J$ will be slightly larger than 
$1\,\%$ for the two-phase flow as listed in table 2.
Table 2. Radial deviation of the 
$J^{\omega }$.
Appendix B. Data validation and resolution test
To validate the accuracy of our simulations, we have simulated two cases with the Taylor number being 
$3.90 \times 10^6$ (
$Re=1600$) and 
$9.52 \times 10^6$ (
$Re=2500$) and validated our results through comparisons with those from Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013) in our previous work (Su et al. Reference Su, Yi, Zhao, Wang, Xu, Wang and Sun2024). Because the 
$Re=6000$ studied in this work is not included in the work of Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013), we additionally simulate a single-phase case at 
$Re = 5600$ and compare the dimensionless conserved quantity 
$Nu_\omega =J^{\omega }/J^{\omega }_{lam}$ with that from Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013), where 
$J^{\omega }_{lam}$ corresponding to the case when the flow is fully laminar. In our work, the minimum flow geometry with a rotational symmetry of six (
$n_{sym} = 6$, i.e. the azimuthal angle of the simulated domain is 
${\rm \pi} /3$) and an aspect ratio of 
$\varGamma =L/d=2{\rm \pi} /3$ is selected to reduce the computational cost while not affecting the results, which has been verified by previous studies (Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2013; Ostilla-Mónico, Verzicco & Lohse Reference Ostilla-Mónico, Verzicco and Lohse2015). We represent the 
$Nu_{\omega }$ from Ostilla et al. (Reference Ostilla, Stevens, Grossmann, Verzicco and Lohse2013) and OpenFOAM as 
$Nu_{\omega,Ot}$ and 
$Nu_{\omega,Op}$, respectively, as presented in table 3. The deviation is obtained by 
$(Nu_{\omega,Op}/Nu_{\omega,Ot})-1$. A deviation of 
$-1.18\,\%$ is found, indicating that our simulation is sufficient to capture the flow field information.
Table 3. Validation of the calculation results for the single-phase flow at 
$Re = 5600$.
To obtain reliable numerical results, the grid's spatial resolutions have to be sufficient. The resolution test is conducted for the single-phase flow and the two-phase flow with 
$5\,\%$ low-viscosity light droplets. A roughly resolved case (
$N_\theta \times N_r \times N_z = 224\times 160\times 144$) with the maximum grid spacing being about 2.02
$\eta _k$, a reasonably resolved case (
$N_\theta \times N_r \times N_z = 336\times 256\times 192$) with the maximum grid spacing being about 1.5
$\eta _k$ and a well-resolved case (
$N_\theta \times N_r \times N_z = 448\times 320\times 288$) with the maximum grid spacing being about 
$1.01\eta _k$ are considered for resolution test as depicted in figure 7. For the grid with 
$448\times 320\times 288$, most of the grid spacing is less than 
$\eta _k$. We have listed the three types of grid resolution in table 4 and compared them with the grid used by Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014) to investigate the bubble dispersion in turbulent TC flow at 
$Re = 5000$. All three grid resolutions have a finer grid spacing in the radial and azimuthal directions than those used by Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014). For the grid spacing in the axial direction, the two finer grid resolutions used in this work have a similar scale to the grid used by Chouippe et al. (Reference Chouippe, Climent, Legendre and Gabillet2014). It is observed that the grid with 
$224\times 160\times 144$ can roughly resolve the single-phase case (see figure 7). For the two-phase case with 
$5\,\%$ low-viscosity light droplets, the grid with 
$224\times 160\times 144$ shows a distinguishable difference of 
$Nu_\omega$ compared with the two finer grid resolutions. For the two finer grids, both the 
$Nu_\omega$ lie within 
$1\,\%$ error bar for the single-phase cases. Because the two-phase interface as well as the disparities in liquid properties between the two phases could introduce additional numerical errors, we provide an error bar of 
$2\,\%$ as a reference. Here, 
$\varDelta _J =1.26\,\%$ for 
$336\times 256\times 192$ and 
$\varDelta _J =1.35\,\%$ for 
$448\times 320\times 288$. In addition, there is a deviation of 
$0.67\,\%$ for 
$\langle J^{\omega }(r) \rangle _r$ when comparing the results from the two finer grids. Therefore, the adopted spatial resolution 
$N_\theta \times N_r \times N_z = 336\times 256\times 192$ is sufficient to obtain reliable results for both single-phase and two-phase cases studied.
Figure 7. Radial dependence of 
$Nu_\omega$ for three different grid resolutions. An error bar indicating a 
$1\,\%$ error is provided for the single-phase case and an error bar indicating a 
$2\,\%$ error is provided for the two-phase flow with low-viscosity light droplets for reference.
Table 4. Grid resolutions. The grid spacing is shown in units of the viscous length scale 
$\delta _{v}$ for single-phase case.
The resolution sensitivity of the statistics discussed in this work is also conducted for the two-phase case with 
$5\,\%$ low-viscosity light droplets. Considering that the distribution of the phase fraction may be dependent on the grid resolution, the distribution of phase fraction for three different grid resolutions is shown in figure 8(a). In addition, given that turbulent statistics are higher-order and more challenging to get converged results compared with the lower-order results such as mean velocity profiles, we additionally show Favre-averaged Reynolds stress, three components of the TKE as well as production term and dissipation term of the TKE transport for three different grid resolutions as depicted in figure 8(b–d). All results are nearly identical for the two finer grid resolutions, indicating that the statistics discussed in this work are qualitatively converged. It is observed that the production term and dissipation term of the TKE transport are visually nearly identical for all three grid resolutions, which is due to the large difference between the maximum and minimum values in the plot.
Figure 8. Radial dependence of (a) phase fraction distribution, (b) Favre-averaged Reynolds stress, (c) three components of the TKE and (d) production and dissipation terms of the TKE transport for three different grid resolutions for the 
$5\,\%$ low-viscosity light droplets case. The terms in (d) are normalised by the 
$|\epsilon _{\varphi =0,i}|$.
To show the effect of the two finer grid resolutions in more detail, we plot the relative differences as a function of the radial position (see figure 9). The relative differences are calculated by 
$Q_1/Q_2-1$, where 
$Q_1$ represents the quantity for the grid with 
$336\times 256\times 192$ and 
$Q_2$ represents the quantity for the grid with 
$448\times 320\times 288$. It is observed that near the two cylinders, the magnitude of relative difference could be larger than 
$10\,\%$ when the quantity 
$Q_2$ is close to zero. This is acceptable considering that the absolute difference would be very small. In the bulk region, the magnitude of the relative difference for the Reynolds stress is within 
$2\,\%$, which has a similar scale as that of the conserved quantity since the conserved quantity in the bulk region is primarily determined by the Reynolds stress. For the three components of TKE, the magnitude of relative differences is within 
$5\,\%$ in the bulk region. Because the production and dissipation terms are high-order quantities, the magnitude of relative difference could be up to nearly 
$10\,\%$ in the bulk region. Since the magnitude of 
$P$ and 
$\epsilon$ in the bulk region are much smaller than those near the inner cylinder, the relative difference is difficult to distinguish visually in figure 8(d). In two-phase flow simulations, the relative differences in figure 9 are acceptable for judging that the grid with 
$336\times 256\times 192$ is sufficient to provide reliable turbulent statistics.
Figure 9. (a) Relative difference of Reynolds stress and three components of TKE for the two finer grid resolutions. (b) Relative difference of the production and dissipation terms of the TKE transport for the two finer grid resolutions.
Appendix C. Mean streamwise velocity profiles
In figure 10, we plot the mean azimuthal velocity profiles normalised by the frictional velocity, 
$u^+=(u_i-\langle u_{\theta } \rangle )/u_\tau$, vs the wall distance 
$y^+=(r-r_i)/\delta _{v}$. Given the differences in fluid properties in the cases with light droplets and low-viscosity light droplets and uneven distribution of the droplets, it is difficult to obtain proper density and viscosity to calculate friction velocity 
$u_\tau$ and viscous length scale 
$\delta _\nu$. Different choices of equivalent density and equivalent viscosity will cause differences in results. For ease of analysis, we have taken the density and viscosity of the continuous phase to calculate the friction velocity and viscous length scale, i.e. 
$u_\tau =\sqrt {\tau _{i}/\rho _f}$ and 
$\delta _\nu = \nu _f/u_\tau$. The shear stress at the inner cylinder is obtained by 
$\tau _{i} = \langle \mu r_i \partial _r \omega \rangle$ for different cases, where 
$\mu$ is the dynamic viscosity of the combined phase. For the single-phase case, 
$u^+$ follows the linear relation 
$u^+=y^+$ well in the viscous sublayer (
$y^+<5$). At 
$y^+>30$, 
$u^+$ does not exhibit a clear logarithmic shape due to the small 
$Re$ in this study (Huisman et al. Reference Huisman, Scharnowski, Cierpka, Kähler, Lohse and Sun2013). Moreover, the velocity profile shows an obvious discrepancy from the classical logarithmic law profile with the typical values 
$\kappa =0.4$ and 
$B=5.2$, which is because the boundary layer is not yet fully developed and the zero-pressure-gradient boundary layer is not satisfied in the TC flow (Ostilla-Mónico et al. Reference Ostilla-Mónico, Van Der Poel, Verzicco, Grossmann and Lohse2014). For the cases with neutral droplets and low-viscosity droplets, 
$u^+$ shifts slightly downwards in the buffer layer and above, which has also been observed in the case of drag enhancement caused by finite-size particles in plane-Couette flow (Wang, Jiang & Sun Reference Wang, Jiang and Sun2023). The downward shift can be interpreted as a thinning of the buffer layer. Differently, 
$u^+$ shifts slightly upwards in the buffer layer and above for the case with light droplets, which is consistent with the case of drag reduction caused by the polymer in turbulent channel flow (Li, Sureshkumar & Khomami Reference Li, Sureshkumar and Khomami2006). The upwards shift can be interpreted as a thickening of the buffer layer. For the case with low-viscosity light droplets, 
$u^+$ shifts upwards again due to their stronger drag reduction effect than that of the light droplets. In addition, due to the effect of the low viscosity of the droplets within the viscous sublayer, 
$u^+=y^+$ is no longer valid.
Figure 10. Mean azimuthal velocity profiles near the inner cylinder. Here 
$u^+=(u_i-\langle u_{\theta } \rangle )/u_\tau$ is the velocity difference from the inner cylinder normalised by the friction velocity 
$u_\tau$ and 
$y^+=(r-r_i)/\delta _{v}$ is the distance from the inner cylinder in units of the viscous length scale. The dotted lines show the linear relation 
$u^+=y^+$ and the logarithmic law 
$u^+=(1/\kappa )\ln y^++B$ with the typical values 
$\kappa =0.4$ and 
$B=5.2$ (Huisman et al. Reference Huisman, Scharnowski, Cierpka, Kähler, Lohse and Sun2013).
Appendix D. The effect of effective Reynolds numbers
In the case of light droplets and low-viscosity light droplets, one may wonder that the drag reduction effect is due to the difference in effective Reynolds numbers. Therefore, it is necessary to correct the Reynolds number to exclude the possibility of drag reduction due to different effective Reynolds numbers. The effective density can be approximated using 
$\rho _\varphi = (1-\varphi )\rho _f+ \varphi \rho _d$. As for the effective viscosity, we note that there is no model available to estimate the effective viscosity of a liquid–liquid two-phase system. Our previous experimental results found that the effective viscosity increases with the volume fraction of the dispersed phase (Yi et al. Reference Yi, Toschi and Sun2021). At low volume fractions and low Reynolds numbers, the effective viscosity is very close to the model of Krieger and Dougherty (KD) (Krieger & Dougherty Reference Krieger and Dougherty1959). In our study, the volume fraction is in the range of 
$0\leq \varphi \leq 10\,\%$ and 
${Re}=6000$. We therefore employ the KD model to estimate the effective dynamic viscosity, i.e. 
$\mu _{\varphi } = \mu _f(1 - \varphi /\varphi _m)^{-2.5\varphi _m}$, where 
$\varphi _m = 0.58$. The Reynolds number corrected by the effective density and effective viscosity can be defined as 
${Re_{eff}}=\rho _\varphi u_i d/\mu _\varphi$, which is 5067 and 4389 for 
$\varphi =5\,\%$ and 
$\varphi =10\,\%$, respectively. The effective Nusselt number can be defined as 
$Nu_{eff}=T/(4 {\rm \pi}\mu _{\varphi } L r_i^2 r_o^2 \omega _i / (r_o^2-r_i^2))$ and the uncorrected drag modulation based on the effective Nusselt number can be written as 
$Nu_{eff}/Nu^{\varphi =0} -1 = (T/T_{\varphi =0})(\mu _f / \mu _\varphi )-1$ as listed in table 5, where 
$Nu^{\varphi =0}=T_{\varphi =0}/(4 {\rm \pi}\mu _f L r_i^2 r_o^2 \omega _i / (r_o^2-r_i^2))$ is the Nusselt number for the single-phase case at 
${Re}=6000$. Since the effective Reynolds number in the two-phase case is different from the single-phase case, we have used 
$T/T_{\varphi =0}-1$ rather than 
$Nu_{eff}/Nu^{\varphi =0} -1$ to characterise the drag modulation as listed in table 1. The corrected drag modulation based on the effective Nusselt number can be obtained by 
$Nu_{eff}/Nu_{eff}^{\varphi =0} -1$, where 
$Nu_{eff}^{\varphi =0}$ is the Nusselt number for the single-phase case at the effective Reynolds number. To obtain the corrected drag modulation, we additionally simulate two single-phase cases at these two effective Reynolds numbers (
${Re_{eff}}=5067$ and 4389) and obtain the corrected drag modulation by 
$Nu_{eff}/Nu_{eff}^{\varphi =0} -1$. See figure 11 for an illustration of the correction procedure, and the drag modulation before and after the correction is listed in table 5. Although the magnitude of drag reduction is reduced when the different Reynolds numbers between the two-phase and single-phase cases are taken into account, the drag reduction effect by the light droplets and low-viscosity light droplets is still significant. Therefore, we conclude that the drag reduction cannot be attributed to the discrepancy in the effective Reynolds number.
Table 5. Drag modulation of light droplets and low-viscosity light droplets.
Figure 11. The correction of the drag modulation characterised by the effective Nusselt number for the case with 
$5\,\%$ low-viscosity light droplets.
