How interfacial tension enhances drag in turbulent Taylor–Couette flow with neutrally buoyant and equally viscous droplets

Abstract

The presence of dispersed-phase droplets can result in a notable increase in a system's drag. However, our understanding of the mechanism underlying this phenomenon remains limited. In this study, we use three-dimensional direct numerical simulations with a modified multi-marker volume-of-fluid method to investigate liquid–liquid two-phase turbulence in a Taylor–Couette geometry. The dispersed phase has the same density and viscosity as the continuous phase. The Reynolds number $Re\equiv r_i\omega _i d/\nu$ is fixed at 5200, the volume fraction of the dispersed phase is up to $40\,\%$, and the Weber number $We\equiv \rho u^2_\tau d/\sigma$ is approximately 8. It is found that the increase in the system's drag originates from the contribution of interfacial tension. Specifically, droplets experience significant deformation and stretching in the streamwise direction due to shear near the inner cylinder. Consequently, the rear end of the droplets lags behind the fore head. This causes opposing interfacial tension effects on the fore head and rear end of the droplets. For the fore head of the droplets, the effect of interfacial tension appears to act against the flow direction. For the rear end, the effect appears to act in the flow direction. The increase in the system's drag is attributed primarily to the effect of interfacial tension on the fore head of the droplets which leads to the hindering effect of the droplets on the surrounding continuous phase. This hindering effect disrupts the formation of high-speed streaks, favouring the formation of low-speed ones, which are generally associated with higher viscous stress and drag of the system. This study provides new insights into the mechanism of drag enhancement reported in our previous experiments.

1. Introduction

Two-phase flows composed of two immiscible and incompressible liquids are ubiquitous in fields such as petroleum, food, pharmaceuticals and cosmetics (Spernath & Aserin 2006; McClements 2007; Wang et al. 2007; Kilpatrick 2012). The last few years have witnessed a renewed interest in two-phase turbulence (Rosti, Brandt & Mitra 2018; Mathai, Lohse & Sun 2020; Rosti & Takagi 2021; Wang et al. 2022; Yi et al. 2023; Girotto et al. 2024; Ni 2024). The presence of a droplet phase can significantly alter the momentum transport within the flow field, leading to changes in the system's drag. However, the underlying mechanism of droplet-induced drag modulation has not been clearly understood and deserves further investigation, especially in the turbulent regime.

For a liquid–liquid two-phase system, the presence of the droplet phase primarily introduces the influence of the two-phase interface. Consequently, existing studies have focused on examining the effect of the two-phase interface on the system's drag, with particular attention paid to interfacial dynamics, including deformation, breakup and coalescence, in shear turbulence. Interface-resolved direct numerical simulations were employed to study the effect of large deformable droplets on drag enhancement in turbulent channel flow using the interface-capturing method (Scarbolo & Soldati 2015; Scarbolo, Bianco & Soldati 2016). These studies found that droplet deformability is crucial for droplet-induced drag enhancement; specifically, the stronger the droplet deformability, the weaker the drag enhancement effect. Interface-capturing simulations of two-phase flow in the Stokes regime were conducted to investigate the effect of coalescence by introducing a short-range repulsive force to prevent droplet merging (de Vita et al. 2019). By comparing the cases allowing droplets to coalesce numerically with the cases using the repulsive force to prevent droplet merging, it was found that droplet coalescence effectively decreases the interfacial surface area, thereby weakening drag enhancement. Conversely, droplet breakup results in higher drag within the system. Interface-capturing simulations of homogeneous isotropic turbulence have also garnered widespread attention for studying the interface's effect on turbulence (Dodd & Ferrante 2016; Maxey 2017; Mukherjee et al. 2019). Dodd & Ferrante (2016) investigated how droplet deformation, breakup and coalescence affect the temporal evolution of turbulent kinetic energy. They showed that droplet coalescence reduces the total interfacial surface area, causing a decrease in surface energy and an increase in local kinetic energy. Recently, Rosti et al. (2019) demonstrated that a statistically stationary state (i.e. a balance between coalescence and breakup rates, and convergence of energy balance) can be reached in homogeneous shear turbulence. The statistically stationary state and the balance between coalescence and breakup rates have also been observed numerically in homogeneous isotropic turbulence (Mukherjee et al. 2019) and wall-bounded turbulence (Soligo, Roccon & Soldati 2019). Therefore, from the perspective of interfacial dynamics, droplet deformation plays a key role in drag enhancement.

Related experiments have also investigated drag modulation in liquid–liquid two-phase turbulence (Bakhuis et al. 2021; Yi, Toschi & Sun 2021; Wang et al. 2022; Yi et al. 2022). In Taylor–Couette turbulence, it has been observed experimentally that when breakup events dominate at low droplet volume fractions, dispersed droplets exhibit a specific size distribution that is well described by the log-normal distribution (Yi et al. 2021). At a fixed Reynolds number, droplets maintain nearly the same average size, while the system's drag increases with the droplet volume fraction (Yi et al. 2022). The effective viscosity as a function of volume fraction shows an increasing trend with the dispersed phase volume fraction, similar to that observed in rigid particle suspensions at moderate volume fractions ($\phi =0\unicode{x2013}40\,\%$) (Krieger & Dougherty 1959), although the effective viscosity for rigid particles consistently exceeds that for dispersed droplets (Yi et al. 2021). Numerical results from de Vita et al. (2019) have recently reported a similar phenomenon. However, the underlying mechanisms remain elusive.

Given the challenge of obtaining sufficient experimental data, interface-resolved simulations of liquid-liquid Taylor–Couette flow at Reynolds number $Re=960$ were conducted by Hori et al. (2023). They reported that the system can be categorized into two regimes based on the Weber number $We$: advection-dominated and interface-dominated regimes for high and low $We$, respectively. In the advection-dominated regime, drag enhancement as a function of volume fraction shows a non-monotonic behaviour. This could be attributed partly to the lower Reynolds number in their study compared to experimental conditions (Yi et al. 2021), and the significant numerical challenges related to coalescence encountered with the adapted interface-capturing method used for the volume fractions studied (Elghobashi 2019). Recent works (Su et al. 2024a,b) investigated numerically the effect of drops with varying density and viscosity in turbulent Taylor–Couette flow. It was found that interfacial tension consistently enhances momentum transport, thereby contributing to drag enhancement. However, the specific mechanism by which interfacial tension induces drag enhancement remains unclear.

In the interface-capturing method, merging of interfaces occurs automatically whenever two interfaces come within one grid cell of each other (Elghobashi 2019). This makes it very difficult to simulate droplets with moderate or high volume fractions. Various methods have been employed to address this issue, such as adaptive grid refinement (Innocenti et al. 2021), film drainage models (Thomas, Esmaeeli & Tryggvason 2010), artificial repulsive forces (de Vita et al. 2019) and multi-marker methods (Coyajee & Boersma 2009). However, these methods are effective primarily when dealing with a small number of droplets. Simulating two-phase turbulence with hundreds or more droplets becomes prohibitively expensive due to significantly increased computational costs. In this work, we investigate liquid–liquid Taylor–Couette turbulence at moderate volume fractions ($\phi =0\unicode{x2013}40\,\%$) using a modified multi-marker volume-of-fluid method. This modified approach allows us to reproduce the experimental drag modulation results reported in Yi et al. (2021). Our objective is to gain an intuitive understanding of how the droplet phase modulates the system's drag.

The paper is organized as follows. In § 2, we describe the numerical method and set-up. In § 3, we discuss the effect of droplets on the angular velocity flux, and analyse the modulation mechanism of viscous stress and angular velocity within the boundary layer. Finally, conclusions are drawn in § 4.

2. Numerical method and setting

The droplet and continuous phases are considered to be immiscible, incompressible and Newtonian. The two-phase immiscible and incompressible flow is governed by the Navier–Stokes equations

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

(2.2)$$\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}_\sigma, \end{gather}$$

where $\boldsymbol {u}$ is the velocity field, $p$ is the pressure, $\boldsymbol {\tau } = \mu (\boldsymbol {\nabla }\boldsymbol {u} + (\boldsymbol {\nabla }{\boldsymbol {u}})^{\rm T})$ is the viscous stress tensor, and $\boldsymbol {f}_\sigma$ represents the interfacial tension. Here, $\rho$ and $\mu$ denote the density and viscosity of the combined phase. In the standard volume-of-fluid (VOF) method, merging of interfaces occurs automatically whenever two interfaces come within one grid cell of each other, a phenomenon known as numerical coalescence (Elghobashi 2019; Soligo, Roccon & Soldati 2021). This makes it challenging to simulate droplets with moderate to high volume fractions. To address this issue, this work employs a modified multi-marker VOF method that was originally proposed by Coyajee & Boersma (2009) to prevent numerical coalescence. Specifically, multiple markers $\alpha _m$ are introduced to represent each group of droplets, where each droplet group is individually marked. The range of $\alpha _m$ is from 0 to 1: $\alpha _m=0$ represents the continuous phase, $\alpha _m=1$ represents the dispersed phase, and $0<\alpha _m<1$ denotes the interface region. This modified multi-marker VOF method enables the simulation of two-phase flow at moderate droplet volume fractions. The evolution of $\alpha _m$ is governed by the transport equation

(2.3)\begin{equation} {\partial _t \alpha_m} +\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha_m \boldsymbol{u})=0, \quad m=1,\ldots,n, \end{equation}

where $n$ is the total number of markers. The interfacial tension is calculated by

(2.4)\begin{equation} \boldsymbol{f}_\sigma=\sum_{m=1}^n\sigma\kappa_m\,\boldsymbol{\nabla}\alpha_m, \end{equation}

where $\sigma$ is the surface tension coefficient, and $\kappa _m = -\boldsymbol {\nabla }\boldsymbol {\cdot }({\boldsymbol {\nabla }\alpha _m}/{|\boldsymbol {\nabla }\alpha _m|})$ is the interface curvature. The density $\rho$ and viscosity $\mu$ of the combined phase are defined as functions of the phase fraction $\alpha$; specifically, $\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 the continuous phase and dispersed phase, respectively. The phase fraction $\alpha$ is defined as $\alpha = \max \{\alpha _1, \ldots, \alpha _n\}$ (Coyajee & Boersma 2009). In this work, we use $\rho =\rho _f=\rho _d$ and $\mu =\mu _f=\mu _d$.

We consider a two-phase flow in Taylor–Couette (TC) turbulence, where the flow is confined between two coaxial cylinders with radii $r_i$ (inner) and $r_o$ (outer). The curvature of the TC system is characterized by the ratio $\eta =r_i/r_o=0.714$. The outer cylinder (OC) is fixed, while the inner cylinder (IC) rotates at a constant angular velocity $\omega _i$. The torque $T$ required to drive the IC is examined to study the droplet-induced drag enhancement. To minimize computational costs without compromising accuracy, we employ a rotational symmetry of order 6 (i.e. the azimuthal angle of the simulated domain is ${\rm \pi} /3$) and aspect ratio $\varGamma =L/d=2{\rm \pi} /3$ in the simulated TC system, where $d$ is the gap width between the cylinders, and $L$ represents the axial length. This choice has been validated for both single-phase and multi-phase TC turbulence (Brauckmann & Eckhardt 2013; Spandan, Verzicco & Lohse 2018; Assen et al. 2022).

In this work, liquid–liquid TC turbulence is simulated with total droplet volume fractions $\phi =0$, $10\,\%$, $20\,\%$, $30\,\%$ and $40\,\%$. The Reynolds number $Re\equiv {\rho } u_i d/\mu$ is fixed at 5200, where $u_i = r_i \omega _i$ is the velocity of the IC. The Taylor number is $Ta \equiv \chi (r_o+r_i)^2(r_o-r_i)^2\omega _i^2/(4\nu ^2) = 4.12\times 10^7$, where $\chi = [(r_i+r_o)/(2\sqrt {r_i r_o})]^4$ and $\nu =\mu /\rho$. The frictional Reynolds number at the IC, $Re_\tau \equiv u_\tau d/\nu$, is 264.86, 281.62, 297.01, 315.37 and 329.76 for $\phi =0$, $10\,\%$, $20\,\%$, $30\,\%$ and $40\,\%$, respectively. Here, $u_\tau$ is the friction velocity, defined as $\sqrt {\tau _{w}/\rho }$, where $\tau _{w}$ represents the shear stress at the IC. The Weber number $We \equiv \rho u_\tau ^2 d/\sigma$ is 7.96, 8.74. 9.85 and 10.77 for $\phi =10\,\%$, $20\,\%$, $30\,\%$ and $40\,\%$, respectively.

No-slip and impermeable boundary conditions are imposed on both cylinder surfaces, while periodic boundary conditions are applied in the axial and azimuthal directions. The IC and OC are subjected to Neumann boundary conditions for the phase fraction, resulting in a default contact angle $90^\circ$. The TC system is discretized using a collocated grid system consisting of $N_\theta \times N_r \times N_z = 336 \times 256 \times 192$ points in the azimuthal, radial and axial directions, respectively. The grids are uniformly distributed 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 /u_\tau$ for single-phase flow. In the radial direction, the grid spacing varies from $0.31\delta _\nu$ near the wall to $2.45\delta _\nu$ at the centre of the gap. In the azimuthal direction, it varies from $2.06\delta _\nu$ near the IC to $2.89\delta _\nu$ near the OC. The grid spacing is uniform in the axial direction, with value $2.89\delta _\nu$. The Kolmogorov scale, denoted as $\eta _k$, is $2.07\delta _\nu$ for single-phase turbulence, determined by the exact dissipation relationships given by $\eta _k/d=(\chi ^{-2}\,Ta\,(Nu_{\omega }-1))^{-1/4}$, where $Nu_\omega =J^{\omega }/J^{\omega }_{lam}$ (Eckhardt, Grossmann & Lohse 2007). Here, $J^{\omega }$ represents the total angular velocity flux, and $J^{\omega }_{lam}$ corresponds to the flux under fully laminar and non-vortical conditions. The maximum grid spacing is $1.39\eta _k$ for single-phase turbulence. To ensure sufficient spatial resolution to resolve the smallest length scales, a resolution test is performed for the single-phase case (see Appendix A).

These simulations utilize the modified multi-marker VOF method with a piecewise-linear interface calculation (PLIC) algorithm implemented in the interFoam solver of the open-source OpenFOAM v8 (Rusche 2003; Chen, Zhao & Wan 2022). Based on the PLIC algorithm, the resolved interface region ($0<\alpha <1$) could be confined within a single layer of grid cells between the two phases (Su et al. 2024b). Therefore, the PLIC algorithm also works in reducing the influence of numerical coalescence. The robustness of OpenFOAM in simulating both single-phase and two-phase TC turbulence has been demonstrated in our previous works (Xu et al. 2022, 2023; Su et al. 2024b). The maximum Courant–Friedrichs–Lewy number is set to 0.2. Temporal discretization employs a blended scheme between the first-order Euler scheme and the second-order Crank–Nicolson scheme, with blending factor 0.9 for robustness and accuracy. Spatial discretization uses a second-order linear upwind scheme for the advection term in the momentum equation. The phase fraction transport equation employs a PLIC scheme to maintain interface sharpness. The PIMPLE algorithm (Holzmann 2016), which is a hybrid version of the PISO algorithm and the SIMPLE algorithm, is used to handle the pressure–velocity coupling to guarantee better stability for two-phase simulations. The pressure equation is solved using the geometric algebraic multigrid solver coupled with the simplified diagonal-based incomplete Cholesky smoother, which is commonly used to speed up the computational efficiency in simulating two-phase flow (Scheufler & Roenby 2019; Chen et al. 2022). Velocity and phase fraction are found using an iterative solver with a symmetric Gauss–Seidel smoother. In the simulations, a residual tolerance $10^{-6}$ is maintained for all variables, except for the phase fraction, which has tolerance $10^{-8}$. The computational accuracy of these settings has been verified by comparing our results with those of Ostilla et al. (2013) in our previous work (Su et al. 2024a,b).

For a fixed volume fraction $\phi$ in the modified multi-marker VOF method, the dispersed phase is divided into $n$ groups, and each group is assigned a marker $\alpha _m$, and the droplet volume fraction is the same for each marker $\phi _{mvf}$. A large value of $\phi _{mvf}$ results in unphysical numerical coalescence, while a too small value of $\phi _{mvf}$ limits the maximum droplet size and is computationally costly as one needs to solve $n$ transport equations. In our study, we choose $\phi _{mvf}=5\,\%$ for all cases – specifically, $n=2$, 4, 6 and 8 for droplet volume fractions $10\,\%$, $20\,\%$, $30\,\%$ and $40\,\%$, respectively. This choice is inspired by previous studies (Rosti et al. 2019; Soligo et al. 2019; Crialesi-Esposito et al. 2022; Mangani et al. 2024), which used the standard VOF method or phase-field method. In these studies, the volume fraction is mainly in the range $3\unicode{x2013}10\,\%$, where the effect of numerical coalescence is considered to be negligible. The choice $\phi _{mvf}=5\,\%$ faithfully reproduces the global drag of the system reported in our previous experimental study (Yi et al. 2021). To ensure that the modified multi-marker VOF method does not qualitatively alter the drag modulation mechanism, two additional simulations using the standard VOF method are conducted and compared with the multi-marker VOF method (see Appendix B). During our simulations, a single-phase case is first simulated to initialize the velocity field. Once a fully developed flow with a pair of Taylor rolls is obtained, droplets of diameter $0.2d$ are uniformly positioned in the simulated domain $\varOmega$. The droplets are then randomly marked by the marker $\alpha _m$ ($m=1,\ldots, n$) in such a way that the total volume of the droplets for each marker is $\phi _{mvf}V_\varOmega$, where $V_\varOmega$ is the volume of the simulated domain. All statistics presented in the paper are collected over at least $3 \times 10^2$ large eddy turnover times, defined as $(r_o-r_i)/(\omega _i r_i)$, after the two-phase flow reaches a statistically steady state.

To demonstrate that the flow reaches a statistically steady state, the temporal evolution of the interfacial surface area $S$ is shown in the inset of figure 1(a). This $S$ fluctuates around a constant value dependent on $\phi$, indicating a balance between droplet breakup and coalescence within the system. In our simulations, the same grid system is employed for both single-phase and two-phase cases, and a resolution test is conducted for both the $\phi =0$ and $\phi =10\,\%$ scenarios (see Appendix A). Due to the significant increase in computational cost with the multi-marker VOF method, conducting a resolution test for $\phi =40\,\%$ becomes impractical. The corresponding grid spacing near the wall is approximately $1.39\eta _k$ for $\phi =0$, and $1.58\eta _k$ for $\phi =40\,\%$, indicating only mild degradation in spatial resolution with increasing $\phi$. In figure 1(a), we present the torque $T$ required to rotate the IC at a constant rate $\omega _i$, superimposed with experimental data from Yi et al. (2021) for comparison purpose. It is evident that $T$ increases with $\phi$ in both numerical and experimental data, indicating drag enhancement due to dispersed interfaces. The datasets agree well with each other at $\phi \leqslant 30\,\%$, with a minor deviation observed at $\phi =40\,\%$. However, the method and grid allocation adopted here suffice for the purpose of this study, i.e. revealing the underlying mechanism of drag enhancement by dispersed droplets.

Figure 1. (a) The torque $T$ needed to drive the IC rotating at a constant rate $\omega _i$. Here, $T$ is normalized by the single-phase torque $T_{\phi =0}$. The simulated results are compared with experimental results from our previous work (Yi et al. 2021). The inset shows the interfacial surface area $S/S_{cyl}$ as a function of time $t/t_0$ at various droplet volume fractions $\phi$, where $S_{cyl}$ is the surface area of the IC, and $t_0$ is the large eddy turnover time. (b) The turbulent stress $J^{\omega }_{T}$, viscous stress $J^{\omega }_{V}$ and interfacial stress $J^{\omega }_\sigma$ contributions to the angular velocity flux. The operator $\sum$ denotes a radially averaged quantity. All the contributions are normalized by the total angular velocity flux for the single-phase flow, $J_{\phi =0}^{\omega }$.

3. Results and discussion

In this work, we fix the OC while sustaining the constant rotational velocity of the IC. The torque $T$ required to drive the IC is examined to study the underlying mechanism of drag enhancement by dispersed droplets. In the TC system, the angular velocity flux $J^{\omega }$ is conserved in the radial direction. Here, $J^\omega$ and the torque at IC $T$ are related through the equation $T=2{\rm \pi} L J^\omega$ (Eckhardt et al. 2007), which provides an efficient way to study the mechanism of the drag enhancement. Inspired by an analogous situation in turbulent channel flows (Picano, Breugem & Brandt 2015; Wang, Jiang & Sun 2023), where it is conventional to decompose the total stress into the viscous part, turbulent part and others, we decompose $J^\omega$ into contributions from the turbulent stress $\tau _{T}$, the viscous stress $\tau _{V}$, and the interfacial tension $\boldsymbol {f}_\sigma$ (Eckhardt et al. 2007; Su et al. 2024b), i.e.

(3.1)\begin{equation} J^\omega\equiv J_{T}^\omega(r) + J_{V}^\omega(r)+J_{\sigma}^\omega(r)= {\rm{const.}}, \end{equation}

where the three terms represent

1. (i) the turbulent stress contribution $J_{T}^\omega (r) =r^3 \langle \rho u^{\prime }_r \omega ^{\prime }\rangle = r^2 \langle \tau _{T}\rangle$,

2. (ii) the viscous stress contribution $J_{V}^\omega (r) = -\mu r^3\,\partial _r \langle \omega \rangle = r^2 \langle \tau _{V}\rangle$,

3. (iii) the interfacial tension contribution $J_{\sigma }^\omega (r) =-\int _{r_i}^{r} \langle {r^2 f_\sigma ^\theta } \rangle \, {\rm {d}}r$.

Here, $u^{\prime }_r=u_r-\langle u_r \rangle$, $\omega ^{\prime }=\omega -\langle \omega \rangle$, $r$ is the radial position, $u_r$ is the radial velocity, $\omega$ is the angular velocity, and $f_\sigma ^\theta$ is the azimuthal component of interfacial tension. The operator $\langle \cdot \rangle$ denotes the average in the axial and azimuthal directions and over time. In this work, all contributions are normalized by the total angular velocity flux for single-phase flow $J_{\phi =0}^{\omega }$. We note that $J_{T}^\omega (r)$, $J_{V}^\omega (r)$ and $J_\sigma ^\omega (r)$ all depend on the radial position. To eliminate the radial dependence, $J_{T}^\omega (r)$, $J_{V}^\omega (r)$ and $J_\sigma ^\omega (r)$ are further averaged in the radial direction, denoted as $\sum J_{T}^\omega$, $\sum J_{V}^\omega$ and $\sum J_\sigma ^\omega$, and are depicted in figure 1(b). This approach has been used widely to study drag modulation in two-phase flows (Ardekani & Brandt 2019; de Vita et al. 2019; Hori et al. 2023; Wang et al. 2023). Both $\sum J_{T}^\omega$ and $\sum J_{V}^\omega$ are virtually not changed by the dispersed droplets. However, $\sum J_\sigma ^\omega (r)$ increases monotonically with $\phi$. Figure 1(b) demonstrates that the drag enhancement is dominated by the interfacial tension contribution, consistent with our previous findings (Su et al. 2024a,b).

In the above analysis, drag enhancement is attributed to the additional contribution of interfacial tension to the angular velocity flux. However, it remains unclear how the introduction of the interface increases the overall drag of the system. Figure 1(b) gives us the impression that viscous stress plays a negligible role in transporting the angular velocity flux. However, this is not the case, especially in the boundary layer, where viscous stress dominates over turbulent stress. It should be noted that averaging in the radial direction masks the important role of the viscous boundary layer, as it represents only a small fraction of the total volume. Taking the droplet volume fraction $\phi =20\,\%$ as an example, $J_{V}^\omega /J^\omega _{\phi =0}$ is much larger than $J_{T}^\omega /J^\omega _{\phi =0}$ and $J_\sigma ^\omega /J^\omega _{\phi =0}$ near the IC, as shown in figure 2(a). Due to the presence of a positive interfacial tension contribution near the OC, a smaller viscous stress contribution ($J_{V}^{\omega }/J_{\phi =0}^{\omega }<1$) is observed, which has also been reported by Hori et al. (2023). In our previous study, we observed that droplets are fragmented by dynamic pressure within the boundary layer (Yi et al. 2022). Given the higher observed viscous stress contribution compared to the single-phase case at the IC surface (figure 2a), further analysis of the boundary layer near the IC may yield insights into the underlying mechanism of drag enhancement by dispersed droplets.

Figure 2. (a) Angular velocity flux and its three contributions as functions of the radial position for the case with $\phi =20\,\%$. The radial position $(r-r_i)/d=0$ corresponds to the IC, and $(r-r_i)/d=1.0$ corresponds to the OC. (b) Viscous stress $\langle \tau _{V} \rangle$ (normalized by the single-phase value at $(r-r_i)/d=0$, i.e. $\langle \tau _{i,\phi =0} \rangle$) as a function of the radial position at various droplet volume fractions. The inset shows a zoom-in of $\langle \tau _{V} \rangle / \langle \tau _{i,\phi =0} \rangle$ near the IC.

Since the viscous stress contribution $J_{V}^{\omega }(r)$ for a specific radial position depends solely on the viscous stress $\tau _{V}$, we will now focus our attention on $\tau _{V}$ for ease of discussion. We show the $\tau _{V}$ distribution at various droplet volume fractions in figure 2(b), where an enlarged view of $\tau _{V}$ within the region $(r-r_i)/d<0.02$ is displayed in the inset. Here, $(r-r_i)/d=0.02$ corresponds to $y^+=5.31$, 5.61, 5.94, 6.31 and 6.60 for cases with $\phi =0$, $10\,\%$, $20\,\%$, $30\,\%$ and $40\,\%$, respectively. Also, $y^+$ is the distance from the IC in units of the viscous length scale at the wall. In the boundary layer, $\tau _{V}$ increases with increasing volume fraction, in line with the drag enhancement observed in figure 1.

Given the definition of viscous stress $\tau _{V} = -\mu r\,\partial _r \omega$, the magnitude of viscous stress is determined by the angular velocity gradient when viscosity and $r$ are fixed. Considering that the IC has a constant angular velocity, it is likely that in the boundary layer, the viscous stress is directly related to the angular velocity. Figures 3(a) and 3(b) show contour plots of the instantaneous angular velocity $\omega$ and viscous stress $\tau _{V}$ on a cylinder surface with $(r_{cut}-r_i)/d=0.0105$ for the $\phi =40\,\%$ case. By carefully examining their spatial distribution, it can be observed that the low-speed streak regions (blue in figure 3a) predominantly correspond to the high viscous stress regions (red in figure 3b). To provide quantitative evidence for this observation, we calculate the joint probability density function (p.d.f.) between $\tau _{V} / \langle \tau _{i,\phi =0}\rangle$ and $\omega / \omega _i$. The joint p.d.f. basically lies along the straight dashed line with a negative slope (see figure 3c), suggesting a negative correlation between $\tau _{V} / \langle \tau _{i,\phi =0}\rangle$ and $\omega / \omega _i$: namely, lower velocity fluid can result in higher viscous stress. As the volume fraction increases, the p.d.f.s of the viscous stress shift rightwards, implying an overall increase in $\tau _{V}$ (see figure 3d). Based on the negative correlation between $\tau _{V} / \langle \tau _{i,\phi =0}\rangle$ and $\omega / \omega _i$, the rightward shift of $\tau _{V}$ means that the angular velocity decreases. Here, we reveal the negative correlation between viscous stress and angular velocity within the boundary layer, where an increase in viscous stress manifests itself as a decrease in angular velocity. Figure 3(d) demonstrates that an overall increase in viscous stress is observed with increasing droplet volume fraction, which would correspond to a decrease in angular velocity. Therefore, understanding the reason for the decrease in angular velocity within the boundary layer becomes crucial for the study of drag enhancement.

Figure 3. Contour plots of (a) the instantaneous angular velocity and (b) viscous stress on the cylinder surface with radius $(r_{cut}-r_i)/d=0.0105$ for two-phase turbulence with $\phi =40\,\%$. Here, $(r_{cut}-r_i)/d=0.0105$ corresponds to $y^+=3.45$. (c) Joint p.d.f. between the angular velocity and the viscous stress for $\phi =40\,\%$. (d) The p.d.f.s of $\tau _{V}$ for different droplet volume fractions. In (c,d), the data are sampled on the same cylinder surface as in (a,b).

To visually compare the angular velocity distribution in single-phase and two-phase turbulence, contour plots of the instantaneous angular velocity for cases $\phi =0$ and $40\,\%$ are illustrated in figures 4(a) and 4(b), respectively. Herring-bone streaks can be observed in single-phase turbulence due to centrifugal instability (Barcilon et al. 1979; Dong 2007). Compared to single-phase turbulence, two-phase turbulence exhibits numerous low-speed streaks. The high-speed streaks are disrupted by the low-speed ones, and their axial size diminishes. Droplet interfaces ($\alpha =0.5$) are depicted as solid lines in figure 4(b). Upon careful inspection, it is observed that droplets occupy predominantly the low-speed streaks. This observation suggests that droplets tend to favour the formation of low-speed streaks. Subsequently, we demonstrate that the presence of droplets indeed impedes the continuous phase, disrupts the high-speed streaks, and promotes the low-speed ones.

Figure 4. Contour plots of the instantaneous angular velocity $\omega /\omega _i$ on a cylinder surface with $(r_{cut}-r_i)/ d=0.02$ for (a) single-phase turbulence with $\phi =0$, and (b) two-phase turbulence with $\phi =40\,\%$. The droplet interfaces ($\alpha =0.5$, solid lines) are also superimposed in (b). (c) The p.d.f.s of angular velocity for the droplet and continuous phases at different cylinder surfaces for the $\phi =40\,\%$ case. The inset shows the average angular velocity of the droplet and continuous phases as a function of wall distance. (d) The p.d.f.s of angular velocity for the droplet and continuous phases at fixed cylinder surface $(r_{cut}-r_i)/d=0.02$ for various droplet volume fractions. The inset shows the average angular velocity of the droplet and continuous phases.

The p.d.f.s of angular velocity for both the continuous and droplet phases are shown in figures 4(c) and 4(d). In figure 4(c), the droplet volume fraction is $\phi =40\,\%$, while the data are sampled on different cylinder surfaces $(r_{cut}-r_i)/d=0.0054$, 0.0105, 0.0159 and 0.02, which correspond to $y^+=1.77$, 3.45, 5.25 and 6.60. The first two surfaces are located inside the viscous sublayer, and the last two are in the buffer layer. In figure 4(d), the data are sampled on the cylinder surface with $(r_{cut}-r_i)/d=0.02$, while two-phase turbulence with different volume fractions is considered. It is observed that the angular velocity of the droplet phase is consistently lower than that of the continuous phase (see figure 4(c) and its inset). The observation persists across various volume fractions investigated here (see figure 4d). Besides, at the fixed radial position $(r_{cut}-r_i)/d=0.02$, the average angular velocity of the droplet phase is nearly invariant at different droplet volume fractions, whereas the average angular velocity of the continuous phase decreases as the droplet volume fraction increases (see inset of figure 4d). The velocity of the droplets being lower than in the continuous phase suggests that the droplet phase will impede the continuous phase surrounding the droplets.

To explain why the droplet phase moves slower than the continuous phase near the IC, we show a cross-section of flow in the $r\unicode{x2013}\theta$ plane for the $\phi =40\,\%$ case in figure 5(a), and ideally divide a single droplet in the vicinity of the IC into two parts along the radial direction, based on the location of the droplet's centre of mass. The part closer to the IC is called fore head of the droplet, and the other part is called the rear end of the droplet. Due to the high shear near the IC, droplets are strongly deformed. However, the droplets do not align with the azimuthal shear direction but show a slight deviation. This deviation from the $\theta$ direction originates from the fact that a fluid parcel closer to the rotating IC will generally have a higher angular velocity. In other words, the fore head of the droplet moves faster than its rear end. To maintain the integrity of the droplet and resist deformation, the rear end will drag the fore head backwards due to the effect of interfacial tension. This process will make the fluid inside the fore head move at a slower velocity than the continuous phase in the same radial position. As a result, the dispersed droplets impede the continuous phase surrounding the fore head of the droplets. On the other hand, the rear end will be dragged forwards by the fore head due to the effect of interfacial tension, thus accelerating the continuous phase surrounding the rear end of the droplets. The accelerating effect of the rear end will be masked by the hindering effect of the fore head, which ultimately causes increased drag.

Figure 5. (a) Contour plot of the instantaneous angular velocity in the $r\unicode{x2013}\theta$ plane for the $\phi =40\,\%$ case. Droplet interfaces ($\alpha =0.5$) are represented by the solid lines. The inset is an enlarged view of the main image marked by the red rectangle. The droplet in the inset is divided into fore head and rear end along the radial direction based on the location of the droplet's centre of mass. (b) The p.d.f. of the sum of the azimuthal interfacial tension acting on the fore head $\sum _{\scriptscriptstyle {Fore}}{f_{\sigma }^{\theta }}$ and the sum of the azimuthal interfacial tension acting on the rear end $\sum _{\scriptscriptstyle {Rear}}{f_{\sigma }^{\theta }}$ is shown for droplets with centres of mass in the range $y^+<24$. The two terms in (b) are normalized by the absolute value of the total azimuthal interfacial tension experienced by all droplets with centres of mass in the range $y^+<24$. Here, the total azimuthal interfacial tension is negative, indicating that the overall effect of these droplets is to impede the surrounding flow.

To demonstrate our argument, we examine the interfacial tension experienced by droplets whose centres of mass fall within the range $y^+<24$. This choice ensures that the majority of the droplets analysed are in the viscous sublayer and buffer layer. For each individual droplet, we compute the sum of the azimuthal interfacial tension acting on the fore head, and plot the p.d.f. of these summed values for different droplets in figure 5(b). Similarly, we perform the analysis on the rear end of each individual droplet. In figure 5(b), the positive values indicate that the effect of interfacial tension appears to act in the flow direction, while the negative values suggest that the effect of interfacial tension appears to act against the flow direction. It can be observed that for the fore head of the droplet, the effect of interfacial tension acts predominantly against the flow direction, which slows down the fluid inside the fore head, thus hindering the surrounding flow. Conversely, for the rear end, the effect of interfacial tension acts predominantly in the flow direction, which drives the fluid within the rear end to move faster, thus accelerating the surrounding flow. To visualize the physical process in which the interfacial tension works, a sketch is shown in figure 6. Here, the effect of interfacial tension acting against the flow direction overwhelms that acting in the flow direction. The droplet-induced drag enhancement is dominated by the effect of interfacial tension on the fore head of the droplet, which leads to a hindering effect of the droplet on the surrounding continuous phase.

Figure 6. A sketch of how the interfacial tension works in a TC system.

Having gained an insight into how interfacial tension works, we are now able to comprehend the role of interfacial tension contribution $J^\omega _\sigma (r)$ to the total angular velocity flux, which is plotted in figure 7(a) for various volume fractions. Although the above discussion focuses on the region near the IC, the analysis of how interfacial tension works can be extended to the entire TC system. Note that in this paper, the fore head always refers to the part of the droplet closer to the IC. As shown in the inset of figure 7(a), the interfacial tension is negative in the region $(r-r_{i})/d<0.02$, i.e. the effect of interfacial tension acting against the flow direction dominates in this region, as just discussed. Specifically, the effect of interfacial tension results in a reduction in average angular velocity with increasing volume fraction in the region $(r-r_{i})/d<0.02$ (see the inset of figure 7b). This leads to an increase in viscous stress with increasing volume fraction as shown in the inset of figure 2(b). According to $J_{\sigma }^\omega (r) =-\int _{r_i}^{r} \langle {r^2 f_\sigma ^\theta } \rangle \, {\rm d}r$, $J^\omega _\sigma (r)$ will increase until $\langle\, f^\theta _\sigma \rangle$ changes its sign from negative to positive, where $J^\omega _\sigma (r)$ attains its local maximum. After the peak, $J^\omega _\sigma (r)$ exhibits a slight decrease. The decrease of $J^\omega _\sigma (r)$ indicates that $\langle\, f^\theta _\sigma \rangle >0$. Near the IC, the droplets are highly deformed. When they are advected away from the IC, the droplets will relax back to a less deformed state. These cylindrical surfaces are filled mainly with the rear end of the highly deformed droplets and the fore head of the less deformed droplets. The slight decrease in $J^\omega _\sigma (r)$ can be attributed to the fact that the effect of interfacial tension experienced by the rear end of the highly deformed droplets is stronger compared to that of the interfacial tension experienced by the fore head of the less deformed droplets.

Figure 7. (a) Interfacial tension contribution $J_\sigma ^\omega (r)$ as a function of radial position at various droplet volume fractions. The inset shows the azimuthal component of the interfacial tension in the near-wall region, and the negative value indicates that the effect of interfacial tension acts against the flow direction. (b) Average angular velocity $\langle \omega \rangle / \omega _i$ as a function of radial position at various droplet volume fractions, with the inset showing the zoom-in of $\langle \omega \rangle / \omega _i$ near the IC.

In the bulk region, $J^\omega _\sigma (r)$ increases slightly with the radial position in the region $0.2\lesssim (r-r_{i})/d\lesssim 0.8$, indicating that the total effect of the interfacial tension in that region is to impede the local flow. Near the OC, droplets are highly deformed due to the high shear near the solid surface. In this case, the rear end of the droplet is closer to the OC where the angular velocity is zero, thus dragging the fore head backwards. This leads to $\langle\, f^\theta _\sigma \rangle <0$, and $J^\omega _\sigma (r)$ obviously increases with the radial position in the region $0.8\lesssim (r-r_{i})/d\lesssim 0.96$. On the other hand, the rear end will be dragged forwards by the fore head, thus resulting in $\langle\, f^\theta _\sigma \rangle >0$, and $J^\omega _\sigma (r)$ shows an obvious decrease with the radial position in the region $0.96\lesssim (r-r_{i})/d < 1$. At the OC with $(r-r_{i})/d = 1$, we have $J^\omega _\sigma (r)>0$, which indicates that the effect of interfacial tension acting against the flow direction dominates the modulation of the flow field within the cylinder gap. As a result, the overall reduction of the average angular velocity with increasing volume fraction is observed within the cylinder gap, as shown in figure 7(b).

Based on the above analysis, we have attributed the drag enhancement to the effect of interfacial tension experienced by the fore head of the droplet, which acts against the flow direction and leads to a hindering effect of the droplet on the surrounding flow. Although the division of the fore head and rear end based on the location of the droplet's centre of mass is somewhat idealized, it is sufficient to reveal the physical mechanism of the droplet-induced drag enhancement. The hindering effect is also evident in the mean azimuthal velocity. In figure 8(a), we present normalized mean azimuthal velocity profiles $u^+ = (u_i - \langle u_{\theta } \rangle ) / u_\tau$ versus wall distance $y^+$. In the single-phase case ($\phi =0$), $u^+$ follows the linear relation $u^+ = y^+$ effectively within the viscous sublayer ($\kern 1.5pt y^+ < 5$), indicating sufficient spatial resolution to resolve the boundary layer. At $y^+ > 30$, $u^+$ does not exhibit a clear logarithmic shape due to the low Reynolds number ($Re$) in this study (Huisman et al. 2013). In the two-phase cases, $u^+$ shifts downwards when $\phi > 10\,\%$. This shift is observed not only in the buffer layer and above, but also in the viscous sublayer, akin to drag enhancement observed with rough walls (Zhu et al. 2016; Xu et al. 2023). In the viscous sublayer, the magnitude of interfacial tension $\langle\, f^\theta _\sigma \rangle$ increases with increasing volume fraction (figure 8(b) and its inset), rendering $u^+ = y^+$ invalid, and indicating the significance of both viscous stress and interfacial tension. By examining the phase fraction distribution, we find that the droplet phase fraction can reach a maximum of $5\,\%$ inside the viscous sublayer. In addition, the magnitude of the azimuthal interfacial tension reaches its maximum at approximately $y^+=3$, indicating that the azimuthal interfacial tension has the greatest effect inside the viscous sublayer. This finding is consistent with previous experimental results of Yi et al. (2022). The authors found that the droplets are fragmented within the boundary layer. Hence the feedback of droplets will attain its maximum within the boundary layer. The downward shift of $u^+$ may be due to the increased effective viscosity by interfacial tension, causing $y^+$ to be overestimated. Future theoretical studies on turbulent boundary layers in two-phase flows should consider the impact of interfacial tension in the viscous sublayer.

Figure 8. (a) Mean azimuthal velocity profiles near the IC, where $u^+=(u_i-\langle u_{\theta } \rangle )/u_\tau$ is the velocity difference from the IC normalized by the friction velocity. The dashed 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. 2013). (b) Profiles for the azimuthal component of the interfacial tension $\langle\, f^\theta _\sigma \rangle$. The inset shows the average phase fraction within the viscous sublayer. The vertical dot-dashed line indicates the location of $y^+=5$.

Our findings emphasize the dominant role of interfacial tension in the drag enhancement caused by droplets. In industrial applications, many existing models (e.g. Euler–Lagrange approach) may not adequately account for the effects of interfacial tension. Our results show that to accurately predict drag and flow characteristics in two-phase turbulent flow systems, models must account for the contribution of interfacial tension, particularly in the boundary layer region. This provides guidance for the improvement and development of two-phase flow models.

The budget analysis adopted in this work is widely used to analyse the drag modulation in two-phase turbulence. In turbulent channel flow laden with dense suspensions (${\phi =0\unicode{x2013}20\,\%}$) of neutrally buoyant spheres, Picano et al. (2015) demonstrated that the drag enhancement is dominated by the particle-induced stress contribution. Interestingly, similar to the profile of $J^\omega _\sigma (r)$ near the IC, the particle-induced stress contribution is positive and increases with wall distance, then decreases further away from the wall (Picano et al. 2015; Wang, Abbas & Climent 2017). In turbulent flow laden with particles, the exact expression for the particle-induced stress remains unknown. Inspired by the interfacial tension contribution induced by the droplets, which can be considered as ‘deformable particles’, we can replace the particle-induced stress effect with the integral of the particle-induced force in the direction normal to the wall. This means that in the near-wall region, the particle-induced stress acts against the flow direction, thus hindering the surrounding continuous phase. Based on the above analysis, we infer that the drag enhancement induced by droplets and particles is dominated by similar mechanisms within the boundary layer. Our results provide one possible explanation for why neutrally buoyant solid particles have a stronger drag enhancement effect compared to droplets as found by de Vita et al. (2019) and Yi et al. (2021). The stronger drag enhancement effect indicates that the hindering effect on the surrounding flow field near the wall is more pronounced for solid particles, which causes larger viscous stress. This may be due to the fact that the deformability of the droplet allows it to adjust its shape when subject to shear flow, thus droplets affect the surrounding flow to a lesser degree. As a result, the drag enhancement caused by droplets is less pronounced than that caused by solid particles.

Given that our work is aimed primarily at reproducing the drag enhancement reported in our previous experimental study and revealing the underlying mechanism, the dimensionless parameters in our simulations mainly correspond to $Re=5.2 \times 10^3$ cases in experiments (Yi et al. 2021). Although we kept many fixed dimensionless parameters ($Re, Ta, \varGamma, \eta, We$, density ratio, viscosity ratio, contact angle) in our simulations, we can propose some new explanations for the effects of different parameters based on our findings. By changing the angular velocity of the IC, it is observed that the effective viscosity decreases with increasing $Re$ ($Ta$ or $We$), suggesting a shear thinning effect (Yi et al. 2021). This can be attributed to the fact that the effect of interfacial tension weakens as the turbulence intensity increases. Specifically, the increase in the near-wall viscous stress due to interfacial tension gradually diminishes relative to the increase in viscous stress due to the increase in turbulence intensity. The effects of the density ratio and viscosity ratio have been studied in our previous work (Su et al. 2024b). We found that decreasing the density and viscosity ratios of the dispersed phase to the continuous phase reduces the contribution of local advection and diffusion terms to the momentum transport, respectively, resulting in drag reduction. The change in contact angle may cause droplets to attach to the IC and form a lubricating layer, thus contributing to the drag reduction. Herein, the effect of the attached droplets on the system drag is dominated by the interfacial tension on the rear end, which appears to act in the flow direction. Specifically, the effect of interfacial tension tends to increase the angular velocity of the fluid near the interface of the attached droplets, and decrease near-wall viscous stress, thus contributing to drag reduction. Therefore, the droplet phase attached to the IC and dispersed within the system will compete to cause drag modulation. For the parameters characterizing the TC device ($\eta$ and $\varGamma$), the choice of aspect ratio $\varGamma =2{\rm \pi} /3$ is to minimize computational cost and corresponds to the cases where the axial length is much larger than the gap width in real experiments. The choice of curvature ratio $\eta =r_i/r_o=0.714$ may limit our findings to only TC systems with curvature ratios close to this value. Further studies are needed to extend the mechanisms of interface-induced drag modulation to TC systems with other curvature ratios and even to channel flow systems.

4. Conclusions

In this study, we investigated the mechanism of drag enhancement by neutrally buoyant droplets in liquid–liquid Taylor–Couette turbulence at Reynolds number 5200 and Weber number approximately 8. Our focus was on the effect of droplet presence within the boundary layer at moderate volume fractions. To achieve this, we employed a modified multi-marker volume-of-fluid method, which accurately reproduces the drag enhancement effect and aligns well with experimental results.

Through angular velocity flux analysis, we identified that interfacial tension plays a pivotal role in drag enhancement. To understand how interfacial tension influences drag enhancement, we examined the effects of droplets on viscous stress and angular velocity near the inner cylinder. Our observations revealed that droplets experience significant deformation and stretching along the streamwise direction due to shear near the inner cylinder. Consequently, the rear end of the droplets lags behind the fore head. To keep the integrity of the droplets and resist deformation, the effect of interfacial tension experienced by the fore head of the droplets appears to act against the flow direction, while the effect of interfacial tension experienced by the rear end appears to act in the flow direction. The effect of the interfacial tension acting against the flow direction overwhelms that of the interfacial tension acting in the flow direction, leading to a hindering effect of the droplets on the surrounding continuous phase. This hindering effect alters high-speed streaks, reducing their size and increasing the occurrence of low-speed streaks, which typically contribute to higher viscous stress and system drag. Furthermore, we observed that the mean streamwise velocity profile no longer follows the linear relation $u^+ = y^+$ in two-phase turbulence, which may be due to the fact that the interfacial tension increases the effective viscosity, resulting in an overestimation of $y^+$. Our findings underscore the necessity of adequately considering interfacial tension in the near-wall region when modelling two-phase turbulence.