Two-dimensional Turbulence in Symmetric Binary-Fluid Mixtures: Coarsening Arrest by the Inverse Cascade

We study two-dimensional (2D) binary-fluid turbulence by carrying out an extensive direct numerical simulation (DNS) of the forced, statistically steady turbulence in the coupled Cahn-Hilliard and Navier-Stokes equations. In the absence of any coupling, we choose parameters that lead (a) to spinodal decomposition and domain growth, which is characterized by the spatiotemporal evolution of the Cahn-Hilliard order parameter ϕ, and (b) the formation of an inverse-energy-cascade regime in the energy spectrum E(k), in which energy cascades towards wave numbers k that are smaller than the energy-injection scale kin j in the turbulent fluid. We show that the Cahn-Hilliard-Navier-Stokes coupling leads to an arrest of phase separation at a length scale Lc, which we evaluate from S(k), the spectrum of the fluctuations of ϕ. We demonstrate that (a) Lc ~ LH, the Hinze scale that follows from balancing inertial and interfacial-tension forces, and (b) Lc is independent, within error bars, of the diffusivity D. We elucidate how this coupling modifies E(k) by blocking the inverse energy cascade at a wavenumber kc, which we show is ≃2π/Lc. We compare our work with earlier studies of this problem.

Binary-fluid mixtures (such as oil and water) have played a pivotal role in the development of the understanding of (a) equilibrium critical phenomena at the consolute point, above which the two fluids mix 1-3 , (b) nucleation 4 , and (c) spinodal decomposition, the process by which a binary-fluid mixture, below the consolute point and below the spinodal curve, separates into the two, constituent liquid phases until, in equilibrium, a single interface separates the two coexisting phases (this phase separation is also known as coarsening) 5,6 . In the presence of flows, the demixing because of spinodal decomposition gets arrested and an emulsion is formed. This process, also known as coarsening arrest, is important in several three-dimensional (3D) and two-dimensional (2D) turbulent flows. The former have been studied recently [7][8][9] . Coarsening arrest in a 2D, turbulent, binary-fluid mixture is also of relevance to problems such as the dynamics of oil slicks on the surface of the ocean, whose understanding is of clear socio-economic and scientific relevance [10][11][12][13] . Oceanic flows have been modelled successfully as 2D, turbulent fluids. Such 2D turbulence is fundamentally different from three-dimensional (3D) fluid turbulence as noted in the pioneering studies of Fjørtoft, Kraichnan, Leith, and Batchelor [14][15][16][17][18] . In particular, the fluid-energy spectrum in 2D turbulence shows (a) a forward cascade of enstrophy (or the mean-square vorticity), from the energy-injection wave number k inj to larger wave numbers, and (b) an inverse cascade of energy to wave numbers smaller than k inj . We elucidate the turbulence-induced arrest of phase separation in a 2D, symmetric, binary-fluid mixture.
Coarsening arrest by 2D turbulence has been studied in ref. 19, where it has been shown that, for length scales smaller than the energy-injection scale π =  k 2 / inj i nj , the typical linear size of domains is controlled by the average shear across the domain. However, the nature of coarsening arrest, for scales larger than  inj , i.e., in the inverse-cascade regime, which is relevant for large-scale oceanic flows, still remains elusive. In particular, it is not clear what happens to the inverse energy transfer, in a 2D binary-fluid, turbulent mixture, in which the mean size of domains provides an additional, important length scale. We resolve these two issues in our study. By combining theoretical arguments with extensive direct numerical simulations (DNSs) we show that the Hinze length scale L H (see refs 8,9) provides a natural estimate for the arrest scale; and the inverse flux of energy also stops at a wave-number scale π  L 2 / H . Coarsening arrest has also been studied in simple shear flows (refs 20-25), which yield coarsening arrest with domains elongated in the direction of shear.
Forced, 2D, statistically steady, Navier-Stokes-fluid turbulence displays a forward cascade of enstrophy, from  inj to smaller length scales, and an inverse cascade of energy to length scales smaller than  inj . In the inverse-cascade regime, on which we concentrate here, E(k) ~ k −5/3 (see, e.g., refs 15,18) and the energy flux Π (k) ~ ε ≡ 〈 ε(t)〉 t assumes a constant value. For the Cahn-Hilliard model, if it is not coupled to the Naiver-Stokes equation, , for large times, where the time-dependent length scale  ∼ t t ( ) 1/3 , in the early Lifshitz-Slyozov [26][27][28][29] regime; if the Cahn-Hilliard model is coupled to the Navier-Stokes equation, then, in the absence of forcing,  ∼ t t ( ) , in the viscous-hydrodynamic regime, first discussed by  , and  ∼ t t ( ) 2/3 , in the very-late-stages in the Furukawa 31 and Kendon 32 regimes. For a discussion of these regimes and a detailed exploration of a universal scaling form for  t ( ) in 3D we refer the reader to ref. 33. We now elucidate how these scaling forms for E(k) and S(k, t) are modified when we study forced 2D turbulence, in the inverse-cascade regime in the coupled Cahn-Hilliard-Navier-Stokes equations.

Results
Cahn-Hilliard-Navier-Stokes equations. We model a symmetric binary-fluid mixture by using the incompressible Navier-Stokes equations coupled to the Cahn-Hilliard or Model-H equations 34,35 . We are interested in 2D incompressible fluids, so we use the following stream-function-vorticity formulation [36][37][38] for the momentum equation: Here u(x, t) ≡ (u x , u y ) is the fluid velocity at the point x and time t, ω = ∇ ×û e ( ) z , φ(x, t) is the Cahn-Hilliard order parameter that is positive in one phase and negative in the other, p(x, t) is the pressure, is the free energy, Λ is the mixing energy density, ξ controls the width of the interface between the two phases of the binary-fluid mixture, ν is the kinematic viscosity, the surface tension ξ = Λ s 2 2/3( / ), the mobility of the binary-fluid mixture is M, and f ω is the external driving force. For simplicity, we study mixtures in which M is independent of φ and both components have the same density and viscosity 33 . We use periodic boundary conditions in our square simulation domain, with each side of length L = 2π. To obtain a substantial inverse-cascade regime, we stir the fluid at an intermediate length scale by forcing in Fourier space in a spherical shell with wave-number where the caret indicates a spatial Fourier transform, ensures that there is a constant enstrophy-injection rate. The higher the Reynolds number Re ∝ 1/ν, the more turbulent is the flow; and the higher the Weber number We ∝ 1/σ, the more the fluctuations in the domains (see Table 1 for definitions of Re, We, and other parameters in our study). To elucidate the , N 2 is the number of collocation points, ν is the kinematic viscosity, M is the mobility parameter, ξ sets the scale of the interface width, the surface tension σ physics of coarsening arrest, we conduct direct numerical simulations (DNSs) of Eqs (1) and (2) (see Methods Section for details).
Coarsening Arrest. In Fig. 1 we show pseudo-gray-scale plots of φ, at late times when coarsening arrest has occurred, for four different values of We at Re = 124; we find that the larger the value of We the smaller is the linear size that can be associated with domains; this size is determined by the competition between turbulence-shear and interfacial-tension forces. This qualitative effect has also been observed in earlier studies of 2D and 3D turbulence of symmetric binary-fluid mixtures [19][20][21][39][40][41][42][43][44] . We calculate the coarsening-arrest length scale c k k we now show that L c is determined by the Hinze scale L H , which we obtain, as in Hinze's pioneering study of droplet break-up 9 , by balancing the surface tension with the inertia as follows: We obtain for 2D, binary-fluid turbulence the intuitively appealing result L c ~ L H (for a similar, recent Lattice-Boltzmann study in 3D see ref. 8). In particular, if we determine L c from Eq. (3), with S(k) from our DNS, we obtain the red points in Fig. 2, which is a log-log plot of σL c versus ε inj /σ 4 ; the black line is the Hinze result (4) for L H , with a constant of proportionality that we find is . 1 6 from a fit to our data. We see from Fig. 2 that the Hinze length scale L H gives an excellent approximation to the arrest scale L c over several orders of magnitude on both vertical and horizontal axes. Note that the Hinze estimate also predicts that, for fixed values of ε inj and σ, the coarsening-arrest scale is independent of D; the plot of L c versus D, in the inset of Fig. 2, shows that our data for L c are consistent (within error bars) with this prediction. Figure 1. Pseudo-gray-scale plots of the order parameter field φ, at late times when coarsening arrest has occured, in 2D symmetric-binary-fluid turbulence with Re = 124. Note that the domain size decreases as we increase the Weber number We from the leftmost to the rightmost panel: We = 1.2 · 10 −2 (R3); We = 5.9 · 10 −2 (R4); We = 1.2 · 10 −1 (R5); and We = 5.9 · 10 −1 (R8). (b) Log-log (base 10) plots of the spectrum S(k), of the phase-field φ, versus k; as We increases (i.e., σ decreases) the low-k part of S(k) decreases and S(k) develops a broad and gentle maximum whose peak moves out to large values of k. (c) Plots versus φ, in the vicinity of the maximum at φ + , of the normalized PDFs P(φ)/P m (φ), where P m (φ) is the maximum of P(φ); the peak position φ + → 1 as We increases (see the inset which suggests that φ − ∼ + We 1 2 1 /2 (black line)).
In Fig. 2(b) we show clearly how the arrest of coarsening manifests itself as a suppression of S(k), at small k (large length scales). This suppression increases as We increases (i.e., σ decreases); and S(k) develops a broad and gentle maximum whose peak moves out to large values of k as We grows. These changes in S(k) are associated with We-dependent modifications in the probability distribution function (PDF) P(φ) of the order parameter φ, which is symmetrical about φ = 0 and has two peaks at φ = φ ± , where φ + = − φ − > 0; we display P(φ)/P m (φ) in Fig. 2(c) in the vicinity of the peak at φ + ; as We increases, φ + decreases; here P m (φ) is the maximum value of P(φ). In particular, our DNS suggests that φ − ∼ + We 1 2 1 /2 , for small We. The modification in P(φ) can be understood qualitatively by making the approximation that the effect of the fluid on the equation for φ can be encapsulated into an eddy diffusivity D e 42,45,46 . The eddy-diffusivity-modified Cahn-Hilliard equation are stable to perturbations. In particular, droplets with linear size < (2π/k d ) decay in the presence of coupling with the velocity field; we expect, therefore, that, in the presence of fluid turbulence, the peak of P(φ) broadens and shifts as it does in our DNS. For a quantitative description of this broadening and the shift of the peak, we must, of course, carry out a full DNS of the Cahn-Hilliard-Navier-Stokes equation as we have done here.
Energy spectrum. We have investigated, so far, the effect of fluid turbulence on the phase-field φ and its statistical properties such as those embodied in S(k) and P(φ). We show next how the turbulence of the fluid is modified by φ, which is an active scalar insofar as it affects the velocity field. In the statistically steady state of our driven, dissipative system, the energy injection must be balanced by both viscous dissipation and dissipation that arises because of the interface, i.e., we must have ε inj = ε ν + ε μ . In Fig. 3(a), we show that ε ν decreases and ε μ increases as we increase We, while keeping ε inj constant, because L c diminishes ( Fig. 1) and, therefore, the interfacial length and ε μ increase. This decrease of L c is mirrored strikingly in plots of the fluid-kinetic-energy spectrum E(k) (Fig. 3(b)), which demonstrate that the inverse cascade of energy is effectively blocked at a wavenumber k c , which we determine below, from the energy flux, and which we find is π  L 2 / c , where L c follows from S(k) (see Fig. 2). The value of k c increases with We; and the inverse cascade is completely blocked for the largest We we use, for which  k k c i nj , the forcing scale. To provide clear evidence that the blocking of the energy flux is closely related to the arrest scale, we show in Fig. 3(c) plots of the energ y f lux is the energy transfer and P(k) is the transverse projector with components P ij (k) ≡ δ ij − k i k j /k 2 . We define k c as the wave-number at which Π E (k) comes within 4% of ε inj . We find that the wave-numer corresponding to the arrest scale 2π/L c (marked by vertical lines for each run) is comparable to k c .
In the presence of the standard viscous term ν∇ 2 u and the Ekman drag αu, it is not possible to see a large range of constant energy flux 47,48 . However, it is possible to attain a large constant energy flux range by carrying out DNSs using hyperviscosity and hypoviscosity 47 (see the Methods Section for details). The plot in Fig. 4(left) shows the energy spectrum and the corresponding energy flux obtained [ Fig. 4(right)] from our runs HR1 and HR2. Consistent with the earlier discussion, we find that the coarsening length L c decreases on increasing We. Furthermore, the formation of arrest-scale domains leads to a blockage of the energy cascade; because we use hypoviscosity, we now see clear evidence of a constant energy flux over a decade for the single-phase Navier-Stokes run. For the binary-fluid case, the energy flux remains constant for a shorter range and then decreases to zero around a wave-number π  L 2 / c . Passive advection. It has been suggested 22,45,46 that coarsening arrest can be studied by using a model in which the field φ is advected passively by the fluid velocity. Such a passive-advection model is clearly inadequate because it cannot lead to the phase-field-induced modifications in the statistical properties of the turbulent fluid (see Fig. 3). The passive-advection case is easily studied by turning off the coupling term φ∇ μ in Eq. (2). We then contrast the results for this case with the ones we have presented above. The parameters we use for the passive-advection DNS are N = 1024, Λ = ξ 2 ,ξ = 0.0176; and we carry out runs for D = 5 · 10 −3 , 1 · 10 −2 , 5 · 10 −2 and 5 · 10 −1 . The evolution of the pseudo-grayscale plots of φ with D, in the left panel of Fig. 5, is qualitatively similar to the evolution shown in Fig. 1. There is also a qualitative similarity in the dependence on D of the scaled PDFs P(φ)/P m (φ); we can see this by comparing the passive-advection result, shown in the middle panel of Fig. 5 for positive values of φ in the vicinity of the peak, with its counterpart in Fig. 2(c). However, there is a qualitative difference in the dependence of L c on D: in the passive-advection case we find L c ~ D 0.27 [Fig. 5 (inset)], which is in stark contrast to the essentially D-independent behavior of L c shown in the inset of Fig. 2(c).

Discussion
In conclusion, our extensive study of two-dimensional (2D) binary-fluid turbulence shows how the Cahn-Hilliard-Navier-Stokes coupling leads to an arrest of phase separation at a length scale L c , which follows from S(k). We demonstrate that L c ~ L H , the Hinze scale that we find by balancing inertial and interfacial-tension forces, and that L c is independent, within error bars, of the diffusivity D. We also elucidate how the coupling between the Cahn-Hilliard and Navier-Stokes equations modifies the properties of fluid turbulence in 2D. In particular, we show that there is a blocking of the inverse energy cascade at a wavenumber k c , which we show is π  L 2 / c . Earlier DNSs of turbulence-induced coarsening arrest in binary-fluid phase separation have concentrated on regimes in which there is a forward cascade of energy in 3D (see ref. 8) and a forward cascade of enstrophy in 2D (see ref. 19). Although studies that use a passive-advection model for φ obtain results that are qualitatively similar to those we obtain for S(k) and the spatiotemporal evolution of φ, they cannot capture the phase-field-induced modification of the statistical properties of fluid turbulence and the correct dependence of L c on D. We find our results to be in qualitative agreement with the earlier studies on the advection of binary-fluid mixtures with synthetic chaotic flows 45,46 ; of course, such studies cannot address the effect of the phase field on the turbulence in the binary fluid.
Some groups have also studied the statistical properties of turbulent, symmetric, binary-fluid mixtures above the consolute point, where the two fluids mix even in the absence of turbulence 40,49,50 . In these studies, there is, of course, neither coarsening nor coarsening arrest.
We hope our study will lead to new experimental studies of turbulence in binary-fluid mixtures, especially in 2D [51][52][53][54] , to test the specific predictions we make for L c and the blocking of the inverse cascade of energy. We also plot the single-phase Navier-Stokes energy spectrum and the energy flux for reference. On increasing We, small domains are formed and these lead to a truncation of energy flux at a wave-number around π  L 2 / c (marked by vertical lines).   , N 2 is the number of collocation points, ν u is the hyperviscosity, ν i is the hypoviscosity, M is the mobility, parameter, ξ sets the scale of the interface width, the surface tension σ = Λξ 2 2/3( / ), 〈 f ω ω〉 is the enstrophy-injection rate, which is related to the energy-injection rate