## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# The importance of fluctuations in fluid mixing

Kadau |

## Supporting Information

#### Files in this Data Supplement:

SI Figure 5SI Text

SI Movie 1

SI Movie 2

SI Movie 3

SI Movie 4

SI Movie 5

SI Figure 6

SI Figure 7

SI Figure 8

SI Figure 9

SI Figure 10

SI Figure 11

SI Figure 12

SI Figure 13

SI Figure 14

SI Figure 15

SI Figure 16

SI Figure 17

SI Figure 18

**Fig. 5.** Schematic representation of experimental 2D RT set-up. The magnet's "Faraday pole pieces" are shown, which provides a vertical gradient in the field such that the magnetic force cÑ(*H*^{2}) is approximately constant over the vertical region around the interface. The cell is placed midway between the pole pieces, and the two magnetically permeable wires of diameter 0.088 cm and relative permeability m_{r} = 2.6 ± 0.5 are attached to the outside are shown. In the presence of the magnetic force, the denser but paramagnetic fluid 1 is on top. Upon switching off the magnet current, the interface becomes unstable, and the heavier fluid falls to the bottom. Note that the magnetically permeable wires play no role once the instability sets in because there no longer is a magnetic field to perturb.

**Fig. 6.** 2D RT mixing described by atomistic MD simulations: penetration velocity for spikes (solid line) and bubbles (dashed line) at Atwood number *A* = 0.75. At *t* = 3,000 *t*_{0} the linear increase of velocity changes. Breakups occur and lead to a Stokes-like constant velocity regime. The simulation consisted of more than 12 million particles interacting via a short-range Lennard-Jones potential and ran for approximately 1 million time steps (5). The exponential growth time t equals 436.9 *t*_{0}. The system has dimensions of (width ´ height ´ depth) 1,400 *r*_{0} ´ 3,063 *r*_{0} ´ 5.25 *r*_{0}, and a gravity of 0.0004 *g*_{0} (for units, see *SI Text*).

**Fig. 7.** The influence of the system width on 2D RT mixing as described by four different atomistic MD simulations with reduced width from left to right at *t* = 4,000 *t*_{0} (top) and *t* = 8,000 *t*_{0} (bottom) (the left frames corresponds to the run shown in Fig. 6). The growth is reduced when merging processes increase the largest horizontal features (to a size comparable to the system).

**Fig. 8.** Influence of the system width *W* on RT mixing (same four runs as shown in Fig. 7).

**Fig. 9.** Chaotic trajectories (time evolution top to bottom) of atomistic MD simulations for two almost identical simulations involving approximately 3.3 million atoms. The only difference is that for the simulation shown on the left, initially the velocity of one light particle near the interface was negatived (*v* ? -*v*). Due to the chaotic nature of the governing equations of motion, the exact outcome is changed; however, the average of macroscopic properties will not be affected (to within fluctuations).

**Fig. 10.** Mass density profile evolution for the atomistic DSMC simulation shown in Fig. 1. Compressibility effects are evident.

**Fig. 11.** Penetration of the falling spike front and the rising bubble front as a function of time for a 2D RT instability as simulated by atomistic simulation (DSMC run shown in Fig. 1). Between *t* = 4 t and 12 t, the penetration of spikes and bubbles shows a quadratic behavior in time as can be seen by the linear region of the penetration-*t*^{2} plot (*Left Inset*) and the penetration velocity plot (*Right Inset*). In this regime, the coefficient a_{B/S} characterizes the growth. Depending on the method of determination, a_{B/S} can vary up to 20% (the a_{B/S} stated in the main graph is calculated by ^{·}*h*^{2}/4*Agh*). After that, spikes start to break-up and are no longer supported by necks to the bulk, and the growth is reduced. The penetration velocity is reduced and fluctuates around a value indicating a Stokes-like behavior. The penetration of individual spikes (S1, S2) and bubbles (B1, B2) shows the variation of movement of individual features.

**Fig. 12.** 2D DSMC simulations of the RT instability for two different Atwood numbers *A* = 0.43 (*Left*) and *A* = 0.91 (*Right*).

**Fig. 13.** The time series of the total velocity spectrum for the 2D DSMC simulation shown in Fig. 1 shows the progression toward larger length scales with time. The velocity field is determined by finding the mass-weighted average particle velocity within each DSMC cell, and the spectrum is then calculated for a horizontal band centered at the initial interface position, which encompasses 1/16 of the area of the domain. Note that the spectra for each of the two velocity components (not shown) are qualitatively similar in shape, although the vertical component yields a spectrum that is approximately five times in magnitude. Power laws corresponding to *k*^{5/3} and *k*^{-3} are shown for reference.

**Fig. 14.** Three-dimensional RT DSMC simulation involving more than 7 billion atoms. The individual development of different spikes (S2, S2) and bubbles (B1, B2), as well as breakups toward the end, is evident. The outer scale Reynolds number at the end is 700 (for dimensions and other specifications of this run, see *SI Text*).

**Fig. 15.** The contour (white) of the interface of the 2D RT mixed interface as simulated by DSMC (same simulation as shown in Fig. 1). The system was divided into 4,096 ´4,096 cells each containing »20 particles. A cell belongs to the contour if |*N*_{H}-*N*_{L}|/(*N*_{H} v+ *N*_{L}) < *C*, where *N*_{H} is the number of heavy particles in the cell, *N*_{L} is the number of light particles in the cell, and the contour parameter *C* = 0.25.

**Fig. 16.** The fractal dimension *D* of the contour shown in Fig. 15 (red) as determined by the slope of the straight line fit of the log-log plot of the number of boxes containing the contour as a function of inverse box size (box counting method). If the contour parameter *C* (see caption Fig. 15) changes, the fractal dimension can change slightly.

**Fig. 17.** Fractal dimension of the same simulation shown in Fig. 15 at two different times. At early times (black symbols), the data exhibits two slopes. For large box sizes (left, blue straight line fit), the contour has a dimension of close to one, i.e., is nonfractal. For smaller box sizes, at scales where mixing already occurred, the fractal dimension is larger than one. At later times (red symbols), the contour shows a fractal dimension, i.e., the data fits a straight line.

SI Figure 18

**Fig. 18.** Fractal dimension of 2D RT mixed interfaces as a function of time for various atomistic DSMC simulations having different gravities *g*, Atwood numbers *A*, and system edge lengths *L*. The dimensions were calculated only for box sizes small enough to be affected by the mixing (solid lines) (as the solid black straight line fit in Fig. 17). The dashed line shows the fractal dimension using all data points (as the dashed black straight line fit in Fig. 17). The exponential growth time t for simulation R1 equals 2,584 *t*_{0} and is the system shown in Fig. 1.

SI Movie 1

**Movie 1.** Atomistic simulation (DSMC run shown in Fig. 2) involving 570 million particles of a 2D RT instability, where the heavy fluid is on top (red) of the light fluid (blue), and the gravitational field pointing downwards. The initially flat interface roughens because of thermal fluctuations, and modes near the most unstable mode according to continuum theory are developing at first. Later, merger processes and complex hydrodynamic flow becomes dominant. At later times, breakups occur that reduce the overall growth of the mixing layer.

SI Movie 2

**Movie 2.** Magnetic levitation RT experiment shows similar features as the atomistic simulated instability in 2D. The technique of magnetic levitation allows for flat initial interface conditions.

SI Movie 3

**Movie 3.** Atomistic simulation (DSMC) involving 7.1 billion particles of a 3D RT instability. As in 2D, thermal fluctuations trigger the onset of the instability, which develops into a complex hydrodynamic flow regime. At later stages, breakups of the falling spikes are visible.

SI Movie 4

**Movie 4.** Molecular dynamics simulation of a 2D RT Instability with a huge initial single mode perturbation. After the spikes (red) almost reached the bottom, the velocities of all of the atoms were negatived (*v*_{i} → -*v*_{i}). With infinite precision (no numerical error in the calculation of the trajectories), the evolution should reverse and the two fluids would be completely separated as in the initial conditions (due to the time reversibility of the underlying equations of motion). The reverse happens for a little while before round-off errors in combination with the chaotic nature of the system pushes toward larger entropy and the mixing-layer growth again.

SI Movie 5

**Movie 5.** Atomistic simulation (DSMC simulation shown in Fig. 4) involving 140 million particles of a 2D RT instability (middle). To the left and right simulations with only a half and a quarter of the system width is shown, respectively. Toward the end, the reduced widths of the smaller systems influences the mixing. However, no systematic differences (except for larger fluctuation in the smaller systems) in the terminal velocity and the time of breakup have been observed. Therefore, the simulation results presented here are not influenced by the boundary conditions (more than 10 l_{max} are initially present).

**SI Text**

**Experiments.** In our experiments we use a glass cell 15 cm in height, 7 cm in width (span-wise), and 0.3 cm thick (SI Fig. 5). Fluid 1 (density r_{1} = 1.394 g cm^{-3}, kinematic viscosity n_{1} = 0.051 cm^{2}/s) is a moderately strong paramagnetic mixture of water, 58.9 wt.-% MnCl_{2}^{.}4H_{2}O, »1 wt.-% surfactant octa(ethylene glycol) dodecyl ether ("C_{12}E_{8}") to reduce surface tension, and a small amount of rhodamine 6G dye; fluid 2 (r_{2} = 0.773 g cm^{-3}, kinematic viscosity n_{2} = 0.043 cm^{2}/s) is weakly diamagnetic hexadecane, which is immiscible with fluid 1 and an interfacial tension of 2.6 erg/cm^{2}. This combination of fluids (with Atwood number A = 0.29) has wetting properties such that the meniscus is virtually absent. For *t* < 0 a vertical magnetic field gradient is applied, and the denser fluid 2 resides above the less dense fluid 1 when c_{1}Ñ(*H _{y}*

^{2}) > (r

_{1}- r

_{2})

*g*. Here, c

_{1}is the (positive) magnetic susceptibility of fluid 1 and we assume that the susceptibility of fluid 2 is negligible. On switching off the magnetic field, the metastable layering becomes unstable in the presence of a perturbation, and the instability occurs in the presence of uniform gravity only (1). Using the technique of planar laser induced fluorescence [PLIF (1)], we illuminated the cell from above with an »1 mm thick "sheet" of light, created by passing a laser beam (frequency doubled Nd:YAG laser at 532 nm) through a cylindrical lens. The sheet passed through the midplane of the cell, causing the dye in fluid 1 to fluoresce. We collected videos of the lower portion of the cell (owing to line-of-sight constraints due to the magnet) as a function of time using a CCD camera at 60 frames per second. On switching off the magnetic force, which has a decay time of 0.065 s, the instability passed first through the linear regime with a growth rate

*n*= (20 ± 3) s

^{-1}and wavelength l

_{max}= (0.7 ± 0.1) cm for the fastest growing mode (as measured). To extract a

_{S}, the spike front

*h*

_{S}is taken as the furthest advance for which at least 5% of the pixels in a given row correspond to the fluorescent fluid 1 (similar procedure as applied for the simulations).

**Atomistic Simulations and the Hierarchy of Methods.** A fluid, like any other material, is composed of atoms and molecules. In the absence of relativistic effects quantum mechanics, is in principle the most general description of matter, and it has been used to model the behavior of individual molecules or small collections of molecules, via numerical solution of the time-dependent Schroedinger equation (TDSE). However, such *ab initio* calculations are extremely computationally intensive, and can currently be used only for very small systems.

For cases where the set of molecules is essentially classical in nature, i.e., when quantum effects may be ignored or accounted for in a simple way, atomistic methods such as molecular dynamics (MD) are the next most general technique (3). In MD, Newton's classical equations of motion are solved numerically for a large number of atoms or molecules, and the resulting trajectories are then sampled to provide statistical estimates of the macroscopic properties and/or flow characteristics of the system. Assuming that we have access to the correct set of inter-particle potentials for the substance in question, MD must by definition yield the correct behavior for these macroscopic properties. Atomistic methods other than MD, such as direct simulation Monte Carlo (DSMC) (4), entail certain additional assumptions or approximations, and so are somewhat less general in their regime of validity than MD.

For flows where the atomistic nature of matter may be ignored, and for which deviations from the ensemble-mean flow due to molecular fluctuations are insignificant, a continuum description of matter is valid. Such a description takes the form of a set of partial differential equations (PDEs) for the evolution of moment fields of the microscopic dynamics, such as density or flow velocity. The most general possible set of continuum PDEs can be derived rigorously from the ordinary differential equations (ODEs) governing MD, or from the associated Liouville equation, and it describes the evolution of the macroscopic moment fields with which we are familiar, e.g., the density and velocity fields, although it is in general not closed. Because this most general system of PDEs is not closed, it cannot be solved directly, and various assumptions and approximations have been proposed to generate closure, although none of the closures commonly used today incorporate fluctuations.

By assuming a particular closure in the form of a simplified expression for the stress tensor, we can obtain the Euler equations or the Navier-Stokes (NS) equations, which are in principle the least general of the possible descriptions of a fluid. Nevertheless, the Euler/NS equations have had enormous success over the decades in describing many of the flow scenarios with which we are familiar in everyday life. This indicates that the approximations inherent in Euler/NS are valid for these familiar regimes, and, in fact, Euler/NS simulations are by far the most popular types of fluid modeling. Furthermore, there is no technique at this time by which MD can be extended to describe the 10^{23} particles of macroscopic fluid samples.

In the context of this work, we have chosen to perform atomistic (i.e., MD and DSMC) simulations, as they are the most general techniques available which can handle a large enough system size to describe a complex flow. This allows a direct comparison with the results of traditional Euler/NS simulations (as well as experiment), and an exploration of when and how the breakdown of the continuum hypothesis happens.

**Molecular Dynamics.** As mentioned above, MD involves the direct numerical integration of Newton's equations of motion for a large number of interacting atoms or molecules. Originally formulated in 1957 by Alder and Wainwright, MD simulations were originally restricted to the consideration of small systems of hard spheres, although they have since been expanded and generalized to describe complex systems of billions of molecules.

Methods for MD simulation are many and varied. Our simulations involved large systems of particles interacting via the continuous, truncated Lennard-Jones (LJ) potential, and were performed on a parallel architecture using the Scalable Parallel Short-range Molecular dynamics (*SPaSM*) code (5,6). The initial state of the Rayleigh-Taylor (RT) system was specified by given values of the global temperature, acceleration due to gravity, and pressure at the interface. Note that since the system is somewhat compressible, it was precompressed according to the known LJ equation of state. The molecular trajectories generated by MD were then used to study the collective behavior of the system via, for example, density and velocity fields.

**Direct Simulation Monte Carlo.** The most time-consuming portion of the MD algorithm is the force evaluation, in which the explicit interactions between each pair of particles is determined. The DSMC scheme was originated by Bird in the 1960s to avoid this computational bottleneck (4). DSMC, which is typically 30-50 times faster than traditional MD, employs a time-stepping algorithm that consists of the following steps:

1. All particles are advected in a straight line (or a parabola in the presence of gravity or other external force) according to their velocity, and for a time equal to the time step.

2. Boundary conditions (e.g., periodic or specular wall) are enforced.

3. Particles are sorted by position into cells.

4. Random pairs of particles are selected from within each cell. Each pair is stochastically accepted or rejected based on its relative velocity, and then the particles are "collided" by adjusting their respective velocity vectors as if they had undergone a hard sphere collision, with uniformly distributed impact parameter. The number of pairs of each species are chosen so that the average collision rates will correspond to the hard sphere collision rates from kinetic theory. Note that the collision rates are proportional to the square of the particle diameter s.

The traditional application of DSMC has been to rarified gas flows, in which the Knudsen number is large, and a cubic mean free path contains a very large number of particles. Usually in such applications, each simulation particle is not taken to represent an actual fluid molecule, but is instead interpreted as an aggregate of a large number of such physical particles. In this way, DSMC has been applied to length scales many orders of magnitude greater than those to which standard MD techniques may currently be applied. However, there is a drawback to the scale-flexible nature of DSMC simulations involving such "effective particles" in that they drastically over-estimate the magnitude of microscale thermal fluctuations. For an unsteady flow such as the RT instability, this can be a serious problem, and we have therefore performed all of our simulations in a regime in which (as in MD) each simulation particle represents a true molecule of the fluid. This allows us to accurately capture the magnitude and behavior of fluctuations. In addition, we are operating in a reasonably dense phase, so that there are only a small number of particles present in a cubic mean free path. This is frowned upon in traditional DSMC, although in reality the only drawback is that the transport coefficients are not given explicitly by simple kinetic-theoretical expressions, as they are in rarified gas DSMC. Instead, they are in principle given by more complex expressions involving the values of the cell size and time step. Although such expressions have not yet been derived in closed form, the transport coefficients themselves may still be measured empirically by standard means, e.g., via Green-Kubo or Einstein relations, as we have done to determine the viscosities of our two fluids. For our non-dimensional simulation parameters (cell size 3.4, time step 0.3, initial interfacial pressure 0.5, and mass densities r_{1} = 2.5, r_{2} = 0.5, for units see main text), the viscosities h_{1} and h_{2} of the light and heavy fluids are 0.65 and 1.46, respectively. Note that this cell size is about 7.5 times the mean free path, and a cubic cell contains about 20 particles. Since DSMC does not allow for an interfacial tension, we increased the collision rate between heavy and light particles by a factor of 10 (as compared to the "regular" collision rate between the same type) to reduce diffusion at the interface.

**Simulation Units.** In the DSMC simulations the diameter *r*_{0} of the colliding hard spheres, the temperature *k*_{B}*T*, and the mass of the light fluid are set to unity. If the light fluid is assumed to be methanol with a molecular diameter *r*_{0} = 0.402 nm and an atomic mass of 32 at 300 K, the dimensions for time and acceleration are *t*_{0} = 1.44 ps and *g*_{0} = 1.94 ´ 10^{14} m/s^{2}, respectively. With this assumption, the quasi 2D DSMC simulation shown in Fig. 1 (3) containing 570 (137) million particles, represents a system with dimensions 5.60 (2.80) mm ´ 5.60 (2.80) mm ´ 2.73 nm, a gravity of 6.13 (24.72) ´ 10^{8} *g*_{Earth}, and a physical simulation time of 90 (35) ns. These simulations including on-the-fly visualization and system dumps were performed by using 65,536 (8,192) CPUs of the Blue Gene/L system for 30 (22) h. The largest 3D DSMC simulation consists of 7.1 billion particles, and dimensions of 1.05 mm ´ 1.05 mm ´ 1.05 mm, a gravity of 7.42 ´ 10^{9} *g*_{Earth}, and a physical simulation time of 13 ns, and was run for 20h on 32768 CPUs of Blue Gene/L (this 3D run had a better load-balance, less communication overhead, and system dumps were performed less frequent as compared to the quasi 2D simulations). Applying a similar procedure for the MD simulations yields *r*_{0} = 0.402 nm, *t*_{0} = 1.11 ps, and *g*_{0} = 3.26 ´ 10^{14} m/s^{2} (6). The biggest quasi 2D MD system has been run on 200 CPUs of an Intel Pentium4 Xeon 2.4Ghz Linux cluster for 8 days. This run consisted of 12 million atoms with the dimensions 0.563 mm ´ 1.23 mm ´ 2.11 nm, had a physical run time of about 10 ns, and a gravity of 1.32 ´ 10^{10} *g*_{Earth} (SI Figs. 6 and 7). It should be noted that, even at this these apparently high gravities, the ratio of the particle interaction force to the gravitational force is still very large, i.e., on the order of 10^{4} to 10^{5}. The largest outer scale Reynolds number achieved in DSMC and MD simulations is 2700 and 1500, respectively.

**Initial Growth.** For the RT instability with a heavy fluid of density r_{1} on top of a light fluid of density r_{2}, linear stability analysis of the Navier-Stokes equations predicts that, for small times, the various interfacial modes will grow exponentially and independently as *e ^{n}(*k)

*t*, with a growth rate dispersion relation

*n*(

*k*). Here,

*k*= 2p/l is the wavenumber of the mode with wavelength l, and

*n*(

*k*) is in general a function of

*g*, r

_{1}, r

_{2}, the kinematic viscosities n

_{1}, n

_{2}, and the interfacial tension g between the two fluids (7). Note that the peak value

*n*

_{max}of

*n*(

*k*) yields a characteristic time scale time scale t= 1/

*n*

_{max}and length scale l

_{max}for the system. In atomistic simulations, the initially flat interface is naturally roughened by thermal fluctuations, which provide excitations in all modes. These modes subsequently grow according to the dispersion relation

*n*(

*k*). Random fluctuations influence the growth rates of individual modes on the small length scales represented by these simulations but, under averaging, there is general agreement with continuum predictions in the initial linear regime (6). The most dominant mode that evolves from an unperturbed interface, as achieved by magnetic levitation techniques as described above (8-9 maxima in a 7-cm-wide system → l

_{max}= 0.78-0.88 cm; see Figs. 1 and 3), is also in good agreement with the most unstable mode as predicted by continuum theory (l

_{max}= 0.83 cm) (see also ref. 1).

**Fractal Dimension of the Interface.** The fractal dimension of the RT interface has been calculated by employing a box-counting scheme (8,9) for our atomistic DSMC simulations of RT instability in 2D and 3D (see SI Figs. 15-18). The fractal dimension for our atomistic simulations are increasing in time and are between 1.6 and 1.8 (2D). This increase in time means that the interface is not self-similar, whereas a constant value is indicative for a self-similar structure. Previous experiments yield an increase of the fractal dimension from 1.5 to 1.8 (8). Continuum simulation result in a non-evolving fractal dimension between 1.5 and 1.8 (8,9). For a 3D simulation we obtained a fractal dimension of 2.6-2.7. The dimension of a slice of the 3D interface is 1.6-1.69, the dimension of an equivalent 2D simulation yields 1.6-1.65.

1. Carlès P, Huang Z, Carbone G, Rosenblatt C (2006) *Phys Rev Lett* 96:104501.

2. Kychakoff G, Howe RD, Hanson RK, Drake MC, Pitz RW, Lapp M, Penney CM (1984) *Science* 224:382-384.

3. Rapaport DC (1995) *The Art of Molecular Dynamics Simulations* (Cambridge Univ Press, Cambridge, UK).

4. Bird GA (1994) *Molecular Gas Dynamics and the Direct Simulation of Gas Flows* (Clarendon, Oxford).

5. Beazley DM, Lomdahl PS (1997) *Comput Phys* 11:230-238.

6. Kadau K, Germann TC, Hadjiconstantinou NG, Lomdahl PS, Dimonte G, Holian BL, Alder BJ (2004) *Proc Natl Acad Sci USA* 101:5851-5855.

7. Chandrasekhar S (1961) *Hydrodynamic and Hydromagnetic Stability* (Oxford Univ Press, Oxford).

8. Dalziel SB, Linden PF, Youngs DL (1999) *J Fluid Mech* 399:1-48.

9. Hasegawa S, Nishihara K, Sakagami H (1996) *Fractals* 4:241-250.