Chaos, Printed
So/irons & Frocralr in Great Britain
Vol
2, No. 5, pp.S31-541.
1992 0
096%0779/92$5.00 + .oO 1992 Pergamon Press Ltd
Growth Pattern Formation Driven by External Fields* GUILLERMO Consejo National
MARSHALL
de Investigaciones
EDUARDO Comision
National
de Energia
Argentina
ARGUIJO
Atomica,
(Received
Cientificas,
12 June
1429 Buenos
Aires,
Argentina
1992)
Abstract-The effects of external flows on diffusion-limited aggregation (DLA) is investigated by a Growth Pattern Formation (GPF) computational model. The GPF introduced consists of a DLA process driven by an external free convection field in a plane viscous incompressible flow. The flow regime is governed by dimensionless parameters which determine the development of different morphologies. It is shown that the GPF introduced is a generalization of previous DLA models and can be used for the simulation of more realistic physical configurations.
INTRODUCTION
Growth pattern formation (GPF), that is, the unstable growth of interfaces or its stable erosion is a common phenomenon in a wide range of problems from physics to biology. They produce complex geometries of fractal or dendritic character and chaotic patterns, see, for instance, Kadanoff [l] and Mandelbrot [2]. A great variety of GPF models driven by linear fields have been proposed for their study, notably the Laplacian Growth Models (LGM) consisting mainly in the solution of the Laplace equation with a random moving interface. The basic LGM is the Diffusion Limited Aggregation (DLA) model introduced by Witten and Sander [3], the Dielectric Breakdown Model (DBM) presented by Pietronero and Wiesmann [4] and the Diffusion Limited Erosion (DLE) or anti-DLA introduced by Patterson [5]. In many real systems the growth mechanism is perturbed by the presence of external fields such as those coming from the inclusion of heat and mass transfer effects. Although many authors are aware of the importance of the role played by the external driving field, apparently, very little is known of its effects on GPF. A pioneering step towards more realistic GPF was accomplished by Meakin [6] with the introduction of a biased random walk (BRW) model to study the effects of particle drift in DLA with a constant external field. There was posed the fundamental question of the relationship between the Hausdorff dimension of the cluster and the Hausdorff dimension of the particle trajectories; it was basically found that, for constant external fields, increasing convection the Hausdorff dimension (HD) of the plane trajectories passed from 2.0 to 1.0 while the corresponding clusters passed from the DLA fractal to 2.0. The effects *Communicated
by Professor
T. Vicsek. 531
G. MARSHALLand E. ARGUIJO
532
of a constant thermal gradient using biased diffusion was also presented by Vicsek [7]. More recently Lam [8] presented extensions of BRW models with notable improvement over previous models. It is not obvious from these BRW models what sort of differential problem was being solved. This subject was treated recently, by Nagatani and Usami [9], with a DBM for electrochemical deposition with nonlinear diffusion and by Marshall [lo] with a DLA driven by a forced convection. With the aim of studying more realistic physical configurations and as a natural extension of previous BRW models, we introduce here following previous work [lo], a DLA model driven by a free convection flow. The long term objective is the study of electrochemical deposition which appears to be a paradigm of GPF driven by external fields. Here, we show that the GPF model includes previous BRW models as particular cases and we study the relationship between the HD of the particle trajectory and that of the cluster as a function of the dimensionless parameters for rain models and rotational flows.
THE BRW AND THE ECDLA MODELS
The basic assumption of the BRW model introduced by Meakin [6] is that, if the external field is constant in space and time and the particles are in a dissipative medium, their trajectories may be described as a random walk with a superimposed drift. In essence, the BRW model consists of the following characteristics: a particle is convected with probability DP (drift probability), while the particle is diffused with probability LY= 1 - DP. More precisely, given a two-dimensional lattice (with unit space step), with P,, a central node and P,, P,, P3 and P4 its four neighbours, east, north, west and south. respectively, and a constant velocity field (north-south), with probability DP a particle moves from its initial position at PO to P,, and with probability LY,it moves to any of its four neighbours. The probability transition can thus be written i=4
where CG+’ is the probability that a particle will be at the position PO at the transition level n + 1, p, = p2 = p3 = p = a/4 and p4 = a/4 + DP = 1 - 3p. The BRW given by expression (1) may come from the discretization of a stochastic differential equation. The question is which one; the difficulty here is that the BRW method consists of a succession of two-steps that are mutually exclusive: convection or diffusion. To answer that question it is easier to pose, first, a differential equation and its discretization and then, to prove that it is identical to the BRW scheme. Hence, first we claim that the BRW model can be described by the two-step split operator C,+AC=O C, -
BC = 0
where C(x, y, t) stands for is a (normalized) velocity, and time coordinates and semi-discretized (time only) (C
n+l
(C”f’
_
-
O
if
a probability measure, A = a,@~, B = a2/3x2 + a2/3y2, u > 0 D is a diffusion coefficient, u + D = 1, x and t are the space 8 is a random variable with uniform distribution in [O. 11. A version of equation (2) yields AC” = 0
if
C”)/k - BC” = 0
if
C”‘)/k
(2)
u=l-D<8<1,
if
+
O
u=l-D
Growth
pattern
533
formation
which can be written U{(Cn+l - C”)/k + AC”}
+ D {(Cn+’
- C”)/k
- BP}
= 0.
(3b)
Equations (3a) or (3b) prescribe a sampling of the convective and diffusive difference operators (in deterministic calculus the split is made in halves). From equation (3b) it follows that (c”+l
-
C")/k
+
uAC” - DBC” = 0.
(3c)
In the limit equation (3~) represents the stochastic differential equation C, + UC, - D(C,,
+ C,,) = 0,
thus showing that the two-step split operator (2) is equivalent (4). A fully discrete version of equation (4) yields
(4) to the integer step operator
i=4 q+l
-
C,” = -ruh(C,”
- C,“) + rD CC:
( i=l where r = k/h2 is the ratio between Equation (5) may be written as
- 4C,”
1
(5)
the time step and the square of the space step. 1=4
c;+l =
cp,c: i=o
where p. = 1 - (4rD + ruh), p1 = pz = p3 = rD and p4 = rD + ruh. Since k and h are the time and space lattice steps, n is the time level and the pi are the same as before, a Monte Carlo solution of equation (6) consists of a succession of random steps in which a particle jumps in the direction of decreasing time from Pi+’ to Py, i = 0 to 4, with probability transitions given by the corresponding pi. With the addition of DLA sticking rules equation (5) constitutes a new scheme that we called the Evolutionary Convective Diffusion Limited Aggregation (ECDLA) model. At this stage, the main difference between the ECDLA given by (6) and BRW described in (1) is the presence of a PO transition in the former, a fact that obviously influences the cluster morphology. But in (6), the PO transition can be eliminated by requiring that 4rD + ruh = 1; this condition determines the ratio r = 1/(4D + uh) and thus the time step, if h is fixed. In this case the ECDLA scheme becomes r=4 co
n+l
=
Epic; i=l
(7)
where p1 = p2 = p3 = rD and p4 = rD + ruh. Making r = 1, D = ~14 and uh = DP, the BRW scheme (1) is obtained. Hence, we proved that the ECDLA model is equivalent to the BRW model for the particular case in which u, D and r are taken as above. For the limiting cases in which D = l/4, uh = DP = 0 (pure diffusion) and D = 0, uh = DP = 1 (pure convection), both models are identical; for intermediate values of uh and D they are similar in an average sense. The equivalence between the stationary version of the ECDLA model, here called Stationary Convective Diffusion Limited Aggregation (SCDLA) and the BRW model is worth analysing. Eliminating the time variable in equation (4) and discretizing it in the usual way, yields a master equation whose solution by Monte Carlo gives
534
G. MARSHALL and E. ARGCIIJO i-4
c;;” = J$,c:
(8)
1-I
where now p , = p2 = p3 = rD, p_, = rD t ruh. r = 1,/h’ and n indicates the iteration level. Again, making r = 1. D = a/4 and uh = DP. the BRW scheme (1) is recovered. Hence, we have proven that the SCDLA model is also equivalent to the BRW model when II. D and r are taken as above. In summary, the BRW model can be associated with a particular case of an evolutionary or with a stationary convective-diffusive differential equation, respectively. The numerical experiments shown below confirm these results. Consider a square box and an aggregation process with a drift, taking place. A uniform velocity field in the 4’ direction is assumed to be an external flow and the ECDLA model given by equation (5) is solved with the Monte Carlo method. Particles are launched at random from a horizontal line inside the box (rain model). Those hitting the lateral or top walls are reflected while the ones hitting the bottom (or the cluster) stick with probability one (the sticking rule is discussed below). Figure 1 shows the results for different values of
(a>
(b)
(d) Fig. 1. ECD1.A
clusters
for U/I = UP=
: (a) 0.5. (h) 0.1. (c) 0.05 and (cl) Il.01
Growth
pattern
535
formation
uh and D (according to Table 1) obtained in a square lattice of 256 x 256 nodes. For comparison, Fig. 2 show the results of the same problem produced with Meakin’s BRW model [6] while Fig. 3 shows the results of both, ECDLA and BRW models, for uh = DP = 1 and 0. Inspection of the corresponding runs in Figs 1 and 2 shows qualitatively their equivalency, as theoretically predicted. This equivalency is supported quantitatively
Table 1. Equivalence
N
1.00 0.99 0.95 0.90 0.75 0.50 0.00
between
uh = DP 0.00 0.01 0.05 0.10 0.25 0.50 1.00
Cc) Fig. 2. BRW clusters
the ECDLA and BRW field uh = DP D = ~~14 0.2500 0.2475 0.2375 0.2250 0.1875 0.1250 0.0000
for different
values
HD (BRW)
of the external
HD (ECDLA)
1.603 1.751 1.813 1.837 1.874 1.902 1.964
1.603 1.744 1.812 1.825 1.885 1.902 1.964
Cd) for DP = (a) 0.5, (b) 0.1, (c) 0.05 and (d) 0.01
G. MARSHALL and E. ARGUUO
536
(a) Fig. 3. ECDLA
@I and BRW clusters
for (a) DP = 1 and (b) DP = 0.
by comparing the estimated HD for both models, shown in the last two columns of Table 1. This is calculated using the box counting method with successive subdivision of II = 0 to 6) and a least squares method. overlaying grids in powers of 2 (h = l/2”, Moreover, the relationship between the HD of the cluster and that of the particle trajectories for increasing convection follow identical trends in both models, as predicted by the theory. The above figures and table also show that the clusters and its associated HD vary substantially in the range uh between 0.05 and 0. This is simply because in that range, for h = l/256 only, convection and diffusion are of the same order. In the classical version of DLA the particle sticks when it hits a nearest neighbour site of the cluster; this variant is called the site-version [ll], [12]. In the bond-version; the particle sticks at the previously visited site, when it hits a site of the cluster. These variants together with absorbing or reflecting boundary conditions influence the cluster and its HD. Figure 3 shows BRW clusters for DP = 1 and 0 obtained with absorbing boundaries (particles hitting lateral and top walls are lost) and with the site-version. Comparison with the corresponding values in Fig. 3 (obtained with reflecting boundaries and with the bondversion) show that the clusters obtained differ. Moreover, for DP = 1, the site-version does have no influence here). Fol not allow for compact clusters with HD = 2 (boundaries
(a) Fig, 4. BRW clusters
for (a) DP = I and (b) DP = 0.05 with absorbing
walls and the site version
Growth pattern formation
537
DP = 0, the site- or bond-versions lose relevance in relation to absorbing or reflecting boundary conditions. To end, we note that the construction of clusters in which the PO transition is allowed, introduces a new dimension into the problem that is analysed elsewhere.
THE GPF MODEL The GPF model introduced here consists in an aggregation process driven by an external Benard flow and is a natural extension of the ECDLA discussed above. The problem is described by the Navier-Stokes equations (9), (10) and the energy equation (11) coupled to a conservation equation (12) for the passive scalar. The governing equations written in transport-vorticity form are au/at
+ u&+x
+ u&&y
= V2w + GraT/ax
V%$ = --Lo
(9) (IO)
aT/at
+ cdT/EJx
+ uaT/EJy = l/PrV’T
(11)
X/at
+ uX/i3x
+ uXZ/3y
(12)
= l/SchV’C
where w is the vorticity, C$is the stream function, u and u are the velocities: u = a$/ay, u = --8$/3x, T is the temperature, C is the concentration of the passive scalar, x, y and t are the space and time coordinates, respectively, Gr , Pr and Sch are the Grashoff, Prandtl and Schmidt dimensionless parameters (for detail see any standard textbook). The computational model assumes an initial configuration and solves (9), (10) and (11) in a lattice using finite differences and deterministic relaxation techniques. The results are then used in the ECDLA model. Assuming for simplicity a steady state flow with known (previously computed) velocity and thermal fields, the sole equation remaining is (12). This can be interpreted also as a stochastic differential equation where the concentration stands for a probability measure. Its solution is obtained in the above-mentioned lattice with a discrete approximation of (12) that can be written as n+l
C,
=
Z,a,Cl
(13)
where J’ represents the nearest-neighbour site of the site i, n is the time level and the summation ranges over all the nearest-neighbour sites. The coefficients a, contain the velocity at the site i together with the dimensionless parameters as well as the space and time step. The master equation (13) is solved in the square lattice with a Monte Carlo method introduced by Marshall [13]. In this process a newly aggregated particle changes the boundary locally and hence the solutions of (9)-(11) that must be recalculated, at every step. This turns out to be intractable in terms of available computer facilities. In the present work a drastic simplification is made by uncoupling (9)-(11) from the ECDLA process, that is, the thermal and velocity fields are obtained first, in a fixed domain, the results then used, unchanged, for the aggregation process. Preliminary numerical results of a parametric study for a rotational flow are presented next. Figure 5 shows the stream function, vorticity and temperatures of the steady state driving field for Gr = 3000 and Pr = 10 in a square plane cavity, obtained with a SOR iterative method and a square lattice with 128 x 128 cells. The boundary values for the temperature are a cold vertical left side at -1 and a hot vertical right side at 1 with a linear variation for the top and bottom sides. Initially velocities and temperature are zero in the domain, with non-slip conditions for the velocities at the boundaries. Figure 6 (left) show, for different SC/Z and constant Pr and Gr numbers, the trajectory of
538
G.
MARSHALL
and E. AKGUIJO
(a)
(b)
Fig. 5. (a) Stream
function.
(b) vorticity
function
and (c) temperature
for C;R = 3000 and P’r = IO
Growth
pattern
539
formation
-1
I
I
!
(a)
I (b)
I: =: =: :: :::: :
I
(d)
Fig. 6. Particle
trajectories
(left)
and GPF
(right)
for Gr = 3000, (d) 100000.
Pr = 10 and Sch =
: (a) 0, (b) 10, (c) 100 and
G. MARSHALL and E. ARC;UIJO
540
Table 2
3000 3000 3000 3000 3000 3000 3000
IO 10 10 10 10 10 IO
0 1 10 100 1000 10 000 100 000
I.590 1.552 1.440 1.1Y5 1.139 1.110 I.070
one particle, launched from the square centre and ending at one of its sides. Increasing the Sch number with Gr fixed is roughly equivalent to increasing Gr while keeping Sdr constant. Increasing Gr (or SC/I), convection is enhanced: the trajectories gradually lose their stochastic motion and become primarily deterministic. However, there always remains a random component resulting in a particle trajectory gradually spiralling away. with back and forth short cycles, until eventually it hits a boundary. Figure 6 (right) shows, for different Sch numbers. the corresponding DLA results. Particles are launched from the centre of the square and stick to the boundaries. The cluster grows until eventually it reaches the particle source. The shape of the clusters is explained by the corresponding trajectories of Fig. 6 (left). Table 2 shows the estimated Hausdorff dimension of the clusters of Fig. 6 (right) as a function of the adimensional numbers. The results indicate that for increasing convection, while the Hausdorff dimension of the trajectories goes from 2.0 to 1.(I, the Hausdorff dimension of the corresponding clusters goes from the DLA fractal to 1.0.
CONCLUSIONS A computational model for the study of the effects of thermal convection in DLA flows has been introduced. It is shown analitically and numerically that the present model includes previous BRW models as particular cases. From a more general point of view. our results show that it is possible to reproduce a wide spectra of GPF and to measure their Hausdorff dimension by playing with the ‘degree’ of stochasticity of the master equations. This, undoubtedly, will permit the simulation of more realistic aggregation processes in general, and in electrochemistry in particular. Arlinowle12Rerne~Ir.v~Part of this work was carried out while the first author was visiting the Department of Bariloche, Argentina. The support of Theoretical Physics at the University of Geneve and the Balseiro Institute. these institutions is gratefully acknowledged. He would also like to thank D. Zanette from the Balseiro Institute for enlightening disc&sions
REFERENCES 1. 2. 3. 4. 5.
L. P. Kadanoff. Simulating hydrodynamics: a pedestrian model, J. .Stut. Phy.r. 39, 267 (lY85) B. B. Mandelbrot. T& Fracral Geomerry of Narurr. Freeman. San Francisco (1982). T. A. Witten and L. M. Sander. Diffusion-limited aggregation. Phy.~. Ret,. B27. 5686 (IY83). L. Pietronero and H. J. Wiesmann. Stochastic model for dielectric breakdown, J. Star. Phy~. 36, UOY(IYW). L. Patterson. Diffusion-limited aggregation and two-tluid displacements in porous media. Pllvv. RCT. Lerr. 52. I621 (1084). P/r>,\. Rcr,. B28. 5221 (lYX3). 6. P. Meakin, Effects of particle drift on diffusion-limited aggregation, 7. T. Vicsek, Pattern formation in diffusion-limited aggregation, Ph!,.r. He\.. Larr. 53, 2281 (1084). 8. L. Lam, R. D. Pochy and V. M. Castillo. Pattern formation in electrodeposits. in L. Lam and H. c‘. Morrlh (eds). ,Yor~lirzectr Stnrcnrres i/l Physical .S~Y~C/II.~. Springer, New York ( IYYO).
Growth
pattern
formation
541
9. T. Nagatani and Y. Usami, Pattern formation in nonlinear diffusion-limited aggregation, Phys. Rev. A39, 2169 (1989). 10. G. Marshall, Growth patterns driven by a nonlinear field. Tenth Int. Conf. on Nonlineur Science, Los Alamos, May 12-25, 1990. 11. C. J. G. Evertsz, Laplacian Fractals, Ph.D. Thesis, University of Groningen (1988). 12. T. Vicsek, Fracful Growth Phenomena. World Scientific, Singapore (1989). 13. G. Marshall, Monte Carlo methods for the solution of nonlinear partial differential equations, Comput. Phys. Comm. 56, 51 (1989).