Diffusion-limited aggregation in laminar flows

International Journal of Modern Physics C Vol. 14, No. 9 (2003) 1171–1182 c World Scientific Publishing Company DIFFUSION-LIMITED AGGREGATION IN LAMINAR FLOWS R. M. H. MERKS∗ , A. G. HOEKSTRA, J. A. KAANDORP, and P. M. A. SLOOT Computational Science Section, University of Amsterdam Kruislaan 403, 1098 SJ Amsterdam, The Netherlands ∗ [email protected] Received 5 March 2003 Revised 16 April 2003 In the diffusion-limited aggregation (DLA) model, pioneered by Witten and Sander (Phys. Rev. Lett. 47, 1400 (1981)), diffusing particles irreversibly attach to a growing cluster which is initiated with a single solid seed. This process generates clusters with a branched morphology. Advection–diffusion-limited aggregation (ADLA) is a straightforward extension to this model, where the transport of the aggregating particles not only depends on diffusion, but also on a fluid flow. The authors studying two-dimensional and three-dimensional ADLA in laminar flows reported that clusters grow preferentially against the flow direction. The internal structure of the clusters was mostly reported to remain unaffected, except by Kaandorp et al. (Phys. Rev. Lett. 77, 2328 (1996)) who found compact clusters “as the flow becomes more important”. In the present paper we present three-dimensional simulations of ADLA. We did not find significant effects of low Reynolds-number advection on the cluster structure. The contradicting results by Kaandorp et al. (1996) were recovered only when the relaxation into equilibrium of the advection–diffusion field was too slow, in combination with the synchronous addition of multiple particles. Keywords: Advection–diffusion-limited aggregation; diffusion-limited aggregation; lattice Boltzmann; moment propagation. 1. Introduction Since the pioneering paper by Witten and Sander,1 diffusion-limited aggregation (DLA) has been of continuous and extensive interest in studies using experimental,2,3 theoretical,4 – 6 and computational1,7 approaches. In computational models of DLA an initial solid seed is placed in the middle of a domain. A particle is then released from a random position at some distance from the seed. The particle carries out a random walk until it hits the seed after which it is added to the developing cluster. This procedure, when iterated, produces aggregates “distinguished by their wispy appearance”.1 More recently, a number of authors8 – 12 have studied ∗ Present address: Biocomplexity Institute, Department of Physics, Swain Hall West 025, 727 East Third Street, Bloomington, IN 47404-7105, USA. 1171 1172 R. M. H. Merks et al. this process under the influence of flow, where Warren et al.10 studied “growth by interception” with a growth probability determined by the fluid velocity. Such advection–diffusion-limited aggregation has been studied not only for understanding abiotic growth processes such as the growth of sedimenting clusters and the crystalization of (falling) snowflakes (see Warren et al.10 and references therein), but also in the context of coral growth modeling.11,13 In models of advection–diffusion-limited aggregation the aggregating particles are transported both by diffusion and by a fluid flow which interacts with the growing aggregate. Such systems should be distinguished from models of ballistic deposition,14 – 17 where particles move straightly or according to a biased random walk and whose governing direction is unaffected by the growing aggregate. In most of these studies the only effect of fluid flow was a preferential growth against the flow. The structure of the clusters as expressed by the fractal dimension was hardly affected, because of the “screening of streamlines from the interior of (. . .) the clusters.”10 These results contradict an experimental study on electrodeposition18 and a simulation study by Kaandorp et al.,11 where it was found that under the influence of a governing flow the aggregates became more compact and the “wispy appearance” was suppressed. To investigate the differences between the work of Kaandorp et al.11 and the other models of advection–diffusion-limited aggregation, in this paper we reinvestigate the effect of low Reynolds-number flows (Re ≈ 0.06–0.3) on advection–diffusion-limited aggregation. Also, we aim to understand why Kaandorp et al.11 found an effect of flow on the fractal structure of the growing aggregates, whereas other authors did not. 2. Model and Simulation Methods The model that we present here is based on the aggregation model by Kaandorp et al.,11 which was originally developed to study the effect of flow on coral morphology. This model is a Meakin19 growth model with added advection. In these models, all possible biased random paths of the particles are ensemble averaged and approximated by a continuous advection–diffusion process where the growing cluster is a particle sink. The aggregation probabilities are then given by the local fluxes into the cluster. The simulation set-up is shown schematically in Fig. 1. We started with an initial solid seed on a solid ground plane. A flow field was calculated using the lattice Boltzmann–Bhatnagar Gross Krook method (LBGK).20 This method is wellsuited for problems of complex geometry21,22 such as DLA clusters. It is resolved on a structured lattice, where we use 18 velocities and a zero velocity (D 3 Q19 ). A no-slip (u = 0) boundary condition was applied at the aggregate and at the solid ground plane. We also applied a no-slip boundary condition at the top plane, following Kovács.12 Hence the clusters were grown in a channel flow. The left, right, front and rear boundaries were periodic. Note that we used periodic boundaries for the flow calculations. Diffusion-Limited Aggregation in Laminar Flows (a) Initial condition (b) Flow 1173 (c) Dispersion (d) Aggregation Fig. 1. Simulation set up. (a) Initial condition, (b) calculation of flow field, (c) advection–diffusion of the growth resource, where the top plane is a source of resource and the ground plane and the solid particles are resource sinks and (d) aggregation, the addition of particles. The sequence (b), (c) and (d) is called a “growth cycle”. All calculations were carried out in three dimensions. When the flow field had sufficiently approached stability (full convergence model, FC) or after a fixed number of time steps (partial convergence model, PC), the dispersion of particles through the fluid was simulated by numerically solving the advection–diffusion equation. The advection–diffusion equation was solved using the moment propagation method.23 – 25 This method redistributes a particle concentration R(x, t), biased by the stable flow field as calculated using the LBGK equations, X  (fi − ∆∗ f eq (u = 0, ρ))R  i + ∆∗ R(x, t) , (1) R(x, t + 1) = ρ x−ci ,t i where the whole quantity inside [. . .] is evaluated at (x − ci , t), fi is the probability that a carrier fluid particle moves into the lattice direction ci and fieq (u = 0, ρ) is the probability that a particle moves along ci in a nonmoving fluid. The diffusion coefficient D is set using the parameter ∆∗ as D= 1 1 ∗ − ∆ 6 6 (2) for the D3 Q19 model we used.25 The ratio between advective and diffusive transport is expressed by the Péclet number Pé = ūL/D, where ū is the mean velocity and L is a characteristic length. Throughout this paper we exclusively use the lattice Péclet number Pélat , in which L = 1 lattice unit, the distance between two nodes of the lattice. Following Kaandorp et al.,11 the top plane was a source of growth resource and both the ground plane and the cluster were resource sinks. We have followed the Meakin19 convention by putting the resource sink at the developing cluster. For the advection–diffusion calculations, the front and rear boundaries were periodic, whereas at the left and right boundaries (i.e., in the flow direction) a no-flux 1174 R. M. H. Merks et al. Table 1. The amount of new particles added to the cluster per growth step in the multi-particle model depending on the amount of available new positions on the cluster (Kaandorp, personal communication). Number of possible positions (n) Number of particles added n < 100 100 ≤ n ≤ 1000 1000 < n ≤ 10000 n > 10000 1 10 50 100 (reflecting) boundary was applied, thus preventing effects of the periodic images of the growth forms. In the Full Convergence model, the moment propagation equation was iterated until the influx of resources at the top plane balanced to certain extent of the outflux of resource at the cluster and the ground plane, i.e., when the change of the total resource mass in the system fell below a threshold, P ∆( x R) P < θAD , ∆t x R (3) where typically θAD = 10−6 . In the Partial Convergence model, the moment propagation equation was iterated for a fixed number of time steps; in this paper we typically applied 50 iterations per growth cycle, thus following Kaandorp et al.11 After a resource dispersion field was obtained, solid nodes were added to the cluster (Fig. 1(d)). Here we study two variants of this aggregation. In the single particle model during each growth cycle a single solid particle is added to the cluster. This follows the Witten and Sander1 model and the original Meakin19 model. In the multiple particle model, however, the number of particles added is increased during the simulation, depending on the size of the cluster (Table 1). In order to compare our results to the results obtained previously by Kaandorp et al.,11 we take diffusion coefficients as used in their simulations. In order to decide which fluid sites should be added to the cluster, we first determined the set of “growth candidates”; this is the set of nearest (face-connected) neighbors of the cluster. In each of these candidates, the growth resource concentration R(xi ) indicates the probability that a particle is present and adheres to the cluster. To enforce the addition of a fixed number of particles, these probabilities were normalized to sum up to 1. Hence the probability P (xi ) that a new particle is added to the cluster at a face-connected neighbor was, where P R(xk ) , Paggr (xk ) = P j R(xj ) j sums over all nearest neighbors. (4) Diffusion-Limited Aggregation in Laminar Flows 1175 The ADLA clusters were characterized using the fractal dimension, which we measured using the radius of gyration Rg ,26 where v ! u 2 u1 X 1 X t (5) xk − xi Rg = n i n k with xi the coordinates of the particles in the cluster and the sums are over all the particles in the cluster. The fractal dimension can be obtained by keeping track of the radius of gyration during the growth of the cluster. For large DLA-clusters, the radius of gyration gets a power-law dependence of the number of particles N , Rg ∼ N β , where the fractal dimension Fr = 1/β. The clusters were further characterized using the compactness C, which we define as the fraction of solid material inside the convex hulla of the cluster, Vobject , (6) C≡ Vhull where Vobject is the volume of the object, and Vhull is the volume of the convex hull, which was determined using the quickhull27 algorithm.b Hypothesis testing was carried out using one-way analyses of variance (ANOVA).c The simulations were carried out on clusters of workstations, the DAS-228 and a Beowulf cluster29 on which a typical simulation presented here took around 15 h (partial convergence) to 50 h (full convergence). 3. Results In Fig. 2 we show two typical clusters developed with the Full Convergence, Single Particle model. In Fig. 2(a) the resource was exclusively transported by means of diffusion, whereas in Fig. 2(b) transport was by diffusion and by advection, where (a) (b) u Fig. 2. Typical three-dimensional clusters developed with the fully converged aggregation model. (a) Diffusion-limited aggregation. ū = 0, D = 0.166. (b) Advection–diffusion-limited aggregation ū = 0.05 (from left to right), D = 0.166; both clusters contain 4500 particles. a Think of the convex hull as the surface given by the tightest gift wrapping around the object. http://www.geom.umn.edu/software/qhull/. c See http://mathworld.wolfram.com/ANOVA.html. b See 1176 R. M. H. Merks et al. Fig. 3. The cross-sections through 3D resource fields and clusters developed with the advection– diffusion-limited aggregation model. The single particle clusters were grown for 4700 (top) and 6150 (middle) growth cycles. The multiple particle clusters were grown for 450 growth cycles. Flows (when applicable) are directed from left to right. Diffusion-Limited Aggregation in Laminar Flows 1177 Table 2. The compactness and fractal dimension. Abbreviations used are: SADLA: Single particle advection–diffusion-limited aggregation; FC: full convergence; PC: partial convergence; MADLA: Multiple particle advection–diffusion-limited aggregation; MADLA II: MADLA with growth function as in Eq. (7); D: diffusion coefficient; Pélat : lattice Péclet-number; Fr : fractal dimension. Mean and standard deviation of the compactness and fractal dimension are calculated over five simulations. Model SADLA, FC SADLA, FC SADLA, PC SADLA, PC MADLA MADLA MADLA II MADLA II Figure (s) ū D Pélat Compactness Fr 2(a), 3(a) 2(b), 3(b) 3(c) 3(d) 3(e) 3(f) 4(a) 4(b) 0.0 0.05 0.01 0.01 0.01 0.01 0.01 0.01 0.167 0.167 0.167 0.00125 0.167 0.00125 0.167 0.00125 0.0 0.3 0.06 8 0.06 8 0.06 8 0.12 ± 0.01 0.12 ± 0.01 0.12 ± 0.01 0.14 ± 0.01 0.10 ± 0.00 0.26 ± 0.00 0.18 ± 0.00 0.38 ± 0.01 2.62 ± 0.36 2.60 ± 0.21 2.61 ± 0.07 2.67 ± 0.34 2.59 ± 0.11 3.03 ± 0.07 2.58 ± 0.11 2.99 ± 0.09 the mean flow velocity was 0.05 in lattice units. The diffusion coefficient was set to 0.166, giving lattice Péclet numbers of 0 and 0.30, respectively. Since the internal structure of these three-dimensional pictures is difficult to interpret, we only show cross-sections in the remainder of the paper. In Fig. 3 we present an overview of the results of the advection–diffusion-limited aggregation models. In these figures the grey scale indicates the particle concentration, and the flow direction is always along the x-axis. In Figs. 3(a) and 3(b) we show the clusters of Fig. 2 in cross-section. The applied flow did not affect the compactness nor the fractal dimension (see Table 2). We also tried to increase the Péclet-number further, by using a smaller diffusion constant. However, this led to extremely high computation times (150 growth cycles with D = 0.083 took four days of computer time on 32 processors of the DAS-2), which made the model computationally intractable. In order to keep the simulation times feasible, we compared the results of the full convergence model to a model where a fixed amount of moment propagation iterations was applied per aggregation step. In Figs. 3(c) and 3(d) we show two results of this model, in which we applied 50 moment propagation iterations and 10 LBGK iterations per aggregation step. In both simulations, the mean flow velocity was set to ū = 0.01. For these parameter settings, the results of the partial convergence model did not differ from the full convergence model: neither the compactness nor the fractal dimension differed significantly. When we increased the Péclet-number, by decreasing the diffusion coefficient to 0.001, we found a “wake” behind the cluster, and the cluster tended to grow into the direction of the flow. The internal structure of the clusters was however hardly affected; the compactness was found to be slightly, though significantly different (Table 2), but the fractal dimension showed no significant difference. In Figs. 3(e) and 3(f) we present results of the multi-particle model. Both simulations were carried out under the influence of a flow of velocity ū = 0.01. When we 1178 R. M. H. Merks et al. used a large diffusion coefficient (D = 0.167) we found open, DLA-like clusters (see Fig. 3(e)). When the diffusion coefficient was decreased (Fig. 3(f)), the clusters became more dense (Table 2). However, when we applied 250 rather than 50 moment propagation per growth cycle, less dense clusters developed (with C = 0.17 ± 0.0), indicating that the advection–diffusion fields had not equilibrated after 50 iterations (data not shown). This suggests that the amount of diffusion (Dt), rather than the Péclet-number, is an important factor in determining the compactness of the cluster. The clusters that developed in the multiple particle model did not become as dense as the clusters reported previously.11,13 In fact, Kaandorp et al.11 used an alternative model for the aggregation probability (Kaandorp, personal communication). In this model, the aggregation probability is the normalized probability that a particle collides with a solid–fluid interface after propagation, P R(xk , t) i fi (xk , t)δ(xk + ci , t) P , Paggr (xk ) = P i fi (xj , t)δ(xj + ci , t) j R(xj , t) (7) P where i sums over the 19 velocities of the lattice and where δ(xi , t) = 1 when xi is a solid node and δ(xi , t) = 0 when it is fluid. Note that this aggregation function favors the addition of nodes enclosed by several solid nodes. The clusters shown in Fig. 4 were developed using this growth function. With a high diffusion coefficient of D = 0.17, these clusters became a bit more dense (Fig. 4(a); Table 2) than the clusters that were developed with the standard growth function (Fig. 3(e)). The clusters were significantly more dense than with the standard growth function when a smaller diffusion coefficient was used (Fig. 4(b)). Note that the clusters are somewhat more dense away from the flow. This effect becomes stronger if the ground plane is removed from the simulation, resulting in a higher velocity around the cluster (data not shown). Multiple particle Partial convergence −3 u=0.01, D=1.25 x 10 u=0.01, D=0.17 Fig. 4. The advection–diffusion-limited aggregation model with the alternative aggregation function as used by Kaandorp et al.11 Flow is directed from left to right, ū = 0.01. Diffusion-Limited Aggregation in Laminar Flows 1179 4. Discussion In this paper we studied three-dimensional advection–diffusion-limited aggregation (ADLA) under the influence of a laminar flow. In this straightforward extension of the diffusion–limited aggregation (DLA) model by Witten and Sander,1 particles are transported by a diffusive process which is biased by a governing flow. We studied two variants of this model, the single particle model and the multiple particle model. In the single particle model one particle is added in a growth cycle. In the multiple particle model, during a growth cycle one or several particles are added, depending on the current size of the cluster. Apart from several complications, which we will discuss below, this model could be seen as a hybrid between the diffusion-limited aggregation model and models of Laplacian growth, where growth occurs over the full surface of the cluster in parallel. We studied two variants of the single particle ADLA model. These variants differ in the way the flow and dispersion fields are relaxed. In the first variant, we relax the advection–diffusion field until the change per unit time of the total mass falls below a threshold (Eq. (3)). Although this full convergence model would be preferred when studying advection–diffusion-limited aggregation in a Meakin model, these simulations quickly became computationally intractable when small diffusion coefficients were used. After adding a new particle to the cluster, it often took hundreds to thousands of time steps before the resource field had sufficiently stabilized. Note that this indicates a major caveat of our simulation methods, because it shows the assumption that the flow and resource fields are in equilibrium cannot be always met. Being aware of this caveat, we also studied an alternative model in which the tracer field was relaxed for a fixed number of time steps (50 per growth cycle). In both cases, the cluster tended to grow slightly into the direction of the governing flow, but the overall structure of the clusters was unaffected by flow, as expressed by the fractal dimension and the compactness. When we used a higher Péclet-number by decreasing the diffusion coefficient (Fig. 3(d)), the cluster became slightly more compact. However, this small effect is most likely unphysical, because we use a fixed advection–diffusion time and because the used Péclet-number is well over the upper-limit above which the moment propagation method no longer produces physical results.25 To keep the computations tractable, we only studied the partial convergence variant of the multi-particle model, which is more closely related to the model studied previously in our group.11 For small Péclet-numbers, the resulting clusters were indistinguishable from the clusters developing with the single particle model, based on the fractal dimension and on the compactness of the clusters. However, when we decreased the diffusion coefficient, keeping the flow velocity constant (as in Kaandorp et al.11 ), the clusters compactified. In the context of accretion models, which are related to models of Laplacian growth, we have observed that growth forms become more compact when the time scales of growth and resource transport come closer together.30 In the ADLA model the growth rate may not be large 1180 R. M. H. Merks et al. enough to reach this region of parameter space. However, when the growth rate is increased as in the multiple particle model, this region will be reached when the resource transport is slowed down by decreasing the diffusion coefficient. In the single particle models the laminar flow did not affect the internal structures of the clusters. In these models, the highest flow velocity is found at a large distance from the clusters, whereas the flow velocity close to the clusters is very low. Thus in close vicinity of the cluster, resource transport is dominated by diffusion, resulting in the DLA-like growth dynamics. Note that the diffusion coefficients for which compact clusters were obtained were extremely low. In fact, the moment propagation method that we have used in this paper becomes invalid for such low values,25 since in that case negative values are propagated through the lattice. The probability of this becomes larger as the flow velocity increases. In our simulations the clusters grow along the flow for such incorrect parameter settings, which may be explained by extremely low or even negative tracer concentrations occurring at the front where the flow velocity is highest (data not shown). With these results, it becomes necessary to reinterpret results of Kaandorp et al. (1996, 2001), who claimed that DLA clusters become more compact under the influence of advection. Our results suggest that this compactification occurs only in the multi-particle model that compactification only results when the net diffusion per growth cycle (Dt) is small and it occurs preferentially at parameter settings where the moment propagation method no longer produces physical results. Applying the methods described in Merks et al. (2002)25 it is easily shown that in the simulations by Kaandorp et al,11 who used a 24-velocity model on a FCHClattice, the lattice Péclet-number may never exceed a value of 2. However, the Péclet-numbers up to 3 as reported by Kaandorp et al.11 were based on the mean velocity. Optimistically estimating umax = 2∗ ū, the lattice Péclet-numbers reported previously should be multiplied by a factor of two to get the maximum Pécletnumbers. Hence, we estimate that negative amounts of tracer may be transported in simulations with reported Pélat > 1.0. A further difference between our results and those presented earlier is the growth direction of the clusters. We found clusters growing mainly towards the source plane, where the clusters may slightly deviate towards the flow, or even away from the flow. Kaandorp et al.11 reported that the clusters had the strong tendency to grow towards the flow; this tendency became stronger when the diffusion coefficient was decreased. This difference may be attributed to the following. In the simulations presented here, the initial condition was a linear concentration field, the stable solution for a system with a source and a sink plane. Kaandorp et al.11 started with all concentrations set to 1.0 (personal commuication). From the analytic, time dependent solution of the diffusion equation, the time to relax from this initial condition to the stable linear solution is estimated as t ≈ L2 /D, where L is the height of the simulation box in lattice units. For the highest diffusion coefficient that they used (D = 0.25), this time already exceeds the total advection–diffusion time available in a simulation of 1000 growth cycles with 50 tracer iterations in total. Diffusion-Limited Aggregation in Laminar Flows 1181 Hence the initial shape of the tracer field affects the growing clusters throughout the time of the simulation. When the diffusion coefficient is decreased, this effect becomes even stronger. In an empty simulation box and with the lowest diffusion coefficient (D = 0.00125), at the end of the total simulation time the concentration would follow a sharp gradient from the sink plane to about one fifth of the simulation box; everywhere else the concentration would be close to 1. The absence of a full topto-bottom resource gradient removes the tendency of the clusters to grow towards the source plane. This may well explain why the growth forms by Kaandorp et al.11 bent more strongly against the flow when the diffusion coefficient was decreased. The moment propagation method that we use in this paper seems not very well-suited for the simulation of ADLA. It is computationally very intensive to reach stable flow and resource fields. Moreover, we can only reach moderate Pécletnumbers.25 Particle tracking methods31 may be better suited for this purpose. Such methods simulate off-lattice random walks biased by the flow field, and would result in a model more closely related to the Witten and Sander1 model. In such a discrete model one would need to track only the single particle to be added to the aggregate, making it computationally much less demanding than the continuous methods which were applied here and in the previous work. In conclusion, we could not confirm the observation that DLA-clusters compactify as “the flow becomes more important (Pé increases)”.11 This finding is in agreement with the previous two-dimensional studies of advection–diffusion-limited aggregation.8– 10,12 Its explanation may be given by the fact that the flow velocity is very low in close vicinity and inside of the cluster, resulting in diffusion-limited resource transport. We expect, however, that higher Reynolds-number, instable and turbulent flows may affect the aggregation process, since these may lead to recirculations inside the cluster. Such effects may well be important in the growth of marine sessile organisms, to which the ADLA-model was previously applied,11,13 and in flow-driven electrodeposition.18 Advective transport was found to affect the compactness in a model of bacterial plaques,32 caused by increased nutrient transport towards the plaque. However, in this model some biomass spreading was allowed, which distributed the extra growth over a small area of the plaque. In diffusionlimited aggregation models, attachment is irreversible and immediately affects the flow and diffusion fields which may explain why we did not find a flow effect. References 1. T. Witten, Jr. and L. Sander, Phys. Rev. Lett. 47, 1400 (1981). 2. P. P. Trigueros, J. Claret, F. Mas, and F. Sagues, J. Electroanal. Chem. 312, 219 (1991). 3. D. D. Chambliss and R. J. Wilson, J. Vac. Sci. Technol. B 9, 928 (1991). 4. M. Muthukumar, Phys. Rev. Lett. 50, 839 (1983). 5. T. C. Halsey, Phys. Today 53, 36 (2000). 6. R. C. Ball and E. Somfai, Phys. Rev. Lett. 89, 135503 (2002). 7. M. Castro, R. Cuerno, A. Sanchez, and F. Dominguez-Adame, Phys. Rev. E 62, 161 (2000). 1182 R. M. H. Merks et al. 8. J. Toussaint, J. M. Debierre, and L. Turban, Phys. Rev. Lett. 68, 2027 (1992). 9. R. Brémond and D. Jeulin, Geostatistical Simulations, eds. M. Armstrong and P. A. Dowd (Kluwer Academic Publishers, 1995), pp. 89–105. 10. P. B. Warren, R. C. Ball, and A. Boelle, Europhys. Lett. 29, 339 (1995). 11. J. A. Kaandorp, C. P. Lowe, D. Frenkel, and P. M. A. Sloot, Phys. Rev. Lett. 77, 2328 (1996). 12. T. Kovács and G. Bárdos, Physica A 247, 59 (1997). 13. J. A. Kaandorp, The Algorithmic Beauty of Seaweeds, Sponges and Corals, eds. J. A. Kaandorp and J. Kübler (Springer-Verlag, Berlin, Heidelberg, New York, 2001), The virtual laboratory, pp. 114–144. 14. M. Vold, J. Colloid Sci. 14, 168 (1959). 15. P. Meakin, P. Ramanlal, L. M. Sander, and R. C. Ball, Phys. Rev. A 34, 5091 (1986). 16. T. Nagatani and F. Sagues, Phys. Rev. A 43, 2970 (1991). 17. T. Nagatani, J. Phys. Soc. Jpn. 60, 2700 (1991). 18. L. López-Tomàs, J. Claret, and F. Sagues, Phys. Rev. Lett. 71, 4373 (1993). 19. P. Meakin, J. Theor. Biol. 118, 101 (1986). 20. Y. H. Qian, D. D’Humières, and P. Lallemand, Europhys. Lett. 17, 479 (1992). 21. S. Chen and G. D. Doolen, Annu. Rev. Fluid Mech. 30, 329 (1998). 22. S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond (Oxford University Press, Oxford, New York, 2001). 23. C. P. Lowe and D. Frenkel, Physica A 220, 251 (1995). 24. P. B. Warren, Int. J. Mod. Phys. C 8, 889 (1997). 25. R. M. H. Merks, A. G. Hoekstra, and P. M. A. Sloot, J. Comp. Phys. 183, 563 (2002). 26. E. Sander, L. M. Sander, and R. M. Ziff, Comput. Phys. 8, 420 (1994). 27. C. B. Barber, D. P. Dobkin, and H. Huhdanpaa, ACM T. Math. Software 22, 469 (1996). 28. H. E. Bal et al., ACM Operating Syst. Rev. 34, 76 (2000). 29. T. L. Sterling, D. Savarese, D. J. Becker, J. E. Dorband, U. A. Ranawake, and C. V. Packer, Proceedings of the 1995 International Conference on Parallel Processing, Vol. I: Architecture, ed. P. Banerjee (CRC Press, 1995), pp. 11–14. 30. R. M. H. Merks, A. G. Hoekstra, J. A. Kaandorp, and P. M. A. Sloot, J. Theor. Biol. 244, 153 (2003). 31. R. S. Maier, D. M. Kroll, H. T. Davis, and R. S. Bernard, Int. J. Mod. Phys. C 9, 1523 (1998). 32. C. Picioreanu, M. C. M. Van Loosdrecht, and J. J. Heijnen, Biotechnol. Bioeng. 69, 504 (2000).