On the fractal dimension of non-Newtonian Hele-Shaw flow subject to Saffman-Taylor instability

The classical Hele-Shaw (H-S) experiment [1] is simple: inject a fluid into an immiscible, more viscous fluid confined between two closely-spaced plates. Such a setup is known to be subject to the Saffman-Taylor (S-T) instability [2], resulting in complex, fractal-like pattern formation. Extensive reviews of experiments [3], and of various aspects of theoretical and computational findings [4, 5, 6, 7], are available. Figure 1(a) shows a typical example of experimental results, obtained in NJIT’s Capstone Laboratory [8]; many such examples can be found in the literature, see [3] and the references therein.

A key reason for the interest in H-S flow and related instabilities, in addition to the easy visualization, is that its mathematical description is similar to Darcy’s formulation modeling porous media flow, which is highly relevant to applications that vary from secondary oil recovery to injection molding [9] and many other natural and man-made setups. Furthermore, there is a close connection between the S-T problem and the Mullins-Sekerka instability relevant to the propagation of a solidification front in an undercooled liquid [10]. Therefore, there are multiple reasons to work towards understanding the nature of pattern formation and emergent interface morphologies.

While a significant body of experimental work has been carried out [11, 12, 13, 14, 15, 16] (some with non-Newtonian fluids), it is difficult to reach precise understanding of the morphology of the emerging patterns, for either Newtonian or non-Newtonian fluids. Such morphology is often characterized by the fractal dimension, $D_{\rm f}$, for which values in the range 1.2 - 2.0 have been reported; see, e.g., [17, 15, 3] for discussion. Most of the (relatively) recent experiments seem to converge to a value close to 1.8 [3], but the differences between results are too great to attempt to identify the influence of non-Newtonian effects in the displaced fluid. This influence is crucial to numerous applications; however, due to inherent limitations of the experimental setups, it is difficult to quantify such differences based on experiments alone [14, 3].

From the theoretical perspective, non-Newtonian free surface flows are nontrivial to model due to often complicated rheological properties. However, the H-S geometry once again comes to the rescue since, in the thin gap limit, for a significant class of non-Newtonian fluids, the most important aspects of their rheology can be reduced to a shear-thinning model involving pressure-gradient dependent viscosity [18, 19, 16, 20]. Such a simplification was carried out systematically some years back [21] showing that, even if elastic effects may be important in general, the thin gap limit reduces the problem to the Darcy-type law ${\bf u}=-b^{2}\nabla p/[12\mu(|\nabla p|^{2})]$, where ${\bf u}$ is the gap (of thickness $b$) averaged fluid velocity, $p$ is the pressure, $\nabla=(\partial_{x},~{}\partial_{y})$, with $(x,y)$ the in-plane coordinates, and $\mu(|\nabla p|^{2})$ is the pressure-gradient dependent viscosity. This simplified viscosity law $\mu(|\nabla p|^{2})$ is obtained by gap averaging the shear-rate dependent viscosity $\mu(|{\bf u}_{z}|^{2})$ (here the $z$-subscript represents a partial derivative in the $z$-direction perpendicular to the plates of the cell); see [18] for details of how to translate between the two viscosity representations. Such modified Darcy formulations can be obtained by systematic asymptotic expansions applicable to a wide class of non-Newtonian fluids, even when the Weissenberg number quantifying the relevance of elastic effects is $O(1)$. When combined with the incompressibility condition, $\nabla\cdot{\bf u}=0$, a nonlinear elliptic problem for the fluid pressure is obtained, which is significantly easier to deal with than the original non-Newtonian formulation involving elastic effects. Despite this simplification, it is still challenging and computationally intensive to solve for the complicated interfaces that ultimately develop due to the S-T instability; the problem is of free-boundary type involving in-plane curvature through the Young-Laplace boundary condition $\bar{p}=-\gamma\kappa$, where $\bar{p}$ is a suitably renormalized pressure, and $\gamma$ the surface tension; see e.g. [7] for details. Thus, since carrying out large simulations is difficult, quantifying the emerging pattern morphology and its dependence on non-Newtonian fluid behavior remains elusive.

An alternative approach is modeling using discrete methods, particularly algorithms based on diffusion-limited aggregation (DLA) [22]. DLA is a version of Monte-Carlo simulation that uses a growing aggregate of particles to simulate solutions of Laplace’s equation with free boundaries. In the context of H-S flow, DLA-type algorithms have been modified to include surface tension effects [23] and used extensively to discuss the pattern formation that emerges in the S-T instability of displaced Newtonian fluids [17]: it was found that surface tension modifies the fractal dimension $D_{\rm f}$ of the emerging patterns from the “pure” DLA value of 1.67 to larger values. Exact values of $D_{\rm f}$ are often unclear, however, and there is a significant body of work discussing the asymptotic value of $D_{\rm f}$ in the limit of large aggregate size, as well as the fractal structure of the emerging patterns [24, 25, 26]. In any case, to the best of our knowledge, whether and to what extent the non-Newtonian character of the displaced fluid influences the value of $D_{\rm f}$, either in experiments or in simulations, is uncertain.

In this paper, we formulate a DLA-type method to model the Hele-Shaw flow of a non-Newtonian fluid. The motivation for this study comes in part from experiments that we briefly report in what follows. The main modeling idea is to incorporate the fluid viscosity’s shear rate dependence via a history-dependent pattern growth. We then discuss how such modification influences the morphology of the emerging patterns, focusing for brevity on the fractal dimension of the emerging patterns only.

The rest of this paper is organized as follows. Section II presents examples of experimental results, which serve as a motivation for the development of the models and simulations that follow. Simulation methods are discussed in Sec. III, with the implementation details relegated to Appendices A - D. The main set of results is presented in Sec. IV, with additional discussion related to reproducibility and the influence of modeling parameters in the Appendix E. Section V provides the summary and further discussion.

Figure 1 shows two examples of experimental patterns obtained by injecting, in a controllable manner, water into an H-S cell initially containing more viscous fluid, that is either Newtonian (glycerol) (a) or non-Newtonian (polyethylene oxide, PEO) (b). The H-S cell consists of two plexiglass plates approximately 1/2 inch thick, placed typically at a distance 100-800 $\mu$m apart, with a 2mm wide hole in the top plate. A syringe is connected to the hole via a plastic tube and water is injected by controlled applied pressure. Both applied pressure and the spacing between the plates were varied in the experiments; the video recordings are available [8]. While the experiments were conducted carefully and consistent results were obtained, we note that they are of limited scope and intended only to illustrate visual clues suggesting that the non-Newtonian response of PEO modifies the morphology of the evolving patterns; the interested reader is directed to the extensive review by McCloud & Maher [3] and numerous references therein showing consistent findings. We note in particular that the surface tension between the considered fluids varies as well as their rheological properties, leading possibly to modification of the emerging lengthscales. We have computed the fractal dimensions, $D_{\rm f}$ of the emerging patterns (using methods discussed later in the text), finding systematically lower values for the non-Newtonian fluid. However, since computing $D_{\rm f}$ involves consideration of different scales that should span orders of magnitude, it is necessary to generate patterns that are much larger than those shown in Fig. 1. For this purpose, we resort to simulations.

The DLA algorithm that provides the basis for this work is relatively simple, particularly the on-lattice version we consider. The Brownian motion required by the DLA model is simulated by a random walk along the four possible directions of the grid. One starts with a seed particle at the center of a lattice, generates a particle (‘walker’) far away, allows the walker to perform a random walk on the lattice until it arrives at a lattice site adjacent to the seed, where it sticks (with unit probability), transforming this particular lattice site from free to occupied (in the algorithm that we use this transformation is permanent; see App. A for more details and references to alternative formulations).

To simulate the S-T instability in the H-S setup, the basic DLA algorithm was modified by Vicsek [23] via the sticking probability. The modification promotes sticking in neighborhoods that contain a large number of occupied sites, since such regions, on average, are characterized by a smaller curvature. Within this model, the sticking probability is defined by

where $n_{\cap}$ denotes the number of occupied sites in some neighborhood of linear dimension $l$ (measured in the number of grid points) centered at the considered sticking location, and $n_{\rm total}=l^{2}$ (Vicsek [23] used $l=17$ and $n_{0}=(l-1)/(2l)$; see Appendix A for further motivation). The original DLA model is recovered by setting $A=0$ and $B=1$; choosing $B\in(0,1)$ and increasing the value of $A$ models a stronger surface tension since this encourages walkers to stick in areas of high local particle density, hence, smaller local curvature.

We wish to simulate both Newtonian and non-Newtonian flows, which requires changes to the algorithm (discussed in some detail in Appendix A); here we provide only an outline. In our Newtonian simulations, we improve Vicsek’s algorithm, modifying it to avoid the formation of ‘holes’ (unoccupied sites surrounded by occupied ones), and by promoting azimuthal symmetry and minimizing grid anisotropy.

Building on this improved Newtonian algorithm, we further adapt it to model interface growth for a non-Newtonian (and in particular, shear-thinning) displaced outer fluid. To this end, we need to modify the sticking probability so as to model shear-thinning behavior. This can be done by realizing that, instead of a shear-rate-dependent viscosity as discussed in the introduction, we may instead think of a viscosity that depends on the velocity magnitude of the displaced fluid, since the two quantities (shear rate and velocity magnitude) are related (fluid velocity is obtained by gap averaging the fluid shear rate). In the context of DLA-based simulations, we can then model shear thinning behavior by making the sticking probability velocity-magnitude dependent. This is the basic idea underpinning the selected approach.

In physical experiments, the particular dependence of viscosity on shear-rate may vary widely. Here, we do not attempt to model any particular fluid; instead (with the above discussion in mind) we assume a simple velocity-dependent non-Newtonian correction to the sticking probability of the form

Here, $C$ and $s$ are two free parameters, chosen in such a way that $p_{\rm nNewt}$ reaches values comparable to $p_{\rm Newt}$ (one could think of these parameters as modeling the shear-thinning properties of the displaced fluid); while $V$ is the velocity magnitude (speed), inversely proportional to the “age” of a particle that has since become part of the aggregate (see App. A for details). This means that younger (more recently attached) particles are more likely to attract walking particles. The local age of the aggregate has two components: a “horizontal” age $T_{H}$ and a “vertical” age $T_{V}$, which are determined pointwise on each lattice site. To illustrate the algorithm, suppose the aggregate comprises $k-1$ particles. Depending on how the $k$th particle attaches, its age is determined differently. If the particle $k$ attaches horizontally, then its vertical age is inherited directly from the occupied lattice site responsible for the attachment (neighbor), and the horizontal age is defined by $T_{H}=k-T_{H}^{\rm neighbor}$. For a vertical attachment, the two ages are calculated analogously. The velocity is then given by the vector $\mathbf{v}=\left[T_{H}^{-1},T_{V}^{-1}\right]$, and the speed by $V=||\mathbf{v}||$. There are additional considerations if a walker attaches to more than one neighbor at once; see App. C for modifications of our algorithm in such cases. The probability of attachment that we use for non-Newtonian simulations includes Newtonian and shear-thinning effects through the following expression,

with $p_{\rm nNewt}$ as defined by Eq. (2).

Figure 2 shows two examples of our Newtonian (a) and shear-thinning (b) simulations. Visually, these figures suggest the formation of skinnier patterns and possibly different morphology for the shear-thinning case. Such a trend is consistent with the experimental results shown in Fig. 1, although the differences between the Newtonian and shear-thinning patterns are less pronounced in the simulations. One reason for this may be the different surface tension between the fluids used in experiments. Returning to the influence of shear-thinning on pattern characteristics in simulations, the results shown in Fig. 2 suggest that the differences are not simple to quantify, and more importantly, we find that they depend on the pattern size, as observed in previous experiments [3] and simulations [25]. Our finding, consistent with Newtonian simulations [25], is that to obtain convergence of the fractal measures quantifying the emerging patterns, millions of walkers are required (the simulation examples in Fig. 2 show $\sim 10^{5}$ walkers, two-to-three orders of magnitude less than required to obtain converged results for $D_{\rm f}$). We therefore proceed to discuss the results of much larger simulations and characterize the morphologies of the resulting aggregates by computing $D_{\rm f}$ using box-counting and correlation algorithms; see App. D for details. We note that additional measures, such as radius of gyration of the aggregates, were considered but did not provide additional insight and are not reported here for brevity.

Figure 3 shows $D_{\rm f}$ as a function of the number of walkers (measured in millions) as the parameters $A,~{}B$ (entering the definition of $p_{\rm Newt}$ in (1)) are modified. First, Figure 3(a) shows the results for $A=0$, without surface tension effects. We observe that large simulations are indeed needed to reach both convergence and an agreement between the two calculations; for fewer than $\sim$10 million walkers, the two measures produce differing results. For sufficiently many walkers, we find the value $D_{\rm f}\approx 1.67$, in agreement with values reported in the literature for on-lattice DLA simulations [4, 25, 27]. We have confirmed that our results are realization-independent; repeating a simulation with different random seed generators produces consistent $D_{\rm f}$ values (for more than 10 million walkers). Regarding the influence of the non-Newtonian correction, recall that within the DLA algorithm, the probability of sticking is unity, and therefore, the non-Newtonian contribution provided by $p_{\rm nNewt}$ is irrelevant in this case (see (3)). Assuming that (Newtonian) DLA represents well the physical S-T problem with vanishing surface tension, the (perhaps not so obvious) conclusion is that $D_{\rm f}$ is not influenced by the non-Newtonian nature of the displaced fluid in this case.

We now consider the effect of changing the surface tension parameter, $A$. Figure 3(b - d) shows the results for $A=1,~{}2,~{}4$. As $A$ increases, $D_{\rm f}$ increases for both Newtonian and non-Newtonian simulations, but differently. We identify a gap between the corresponding Newtonian and non-Newtonian fractal dimensions, $\Delta D_{\rm f}=D_{\rm f}^{\rm Newt}-D_{\rm f}^{\rm nNewt}$ and find that the size of $\Delta D_{\rm f}$ depends non-monotonously on the value of $A$, with the maximum value observed for $A=2$, see Fig. 3(c). For yet larger values of $A$, $\Delta D_{\rm f}$ is found to disappear completely (figures not included for brevity).

The results of Fig. 3 lead naturally to the following question: for a given surface tension parameter $A$, is $\Delta D_{\rm f}$ non-negative for all choices of shear-thinning parameters $(C,s)$? To investigate this question, we carried out additional simulations with several different $(C,s)$ pairs, simulating fluids with different shear-thinning properties. The results, given in Table III in App. E.2, suggest that the answer to this question is affirmative: while the size of $\Delta D_{\rm f}$ depends on $(C,s)$ and on $A$, our numerical experiments always find $\Delta D_{\rm f}\geq 0$.

In this work, we discuss the fractal dimension, $D_{\rm f}$, of the patterns that form due to the Saffman-Taylor instability of Hele-Shaw flow. We find that the shear thinning property of the more viscous fluid leads to an overall decrease of $D_{\rm f}$, compared to its Newtonian counterpart. This result is obtained from large-scale simulations of our proposed model, which modifies the well-known approach based on the DLA algorithm for simulating unstable Hele-Shaw two-fluid flow. Such a result would be difficult to obtain either experimentally or via continuous PDE-based simulations since extremely large patterns are needed to produce accurate and converged results.

The manner in which $D_{\rm f}$ changes due to non-Newtonian behavior is nontrivial. First, for vanishing surface tension, non-Newtonian behavior is not relevant since our model reduces to classical DLA simulations. For very large surface tension, the difference in $D_{\rm f}$ between Newtonian and non-Newtonian simulations becomes small, which may be expected in hindsight since the local curvature effects are dominant. For intermediate surface tension values, however, we find consistently smaller values of $D_{\rm f}$ for non-Newtonian fluids, compared to Newtonian ones. We hope that our findings will encourage new research, both experimental and theoretical, that will make progress in identifying the general connection between fluid rheological properties and the fractal dimension of the emerging patterns.