We simulate Hele-Shaw flow of both Newtonian and non-Newtonian fluids using diffusion-limited aggregation (DLA) subject to a surface tension condition, and observe the Saffman-Taylor instability as in experiments. Our program generates a DLA by simulating a brownian motion and calculates the fractal dimension whenever a particle (“walker”) sticks to the aggregate. We then compare these dimensions to experimental values taken from images of a true thin plate experiment.

Diffusion limited aggregation(DLA) is a process where random walkers walk until they collide with and stick to an aggregate. As they walk in all directions with uniform probability, it can be shown that they find a solution to the Laplace equation. Hele-Shaw flow between two plates follows Darcy’s Law, which for incompressible liquids results in the Laplace equation. Therefore, we can use DLA to model the behavior of Hele-Shaw flow. Using this model, we hope to find a connection between DLA simulations and experimental results through use of the fractal dimension, as well as other measures.

An aggregate is stored in a two dimensional array of integers that represents the walking space of the walkers. Each walker represents a particle undergoing Brownian motion which then places an integer representing its walk number at the point in the walking space at which it sticks to the aggregate. At the beginning of a simulation, the space is initialized by placing a 2x2 seed of 1’s at its center. Each walk starts by calculating a random position around a circle of radius R + b (R = maximum radius of the aggregate, b = some constant buffer around the maximum radius). Once this position is calculated, the walker then begins to choose one of 4 directions randomly {up, down, left, right} making sure the direction it picks will leave it within the bounding circle. It then implicitly updates its location or locale checking if its projected location would be adjacent to the aggregate. If it is, the walker then calculates its probability of sticking. From this the walker can decided if it sticks to the aggregate or not. If it does decide to stick, it checks to ensure it would not cause any holes in the aggregate and places the current walk number on the space at that point, ending the walk. The maximum radius of the boundary is then re-calculated, which may increase the size of the bounding circle as well. After a successful walk, the same walker object is re-originated to a point along the bounding circle and the process is repeated until the desired number of walks is attained.

Though standard DLA works well for crystal growth and other solid formations, we have to take into account another factor when modeling a liquid: surface tension. Standard DLA creates long, thin, wispy tendrils that would collapse in on themselves under the force of surface tension if made from real liquid. Though it is clear how this force acts on normal particles, how can we model this effect using a static aggregation?

The force of surface tension at a point is correlated to the curvature at that point. As the liquid becomes more and more curved, particles move from zones of high curvature to zones of low curvature to minimize the surface area. We can model this by modifying the probability that a walker will stick to the aggregate; rather than sticking 100% of the time, the stick rate increases as the curvature does. We do so according to the following formula:

\[ P=A*(\frac{N}{L^2}-\frac{L-1}{2L})+B\]

where N is the number of stuck walkers within a box of side length L centered at the point the walker would stick, and A and B are tuning parameters of the model. The term within the parentheses is a discrete version of the curvature; as it increases, so does the probability of sticking.

A (shear thinning) non-Newtonian fluid differs from a Newtonian one in that as the fluid flows faster, the viscosity decreases, and with it, the resistance against an ever increasing flow rate. In our model, we incorporate this as a probability adjustment which is added to the probability obtained from the surface tension model. The formula is:

Within the simulation, however, we do not have an exact velocity, as we instead have a state based aggregate that changes from step to step. Instead, we calculate an artificial velocity, which is simply the inverse of the age of adjacent stuck walkers. We find both components and then find the magnitude, which we use as V. This is correlated with the velocity, as a faster moving surface will have relatively young walkers at its edge which are replaced quickly. This means that limb formations with larger artificial velocities will grow more rapidly than older branches and influence the overall direction of growth of the aggregate. Using the Newtonian simulation as a basis we derived the average growth of the aggregate over 106 walks:

To prevent unrealistic holes in the aggregate, before sticking, each walker checks to see if it would make a hole. This is done by examining the connectedness of the region around the potential sticking point. If the region was previously disconnected, connecting it will create a hole, as doing so will disconnect the complement.

Implementation wise, this is done by finding the number of flips when checking around the point clockwise, where a flip is a transition from occupied to empty or vice versa, and rejecting if the number is greater than 2.

To calculate the fractal dimension of the aggregate, we use a modified version of the standard box counting method. Box counting measures the count of filled or partially filled boxes using progressively smaller grids, and then finds a power law correlation using a log-log plot, with the slope of that plot (or rather of the least squares fit line) being the fractal dimension. Normally, box counting begins with the coarsest possible grid i.e. a single square containing the entire simulation. However, we found that this gives biased results as the aggregate changes in size, so instead, we only fit the line using the 5 data points with the finest grid. The choice of 5 is arbitrary here; what is important is that it stays consistent between all measurements.

Running large scale newtonian simulations for aggregates of 106 yield an apparent fractal dimension of 1.88. Performing a linear fit on the resulting data shows a very small correction to this number.

The average newtonian fractal dimension of 5 different randomly seeded simulations using A = 2 and B = .15 was 1.881. Using the average aggregate boundary speed at 10000 and 99000 walks to calculate Cnn = 371.6 and α = 1.379, we then found the average fractal dimension of 5 randomly seeded non-newtonian simulations to be 1.717.

Newtonian Non-Newtonian 1.880328 1.721643 1.883453 1.713094 1.881879 1.713733 1.880099 1.718442 1.877978 1.719355 Avg Avg 1.880747 1.717253

Next, using the slightly changed values of A and B we calculated Cnn = 1186.2 and α = 1.859.

Newtonian Non-Newtonian 1.864648 1.7725 1.863674 1.7837 1.860703 1.7743 1.855549 1.7732 1.856701 1.7751 Avg Avg 1.860255 1.7758

We can see a higher non-newtonian fractal dimension and a lower newtonian fractal dimension. This means that the newtonian and non-newtonian form a more “similar” image with the changed variables.

Running the non-newtonian version with the same base parameters and choosing newtonian tuning variables Cnn = .85 and α = .233 yields a fractal dimension around 1.82. The non-newtonian simulation shows a small positive correction when held to a linear best fit.

By varying the tuning variables we can see a fairly large variation in the fractal dimension, as can be seen in the following graphs of aggregate particle number (Walk #) vs. fractal dimension:

The following images represent the Newtonian and Non-Newtonian Diffusion Limited Aggregation with 106 particles each. It is easily apparent that the more viscous Non-Newtonian model has thin spiny fingers when compared to its Newtonian counterpart. The change in speed on the larger limbs ensures that long branches grow longer more quickly.

For the purpose of data analysis, we needed frame-by-frame images of the DLA at each walk number. We decided that instead of creating and storing m png image for each of n walks, we would create one png image from which we could later derive each individual walk. We did this by using the integer of the walk at a given point on the aggregate as the color of a pixel on the image. This allows each pixel to have a unique color ranging from #000000 to #FFFFFF (in hexadecimal) or 16777215 - 1 in decimal (the -1 is to allow for background color). From this an image processor can simply shear off any colors above or below a certain desired walk number to get the image of the DLA at that point in the simulation. Example of an encoded Newtonian image:

Newtonian

Non-Newtonian

Run it yourself - Code for the simulations is available at: https://bitbucket.org/NJHale/math451h/src/60babc48cb2388a4247145cce0598343dd582976?at=java-dev

The simulation was done in Java, using the following class structure:

To store the aggregate as well as the times at which each walker stuck to it, we create a Space, which stores this information in 2D arrays of integers. This class also contains a method for rendering and displaying the space as an image.

We then create a DLAUtil class, which contains all of the methods required to find the fractal dimension of a provided Space, including box counting and linear regression.

Next we create a RadialWalker, which includes logic for both Newtonian and non-Newtonian walkers that spawn radially around a seed. In order to avoid memory overhead, we only create one RadialWalker, and reset it to a new spawn location whenever we would create another.

After, the DLASimulation runs a simulation using a particular Walker and Space and the number of walks to run. Depending on parameters, it also writes the simulation to image files every thousand walks, and displays the simulation on screen in real time. This class also contains the method for calculating the average speed of the boundary, which was used to find appropriate values of Cnn and alpha.

Finally, the Main class accepts command line arguments and runs the simulation according to the specified parameters, initializing an ArraySpace and a RadialWalker and feeding them into DLASimulation. The parameters include the size of the Space, the number of walks to perform, Newtonian and non-Newtonian probability formula constants, and seed values for randomization.