Qing
Zhang
a,
Amin
Amooie
b,
Martin Z.
Bazant
bc and
Irmgard
Bischofberger
*a
aDepartment of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA. E-mail: irmgard@mit.edu
bDepartment of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
cDepartment of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
First published on 11th January 2021
The displacement of a fluid by another less viscous one in a quasi-two dimensional geometry typically leads to complex fingering patterns. In an isotropic system, dense-branching growth arises, which is characterized by repeated tip-splitting of evolving fingers. When anisotropy is present in the interfacial dynamics, the growth morphology changes to dendritic growth characterized by regular structures. We introduce anisotropy by engraving a six-fold symmetric lattice of channels on a Hele-Shaw cell. We show that the morphology transition in miscible fluids depends not only on the previously reported degree of anisotropy set by the lattice topography, but also on the viscosity ratio between the two fluids, ηin/ηout. Remarkably, ηin/ηout and the degree of anisotropy also govern the global features of the dendritic patterns, inducing a systematic change from six-fold towards twelve-fold symmetric dendrites. Varying either control parameter provides a new method to tune the symmetry of complex patterns, which may also have relevance for analogous phenomena of gradient-driven interfacial dynamics, such as directional solidification or electrodeposition.
The phenomenon of viscous fingering has played an important role in elucidating the basic principles of these two types of growth,12–15 as well as methods to control the resulting patterns.16–27 Viscous fingers result from the Saffman–Taylor instability, when one fluid is displaced by another less viscous one in the quasi-two dimensional geometry of a Hele-Shaw cell.28,29 It has been shown that dendritic growth requires anisotropy in the interfacial dynamics.7,30–32 In its absence, dense branching is instead the generic mode of growth.13 Anisotropy fixes the tip of an advancing interface into a stable parabolic shape that prevents it from splitting7,33–35 and introduces global symmetries along preferred growth directions, which are also seen in discrete models of diffusion-limited aggregation on crystal lattices.36–38 Experimentally, anisotropy can be introduced either externally in the growth environment or internally in one of the fluids. External anisotropy can be imposed by engraving ordered channels on one of the plates of a Hele-Shaw cell, by using channels confined with elastic membranes or by placing air bubbles at the tips of growing fingers.30,33,39–43 Internal anisotropy can be induced by replacing one of the fluids with a liquid crystal in the nematic phase.44
Previous studies have considered a particular limit of the viscous-fingering instability; the limit where the viscosity ratio between the less-viscous inner fluid and the more-viscous outer fluid, ηin/ηout, is very low, which is typically the case when air or water displace a viscous liquid. The patterns are then characterized by one single growing length scale, the finger length. Under these conditions, experiments using a Hele-Shaw cell with engraved ordered channels have identified the degree of anisotropy, defined as the ratio between the channel height h and the plate spacing b, h/b, as a control parameter for the morphology transition from dense-branching to dendritic growth;30 dendritic structures form beyond a critical value of h/b. When the two fluids are miscible, the degree of anisotropy is the only control parameter for the morphology transition. In the case of two immiscible fluids, the capillary number sets the critical h/b for the transition.31,45 For miscible fluids and for immiscible fluids at high capillary number, the dendrites directly reflect the underlying symmetry of the lattice; four-fold symmetric dendrites grow in a four-fold symmetric lattice, six-fold symmetric dendrites grow in a six-fold symmetric lattice.10 Dendrites grow in the direction of the channels, which are the regions of largest effective plate spacing within which the flow velocity is highest.30,31
We here reveal how a previously unexplored control parameter, the ratio of the viscosities of the inner and the outer fluid, ηin/ηout, modifies both the morphology transition and, remarkably, the symmetry of the dendritic structures in miscible fluids in anisotropic environments. Recent studies in isotropic environments have identified the viscosity ratio as an important control parameter that governs not only the onset of the instability,46–50 but also the global features of the patterns introducing a second length scale, the radius of a central region of complete outer-fluid displacement that grows concomitantly with the fingers.49,51–53 This central region becomes increasingly larger, and therefore the relative length of the fingers increasingly smaller, as the viscosity ratio between the two fluids increases. Here we show that a morphology transition from dense-branching to dendritic growth can occur over a large range of viscosity ratios. We engrave channels creating a six-fold symmetric lattice on one of the plates and show that the critical degree of anisotropy, h/b, required for the transition to dendritic growth depends on the viscosity ratio between the two liquids. Remarkably, the dendrites can adopt a rich variety of emergent structures: they exhibit six-fold symmetric growth far from the morphology boundary and systematically transition towards twelve-fold symmetric structures as the boundary is approached. Our study reveals novel ways to tune both the morphology transition and the symmetry of dendritic patterns by either controlling the viscosity ratio between the two fluids or the geometric features of the growth environment.
![]() | ||
Fig. 1 (a) Schematic of the modified Hele-Shaw cell. Top image: Top view of the bottom plate of the Hele-Shaw cell with an engraved six-fold symmetric lattice, with width of the lattice channels w and distance between the edges of two channels d. Bottom image: Side view of the modified Hele-Shaw cell, denoting the plate spacing b and the channel height h. (b) Examples of dendritic growth (top, for h/b = 0.5, h = 50 μm, b = 100 μm) and dense-branching growth (bottom, for h/b = 0.04, h = 10 μm, b = 254 μm) at low viscosity ratio ηin/ηout = 0.0013. The scale bar is 1 cm. (c) Morphology diagram controlled by the viscosity ratio ηin/ηout and the degree of anisotropy h/b. Blue symbols denote dense-branching growth, black symbols denote dendritic growth. Experiments are performed with engraved plates with different channel heights h and plate spacings b and at different volumetric flow rates q. (▽) h = 10 μm, q = 1 ml min−1; (×) h = 28 μm, q = 1 ml min−1; (□) h = 28 μm, q = 10 ml min−1; (⋄) h = 50 μm, q = 1 ml min−1; (△) h = 50 μm, q = 10 ml min−1; (○) h = 250 μm, q = 1 ml min−1; (+) h = 250 μm, q = 10 ml min−1. The value of b for each experiment is listed in Table S2 of the ESI.† The solid line denotes a fit to (h/b − (h/b)*)/(ηin/ηout) = A (A = 3 and (h/b)* = 0.04 are best-fit parameters). |
The miscible fluids used in our study are glycerol (PTI Process Chemicals) and water (VWR). We tune the viscosity of the inner fluid by mixing glycerol and water in different proportions and we use pure glycerol as the outer fluid. Details on the composition of the water–glycerol mixtures and their viscosities are reported in Table S1 of the ESI.† The fluids are injected through a 2 mm diameter hole in the center of one of the plates at a precise volumetric flow rate set by a syringe pump (Harvard PHD 2000). We use flow rates of 1 ml min−1 and 10 ml min−1, which allows us to probe an order of magnitude difference in flow rate while staying in the high Péclet number regime (Pe = Ub/D12, where U is the fingertip velocity and D12 is the inter-diffusion coefficient20), here ranging between 2100–45240. Within this regime, the inter-diffusion of the fluids is negligible so that the fluids remain separated by a well-defined interface.20 The patterns are recorded with either a Point Grey camera (Grasshopper 3 GS3-U3-91S6M) at frame rates up to 9 fps or a LUMIX GH5 camera at frame rates up to 60 fps.
We complement the experiments with two-dimensional (2D) high-resolution numerical simulations using the finite element software COMSOL Multiphysics (v5.4), which allows us to access the pressure distribution in the fluids. Our model replicates the geometry of the Hele-Shaw cell in terms of the cell diameter and the six-fold symmetric lattice dimensions. The lowest viscosity ratio we can access in our simulations (ηin/ηout = 0.006) is slightly higher than that probed in experiments (ηin/ηout = 0.0011), as for very low ηin/ηout the fingertip velocity becomes too fast compared to the mean flow velocity, which results in numerically unstable solutions. This is a known issue for a number of numerical approaches,54 but one that could be overcome by designing numerical schemes suited for low viscosity ratios.55 It, however, does not prevent us from accessing the full range of patterns observed in the experiments.
We employ the finite element method to solve the partial differential equations. We couple the convection–diffusion mass-transport equation from the Transport of Diluted Species Module with the continuity equation for the single-phase, incompressible flow velocity from the Darcy's Law Module. The governing equations are:
![]() | (1) |
![]() | (2) |
∇·u = 0 | (3) |
The flow in a Hele-Shaw cell can be approximated as quasi-2D as the plate spacing, b, is much smaller than the radial dimension. The gap-averaged velocity of the fluids is then with k = bi2/12, where bi describes the gap thickness at any point of the textured surface. The spatial variability in the plate spacing is incorporated in the numerical model by defining a binary spatial distribution of permeability [L2], consisting of a permeability value for the obstacles (assigned to the triangles forming the lattice cells), denoted as k1, and a permeability value for the channels (assigned to the background domain), denoted as k2. The ratio between the two permeabilities, k2/k1, is (1 + h/b)2. We use an exponential mixing rule for the mixture viscosity η and η = ηoute−Mc, where ηout is the viscosity of the outer fluid and ηin/ηout = e−M. For miscible fluids, both the pressure and the normal velocity are continuous at the interface. We define a small circular inlet region around the cell center which provides a smooth boundary, to avoid a point-source injection that could lead to a singularity in the domain. A normal inflow velocity for flow and a Dirichlet boundary condition (c = 1) for transport are applied at the perimeter of the circular inlet region, and atmospheric pressure (open-flow) condition for flow and an outflow condition (n·D∇c = 0) for transport are imposed on the outer cell boundary. The initial conditions in the entire domain are c = 0 and p = 0. The absolute values of the injection velocity and the permeability within the computational domain differ from those in experiments. This does not affect the resulting patterns, as only the ratio of the permeabilities governs the pattern morphology.
We solve for pressure and concentration fields in a fully coupled approach using the Parallel Direct Sparse Solver Interface (PARDISO) and Newton's method with dynamic damping for highly nonlinear systems. The implicit Generalized-α Method is used for the time stepping scheme.56,57 We use the default discretization settings that govern the order of discretization in the shape functions for the dependent variables of each module: first-order discretization for the convection–diffusion equation and second-order discretization for Darcy's law, as these settings work efficiently and robustly. The optimal mesh resolution is found with these discretization orders. The annular mesh area used is 0.00606 m2 discretized by 222162 triangular elements. We have confirmed the numerical validity and convergence of our simulations (see ESI†). The discretization by a triangular mesh provides a source of perturbation sufficient for the instability to occur; the apparent slight asymmetry of the computed patterns is mesh driven due to the spatial non-uniformity of the perturbation and the triangularization of the domain.
In agreement with previous studies at very low viscosity ratios and high capillary numbers, we find that the morphology transition from dense-branching growth to dendritic growth occurs above a value of h/b ≈ 0.0530,45 for our lowest ηin/ηout. Below this value, fingers grow by repeated tip-splitting which results in dense-branching growth, above this value, the fingertip is stabilized which results in dendritic growth, as shown in Fig. 1b. Remarkably though, this critical h/b depends strongly on the viscosity ratio: as ηin/ηout increases, a larger h/b is needed to transition from dense-branching growth to dendritic growth, as shown in Fig. 1c. We find that neither the absolute values of the channel height h and the plate spacing b, nor the volumetric flow rate are control parameters for the morphology transition, as shown by the different symbols in Fig. 1c, which denote experiments performed with plates of various channel heights h ranging from 10 μm to 250 μm, various plate spacings b ranging from 125 μm to 1350 μm, and at two volumetric flow rates of 1 ml min−1 and 10 ml min−1. For a given viscosity ratio, any combination of h and b yielding a certain value of h/b leads to the same growth morphology.
To quantify the change from six- towards twelve-fold symmetry, we define the length of the main dendrites, Rm, corresponding to the structures growing in the direction of the six straight channels, and the length of the sub dendrites, Rs, corresponding to the structures growing at an angle of 30° with respect to the six straight channels, as shown in the inset of Fig. 3b. The ratio Rs/Rm exhibits a transient regime at early times and then remains almost constant in time for fully developed patterns. Moreover, Rs/Rm is independent of the interfacial velocity for the range of flow rates investigated, as shown in Fig. 3a, where we normalize the time by t40mm denoting the time when Rm = 40 mm. To compare the patterns formed at different viscosity ratios, we measure Rs/Rm when Rm = 40 mm, which is well within the fully developed regime. For a fixed h/b, the ratio Rs/Rm monotonically increases with viscosity ratio. In addition, a decrease in h/b leads to an increase in Rs/Rm, as shown in Fig. 3b. We can rescale all data by normalizing the viscosity ratio with (h/b − (h/b)*), as shown in Fig. 3c. The factor (h/b)* will become evident in the discussion of the morphology boundary. The numerical results are in good qualitative agreement with the experiments and exhibit the same scaling with h/b, but yield slightly lower values of Rs/Rm compared to the experimental results. This is likely due to the 2D nature of the simulations (where we average the flow in the third dimension across the gap and assume a parabolic velocity profile in the gap direction61–63), which do not capture effects related to the partial displacement of the outer fluid or to the three-dimensional tongue-like structures that form between miscible fluids in a Hele-Shaw cell.49,64,65 Exploring further improvements to the model, e.g., solving Stokes flow in the full 3D domain, and a deeper investigation into quantitative comparisons with experiments are interesting topics for future work.
![]() | (4) |
Why do six-fold dendritic patterns only form far from the morphology boundary, and what leads to the growth of an additional generation of dendrites as we approach the boundary? The importance of the viscosity ratio and the degree of anisotropy for determining Rs/Rm can be seen in a simplified analysis taking into account the effective permeability at different locations corresponding to the growth of sub dendrites or main dendrites, as detailed in the ESI.† Note that the effective permeability in our system is isotropic and lacks a macroscopic preferred direction for single-phase flow. In general, the permeability tensor must be symmetric (by Onsager reciprocity for Stokes flow) and positive definite (by the Second Law of Thermodynamics) and thus represented by an orthogonal matrix,66 so its eigenvectors, corresponding to the fastest and slowest directions, must be mutually perpendicular.67 This orthogonality is incompatible with triangular symmetry, thus the permeability eigenvalues in our textured Hele-Shaw cell must be degenerate, implying isotropic single-phase flow.
For two-phase flow, however, the gradient of viscosity at the interface between the two fluids can locally break the symmetry and induce an anisotropic effective permeability near the interface. Using concepts derived for the hydrodynamics of slippage on textured surfaces for two-phase flows over hydrophobic surfaces,68,69 we consider that the more-viscous outer fluid is partially trapped in the texture as the tip of the less-viscous fluid passes over the texture in the middle of the channel along the “path of least resistance”. For small textures, the trapped fluid leads to a local effective slip length tensor,70,71bslip, which causes the effective permeability tensor to become anisotropic and orthogonal in the vicinity of the interface,67 leading to the appearance of sub dendrites that impart this square symmetry to the pattern. In the limit of “weak anisotropy” in the slip tensor, Tr(bslip) ≪ b, as is the case for our experiments, we find that the interface velocities of the sub dendrites and the main dendrites, and therefore Rs/Rm, are indeed governed by ηin/ηout and h/b.
To get further insight into the growth of the dendrites, we consider their macroscopic path selection. The main dendrites Rm grow along the six straight channels. The sub dendrites Rs select a path at a 30° angle from these straight channels. At early stage, two fingers form between each pair of neighboring main dendrites on each side of the 30° direction, due to the anisotropy of the lattice. This is observed at any viscosity ratio, as shown in Fig. 4a. Whether these fingers will merge towards each other and grow into a sub dendrite or merge with the main dendrites resulting in a six-fold symmetric pattern depends on the pressure distribution imposed both globally by the main dendrites and locally at the tip of the sub dendrites. At low ηin/ηout and high h/b, the rapid growth of the main dendrites sets up a large pressure gradient at their tip which in turn induces a small pressure gradient in the 30° direction, as shown in Fig. 4a, which prevents the sub dendrites from growing. With increasing ηin/ηout and decreasing h/b, however, the sub dendrites themselves build locally a high pressure gradient at their tips which amplifies their growth. We provide further details on the growth of the sub dendrites in Fig. S5 and S6 of the ESI.†
Once the sub dendrites have emerged, they continue to grow along the 30° direction following a zig-zag path, as illustrated in Fig. 4b. As the tip of the sub dendrite reaches a lattice junction, indicated by a red dot, the path towards the 30° direction (red arrow) is selected, rather than the straight path (blue arrow). This selection results from the pressure profile induced in the outer fluid by the main dendrites, which effectively shields the sub dendrites from growing towards the main dendrites and pushes them towards the 30° direction. Indeed, when the tip of a sub dendrite reaches the entrance of a lattice junction, as schematically shown in the zoomed-in region in Fig. 4c, it does not grow straight towards channel 2, but is deviated towards the 30° direction as a result of the global pressure distribution built up by the neighboring main dendrites. The local pressure distribution at the tip of the sub dendrite then induces a maximum pressure gradient towards channel 1, and most of the flow goes into channel 1. It is this combination of the global pressure distribution from the main dendrites and the local pressure distribution from the tip of the sub dendrites that leads to the rich pattern selection in dendritic growth.
These different paths selected by the main dendrites and sub dendrites also reveal the origin of the maximum value of Rs/Rm ≈ 0.85. It reflects the condition where the velocity of the main and sub dendrites becomes approximately equal. As the path selected by the sub dendrites deviates from the radial direction at each junction, the total path is times longer than that of the main dendrites in the straight radial channels. The length of the main dendrite, Rm, is therefore
, i.e.,
. Interestingly, our experiments show that once this condition is reached, a further increase in ηin/ηout or decrease in h/b induces the morphology transition to dense-branching growth. This suggests that the morphology transition occurs when the difference between the pressure gradient in the straight channels and the 30° direction becomes negligible, and therefore the role of anisotropy becomes negligible, such that the parabolic tips can no longer be stabilized allowing for tip-splitting to occur. Reaching a full understanding of this morphology transition can be topic of further research.
This diversity of different dendritic patterns provides novel opportunities for tuning the growth of complex structures, not only in viscous fingering, but perhaps also in other cases of interfacial motion limited by gradient-driven transport processes, which lie in the same universality class.72 In general, we expect that dendritic growth following the preferred directions of an anisotropic environment will tend to acquire orthogonal symmetry for “weak anisotropy”, whenever transport near the interface is governed by a local effective conductance tensor, which must be orthogonal like the effective slip tensor in a weakly textured Hele-Shaw cell.67 For example, in template-assisted directional solidification,73,74 a similar morphological transition may arise, controlled by the ratio of thermal diffusivities (analogous to the ratio of inverse viscosities here), whenever the pattern is controlled by the conduction of latent heat away from the interface in the liquid phase. Similarly, in template-assisted electrodeposition,75–78 it may be possible to tune the symmetry of dendritic patterns by varying the strength of diffusion anisotropy in the electrolyte domain. Active control of anisotropic dendritic growth may also be achieved, for example, by applying electric fields to control viscous fingering18,19 over patterned, charged surfaces79 having anisotropic electro-osmotic slip tensors.80
Footnote |
† Electronic supplementary information (ESI) available: Information about experimental parameters. The role of diffusion in the numerical simulations. The convergence of the numerical simulations. Movies of dendritic growth for viscosity ratios ηin/ηout = 0.0013 and ηin/ηout = 0.05. A simplified model to account for the effect of h/b and the viscosity ratio on the symmetry selection. Additional details on the growth of sub dendrites. See DOI: 10.1039/d0sm01706j |
This journal is © The Royal Society of Chemistry 2021 |