# 3D Spinodal Decomposition in the Inertial Regime

\[

## Abstract

We simulate late-stage coarsening of a 3D symmetric binary fluid using a lattice Boltzmann method. With reduced lengths and times, $`l`$ and $`t`$ respectively (with scales set by viscosity, density and surface tension) our data sets cover $`1\stackrel{<}{}l\stackrel{<}{}10^5`$, $`10\stackrel{<}{}t\stackrel{<}{}10^8`$. We achieve Reynolds numbers approaching $`350`$. At Re $`\stackrel{>}{}100`$ we find clear evidence of Furukawa’s inertial scaling ($`lt^{2/3}`$), although the crossover from the viscous regime ($`lt`$) is very broad. Though it cannot be ruled out, we find no indication that Re is self-limiting ($`lt^{1/2}`$) as proposed by M. Grant and K. R. Elder$`[`$Phys. Rev. Lett. 82, 14 (1999)$`]`$.  

(Received 25 February 1999)

PACS numbers: 64.75+g, 07.05.Tp, 82.20.Wt

\] When an incompressible binary fluid mixture is quenched far below its spinodal temperature, it will phase separate into domains of different composition. Here we consider only fully symmetric 50/50 mixtures in three dimensions, for which these domains will, at late times, form a bicontinuous structure, with sharp, well-developed interfaces. The late-time evolution of this structure remains incompletely understood despite theoretical , experimental  and simulation  work.

As emphasized by Siggia  and Furukawa , the physics of spinodal decomposition involves capillary forces, viscous dissipation, and fluid inertia. Thus, assuming that no other physics enters, the control parameters are interfacial tension $`\sigma `$, fluid mass density $`\rho `$, and shear viscosity $`\eta `$. From these can be constructed only one length, $`L_0=\eta ^2/\rho \sigma `$ and one time $`T_0=\eta ^3/\rho \sigma ^2`$. We define the lengthscale $`L(T)`$ of the domain structure at time $`T`$ via the structure factor $`S(k)`$ as $`L=2\pi S(k)𝑑k/kS(k)𝑑k`$. The exclusion of other physics in late stage growth then leads us to the dynamical scaling hypothesis : $`l=l(t)`$, where we use reduced time and length variables, $`lL/L_0`$ and $`t(TT_{int})/T_0`$. Since dynamical scaling should hold only after interfaces have become sharp, and transport by molecular diffusion ignorable, we have allowed for a nonuniversal offset $`T_{int}`$; thereafter the scaling function $`l(t)`$ should approach a universal form, the same for all (fully symmetric, deep-quenched, incompressible) binary fluid mixtures.

It was argued further by Furukawa  that, for small enough $`t`$, fluid inertia is negligible compared to viscosity, whereas for large enough $`t`$ the reverse is true. Dimensional analysis then requires the following asymptotes:

$$lbt\text{}tt^{}$$
(1)   
$$lct^{2/3}\text{}tt^{}$$
(2)   
where, if dynamical scaling holds, amplitudes $`b,c`$ (and the crossover time $`t^{}`$) are universal. The Reynolds number, conventionally defined as, $`\mathrm{R}e=\rho /\eta LdL/dT=l\dot{l}`$, becomes indefinitely large in the inertial regime, Eq. (2).

In a recent paper, Grant and Elder have argued  that the Reynolds number cannot, in fact, continue to grow indefinitely. If so, Eq. (2) is not truly the large $`t`$ asymptote, which must instead have $`lt^\alpha `$ with $`\alpha \frac{1}{2}`$. Grant and Elder argue that at large enough Re, turbulent remixing of the interface will limit the coarsening rate , so that Re stays bounded. A saturating Re (which they estimate as Re $`10100`$) would require any $`t^{2/3}`$ regime to eventually cross over to a limiting $`t^{1/2}`$ law. But if a single length scale $`lt^{1/2}`$ is involved, a saturating Re implies balance between viscous and inertial terms ($`t^{3/2}`$), while the driving term (interfacial tension) remains much larger than either ($`t^1`$). This suggests a failure of scaling altogether, with at least two length scales relevant at late times. In any case, the arguments of Grant and Elder are far from rigorous; the coarsening interfaces could, remain one step ahead of the remixing despite an ever-increasing Re which, if applied to a static interfacial structure, would break it up. Thus Eq. (2) cannot yet be ruled out as a limiting law.

In what follows we present the first large-scale simulations of 3D spinodal decomposition to unambiguously attain a regime in which inertial forces dominate over viscous ones. We find direct evidence for Furukawa’s $`lt^{2/3}`$ scaling, Eq. (2). Although a further crossover to a regime of saturating Re cannot be ruled out, we find no evidence for this up to Re $`350`$. Our work, which is of unprecedented scope, also probes the viscous scaling regime $`[`$Eq. (1)$`]`$, and the nature of the crossover between this and Eq. (2). Full details of our results  and of the simulation algorithm  will be published elsewhere.

Our simulations use a lattice Boltzmann (LB) method  with the following model free energy:

$$F=𝑑𝐫\left\{\frac{A}{2}\varphi ^2+\frac{B}{4}\varphi ^4+\stackrel{~}{\rho }\mathrm{ln}\stackrel{~}{\rho }+\frac{\kappa }{2}|\varphi |^2\right\},$$
(3)   
in which $`A`$, $`B`$ and $`\kappa `$ are parameters that determine quench-depth ($`A/B1`$ for a deep quench) and interfacial tension ($`\sigma =\sqrt{8\kappa A^3/9B^2}`$); $`\varphi `$ is the usual order parameter (the normalized difference in number density of the two fluid species); $`\stackrel{~}{\rho }`$ is the total fluid density, which remains (virtually) constant throughout .

The simulation code follows closely that of  (for details see ) and uses a cubic lattice with nearest and next-nearest neighbor interactions (D3Q15). It was run on Cray T3D and Hitachi SR-2201 parallel machines with system sizes up to $`256^3`$. The LB method allows the user to choose $`\eta ,\sigma ,\rho `$ (we set $`\rho =1`$ without loss of generality), along with the order-parameter mobility $`M`$ defined via $`\dot{\varphi }=M(\delta F/\delta \varphi )`$. Although it plays no role in the arguments leading to Eqs. (1) and (2), $`M`$ must be chosen with some care to ensure that at late times (a) the algorithm remains stable, (b) the local interfacial profiles remain close to equilibrium (so that $`\sigma `$ is well-defined), and (c) the direct contribution of diffusion to coarsening is negligibly small. Table I shows the parameters used for our eight $`256^3`$ runs.

In all runs, the interface width is $`\xi 5\sqrt{\kappa /2A}3`$ in lattice units. This was found  to be the minimum acceptable to obtain an accurately isotropic surface tension. To minimise diffusive effects, data for which the diffusive contribution to the growth rate was greater than 2% was discarded ; this corresponded to a minimum value of $`L`$ of $`15<L_{min}<24`$, depending on the run parameters. The large size of our runs allowed a ruthless attitude to finite size effects: we use no data with $`L>\mathrm{\Lambda }/4`$, with $`\mathrm{\Lambda }`$ the linear system size. In our $`256^3`$ runs, these filters mean that the good data from any single run lies within $`20\stackrel{<}{}L64`$, a comparable range to previous studies . Datasets of high and low $`L_0`$ are well fit respectively by $`\alpha =1`$ and $`\alpha =2/3`$ (see Fig.1).

However, as emphasized by Jury et al. , meaningful tests of scaling are best made not by looking at single data sets but by combining those of different parameter values. To this end, the good data from each run were fit to $`L=B(TT_{int})^\alpha `$, so as to extract an intercept $`T_{int}`$; we then transformed the data to reduced physical units $`l`$ and $`t`$ defined above. The exponent $`\alpha `$ was first allowed to float freely; this gave reproducible values at large $`l`$ and $`t`$ (e.g., $`\alpha =0.69`$ and $`0.67`$ for the last two data sets in Table 1), but more scattered ones at small $`l`$ and $`t`$ ($`\alpha =0.88`$, $`0.86`$ and $`1.16`$ for the first three data sets). In the latter region the floating fit is relatively poorly conditioned; it also gives large relative errors in $`T_{int}`$ (see Fig.1). In contrast, fits to $`\alpha =1`$ for these three data sets gave much better data collapse with consistent values of $`b`$ ($`b=0.073`$, $`0.073`$ and $`0.072\pm 0.01`$). Thus we are confident of $`\alpha =1`$ in this region. For the remaining data sets we estimate errors in individual exponent values at around 10% and in reduced time $`t`$ around 3% to 10%. Figure 2 shows all our data sets on a single plot using reduced variables $`l`$ and $`t`$. Such a plot is necessarily log-log, since our data sets span seven decades in $`t`$ and five in $`l`$, a range which exceeds all previous studies combined.

These LB results are fully consistent with the existence of a single underlying scaling curve $`l=l(t)`$, in which viscous ($`l=bt`$) and inertial ($`l=ct^{2/3}`$) asymptotes are connected by a long crossover whose breadth justifies our use of a single floating exponent $`\alpha `$ in the fits used above to extract $`T_{int}`$ for each run. Although we cannot rule it out for still larger times $`t`$, we see no evidence for a further crossover to a regime with asymptotic exponent $`\alpha 1/2`$ as demanded by Grant and Elder .

Before considering our results in more detail, we discuss their relation to others previously published. We restrict attention to those 3D data sets for which reliable estimates of $`L_0`$ and $`T_0`$ exist . Datasets of Laradji et al.  and of Bastea and Lebowitz  are shown on Fig.2 (fitted to $`\alpha =1`$ ). These lie in an $`l,t`$ range ($`1\stackrel{<}{}l\stackrel{<}{}20`$) in which our own data shows viscous (linear) scaling $`[`$Eq. (1)$`]`$; both data sets were claimed to confirm the linear law by their authors, but with differing values of $`b=0.13`$, 0.3. Our own $`b`$ values are lower than either (see above and Fig.2). As noted above, we took special care to ensure that the diffusive contribution to coarsening was small; we have found that, for matching $`L_0,T_0`$ values, LB data sets similar to those of Refs. can be generated using too large a mobility $`M`$. We hypothesize therefore that both data sets have strong residual diffusion, leading to an overestimate of $`b`$. Likewise the data of Appert et al. , which lies in the crossover regime of our scaling plot, asymptotes to our data from above; this suggests that their fitted exponent $`\alpha 2/3`$ is too low because of diffusion.

A different explanation, based on a possible nonuniversality of the physics of topological reconnection of domains, was suggested by Jury et al. , whose dissipative particle dynamics (DPD) results also appear in Fig.2 (inset) . These authors found that each data set was well fit by a linear scaling, Eq. (1), but with a systematic increase of the $`b`$ coefficient upon moving from upper right to lower left in the scaling plot. Their alternative suggestion was that their own data, and that of Refs., were part of an extremely broad crossover region, $`1\stackrel{<}{}t\stackrel{<}{}10^4`$ in reduced time. Our LB data support the idea of a broad crossover, but instead places it at $`10^2\stackrel{<}{}t\stackrel{<}{}10^6`$. Note that, unlike those of Refs., all the data sets of Jury et al. do lie very close to our own (Fig.2 inset). Since the two simulation methods are entirely different, this lends support to the idea of a universal scaling, although the fact that each DPD run is best fit by a locally linear growth law does not . The latter could be partly due to finite size effects; to obtain enough data, Jury et al. included results up to $`L=\mathrm{\Lambda }/2`$, whereas we reject all data with $`L>\mathrm{\Lambda }/4`$.

The arguments of Ref.  involve the intrusion of a second length scale, alongside $`L_0`$, which in the LB context is the interfacial width $`\xi `$ (or more generally, a molecular scale). The ratio $`h=\xi /L_0`$ for real fluids is in the range $`0.05`$ (water) to $`10^7`$ (glycerol). In simulations, $`\xi `$ cannot be smaller than the lattice spacing, and the inertial region is achieved by setting $`L_01`$, so $`h1`$. In this sense our interface is “unnaturally thick”: simulation runs that enter the inertial regime do so directly from a diffusive one, without an intervening viscous regime. However, this should not matter if $`l(t)`$ follows a universal curve, as our results (in contrast to Ref.), in fact suggest. But the microscopic length still plays an interesting role, as follows. As a fluid neck stretches thinner and thinner before breaking, it shrinks laterally to the scale $`\xi `$; diffusion then takes over to finish the job of reconnection. So, although our work involves length scales where the direct contribution of diffusion to domain growth is negligible, we must ensure that it is handled correctly at smaller scales. This factor limits the accessible range of $`l`$ and $`t`$, not only at the lower  but also at the upper end .

The breadth of the viscous-inertial crossover is somewhat less extreme when expressed in terms of Re (see above); our data span $`0.1\stackrel{<}{}`$Re$`\stackrel{<}{}350`$ and the crossover region is roughly $`1\stackrel{<}{}`$Re$`\stackrel{<}{}50`$. Re values (at $`L50`$) for each run are shown in Fig.3 against reduced time $`t`$. Data are consistent with Re $`t^{1/3}`$ as predicted from Eq. (2). Note that, in simulating high Re flows, one should strive to ensure that the dissipation scale  (defined as $`\lambda _d=(\eta ^3/ϵ\rho ^3)^{1/4}`$, with $`ϵ`$ the energy dissipation per unit volume) always remains larger than the lattice spacing. This ensures that any turbulent cascade (whose shortest scale is $`\lambda _d`$) remains fully resolved by the grid. Equating dissipation with the loss of interfacial energy, one has $`ϵd(\sigma /L)/dT`$ and so, in reduced units, $`\lambda _d(l^2/\dot{l})^{1/4}`$. Comparable $`ϵ`$ values are found directly from our simulated velocity data; and $`\lambda _d`$ remains larger than the grid size for all our runs .

A decisive check that we really are simulating a regime where inertial forces dominate over viscous ones, is based directly on the velocity fields found in our simulations . From these we calculated rms values of the individual terms in the Navier-Stokes equation ($`\rho =1`$), $`\left(𝐯/t+𝐯𝐯\right)=\eta ^2𝐯𝐏.`$ Here $`𝐏`$, the pressure tensor, contains the driving terms arising from interfacial tension. Ratios $`R_1=𝐯/t_{rms}/\eta ^2𝐯_{rms}`$ and $`R_2=𝐯𝐯_{rms}/\eta ^2𝐯_{rms}`$, were then computed; these can be seen in Fig.3.

The ratio $`R_2`$ is closely related to the Reynolds number Re: it differs in representing length and velocity measures based on the rms fluid flow rather than on the interface dynamics and, because the length scales associated with the velocity gradients are smaller than the domain size, is significantly smaller than Re. The dominance (by a factor ten) of inertial over viscous forces is, at late times, nonetheless clear (Fig.3).

We finally ask whether, at the largest Re values we can reach, there is in fact significant turbulence in the fluid flow. One quantitative signature of turbulence is the skewness $`S`$ of the longitudinal velocity derivatives; this is close to zero in laminar flow but approaches $`S=0.5`$ in fully developed turbulence . We do detect increasingly negative $`S`$ as Re is increased but reach only $`S0.3`$ for Re $`350`$ . This suggests that at our highest Re’s, turbulence is at most partially developed – a view confirmed by visual inspection of velocity maps . Grant and Elder’s suggestion of an eventual transition to turbulent remixing thus remains open.

In conclusion, we have presented LB simulation data for 3D spinodal decomposition which spans an unprecedented range of reduced time and length scales. At $`t\stackrel{<}{}10^2`$ (Re $`\stackrel{<}{}1`$) we observe linear scaling, as announced in the previous literature . This is followed by a long crossover ($`10^2\stackrel{<}{}t\stackrel{<}{}10^6`$, or $`1\stackrel{<}{}`$Re$`\stackrel{<}{}50`$) connecting to a regime in which inertial forces clearly dominate over viscous ones (see Fig.3); our work is the first to unambiguously probe this regime in 3D. In the region so far accessible ($`10^6\stackrel{<}{}t\stackrel{<}{}10^8`$, or $`50\stackrel{<}{}`$Re$`\stackrel{<}{}350`$) Furukawa’s prediction of $`t^{2/3}`$ scaling is obeyed, to within simulation error. An open issue is whether this regime marks the final asymptote or whether a further crossover occurs to a turbulent remixing regime (saturating Re) as proposed by Grant and Elder . If it does, we have shown that any limiting value of Re must significantly exceed their estimate of $`10100`$.

We thank Craig Johnston, Simon Jury, David McComb, Patrick Warren and Julia Yeomans for valuable discussions. Work funded in part under the EPSRC E7 Grand Challenge.