2892
Ind. Eng. Chem. Res. 2006, 45, 2892-2896
Cellular Automata Model of Phase Transition in Binary Mixtures† Carolina Vannozzi,‡ Diego Fiorentino,§ Matteo D’Amore,§ David S. Rumshitzki,⊥ Andreas Dress,# and Roberto Mauri*,‡ Department of Chemical Engineering, UniVersity of Pisa, 56126 Pisa, Italy, Department of Chemical and Food Engineering, UniVersity of Salerno, Italy, Department of Chemical Engineering, City College of CUNY, New York, New York 10031, and Institute for Computational Biology, Shanghai Institute for Biological Sciences, Shanghai, People’s Republic of China
We investigate the dynamical behavior of spinodal decomposition in binary mixtures using a variation on the Rothman-Keller cellular automaton, in which particles of type A(B) move toward domains of greater concentration of A(B). Domain growth and system morphologies are determined via the pair correlation function, revealing that (i) the characteristic domain size grows as R(t) ≈ t1/3 and (ii) the configurations of the mixtures at different times are self-similar. These results did not change when different game rules were adopted, indicating that self-similarity and the 1/3 scaling law constitute fundamental properties of any diffusiondriven phase separation process. The same model also was applied to describe the mixing process, and, as expected, it was determined that the characteristic time was dependent on the square of a characteristic linear dimension of the system. I. Introduction The dynamics of the phase separation of a binary mixture is a challenging problem that has attracted recent as well as past interest (reviews on phase separation can be found in Hohenberg and Halperin1 and Gunton et al.2). Traditional theoretical approaches to phase separation have either taken a continuum mechanics, mean-field approach (van der Waals,3 Cahn and Hilliard4) or, more recently, have used molecular dynamics, i.e., integration of Newton’s equations (see the review in Langer5), to understand the behavior of the phase separation of binary mixtures. Both approaches have limitations. Continuum mechanics breaks down when the interface between the two phases becomes only a few molecules thick; in addition, the resulting equations are difficult to solve for all but the simplest geometries, and they are subject to strong instabilities when linearized. More importantly, still at issue are the correct constitutive laws for these theories. On the other hand, molecular dynamics simulations are limited by the number of molecules that can be followed and by the fact that the computational complexity increases significantly for nonspherical structures. Cellular automaton models avoid each of these problems (a classic collection of papers on cellular automata intended as a companion to von Neumann’s theory of self-reproducing automata can be found in the work by Burks6). They use welldefined rules, uncomplicated by instabilities, with a simplicity that allows a simulation to follow a large number of molecules (105-106), even when a “structure irregularity” is included and a fairly small computer is used. As shown by Rothman and * To whom correspondence should be addressed. Tel.: ++39-050511248. Fax: ++39-050-511248. E-mail:
[email protected]. † This paper was originally submitted for publication in the Reuel Shinnar Festschrift special issue of Industrial and Engineering Chemistry Research (published in 2004, Volume 43, Issue 2). It was received in its original form May 2, 2003 and accepted on November 4, 2003; however, because of a series of communication lapses, it was not available for print until now. ‡ University of Pisa. § University of Salerno. ⊥ City College of New York. # Shanghai Institute for Biological Sciences.
Figure 1. Initial state of an 100 × 100 system with φ ) 0.5.
Figure 2. Particle pair with its 10 neighbors.
Keller7 and by many other authors (see the review article by Rothman and Zaleski8), automata are particularly valuable in modeling phase separation, because they can include any microscopic feature within the model providing, in addition, testable macroscopic predictions that are valid on the length and time scales of the physical experiments. The purpose of this paper is to take a look at the diffusioncontrolled dendritic growth of the bicontinuous structures that form during the phase separation of a critical binary mixture. This process, which is generally referred to as spinodal decomposition, was studied analytically by Lifshitz and Slyozov (see the work of Lifshitz and Pitaevskii9), who found that, during diffusion-controlled phase separation, the mean size of the
10.1021/ie051240w CCC: $33.50 © 2006 American Chemical Society Published on Web 03/17/2006
Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 2893
mixtures is studied using cellular automaton techniques, modifying the original Rothman-Keller approach by introducing a parameter that endows the system with a stochastic evolution. A similar approach was followed by Alexander et al.,19 although using different “rules of the game”. After the characteristics of the cellular automaton are described in section II, in section III, the results of the two-dimensional (2D) simulations are reported, yielding the scaling laws of the binary phase separation, which are then discussed in the final section (section IV). II. Cellular Automaton Figure 3. Rules of the game.
single-phase domains grows in time as t1/3. This result was confirmed experimentally for critical mixtures by Goldburg and co-workers (see the work of Chou and Goldburg10 and references therein), Knobler and co-workers (see the work of Wong and Knobler11 and references therein), and Beysens and coworkers (see the work of Geunoun et al.);12 later, the same 1/3 growth law was found to hold for off-critical quenches as well (see the work of Cummings et al.13 and White and Wiltzius14) and even for isolated drops in secondary emulsions (see the work of Califano and Mauri15). Similar results were found in numerical simulations (see the work of Tanaka and Araki16 and Vladimirova et al.17), indicating that phase separation is a selfsimilar process, similar to any other Landau-type phase transition.18 In this article, the evolution of phase-separating critical
Our automaton uses a regular 2D square lattice with size n and toroidal symmetry (see Binder20,21). Lattice sites can be occupied by two different types of “particles”, A or B (we can think of them as molecules or a combination thereof), and no holes are permitted in the grid, which corresponds to assuming that the mixture is incompressible. Graphically, A particles are represented as black cells, and B molecules as white ones. Initially, the system is well-mixed, which means that its initial state is a random configuration of A and B particles with fixed overall composition φ. In Figure 1, a typical initial configuration is shown with n ) 100 and φ ) 0.5, i.e., with the same number of A and B particles. Physically, in a phase-separating binary mixture, the attractive force between two particles of different types (i.e., a black and a white particle), F12, generally is different from the forces between two particles of the same type (i.e., two black or two
Figure 4. Snapshots of the morphology of a 100 × 100 system at (a) t ) 1000, (b) t ) 5000, (c) t ) 10 000, and (d) t ) 30 000.
2894
Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006
white particles), F11. From a thermodynamic viewpoint, the difference between these two forces, F11 - F12, is proportional to the excess free energy (or the enthalpy of mixing), which, in most cases (i.e., for regular mixtures) is independent of temperature; this so-called enthalpic force can be either positive or negative, helping the mixture to either separate or mix, respectively (when the enthalpic force is zero, the mixture is called ideal, because it behaves as an ideal gas). On the other hand, there is an entropic force as well, derived from the ideal part of the free energy, which is proportional to the temperature and is always helping the mixture to mix. Accordingly, a binary mixture phase-separates when the enthalpic driving force is positive and prevails over its entropic counterpart. We translated this idea into the following rule of the automaton: two adjacent particles (i,j) will swap positions when Nij g 5, and they will not when Nij < 5. Here, Nij is the sum of the unlike neighboring particles (i.e., A particles surrounding B particles and B particles surrounding A particles) of a generic particle pair (ij), where the 10 neighboring sites of a generic particle pair are indicated in Figure 2 as the five light gray particles surrounding the black A particle and the five dark gray particles surrounding the white B particle. In the following, we will refer to Nij as the local driving force, as it induces, when large enough, the local antidiffusive transport process that is eventually responsible for the phase separation. Obviously, this rule (and the concept of local driving force as well) is applicable only when the two adjacent particles are of different types; otherwise, the switch would leave the system unchanged. In addition, the particle pairs can be located both horizontally and vertically, because we must account for switches in both directions. Referring to Figure 3, in the top configuration, the A particle of the pair is surrounded by two B particles, whereas the B particle is surrounded by four A particles, and therefore, the switch will occur; on the other hand, in the bottom configuration, the adjacent unlike particles are two and one, respectively, and therefore, they will not swap. Clearly, this rule allows only the switch between two adjacent particles and therefore corresponds, macroscopically, to a diffusive transport mechanism. Convection, on the other hand, corresponding to a net translation of the single-phase domains, is not modeled by our automaton. Basically, this rule means that a particle pair will swap when the new configuration is more “convenient” than the old one, that is, when the swap increases the number of surrounding particles of the same type, thereby decreasing the potential energy (i.e., the enthalpy) of the system. When Nij ) 5, i.e., in the absence of any enthalpic driving force, it is irrelevant, from a purely energetic point of view, whether the pair swaps, and therefore, in the absence of thermal fluctuations (i.e., at 0 K), there should be no swaps. Therefore, forcing adjacent particles to swap when Nij ) 5, as we do here, accounts for a nonvanishing temperature (i.e., an entropic force), which however, because of the extreme simplicity of the model, we cannot indicate quantitatively. During the simulation, at each time step, the automaton first chooses a lattice site (i,j) at random, selects as the adjacent site either (i + 1,j) or (i,j + 1), computes Nij, and finally swaps the particles of the pair when Nij g 5. III. Results Several simulations were run to study the phase separation of an incompressible binary mixture. We considered only the critical case, corresponding here to φ ) 0.5, where half of the particles are of type A and the other half are of type B. As the
Figure 5. (a) Radial pair correlation function, C(r,t), after an instantaneous critical quench as a function of r for different times t and (b) as a function of the self-similar parameter z ) rt-1/3.
initial state of the system is completely random, it corresponds to that following an instantaneous quench from a state above the spinodal point, where A and B are miscible, to a state below, where they are immiscible. All simulations were performed on a 100 × 100 2D lattice with toroidal boundary conditions. The number of lattice sites (and therefore the number of particles) is large enough to ensure the statistical accuracy of the data obtained during these simulations. All simulations were allowed to evolve for 30 000 time steps, and the state of the system was recorded at constant intervals. Figure 4 shows the state of a single system at four different times, with the black and white domains denoting the zones occupied by A and B particles, respectively. As expected, during the early stages of the process particle switches are very frequent, reflecting the fact that the driving force is very large (see Figure 4a-c). On the other hand, the late stages (see Figure 4d) are characterized by domains separated by relatively sharp interfaces, whose size increases monotonically with time. In addition, the configurations exhibit an interwoven, bicontinuous topology that is common in spinodal decomposition. To make a quantitative analysis of the domain growth, it is useful to evaluate the time-dependent domain size, which
Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006 2895
Figure 6. Snapshots of the morphology of a 100 × 100 mixing system at t ) 1, 500, 1000, 5000, 8000, and 10000.
involves considering the product [Θ(i,j) Θ(i + m,j + n)], where Θ(i,j) ) 1 or -1 depending on whether the lattice site (i,j) is occupied by an A or B particle, respectively.22 At this point, the pair correlation function can be defined as the average of this product over all pairs of lattice sites whose distance r ) (k2 + l2)1/2 is kept constant, i.e., with
C(r,t) )
1 4(n - 2r)
n
n
∑ ∑ Θ(i,j)‚Θ(i + k,j + l) 2 i)1 j)1
where the factor 4 in the denominator accounts for the four neighbors of each particle, whereas the normalizing factor (n - 2r)2 ensures that each pair is counted once. In Figure 5, the
correlation function C(r,t) is plotted as a function of r at different times t, showing that the different curves collapse into a selfsimilar master curve, C(z), with z ) rt-1/3, although it should be stressed that the master curve is rather murky, probably due to the fact that a 100 × 100 domain is far too small. In any case, this result reveals that, defining the single-phase domain size R as the distance at which the correlation function C(r,t) crosses the horizontal axis [i.e., C(R,t) ) 0], satisfies the scaling law R(t) ∝ t1/3, in agreement with the theoretical predictions by Lifshitz and Slyozov9 and the numerical results based on the solution of the mean-field equations (see the work of Tanaka and Araki16 and Vladimirova et al.17). Both the self-similarity of the mixture morphology and the 1/3 power law of the droplet
2896
Ind. Eng. Chem. Res., Vol. 45, No. 8, 2006
growth rate agree well with all the experimental results mentioned in the Introduction, indicating that they must be fundamental properties of the physical process and do not depend on the particular type of binary mixture considered. At the end, we checked that in 50 × 50 and 150 × 150 domains the morphology of the mixture was still bicontinuous and followed the 1/3 power law. The simulations were also repeated introducing one of the following changes in the swapping rule: (a) when Nij ) 5, the swapping of neighboring particles takes place with a probability P ) 2/3 (instead of 1), thereby describing phase separation at a different temperature, or (b) the number of neighboring particles that are considered is only 6 (instead of 10), and therefore, the pair will swap when Nij g 3. In both cases, the results were simply shifted in time with respect to the “old” ones, therefore reinforcing the conclusion that self-similarity and the 1/3 scaling law are fundamental properties of any diffusion-driven phase separation process. Finally, we also modeled the mixing process of two perfectly mixable fluids. As we did previously, we translated the physics of mixing into the following simple rule of the automaton: two adjacent particles A and B will swap when Nij < 5, and they will not swap when Nij g 5. The results are displayed in Figure 6, which shows that the correlation fucntion C(r,t) decays to zero within a time that is dependent on the sqaure of the size of the domain, as expected. IV. Conclusions We constructed a cellular automaton that simulates the phase separation undergone by a binary mixture that is brought from its single-phase, homogeneous initial state to the region of instability of its phase diagram. Using a very simple “game rule”, we found that the characteristic size R of the single-phase domains grows in time as tn, where n = 1/3, in agreement with both the theoretical predictions and the mean-field results. More importantly, we verified that the configurations of the mixtures at different times are self-similar, again in agreement with both experimental findings and numerical simulations. These results are far from obvious, as the only “physics” built into the automaton is that the A and B particles composing the mixture are attracted to particles of the same type more than to particles of different type (i.e., A is attracted to A and B to B more than A is attracted to B). On the other hand, both the integration of the mean-field equations of motion and the implementation of molecular dynamics techniques are much more sophisticated approaches, with many physical assumptions on the system. Therefore, the fact that we basically obtain the same results indicates that self-similarity and the 1/3 scaling law are fundamental properties of any diffusion-driven phase separation process, and, in fact, our results did not change when we adopted different game rules. The same model also was applied to describe the mixing process, and, as expected, it was determined that the characteristic time was dependent on the square of a characteristic linear dimension of the system.
Literature Cited (1) Hohenberg, P. C.; Halperin, B. I. Theory of Dynamic Critical Phenomena. ReV. Mod. Phys. 1977, 49, 435. (2) Gunton, J. D.; San Miguel, M.; Sahni, P. S. Spinodal Decomposition of Binary Mixtures. In Phase Transition and Critical Phenomena; Domb, C., Lebowitz, J. L., Eds.; Academic Press: London, 1983; Vol. 8. (3) van der Waals, J. D. The Thermodynamic Theory of Capillarity under the Hypothesis of a Continuous Variation of Density. Reprinted in J. Stat. Phys. 1979, 20, 200. (4) Cahn, J. W.; Hilliard, J. E. Free Energy of a Nonuniform System. III. Nucleation in a Two-Component Incompressible Fluid. J. Chem. Phys. 1959, 31, 688. (5) Langer, J. S. Spinodal Decomposition. In Systems Far from Equilibrium; Garrido, L., Ed.; Lecture Notes on Physics No. 132; SpringerVerlag: Berlin, 1980. (6) Burks, A. W. Essays on Cellular Automata; University of Illinois Press: Urbana, IL, 1970. (7) Rothman, D. H.; Keller, J. M. Imiscible Cellular-Automaton Fluids. J. Stat. Phys. 1988, 52, 1119. (8) Rothman, D. H.; Zaleski, S. Lattice-Gas Models of Phase Separation: Interfaces, Phase Transitions and Multiphase Flows. ReV. Mod. Phys. 1994, 66, 1417. (9) Lifshitz, E. M.; Pitaevskii, L. P. Physical Kinetics; Pergamon Press: New York, 1984; Chapter 12. (10) Chou, Y. C.; Goldburg, W. I. Phase Separation and Coalescence in Critically Quenched Isobutyric-Acid-Water and 2-6-Lutidine-Water Mixtures. Phys. ReV. A 1979, 20, 2105 and references therein. (11) Wong, N. C.; Knobler, C. Light-Scattering Studies of Phase Separation in Isobutyric Acid + Water Mixtures: Hydrodynamic Effects. Phys. ReV. A 1981, 24, 3205 and references therein. (12) Guenoun, P.; Gastaud, R.; Perrot, F.; Beysens, D. Spinodal Decomposition Patterns in an Isodensity Critical Binary Fluid: DirectVisualization and Light-Scattering Analyses. Phys. ReV. A 1987, 36, 4876. (13) Cumming, A.; Wiltzius, P.; Bates, F. S.; Rosedale, J. H. LightScattering Experiments on Phase Separation Dynamics in Binary Fluid Mixtures. Phys. ReV. A 1992, 45, 885. (14) White, W. R.; Wiltzius, P. Real Space Measurement of Structure in Phase Separating Binary Fluid Mixtures. Phys. ReV. Lett. 1995, 75, 3012. (15) Califano, F.; Mauri, R. Drop Size Evolution During the Phase Separation of Liquid Mixtures. Ind. Eng. Chem. Res. 2004, 43, 349. (16) Tanaka, H.; Araki, T. Spontaneous Double Phase Separation Induced by Rapid Hydrodynamic Coarsening in Two-Dimensional Fluid Mixtures. Phys. ReV. Lett. 1998, 81, 389. (17) Vladimirova, N.; Malagoli, A.; Mauri, R. Diffusion-Driven Phase Separation of Deeply Quenched Mixtures. Phys. ReV. E 1998, 58, 7691 and references therein. (18) Le Bellac, M. Quantum and Statistical Field Theory; Clarendon Press: Oxford, U.K., 1991; Chapter 2. (19) Alexander, F. J.; Edrei, I.; Garrido, P. L.; Lebowitz, J. L. Phase Transitions in a Probabilistic Cellular Automaton: Growth Kinetics and Critical Properties. J. Stat. Phys. 1992, 68, 497. (20) Binder, P. M. Parametric Ordering of Complex Systems. Phys. ReV. E 1994, 49, 2023. (21) Binder, P. M. Domains and Synchronization in High-Dimensional Cellular Automata. Phys. ReV. E 1995, 51, R839. (22) Bruschi, M.; Santini, M. Cellular Automata in 1+1, 2+1 and 3+1 Dimensions, Constants of Motion and Coherent Structures. Physica D 1994, 70, 185. (23) Vannozzi, C. Phase Separation by Cellular Automata. M.S. Thesis, Pisa University, Pisa, Italy, 2004.
ReceiVed for reView November 8, 2005 ReVised manuscript receiVed January 27, 2006 Accepted February 28, 2006 IE051240W
