3.1. Vorticity generation and the instability in the core region

Well within the core, we have $r\ll r_c$, so $q \ll 1$. Taylor expanding in this limit, after some algebra and retaining the first surviving term in the expansion, (3.3) reduces to

(3.4)\begin{equation} \frac{\Delta U}{U_c} \sim{-}\frac{D}{ We } t_n^2,\end{equation}

where $D=4.5[1-\exp (-1)]$. The vorticity produced, and thus the jump in velocity, are quadratic in time and proportional to the surface tension. In this limit, the interface can be closely approximated by a straight line rotating at a constant rate. We may then write the relevant dispersion relation (Chandrasekhar Reference Chandrasekhar1981), in the case where there is no density contrast ($At=0$), in non-dimensional form as

(3.5)\begin{equation} \frac{\omega}{k {U}_c} ={-}\frac{\Delta U} {2 U_c} \pm \left[\frac{k r_c}{2 We} - \frac{1}{4}\left\{\frac{\Delta U}{U_c}\right\}^2 \right]^{1/2}, \end{equation}

where the real and imaginary parts of $\omega$, respectively, give the frequency and growth rate of a perturbation of wavenumber $k$. Substituting (3.4) in (3.5) we get

(3.6)\begin{equation} \frac{\omega}{k {U}_c} = \frac{D} {2 We} t_n^2 \pm \left[\frac{k r_c}{2 We} - \frac{1}{4}\left\{\frac{D}{We}\right\}^2 t_n^4 \right]^{1/2}.\end{equation}

Note that while we derive this expression under the inviscid approximation, since there are no singular effects of viscosity in this problem, it is applicable in viscous flow at high Reynolds numbers. Since viscosity multiplies the highest derivative in the Navier–Stokes equation, in certain high Reynolds number flows this term can provide a singular perturbation to the problem, and qualitatively change the nature of the solution. These effects, however, typically occur at walls. In the absence of walls, as in the present case, we may expect that high Reynolds number viscous results will agree well with inviscid results. The first term within the square bracket stands for the standard stabilising action of surface tension, increasing with wavenumber. The second term, on the other hand, indicates the destabilising effect of surface tension. As seen from (3.6), the destabilising term is quadratic in the surface tension $\propto 1/We^2$, while the stabilising term is only linear in this quantity.

Moreover, it is evident from (3.6) that the destabilising action of surface tension must win over its stabilising action at some time for any Weber number. In other words, the interface within the core becomes KH unstable when

(3.7)\begin{equation} t_n > \left(\frac{2 k r_c We }{D^2}\right)^{1/4}.\end{equation}

In figure 2 we plot the growth rate of perturbations given by the imaginary part of $\omega$ in (3.6) as a function of the wavenumber $k r_c$. It is seen that the higher the surface tension, the faster the instability grows. But at any given time, for any finite $We$, there is a cutoff wavenumber beyond which flow is stable. Equation (3.6) shows that the cutoff wavenumber depends only on $We$ and $t_n$. We observe that at $t_n =1$, the lowest surface tension ($We=10\ 040$) is stable for all wavenumbers. However, with time, the vorticity on the interface increases in all cases, destabilising longer wavelengths even at higher $We$, as seen in figure 2(b) at $t_n=5$.

Figure 2.

Non-dimensional growth rate of perturbations ( $r_c{\rm Im}(\omega )/U_c$) on the interface vs the wavenumber $kr_c$; (a) $t_n = 1$ and (b) $t_n = 5$. Here, $r_c = 0.1$ and $At = 0$. As $We$ decreases, the growth rate increases, showing the destabilising effect of surface tension.

3.2. Vorticity generation and the instability in the region outside the core

Well outside the vortex core (refer Appendix A for details), $q \gg 1$, using (3.3) we have

(3.8)\begin{equation} \frac{\Delta U}{U_c} \sim{-}\frac{40 [1-\exp({-}1)] t_n^2 }{q^2 We } \pm \frac{t_n At}{[1-\exp({-}1)]^2 q^{3/2}}. \end{equation}

The first term on the right-hand side of (3.8) corresponds to the contribution of surface tension to the velocity jump, which is negative. The second term on the right-hand side is the contribution due to the density contrast, it is negative when the inner fluid is heavier than the outer fluid and positive otherwise. The surface tension contribution to the vorticity on the interface decrease as $r^{-4}$ and the buoyancy contribution to the vorticity decreases as $r^{-3}$. To estimate the instability, by following a procedure analogous to Dixit & Govindarajan (Reference Dixit and Govindarajan2010), we approximate the interface as a circle of radius $r$ and perturb it at an azimuthal wavenumber $m$, to get (refer to Appendix B for details)

(3.9)\begin{align} \frac{\omega r_c}{mU_c} &= (q)^{{-}1}Q+ \left( \frac{1-At}{2}\right)\left( \frac{\Delta U}{U_c}\right)q^{{-}1/2} \nonumber\\ &\quad \pm q^{{-}1/2} \left[-\left\{\frac{(1-At) }{2m } + \left( \frac{1-At^2}{4}\right)\right\} \left(\frac{\Delta U}{U_c}\right)^2\right. \nonumber\\ &\quad - \left.\frac{(1-At) Q}{m q^{1/2}} \left(\frac{\Delta U}{U_c}\right) -\frac{At Q^2}{m q} + \frac{m}{2 q^{1/2}We} \right]^{1/2}, \end{align}

where $Q = [1-\exp (-q)]/[1-\exp (-1)]$. On examining this expression in the light of (3.8), we see that, at $At=0$, the largest destabilising term is $O(q^{-2})$ smaller than the stabilising term, so instability is not expected. When $At>0$, however, the spiralling interface can be unstable, as found by Dixit & Govindarajan (Reference Dixit and Govindarajan2010), but now surface tension stabilises the flow at high azimuthal wavenumber. Thus, the instability within a vortex core is driven by surface tension and that outside by density differences.

4.1. Simulation details

We conduct DNS using an open-source volume of fluid code Basilisk to solve (2.1a,b) and (2.2). The surface tension force given by (2.3) is modelled as a continuum surface force in the cells containing the interface. Curvature is calculated using a heights function method. The details of the formulation are given in Popinet (Reference Popinet2018). We place a Lamb–Oseen vortex at the interface in the centre of the domain, as shown in figure 1(a,b) and allow it to evolve with time. Since our computational domain, of length $L=5{\rm \pi} r_c$, is much larger than our vortex, the far-field boundary conditions do not affect the results. We use free-slip conditions on all sides, and have checked that periodic boundary conditions give practically indistinguishable answers. Note that our simulations include centrifugal effects of density contrast when $At \ne 0$. All simulations shown here are viscous, i.e. $\mu \ne 0$. We discretise the domain with $2048^2$ collocation points, and vary the Weber number $(We): 10\ 040, 100.4, 33.3$ and $12$. The units of $\varGamma$, $r_c$, $\bar {\rho }$, $\mu$ and $\sigma$ are arbitrary in the simulations, and care is taken to define non-dimensional numbers as done in the theory for purposes of comparison. We restrict our simulations to times over which the sum of the interfacial and kinetic energies remains constant.

4.2. Vorticity generation and the resulting instabilities

Up to a time $t_n \sim We^{1/2}$ we expect the vorticity generated on the interface in the numerical simulations to closely follow (3.3). This is shown to be the case in figure 3(a), where $t_n=1$ and $t_n = 2$ for $We=12$. To calculate the velocity jump across the interface at the point $r_s$, we numerically integrate the vorticity along the normal. We have checked that $\Delta U$ is insensitive to grid resolution. As predicted, $\Delta U \sim t_n^2$. Interestingly, the velocity jump changes sign, going from negative to positive below $r_c$. There is one more sign change at large $r$ (not shown here), and the velocity jump in the distant spiral arms is negative again, although very small. There is a small deviation of the simulated data points from (3.3). This is to be expected since surface tension and viscous effects are neglected in the derivation of (2.5). Viscous effects in the simulation will bring about a small reduction in the velocity jump as compared with theory. Secondly, from (2.7) we expect surface tension to affect the basic interface shape significantly after a time ${\sim}We^{1/2}$. Figures 3(b) and 3(c) show the perturbation vorticity at two times, calculated by subtracting the initial vorticity from the instantaneous field. In an inviscid flow, by Kelvin's circulation theorem, the vorticity created would move with the interface for eternity. Viscous flow is quite different in this aspect and this difference is an important feature of viscous multiphase flow. In our case, because the Reynolds number is high, the vorticity stays attached to the interface for quite a while. However, later, we observe, especially in the high surface tension cases such as in figure 3(c), a peeling off of vortices from the interface. This is because the new vorticity tends to roll up, as is normal for any shear layer. To stay attached, the interface would need to roll up along with the vorticity, causing a large increase in its length. This is resisted by the interface, which becomes shorter locally and the vorticity peels off. It is seen that at the larger time, instability has set in, which is consistent with the prediction of instability when $t_n > 2.08$ for $kr_c = 2 {\rm \pi}$ by (3.4). The instability becomes visible in the simulations to the naked eye at around $t_n=5$ (see movies in the supplementary material available at https://doi.org/10.1017/jfm.2022.97). For the case where surface tension ($We = 10040$) is much lower, the time above which instability can occur, from (3.4), is $11.19$ and for $We=12$ it is $2.08$, but because the growth rate is minuscule as predicted by equation (3.5) and shown in figure 2, an instability does not become visible during the time of our simulation. Note that in simulations $\mu \ne 0$ whereas we do the stability analysis for the inviscid case. The viscosity in the numerical simulations reduces instability growth rate somewhat.

Figure 3.

Vorticity generation with $We=12, At = 0, {Re} = 10114$. (a) Velocity jump measured across the interface, as a function of interface location within the vortex core. Symbols: DNS, line: (3.3). The perturbation vorticity from the DNS is shown for time $t_n=1$ in (b) and for $t_n = 15$ in (c). The instability is well developed by $t_n=15.$ All lengths are scaled by $r_c$.

The perturbation vorticity for $We=10\ 040$ and $12$ is plotted in figure 4. In figure 4(a), where $At = 0$, we find, in accordance with our expectations from the discussion in §§ 3 and 4.2, that a minuscule amount of vorticity is generated on the interface but no instability happens. Further, at $r \sim r_c$, this vorticity is positive. However, for $We=12$, with $At = 0$, strong negative vorticity is generated within the core, and significant positive vorticity occurs in the vicinity of $r \sim r_c$. As time progresses, vorticity on the interface within the core causes a KH instability, and into a breakdown into small patches of negative vorticity within the core. At even later times, as is evident from figure 4(b) positive and negative perturbation vorticities are interspersed, and dynamics in the core region appears to be chaotic. This state causes the destruction of the initial Lamb–Oseen vortex. We notice a small level of asymmetry in the vorticity contours shown in figure 4(b), and believe this to be a numerical artefact. We now examine, in figure 4(c), what happens when there is a density difference in the two fluids. It is clear that when the surface tension is low and the density contrast is high, there is a generation of alternating spirals of positive and negative vorticity in the region outside the core. The blue arms of the spiral are unstable to centrifugal Rayleigh–Taylor (CRT) modes, while the red arms are stable (Dixit & Govindarajan Reference Dixit and Govindarajan2010). At later time we see a nonlinear breakdown of the flow, as in the figure. When surface tension and density contrast are both high, just outside the vortex core, the CRT mode is stabilised by surface tension but displayed at long distances from the core, as is to be expected from our simplified theory. And the vortex gets disrupted due to the rapidly growing KH instability. We observe from figure 4(d) that the heavy fluid moves outward from the core of the vortex due to centrifugal effects. Note that the different panels of figure 4 are shown at different times. The times were chosen so as to make the instability and its nonlinear development evident.

Figure 4.

Perturbation vorticity for $We=10040$ and $12$ with and without density contrast for ${Re} = 10\ 114$: (a) $We=10\ 040$, $At = 0$ at $t_n=25$. Perturbation vorticity of small magnitude, in the neighbourhood of $r \sim r_c$, is generated. Vorticity generated within the core is small; (b) $We=12$, $At = 0$ at $t_n = 25$. The negative vorticity is mainly from the interior of the vortex; (c) $We=10\ 040$, $At = 0.05$ at $t_n = 40$. The alternating signs of density jump on the spiral arms cause alternating positive and negative vorticity; (d) $We=12$, $At = 0.5$ at $t_n = 14$. Due to the effect of either surface tension or density contrast or both, the original flow structure is disrupted in (b),(c) and (d). The black circle represents the vortex core.

4.3. Evolution of the interface

The interfaces for the cases in figure 4 are shown in figure 5. At low surface tension and zero density difference, the spiralling interface of figure 5(a) is indistinguishable from that at zero surface tension seen earlier in figure 1(c). This is because the time $T_s$, given by (2.7), beyond which surface tension effects will be significant, is larger than our simulation time. At high surface tension, we observe from figures 4(b) and 5(b) that the interface shape is completely different from the vorticity distribution. Although vorticity is generated only at the interface, it rolls up into small-scale structures independent of the interface due to the KH instability, and with time spreads through the vortex core. The interface meanwhile adopts a relatively short and straight shape in the central region, in deference to the high surface tension. In contrast, due to the low surface tension in figure 5(c), the interface closely mimics vorticity contours of figure 4(c). Here, the CRT and KH instability are on display due to the density contrast. Far away from the core, alternating positive and negative vorticity appears on the spiral arms, indicating, from (3.8), that in this case density contrast is dominant over surface tension. Meanwhile, in the core, the vortices peel off from the interface, as a consequence of the KH instability. An important consequence of centrifugal forces is seen in figure 5(d). We started with a symmetry in the flow when rotated by the angle ${\rm \pi}$. In this image we see that this symmetry is completely broken, with the light fluid occupying the initial vortex core, and the heavy fluid having been expelled by centrifugal forces. We thus have a competition between the response to the vortex, which increases the length of the interface, and surface tension, which acts to reduce it. An independent estimate of the interface length can be obtained from energy balance. As there is no external forcing and small viscous dissipation here, the total kinetic energy is given by

(4.1)\begin{equation} \partial_t E= \langle \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{F}_{\sigma} \rangle +\epsilon_{\mu},\end{equation}

where the angled brackets refer to an average per unit area taken over the whole domain, $E =\langle \rho u^2/2 \rangle$ and $\epsilon _{\mu }$ is the viscous dissipation rate. We also have $\langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {F}_{\sigma } \rangle = - ({1}/{A})\partial _t \int {\sigma dl}$, where $dl$ is an interfacial line element and $A$ is the total area of the domain (Joseph Reference Joseph1976). At high Reynolds number and significant levels of surface tension, we may neglect the viscous dissipation, to obtain

(4.2)\begin{equation} \partial_t \left(E + \frac{\sigma S}{A} \right) = 0, \quad {\rm i.e.} \ - \frac{E-E_0}{ \sigma} = \frac{S}{A}, \end{equation}

where $S = \int \, {\rm d} l$, and $E_0$ the initial kinetic energy. We obtain the length $S$ of the interface numerically as a function of time from all four simulations, by tracking jumps in $c$ from $0$ to $1$. Obtaining the interface length correctly is a hard test for simulations. We perform this test by defining a simulation Weber number $We_c \simeq \rho U_c^2 r_c S/[A (E_0-E)]$, where the kinetic energy is obtained from the simulations. This quantity is plotted against time in figure 6(a), and is reassuringly close to the prescribed Weber number in all cases except for the highest Weber number. In this case, we do not expect an agreement, since surface tension effects are small, becoming comparable to viscous effects, which we have neglected in estimating $We_c$. We have checked in inviscid simulations (not shown) of all four cases that $We_c$ at all times is in excellent agreement with $We$. This confirms that the discrepancy at high $We$ is due to viscous effects. The test verifies the reliability of the numerics. The interface length is plotted in figure 6(b). At the highest Weber number is the practically unbridled lengthening of interface by its winding up into a spiral. As the surface tension is made higher, it competes better against the dynamics as dictated by the central vortex and suppresses the increase in interface length.

Figure 5.

Regime occupied by the two fluids, one shown in black and the other in white. (a) A spiral interface is seen at $We=10\ 040$ with $At = 0$ at $t_n=25$; (b) $We=12$, $At=0$, $t_n=25$. The fine structure is erased, and a straight interface is seen in the central portion; (c) $We=10\ 040$ with $At = 0.05$ at $t_n=40$. The CRT instability in the nonlinear regime is visible; (d) $We=12$, $At = 0.5$ at $t_n = 14$. The interface is displaced from the core of the vortex when both surface tension and density contrast are present. Here, ${Re} = 10\ 114$.

Figure 6.

(a) Weber number estimated by the energy balance (symbols) compared against the actual Weber number (lines). (b) Increase in interface length (scaled by total area) with time.

4.4. Two-dimensional turbulence

To evaluate whether the mechanism we propose, of destabilisation of a flow by surface tension, is of relevance in a general 2-D turbulent flow of two immiscible fluids, we conduct DNS – at a Reynolds number (${Re}_T = \bar {\rho }u_{rms} L/\mu$, where $u_{rms}$ is the root mean square velocity at $t=t_i$) of $5.3 \times 10^5$ unless otherwise specified – and $At=0$ on a doubly periodic box of length $L=2{\rm \pi}$ and discretise it with $2048^2$ collocation points. We generate the initial vorticity field by first performing a single-fluid simulation. Here, we start with white noise for vorticity, and run the simulation until we attain a $k^{-3}$ scaling in the energy spectrum. Once this energy spectrum has been attained, we set the time to zero, and restart the simulation, this time with two fluids, separated by a flower-shaped interface. Such an interface shape provides some irregularity, and is obtained by superimposing a sinusoidal perturbation on a circle. The volume fraction of the minority phase is $0.225$, and simulations are conducted at four different $We_T$ ($= \rho u^2_{rms} L/\sigma$): $\infty$, $5 \times 10^5$, $5 \times 10^4$ and $5 \times 10^3$ in runs $TR1$, $TR2$, $TR3$ and $TR4$, respectively. The initial conditions for the two-fluid simulation are shown in figure 7(a,d). The relevance to 2-D turbulence of the model flow we studied in previous sections becomes apparent when we examine figure 7(d) in some detail. We started out our single-fluid simulation with white noise, and, as expected, the vorticity organised itself primarily into concentrated patches. While some of these patches, due to the very nature of 2-D turbulence, have deviated from a circular shape, most of the vortex patches are fairly circular. Moreover, as shown in Appendix C, the vorticity profiles within the patches are well described by the Lamb–Oseen. This is not too surprising, given that Lamb–Oseen vortices are the natural steady state for a patch of vorticity (in the absence of forcing) in viscous flow, but it is reassuring to check this. Thus, studying an interface near a single Lamb–Oseen vortex is a good model for studying two-fluid turbulence. It remains to be argued why we chose a straight interface going through the centre of the vortex. This was for simplicity of algebra, and it provides the physics of the instability. It can be checked that a displaced interface, or one that is initially of a different shape, will be wound into a spiral-like shape in the vicinity of the vortex.

Figure 7.

Density (a–c) and vorticity (d–f) contours for viscous simulations at $t_i/\tau _e=9.4$ with ${Re} = 5.3 \times 10^5$. (a,d) Initial profiles, (b,e) profiles at $t^{*}/\tau _e =14.6$ for $We_T=5\times 10^5$, (c, f) profiles at $t^{*}/\tau _e =14.6$ for $We_T=5\times 10^3$ where $\tau _e = L/u_{rms} = 0.745$ is the eddy turn over time. Small-scale vortical structures occur in greater abundance as the surface tension increases. The colour scale in (d–f) has been adjusted for a better view of small-scale vorticity.

The flow consists of vortices of several scales (with vortex core radii ranging from approximately $0.001L$ to $0.03L$), encompassing a range of Weber numbers based on each, with the largest being an order of magnitude smaller than the smallest $We_T$. Moreover, the interface is not particularly designed to pass through the centre of any vortex. So, we do not expect an exact agreement with theory, but do expect vorticity generation due to surface tension, and instabilities resulting in small-scale vortices. The regions occupied by the two fluids at $t/\tau _e=14.6$ for high and low Weber number are shown in figures 7(b) and 7(c) respectively, and the corresponding vorticity distributions are shown in figures 7(e) and 7( f) respectively. At $We_T=5 \times 10^5$, the two fluids, which displayed many spiralling interfaces at short times (see movies in the supplementary material), are mingled intricately by this time, with elongated fine structures of the inner fluid. This picture is already qualitatively different from flow at zero surface tension, in that we see vorticity generation due to surface tension along the interfaces, as predicted. But the instability, and the vorticity being peeled off from the interfaces is not yet visible. The picture is starkly different at high surface tension. There is an explosion of small-scale vorticity everywhere in the flow, and the initial vortices have been disrupted completely. Because of these small-scale vortices, there is an enhancement of energy at large wavenumbers with increase in surface tension.

The total kinetic energy in the system in shown as a function of time in figure 8(a). The two-phase simulation is started at $t=t_i$. During the time $t< t_i$, the initial condition is prepared by running the single-phase simulation, and so this portion of the plot is not relevant to our discussion. We see that kinetic energy goes down with time. A very small part of this may be attributed to viscous dissipation whereas, especially at high surface tension, a significant part of the kinetic energy is converted to surface energy, manifested as an increase in the total length of the interface. Turbulent kinetic energy spectra (given by $E(k)=\varSigma | u_k^2|$, where $u_k$ is the velocity in spectral space) are shown in figure 8(b). With increase in surface tension, there is a significant enhancement of kinetic energy at large wavenumbers. The small-scale vorticity generated by the action of surface tension is the reason for this increase. Incidentally, earlier studies (Li & Jaberi Reference Li and Jaberi2009; Trontin et al. Reference Trontin, Vincent, Estivalezes and Caltagirone2010) have obtained a similar increase in kinetic energy with decrease in Weber number, but have not explored the underlying mechanism. The details of the instability might be different in the case of Trontin et al. (Reference Trontin, Vincent, Estivalezes and Caltagirone2010) due to vortex stretching in their case.

Figure 8.

(a) The variation in time of the total kinetic energy ($\rho u^2/2$). At $t_i/\tau _e$, we start our multiphase simulations with the flower-shaped interface and non-zero surface tension (figure 7a,d). (b) Energy spectra ($E(k)$) versus wavenumber ($k$) at $t^*/\tau _e$, which is shown in (a). Note that $E(k)$ is normalised in each case by $E(1)$. The fraction of energy in small scales is higher at higher surface tension. The vertical dashed lines correspond to the characteristic length scale $L_{\sigma }$ in each case. This is compared in the inset at different $We$ with the estimated length scale ($L_{est}$). The black dash-dotted line indicates a slope of $-3$.

In figure 8(b), which is shown at $t/t_e=14.6$, we had seen a significant component of energy in the small scales, and the small-scale vortices were visible in figure 7( f). We perform a consistency check now to verify that the instability due to surface tension did indeed set in well before this time. The theory is for a single vortex, so we choose the largest positive vortex (vortex 2 in figure 10b) in Appendix C. We first fit this vortex to a Lamb–Oseen vortex of appropriate circulation and width, yielding $r_c \simeq 0.2$, as shown in figure 11(b) of Appendix C. Using this and the fact that the maximum vorticity is $\varGamma /({\rm \pi} r_c^2)$ for a Lamb–Oseen vortex, we obtain $\varGamma =38.96$ and $U_c=19.6$. Thus, for this single vortex, the effective $We=768.32$ (run $TR4$). Using (3.7), and for a typical wavenumber of $kr_c=2$, we estimate the time ($t_n= t/T_c$) above which the instability should take place as $t_n= 4.4 + t_i/\tau _e$ (we add $t_i/\tau _e$ because the two-phase simulation starts at time $t_i$). Since the inertial time scale $T_c= 0.0065$, we expect the instability should occur after $t/\tau _e = 9.783$. Indeed, in figures 7( f) and 8(b), we observe significant small-scale vorticity at $t/\tau _e = 14.6$, clearly indicating the development of an instability. The interface near smaller vortices will go unstable at earlier times than near the largest vortex. Due to the nonlinear nature of the flow, our theory cannot predict the size distribution of vortices at later times.

Another check we perform is to calculate the characteristic length scale, $L_{\sigma }$, of a blob of the minority fluid, as a function of time, and compare it against the Hinze scale, $L_{est}$, obtained by balancing the inertial forces with the surface tension (Hinze Reference Hinze1955; Perlekar et al. Reference Perlekar, Benzi, Clercx, Nelson and Toschi2014; Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Perlekar Reference Perlekar2019). The two are defined as

(4.3a,b)\begin{equation} L_{\sigma} = 2 {\rm \pi}\frac{\varSigma_{k} C(k)}{\varSigma_{k} kC(k)} \quad {\rm and} \quad L_{est} \sim \sigma^{1/3} \beta^{{-}2/9}, \end{equation}

where the indicator function spectrum (Perlekar, Pal & Pandit Reference Perlekar, Pal and Pandit2017) is $C(k)\equiv \sum |c_k|^2$ ($c_k$ is the Fourier transform of indicator function $c$) and $\beta = \nu \langle |\nabla \varOmega |^2 \rangle$ is the enstrophy dissipation rate. The inset in figure 8(b) shows very good agreement between the estimated and the calculated length scales. We notice in figure 8(b) that the range of small-scale turbulence affected by the surface tension increases with surface tension, and is of the same order of magnitude as the characteristic length scale $L_\sigma$.

To quantify the chaotic mixing process and the effect of surface tension on it, we divide the domain into $n_b=256$ equal area square boxes, $16$ in each dimension, and calculate the fraction of the minority fluid ($a$) occupying each box. The volume fraction ($A$) of the minority phase averaged over the entire computational domain is $0.225$, and this remains invariant in time. The probability density function (p.d.f.) $P(a/A)$ for $TR2$ and $TR4$, at low and high surface tension respectively, are shown in figure 9. It is instructive to compare these with the p.d.f. of $(a/A)$ for a circular patch of minority fluid embedded in a square patch of majority fluid, and occupying the same volume fraction of $0.225$ as in our simulations. When this square patch including both fluids is split into $n_b$ boxes, many boxes do not have minority fluid ($a=0$) and some boxes have $a=1$, i.e. $a/A=1/0.225$, giving peaks at $a/A=0$ and $a/A=4.44$. Any entries between $a/A=0$ and $4.44$ indicate a box occupied in part by both fluids. It is clear that significant mixing has taken place due to turbulence in both $TR2$ and $TR4$. When ideally mixed, every box would have $a=A$, giving a single peak at $1$ in the p.d.f. of $a/A$, and no entries elsewhere. But we do not expect ideal mixing in our simulation, since surface tension would compete with the mixing tendency of turbulence, to minimise interface length. Therefore, it would impart to the minority fluid a tendency to form blobs. Consistent with this expectation, we see in figure 9 weaker mixing at higher surface tension. This is also consistent with the size distribution shown in Trontin et al. (Reference Trontin, Vincent, Estivalezes and Caltagirone2010), where the number of large droplets increase with the decrease in $We$.

Figure 9.

Value of $P(a/A)$ versus $a/A$ for two runs $TR2, We = 5 \times 10^5$ and $TR4, We= 5 \times 10^3$ time averaged from $14.1\tau _e$ to $14.6\tau _e$. Black dashed line corresponds to $a/A=1$. The black solid line with circles represents the p.d.f. for a circular patch of minority fluid of the same volume fraction (=0.225) as that of the minority fluid in the simulations. As the $We$ increases, most boxes have the same volume fraction as that of the total domain.