Numerical investigation of controlling interfacial instabilities in non-standard Hele-Shaw configurations
TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Numerical investigation of controllinginterfacial instabilities in non-standardHele-Shaw configurations
Liam C. Morrow, Timothy J. Moroney, and Scott W. McCue † School of Mathematical Sciences, Queensland University of Technology, Brisbane QLD 4001,Australia(Draft manuscript as of 4 July 2019)
Viscous fingering experiments in Hele-Shaw cells lead to striking pattern formationswhich have been the subject of intense focus among the physics and applied mathematicscommunity for many years. In recent times, much attention has been devoted to devisingstrategies for controlling such patterns and reducing the growth of the interfacial fingers.We continue this research by reporting on numerical simulations, based on the level setmethod, of a generalised Hele-Shaw model for which the geometry of the Hele-Shawcell is altered. First, we investigate how imposing constant and time-dependent injectionrates in a Hele-Shaw cell that is either standard, tapered or rotating can be used toreduce the development of viscous fingering when an inviscid fluid is injected into aviscous fluid over a finite time period. We perform a series of numerical experimentscomparing the effectiveness of each strategy to determine how these non-standard Hele-Shaw configurations influence the morphological features of the inviscid-viscous fluidinterface. Surprisingly, a converging or diverging taper of the plates leads to reducedmetrics of viscous fingering at the final time when compared to the standard parallelconfiguration, especially with carefully chosen injection rates; for the rotating platecase, the effect is even more dramatic, with sufficiently large rotation rates completelystabilising the interface. Next, we illustrate how the number of non-splitting fingers canbe controlled by injecting the inviscid fluid at a time-dependent rate while increasingthe gap between the plates. Our simulations compare well with previous experimentalresults for various injection rates and geometric configurations. We demonstrate how thenumber of non-splitting fingers agrees with that predicted from linear stability theoryup to some finger number; for larger values of our control parameter, the fully nonlineardynamics of the problem lead to slightly fewer fingers than this linear prediction.
1. Introduction
A standard Hele-Shaw cell is an experimental device (figure 1) consisting of two parallelplates separated by a small gap filled with a viscous fluid. Fluid flow in this device hasreceived significant attention largely due to the interfacial patterns that form when aninviscid fluid is injected into the viscous fluid. These viscous fingering patterns form dueto the Saffman-Taylor instability (Saffman & Taylor 1958), and are characterised by theirdistinctive branching and tip-splitting behaviour. Closely related interfacial instabilitiesappear in a wide variety of phenomena, including saturated flow in porous media (Homsy1987), the growth of bacterial colonies (Ben-Jacob et al. † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] J u l L. C. Morrow, T. J. Moroney, and S. W. McCue
Q b
Figure 1.
Illustration of the standard Hele-Shaw cell experiment with two parallel platesseparated by a small gap filled with viscous fluid (shaded). An inviscid fluid (white) is injectedinto the viscous fluid at a rate Q . The immiscible fluids are separated by a sharp interface,which becomes increasingly unstable as it expands, forming distinct viscous fingering patternscharacterised by their branching and tip-splitting morphology. describe these processes (Ben-Jacob & Garik 1990; Liang 1986; Li et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. ontrolling interfacial instabilities in Hele-Shaw flow et al. (2012) and determinehow effective this strategy is over a range of parameter values, including a number ofexamples in which there is significant fingering (the only numerical example of thisstrategy provided by Dias et al. (2012) involved a near-circular interface). For the caseof tapered plates, we extend the work of Al-Housseiny & Stone (2013); Bongrand & Tsai(2018), which involved experiments and linear stability theory, by performing a series ofnumerical simulations over a wider range of injection rates and taper angles. We find thatour new optimal injection rate appears to noticeably reduce the fingering pattern (via areduction in the isoperimetric and circularity metrics) for the converging case, producingan atypical fingering pattern with short and stubby fingers (which appear similar to thoseobserved by Pihler-Puzovi´c et al. (2012, 2013)). On the other hand, for the diverginggeometry, this optimal injection rate also appears to reduce the instability, althoughwith a much less dramatic effect. Finally, the case of rotating plates with injection ofinviscid fluid has not been considered previously in the literature. Here, we explore casesin which fingers initially develop in the usual way; however, the centrifugal force acts tostabilise the interface so that the bubble ends up tending to a circle in shape in the longtime limit. Physically speaking, this effect is due to the centrifugal force propelling thedense fluid outward which stabilises the interface.The second objective with which we shall be concerned involves adjusting the flow rateand the geometry of the experimental apparatus in an attempt to prevent ongoing tipsplitting so that the pattern evolves with a predetermined number of fingers. For example,numerical, weakly nonlinear and experimental studies indicate that for a standard Hele-Shaw cell, an injection rate with the scaling Q ∼ t − / can produce N -fold symmetricbubbles whose shape is independent of the initial condition, and can be controlled by thestrength of the injection rate (Brener et al. et al. et al. (2015), who applylinear stability analysis and experimental results to provide evidence that a constantnumber of fingers should develop if the parallel plates are separated via the scaling b ∼ t / . We extend this work by providing numerical evidence confirming that the numberof non-splitting fingers can be controlled by implementing a more complicated time-dependent injection rate at the same time as separating the plates, as proposed by Zheng et al. (2015). Further, our simulations provide insight into how interactions betweenneighbouring fingers can influence the evolution of the interface extending beyond linearstability analysis; indeed, the number of non-splitting fingers is observed to be less thanthat predicted by linear stability analysis for a sufficiently large control parameter. Ournumerical results here are consistent with the very recent findings of Vaquero-Stainer et al. (2019), who use a finite element scheme to also explore the b ∼ t / scaling.The outline of this paper is as follows. In §
2, we summarise a generalised model forHele-Shaw flow in a non-standard geometry, for which the gap between the plates dependson both time and space and the plates are allowed to rotate. In §
3, we compare numerical
L. C. Morrow, T. J. Moroney, and S. W. McCue simulations for the standard Hele-Shaw configuration (parallel stationary plates withconstant injection rate) with experimental results and predictions from linear stabilityanalysis. We show that the numerical simulations agree with predictions from linearstability analysis for small time, and reproduce key morphological features observed inexperiments for large time. In §
4, we consider the objective of injecting a prescribedamount of fluid over fixed period of time for the different geometric configurations ofparallel, tapered or rotating plates. Subsequently, in §
5, we study the objective ofpreventing tip-splitting and controlling the number of viscous fingers by carefully alteringthe time-dependent gap and/or injection rate. Finally, in §
2. Mathematical model
Governing equations
The geometry we consider involves the injection of an incompressible inviscid fluid(with flow rate Q ) through a small orifice in the centre of a Hele-Shaw cell otherwisefilled with a viscous fluid (figure 1). We assume the two fluids are immiscible and denotethe simply connected domain of inviscid fluid by Ω( t ) and the interface between the twofluids by ∂ Ω. In our model, the viscous fluid is infinite in its extent, and so while theproblem is driven by injection of an inviscid fluid at a point, we can also interpret theflow as being driven by a suction of viscous fluid from infinity. A feature of our model isthat we allow the small gap between the plates, b , to depend on both space and time.We use a two-dimensional model of Hele-Shaw flow in a rotating frame that is derivedby averaging Stokes flow over the small gap between the plates. Denoting ˆ p , v , µ , and ρ as the pressure, velocity, viscosity, and density of the viscous fluid, the governing fieldequations modified to incorporate Hele-Shaw plates rotating at angular velocity ˆ ω are(Carrillo et al. v = − b µ ` ∇ ˆ p − ˆ ω ρr e r ˘ , x ∈ R \ Ω( t ) . (2.1)We note that the only effect of the rotation considered is the centrifugal force, andthe Coriolis force is neglected. The centrifugal term in (2.1) is removed by introducing p = ˆ p − ˆ ω ρr /
2, and thus we have v = − b µ ∇ p, x ∈ R \ Ω( t ) , (2.2) ∇ · p b v q = − ∂b∂t . x ∈ R \ Ω( t ) , (2.3)noting that for all the cases we consider ∂b/∂t is spatially uniform. Equation (2.2) isanalogous to Darcy’s law, which provides an intimate connection between Hele-Shawflow and porous media flow (Homsy 1987). Equation (2.3) ensures that the fluid’s volumeis conserved, and reduces to the traditional divergence free condition, ∇ · v = 0, in thestandard configuration for which the plates are parallel and stationary. Note that we shallignore the pressure gradients in the inviscid domain Ω( t ), which makes this a one-phaseHele-Shaw model. The pressure of the bubble is taken as the reference pressure, so p = 0at all times in the bubble.By substituting (2.2) into (2.3), we have the Reynolds lubrication equation ∇ · ˆ b µ ∇ p ˙ = ∂b∂t , x ∈ R \ Ω( t ) . (2.4) ontrolling interfacial instabilities in Hele-Shaw flow ∇ · (cid:18) b µ ∇ p (cid:19) = ∂b∂t b µ ∂p∂r ∼ − Q πr + r ∂b∂t as r → ∞ p = − σ (cid:0) κ + b (cid:1) − ωr v n = − b µ ∂p∂n Figure 2.
A schematic of our generalised Hele-Shaw model with spatially and/or temporallydependent plate gap thickness and rotating plates (2.2)-(2.7). The viscous fluid is representedby the shaded (blue) region and the inviscid bubble is represented by the white region.
The boundary conditions on the interface are p = − σ ˆ κ + 2 b ˙ − ωr , x ∈ ∂ Ω( t ) , (2.5) v n = − b µ ∂p∂n , x ∈ ∂ Ω( t ) , (2.6)where the centrifugal parameter ω = ρ ˆ ω /
2. The dynamic boundary condition (2.5)incorporates the effects of surface tension via the Young-Laplace equation, where σ isthe surface tension parameter and κ is the signed curvature of the interface in the lateraldirection. The term 2 /b in (2.5) represents the curvature in the transverse direction forthe case where the fluid is perfectly wetting (McLean & Saffman 1981). The kinematicboundary condition (2.6) equates the velocity of the interface to the velocity of the viscousfluid on the interface. We note that both viscous stresses and the effect of a thin wettingfilm left behind by the viscous fluid are ignored. The far-field boundary condition is b µ ∂p∂r ∼ − Q πr + 12 r ∂b∂t r → ∞ , (2.7)where Q is a time-dependent flow-rate at which the inviscid fluid is injected. This formof the far-field boundary condition ensures the rate of change of volume of the inviscidbubble is indeed given by Q . Our model is summarised by a schematic in figure 2.We note in passing that the complementary geometry with viscous fluid in Ω( t )and inviscid fluid in R \ Ω( t ) (that is, the opposite case with the fluids swapped) hasattracted interest in the literature. For that scenario, both the lifting and centrifugalconfigurations produce fingers which appear to be distinct from traditional Saffman-Taylor fingers. These problems with the complementary geometry have been studiedthrough a combination of experimental, analytical, and numerical techniques (Alvarez-Lacalle et al. et al. et al. et al. et al. et al. et al. L. C. Morrow, T. J. Moroney, and S. W. McCue
Numerical scheme
Many numerical schemes used to study viscous fingering in a standard Hele-Shawcell, where the governing equation for pressure (2.4) reduces to Laplace’s equation ∇ p = 0, implement a boundary integral method (Dai & Shelley 1993; DeGregoria &Schwartz 1986; Li et al. et al. et al. et al. et al.
3. Review of standard configuration
Most mathematical studies investigating the influence of manipulating the geometryof the Hele-Shaw cell on viscous fingering are performed using linear stability analysis.While this approach provides a useful tool for understanding the qualitative behaviour ofsolutions, as well as for deriving strategies for controlling viscous finger development, itis only accurate for small time and, as such, does not capture the full nonlinear dynamicsof the problem. In this section, we review linear stability analysis for the standard Hele-Shaw problem where the plates are parallel and stationary ( b constant) and the inviscidfluid is injected at a constant rate Q . Further, we show that our numerical simulationsare consistent with predictions made by linear stability analysis when time is sufficientlysmall, and can accurately reproduce experimental results for longer times.Considering (2.2)-(2.7) in polar coordinates ( r, θ ) with p = p ( r, θ, t ) and the interface ∂ Ω denoted by r = s ( θ, t ), we assume a perturbed circular solution p ( r, θ, t ) = p ( r, t ) + ε ∞ (cid:88) n =2 P n ( t ) r − n cos nθ + O ( ε ) , (3.1) s ( θ, t ) = s ( t ) + ε ∞ (cid:88) n =2 γ n ( t ) cos nθ + O ( ε ) , (3.2)where ε (cid:28)
1. The leading order radius of the interface becomes s = ˆ s (0) + Qtπb ˙ / . The resulting differential equation for the n th mode of perturbation is (Paterson 1981) γ n γ n = n − s ˆ Q πbs − n ( n + 1) b σ µs ˙ , (3.3) ontrolling interfacial instabilities in Hele-Shaw flow Figure 3. ( a - c ) Experimental results from Chen (1987), reproduced with permission fromSpringer Nature, comparing the development of viscous fingers for different injection rates and( d - f ) the corresponding numerical simulations. The gap thickness between plates is 7 . × − cm, and the injection rate and final times (left to right) are: Q = 2 . × − mL/s and t f = 65s; Q = 4 . × − mL/s and t f = 490 s; and Q = 1 . × − mL/s and t f = 1650 s. Numericalsolutions are plotted in time intervals of t f /
10. Additionally, σ = 20 g/s , µ = 10 . · s),and ω = 0 g/(s · mL). The initial condition is of the form (4.3) with R = 0 .
25 cm. Simulationsare performed on the domain 0 ≤ r ≤ ≤ θ < π using 750 ×
942 equally spaced nodes. where the most unstable mode of perturbation, n max , is predicted to be n max = d ˆ µQs πσb ˙ . (3.4)Equation (3.4) comes from setting ∂ ( γ n /γ n ) /∂n = 0 and solving for n . As such, n max isnot an integer, and so in practice the most unstable mode is the closest integer to n max .Note that, given s is an increasing function of time, then n max also increases in time,which means the most unstable mode predicted by linear stability is a dynamic property(while not strictly relevant for fully nonlinear pattern formation, this observation isclosely related to the ongoing tip-splitting that occurs for longer times).In figure 3 we compare experimental results obtained by Chen (1987) with our nu-merical simulations, and test some predictions made by linear stability analysis. Thisfigure illustrates the classic pattern formation in a standard Hele-Shaw configuration forthree different injection rates in decreasing order. For early times, we can apply (3.4)to predict the number of fingers that are produced. Using figure 3( e ) as an example,equation (3.4) (with the appropriate parameter values) predicts that n max ≈ e ) shows this prediction is consistent with the numerical simulation.Another straightforward result from (3.3) is that, for fixed σ , b and µ , increasing theflow rate Q results in a positive contribution to γ n /γ n , which has a destabilising effecton each mode. Further, we see from (3.4) that increasing Q increases the most unstablewave number. These observations are consistent with the experimental measurementsperformed by Chen (1987) (and many others), which appear to show that increasing L. C. Morrow, T. J. Moroney, and S. W. McCue the injection rate results in larger wave numbers becoming more unstable, leading tobranching and tip-splitting. Of course, for later times, nonlinear effects become significantand linear stability analysis no longer provides an accurate description of the solution.In this nonlinear regime, our numerical simulations are able to reproduce the mainmorphological features of these experiments for the different injection rates considered,as we can see by comparing images in each column of figure 3. We view this comparisonas a preliminary test of our numerical method.
4. Reducing growth of viscous fingering pattern
In this section, we investigate strategies for controlling viscous fingering when aprescribed amount of the inviscid fluid is injected over a finite period of time. To begin, in § § et al. (2012). Following these parts, in § I ( t ) = L πA , (4.1)where L and A are the arc length and area enclosed by the interface, and the ratio of thetip to base radii, which we refer to as the circularity ratio (the “roundness”), defined as C ( t ) = R outer R inner , (4.2)where R outer is the radius of the smallest circle (centred at the origin) that completelyencloses the bubble and R inner is the radius of the largest circle (centred at the origin)that contains only inviscid fluid. Both I and C will be unity when the interface is circularand increase as the instabilities cause the interface to deform away from a circle.In this section, the initial condition of the interface is s ( θ,
0) = R ˜ . (cid:88) n =2 cos p n p θ − πθ n qq ¸ , (4.3)where θ n is a uniformly random number between 0 and 1. Simulations are performed onthe domain 0 ≤ r ≤ . ≤ θ < π using 750 ×
628 equally spaced nodes. For eachparameter combination considered, 10 simulations are performed, and I and C are bothaveraged over these simulations. For simulations in this section we use σ = 3 g/s and µ = 1 .
58 g/(cm · s). ontrolling interfacial instabilities in Hele-Shaw flow Time-dependent injection rate
The first strategy we consider is proposed by Dias et al. (2012), who, using linear sta-bility analysis and optimal control theory, derived the optimal time-dependent injectionrate when the plate gap thickness is uniform. By seeking solutions of the form (3.1) and(3.2), Dias et al . showed that the growth rate of the most unstable perturbations to thecircular solution, s , when the inviscid fluid is injected over the time interval 0 ≤ t ≤ t f ,are minimised when Q ( t ) = 2 πb ( R f − R ) t f ˆ R + R f − R t f t ˙ , (4.4)where s (0) = R and s ( t f ) = R f . The average of (4.4) over this time period is¯ Q = πb ( R f − R ) t f . (4.5)We are interested in comparing results from the linear injection rate (4.4) with a constantinjection rate where Q = ¯ Q , so that in both cases the same amount of fluid is injectedover the fixed time period. Using both experiments and numerical simulations, Dias et al. (2012) showed that (4.4) does suppress the growth of viscous fingers comparedto (4.5); however, only cases for which the injection rate is sufficiently low that viscousfingers were completely suppressed were considered. We extend this work by performingsimulations over a much wider range of injection rates to better compare the developmentof viscous fingers between the injection rates of the forms (4.4) and (4.5). Figure 4presents numerical solutions for the constant (top row) and linear (second row) injectionrates. The columns from left to right are for increasing values of ¯ Q . We observe that thelinear injection rate appears to inhibit viscous fingering, and in particular, tip-splittingis delayed resulting in shorter fingers than the corresponding constant injection case.These are only visual observations. Focusing on the representative case ¯ Q = 2 mL/s(fourth column of figure 4), a more quantitative measure is provided in figure 5, wherethe isoperimetric ratio I is plotted against time for both injection schemes. Initially, theisoperimetric ratio of the linear injection rate case (solid blue curve in figure 5( b )) growsmuch more slowly than the constant injection rate case (solid blue curve in figure 5( a )),which is simply because the linear injection rate is lower than the constant injectionrate for the first half of the simulation, and so the interface is less unstable. For latertimes, the isoperimetric ratio of the linear injection rate case grows faster correspondingto times for which the linear injection rate is faster. Despite this switch in behaviour, theoverall effect of the linear injection rate (4.4) is to noticeably reduce the isoperimetricratio at the final time t f when compared to the constant injection rate. These numericalresults provide new quantitative evidence for how well Dias et al. (2012)’s ‘optimal’ flowrate works in practice (we return to figure 5 below).To investigate the robustness of the linear injection strategy as ¯ Q is varied, we computeboth the isoperimetric ratio I and the circularity ratio C according to (4.1) and (4.2) at t = t f for both the constant and linear injection rates, shown in figure 6. The results forthe constant injection rates are denoted by (navy blue) • , while the linear injection rateis indicated by (red) IJ . Recall that for each data point, 10 simulations are performed andthe error bars indicate plus or minus one standard deviation. Over the range of ¯ Q valuesconsidered, both of the measures I and C are considerably lower for the linear injectionrate case compared to the constant injection rate. Thus we conclude the linear injectionscheme is successful in reducing the fingering pattern, regardless of ¯ Q .0 L. C. Morrow, T. J. Moroney, and S. W. McCue
Parallel plates ( b = 0 . b = 0 . α = 5 . × − r = 7 cm b = 5 × − cm) with constant injectionConverging plates ( α = 5 . × − r = 7 cm b = 5 × − cm) with optimal injectionDiverging plates, ( α = − . × − r = 7 cm b = 5 × − cm) with constant injectionDiverging plates ( α = − . × − r = 7 cm b = 5 × − cm) with optimal injection Figure 4.
Comparison of numerical solutions with different injection schemes and plate gapthicknesses. For each column, the average injection rate is ¯ Q = 0 .
6, 1.2, 1.6, and 2 mL/s and thefinal time is t f = 25 .
62, 12.81, 9.61, and 7.78 s. Rows 1, 3, and 5 have constant injection rate¯ Q . For row 2, the injection rate is (4.4), while for rows 4 and 6, the injection rate is determinedfrom the solution to (4.13). We use R = 0 . R = 0 .
37 cm for rows 3and 4, and R = 0 .
94 cm for rows 5 and 6. Additionally, R f = 5 cm for rows 2, 4, and 6. Forall simulations, the initial volume of the bubble is approximately 0.157 mL, and the volume at t = t f is 15.7 mL. Profiles are plotted in time intervals of t f /
10. The scale bar represents alength of 2 cm. ontrolling interfacial instabilities in Hele-Shaw flow Figure 5.
Isoperimetric ratio for different plate tapering configurations with ( a ) constant and( b ) time-dependent injection rates, where the average injection rate of all configurations is ¯ Q = 2mL/s and final time t f = 7 .
68 s. Solid (blue) curve denotes the configuration with parallel plates b = 0 . α = − . × − , r = 7 cm, and b = 0 .
395 cm, and dashed (red) curve is for tapered plate with α = 5 . × − , r = 7 cm, and b = 5 × − cm. For ( b ) the injection rate for the solid curve is given by (4.4)with R = 0 . R f = 5 cm, while for the dashed and dotted curves, the injection rate isdetermined from the solution to (4.13) with R = 0 .
37 cm and R f = 5 cm, and R = 0 .
94 cmand R f = 5 cm respectively. Shaded (blue) region represents one standard deviation above andbelow the mean. Figure 6. ( a ) Isoperimetric (4.1) and ( b ) circularity (4.2) ratios at the final time t = t f fordifferent injection schemes and plate tapering. The parameter ¯ Q (mL/s) is the average injectionrate over a simulation, and the final time is t f = ((99 π/
20) mL) / ¯ Q . Navy blue ( • ) and redcurves ( IJ ) denote cases with parallel plates b = 0 . İ ) and purple ( § ) denotes caseswith tapered plates where α = 5 . × − , r = 7 cm, and b = 5 × − cm, and green ( ‚ ) andlight blue ( đ ) is tapered plates where α = − . × − , r = 7 cm, and b = 0 .
395 cm. For • , İ , and ‚ , the inviscid bubble is injected at a constant rate, ¯ Q . For IJ , injection rate is (4.4) with R = 0 . R f = 5 cm. For § and đ , injection rate is computed from the solution to (4.13)with R = 0 .
37 cm and R f = 5 cm, and R = 0 .
94 cm and R f = 5 cm respectively. Tapered Hele-Shaw geometry
We now turn our attention to the configuration where the gap between the plates islinearly tapered in the direction of the flow such that b ( r ) = (cid:40) b − α ( r − r ) if r ≤ r ,b if r > r , (4.6)together with a constant injection rate. The parameter α = ( b (0) − b ) /r controls thegradient of the taper. The influence of tapering the plates of a Hele-Shaw cell has beenstudied using linear stability analysis in both channel and radial geometry (Al-Housseiny et al. | α |(cid:28)
1, Al-Housseiny & Stone2
L. C. Morrow, T. J. Moroney, and S. W. McCue
Figure 7. ( a - c ) Experimental results from Bongrand & Tsai (2018), reproduced with permissionfrom the American Physical Society, comparing the viscous fingering pattern for differentinjection rates and plate configurations, and ( d - f ) the corresponding numerical simulations.The gap thickness is of the form (4.6) with ( d ) Q = 2 / b = 0 .
12 cm, α = 0, and t f = 14s, ( e ) Q = 2 / b = 1 . × − cm, α = 6 . × − , r = 7 cm, and t f = 30 s, and ( f ) Q = 11 / b = 5 × − cm, α = 4 . × − , r = 7 cm, and t f = 4 . R = 0 . (2013) derived an ordinary differential equation for γ n γ n γ n = n − s ˆ s − n ( n + 1) b ( s ) σ µs ˙ − ˆ nσ µs − s b ( s ) ˙ α, (4.7)where s = Q/ πb ( s ) s , which suggests that the most unstable mode of perturbation is n max = d
13 + 4 µs s b ( s ) σ − s α b ( s ) . (4.8)Of course, by setting α = 0, (4.7) and (4.8) reduce to (3.3) and (3.4). As noted by Al-Housseiny & Stone (2013), diverging plates ( α <
0) introduces a negative offset to γ n /γ n in the form of s α/b ( s ), while for converging plates ( α > § α < α > s , is initially higher for the divergingcase ( α <
0) than it is for the converging case ( α > α > et al. (2018) and numerically by Jackson et al. (2017) (who only considered the evolutionof 8-fold symmetric bubbles). We extend these studies by providing insight into howeffective tapering the plate gap is at reducing the development of viscous fingering bycomparing simulations over a range of values of α and Q to the corresponding parallel ontrolling interfacial instabilities in Hele-Shaw flow α , b , and Q , shown in figure 7. For parallel plates ( α = 0)our simulations are able to reproduce the classic viscous fingering patterns observedexperimentally (figure 7( a ),( d )). When the plates are converging and the injection rateis sufficiently low, experimentally it is observed that the interface is stabilised, whichis reproduced by our numerical simulations (figure 7( b ),( e )). For a faster injection rate,the interface is unstable and develops fingers that appear slightly different to traditionalviscous fingers (although the mechanism is presumably the same); our numerical solutionis able to reproduce this morphology (figure 7( c ),( f )).Returning to our control objective that involves injecting the same volume of inviscidfluid over a fixed period of time, we compare numerical simulations in figure 4 for differentvalues of α over four different constant injection rates. We see that tapering the platesin the direction of flow ( α = 5 . × − ; third row of figure 4) delays tip-splitting andproduces shorter fingers compared to the corresponding parallel plate configuration (firstrow of figure 4). Furthermore, for the lowest of the four injection rates, ¯ Q = 0 . t f , the interface develops numerous shortfingers (Bongrand & Tsai (2018) refer to these as “wavy” fingers) which does not occurwhen α = 0. On the other hand, tapering the plates so they are diverging in the directionof flow has a qualitatively different effect. Here ( α = − . × − ; fifth row of figure 4),simulations indicate that the interface develops numerous long fingers and tip-splitting isreduced compared to the parallel plate case. Interestingly, increasing ¯ Q does not appearto significantly increase either the number or length of fingers that develop compared tothe other configurations considered.To better quantify how the interfacial instabilities develop when the plates are taperedover the duration of a simulation, we compare the isoperimetric ratio I for the taperedand parallel plate configurations for a particular flow rate, shown in figure 5( a ). Whenthe plates are converging (dashed red), the isoperimetric ratio initially grows much slowercompared to when the plates are parallel (solid blue), while for later times, it increasesat a faster rate. This behaviour can be explained by noting that, in order to injectthe required volume of fluid over the time period, the interface must be slower in thetapered case for small times and faster for later times (analogous to the linear injectionrate (4.4)). For the diverging case (dotted yellow), the opposite trend is observed; here,the isoperimetric ratio initially grows more rapidly than the parallel plate case as thevelocity of the interface is initially higher. However, as the bubble expands, we find thatthe growth of I slows for times leading up to t f . The leading order effects in (4.7) suggestthe two mechanisms responsible for this reduction in the growth of I are the decrease inthe interface’s speed and an increase in the stabilising effect of surface tension. Despitethe differences in stabilising and destabilising effects for the converging and divergingcases, both result in a reduction in I at t = t f .Both the isoperimetric ratio I and circularity ratio C at the final time t = t f areshown in figure 6 for various values of ¯ Q . For ¯ Q ≤ İ ) produces a more circular interface than the parallel platecase (navy blue, • ). Furthermore, figure 6( a ) indicates that I may increase above theother configurations if ¯ Q >
L. C. Morrow, T. J. Moroney, and S. W. McCue -0.06 -0.04 -0.02 0 0.02 0.04 0.06051015
Figure 8.
Isoperimetric ratio, I , at final time t f as a function of the gradient of the taper, α ,for injection rate Q = 2 (blue, • ), 1.6 (red , ‚ ), and 1.2 (yellow, ˛ ) mL/s with t f = 7 .
78, 9.61, and12.81 s. For all simulations, r = 7 cm, and b and R are chosen such that for each simulation,the volume of the inviscid bubble is 0.157 mL at t = 0 and 15.7 mL at t = t f . The value of I iscomputed by performing 10 numerical simulations with initial condition (4.3), and averaging I between the simulations. the maximum radius of the interface can increase above 7 cm, which, according to (4.6),is where the plates are no longer tapered. While this increase in normal velocity results ina sharp increase in I as ¯ Q increases, it does not appear to significantly increase the lengthof the fingers as a corresponding sharp increase in C is not observed. For the divergingplate case (green, ‚ ), both the isoperimetric and circularity ratios are larger than thatof the parallel plate case when the injection rate is slow. However, as ¯ Q is increased,both of these quantities become smaller compared to the parallel case. Further, over therange of values of ¯ Q considered, there is relatively little variation in I and C compared tothe other configurations. Thus our results indicate that compared to the correspondingparallel configuration, we can produce a more circular interface for both slower and fasterinjection rates by imposing linearly converging and diverging plates, respectively.To clarify the influence of α on the development of viscous fingering, we compute theisoperimetric ratio from numerical simulations for values of α between − . × − and5 . × − with Q = 2 (blue, • ), 1.6 (red , ‚ ), and 1.2 (yellow, ˛ ) mL/s, and show theresults in Figure 8. This figure indicates that for each value of Q considered, I is a non-monotonic function of α , and has a maximum at α ≈ − .
02. Furthermore, for Q = 1 . α = − . × − and 5 . × − result in a reduction of I compared to α = 0 (that is, tapering either way reduces the fingering pattern). However, choosing α = 5 . × − results in the smallest value of I for each of the injection rates considered(for these injection rates, converging plates has a greater effect of reducing the fingeringpattern than diverging plates).4.3. Tapered plates with time-dependent injection In § § et al. (2012), deriving an optimal injection rate that attempts tominimise viscous fingering when the gap thickness is of the form (4.6). We also performnumerical simulations to investigate the effectiveness of this configuration. In contrast to ontrolling interfacial instabilities in Hele-Shaw flow Figure 9.
Sketch of different injection rates as a function of time. Solid (blue) curve is constantinjection rate ( Q = 1 . b = 0 . R = 0 . b = 0 . α = 5 . × − , b = 5 × − cm, r = 7 cm, and R = 0 .
37 cm. Dash-dotted (purple) curve is optimal injection rate computedfrom the solution to (4.13) with α = − . × − , b = 0 .
395 cm, r = 7 cm, and R = 0 . , t f ] is the same for allinjection rates. Additional parameters are t f = 8 .
54 s and R f = 5 cm. the strategies discussed in § § µ s s /σb (cid:29) n max ≈ d µ s s σb ( s ) − s α b ( s ) . (4.9)Furthermore, (4.7) evaluated at n = n max reduces to λ ( s , s ) = γ n max γ n max ≈ c µb ( s ) σ ˆ s − σα µ ˙ / − s s + s αb ( s ) . (4.10)The idea presented by Dias et al. (2012) is to determine an injection rate that minimisesthe integral (cid:90) t f λ d t, (4.11)which is found from the solution to the Euler-Lagrange equationdd t ˆ ∂λ∂ s ˙ = ∂λ∂s . (4.12)By substituting (4.10) into (4.12), we arrive at the second order nonlinear differentialequation b ( s ) : s + 2 α s α σ s µ = α σ µ , (4.13)with boundary conditions s (0) = R and s ( t f ) = R f . When α = 0, (4.13) reducesto : s = 0, the same equation derived by Dias et al. (2012) as expected. We solve (4.13)numerically and compute the optimal injection rate as Q ( t ) = 2 π s s b ( s ).6 L. C. Morrow, T. J. Moroney, and S. W. McCue
The result of this computation is presented in figure 9, where we compare the time-dependent injection rate for the tapered and parallel plate cases. This figure illustratesthat the optimal flow rate for the tapered geometry when the plates converge in thedirection of the flow (dotted green curve) is non-monotone in such a way that it is lowerthan the corresponding constant (solid blue curve) rate both when t = 0 and t = t f . Thischoice of injection rate acts to slow the speed of the interface for small time (analogous tothe linear injection rate (4.4)), while for later times it acts to prevent the rapid increasein speed that occurs when the interface reaches the region in which the gap betweenthe plates is smaller. For the diverging plates, figure 9 shows that the optimal injection(dash-dotted purple curve) rate is monotonically increasing such that the normal velocityof the interface is slowed down when the gap between plates is smallest in exchange fora faster injection rate for later times when the gap becomes larger.To illustrate the effect of implementing our new optimal injection rate when the platesare tapered, we include numerical results in rows four and six of figure 4. To ensure thecomparison with the previous three configurations is fair, we compute the solutions for thesame average flow rates (we use ¯ Q = 0 .
6, 1 .
2, 1 .
6, and 2 mL/s in columns 1-5, respectively)over the same period of time. When imposing a time-dependent injection rate while theplates linearly converge in the direction of the flow (row four), the interface appears toremain stable until near the very end of the simulation, where numerous stubby fingersform over the final time interval. It is interesting to note these fingers appear significantlyshorter than the fingers that develop for the other configurations, as we discuss below.In comparison, when the plates diverge in the direction of the flow, row 6 of figure 4indicates that the optimal injection rate appears to have less impact on the morphologyof the interface, both in regards to the number and length of fingers that develop.For the representative case ¯ Q = 2 mL/s (third column of figure 4), the isoperimetricratio as a function of time is shown in figure 5(b). Compared to the corresponding casewhere the injection rate is constant shown in figure 5( a ), we see that the optimal injectionrate acts to reduce the growth rate of I for small time for all three values of α . We seethat when the plates are converging and the injection rate is optimal (dashed red), theinterface remains essentially circular for almost all of the simulation. For late times,the isoperimetric ratio sharply increases as the interface appears to ‘switch’ from stableto unstable and short fingers (observed in the fourth row of figure 4) begin to develop.When compared to the corresponding constant injection case, imposing a time-dependentinjection rate decreases I ( t f ). For diverging plates (yellow dotted), we observe that theisoperimetric ratio appears to grow slowly over the first half of the simulation, and fasterover the second half, which is the opposite behaviour to the constant injection case. Whilethe time-dependent injection rate has resulted in a decrease in I ( t f ), this decrease is lesssubstantial compared to the converging plate configuration.In figure 6, we compare the isoperimetric ratio I and the circularity ratio C at t = t f asa function of ¯ Q for both the constant and optimal injection rates with α < đ ) and α > § ). As was observed for the constant injection rate case (yellow, İ )discussed in § I is less than the other configurations for¯ Q < . Q becomes large. This increase in I correspondsto the large number of fingers that develop for late times, as seen in fourth row of figure4. While the number of fingers significantly increases for large values of ¯ Q , the lengths ofthese fingers are short compared to the diverging case, and this is reflected in the valueof C , which is the lowest over the range of ¯ Q considered for each of the configurations. Aswe noted for the constant injection rate case in § ontrolling interfacial instabilities in Hele-Shaw flow Constant injectionOptimal injection
Figure 10.
Numerical solution to (2.2)-(2.7) with (left to right) ω = 0, 10, 20, 30, and 40g/(s · mL). Top row is for rotating plates with constant injection rate Q = 1 . R = 0 . R f = 5 cm. For all simulations, t f = 9 .
61 and b = 0 . sixth rows of figure 4, imposing the time-dependent injection rate does not appear tosignificantly impact either the number or length of fingers, and thus we see a relativelysmall reduction in I and C . We conclude that when compared to the correspondingconstant injection case, imposing a carefully chosen injection rate does result in a morecircular interface for both converging and diverging plates depending on the choice of ¯ Q ,but the interface can exhibit very different morphological features dependent on α . Wediscuss this issue further in §
6. 4.4.
Rotating plates
In this subsection, we now turn our attention to the case for which the gap between theplates is constant and the Hele-Shaw cell is rotated while the inviscid bubble is injected.By seeking solutions to (2.2)-(2.7) of the form (3.1) and (3.2) when b is a constant and ω >
0, we find that γ n γ n = ( n − s ˆ Q πbs − n ( n + 1) b σ µs ˙ − b nω µ , (4.14)and n max = d ˆ Qµs b πσ − ωs σ ˙ . (4.15)Thus, the centrifugal force acts as a stabilising term as it contributes a negative offsetto γ n /γ n and decreases the most unstable mode of perturbation. Furthermore, by notingthat γ n ( s ) = γ n (0) s n − exp ˆ b nπ (( n − σ − s ω )6 Qµs ˙ , (4.16)it follows that when n ≥ Q >
0, and ω > γ n ( s ) /s → s → ∞ . By comparison,when ω = 0, γ n ( s ) /s → ∞ as s → ∞ . This suggests that when the plates arerotating, there exists a critical radius where all modes of perturbation will be stable and8 L. C. Morrow, T. J. Moroney, and S. W. McCue
Figure 11.
Isoperimetric ratio of the numerical solution to (2.2)-(2.7) (top to bottom) ω = 0,10, 20, 30, and 40 g/(s · mL) for ( a ) constant and ( b ) time-dependent injection rate (4.17). Forall simulations, the same amount of is injected over time period where t f = 9 .
61 and the averageinjection rate is ¯ Q = 1 . R = 0 . R f = 5 cm. Shaded (blue) region represents one standard deviation above and below the mean.For all simulations, t f = 9 .
61 and b = 0 . the interface will become circular. Interestingly, this result is similar to case where theplates are stationary and the bubble is contracting ( ω = 0 and Q < n ≥ γ n ( s ) /s → s → + (Dallaston & McCue 2013).To illustrate the nonlinear behaviour of solutions to (2.2)-(2.7) when ω >
0, we performnumerical simulations for different values of ω , shown in the first row of figure 10. For ω = 10 g/(s · mL) (second column), we see that both the number and length of fingersthat develop is less than that for the case in which the plates are stationary (first column).For larger values of ω , we find that fingers initially develop; however, as time increases,the base of these fingers appear to be ‘pulled’ towards the finger tips, and in the case of ω = 40 g/(s · mL) (fifth column), the interface appears to be essentially circular at t = t f .We compute the corresponding isoperimetric ratio of these simulations, shown in figure11( a ). For cases where ω >
0, we find that while I initially increases as it does for when ω = 0, there exists a turning point (denoted by red dots) after which I monotonicallydecreases. This is consistent with the behaviour predicted by linear stability analysis,and thus our results suggest that when ω >
0, the interface will become circular after asufficient amount of time has passed.The explanation for why the centrifugal force causes the interface to become circularfor sufficiently large times relates to the dynamic boundary condition (2.5). This equationindicates that as the interface expands and fingers develop, the centrifugal force creates apressure differential between the base and tip of the fingers. From the kinematic boundarycondition (2.6), we see that this pressure differential acts to increase the normal velocityof the interface at the base of the finger. By comparison, when ω = 0 this pressuredifferential is absent and the normal velocity of the base fingers tends to be slower thanthat of the tips (see row one figure 4 for example). Furthermore, the normal velocityof the base of the finger increases linearly in r , suggesting that as the interface grows,the effect of the centrifugal force becomes stronger. Thus, while surface tension can bethought of as ‘penalising’ regions where curvature is high, the centrifugal term acts topenalise longer fingers.Regarding our objective of reducing the growth of viscous fingers when a prescribedamount of fluid is injected over a finite period of time, figure 11( a ) suggests that for aparticular Q , the interface will be essentially circular at t = t f for sufficiently large ω .Thus, we wish to determine the minimum value of ω that ensures the interface will becircular at t = t f . We perform a parameter sweep of ω for a particular value of Q , and ontrolling interfacial instabilities in Hele-Shaw flow Figure 12.
Stability diagram denoting the critical centrifugal parameter, ω c , such that I ( t f ) = 1 .
01 for constant ( • , blue) and time-dependent ( ◦ , red) injection rates. Theparameter ¯ Q (mL/s) is the average injection rate over a simulation and the final time is t f = ((99 π/
20) mL) / ¯ Q . The initial condition is of the form (4.3) with R = 0 . Q is (4.17) with R f = 5 cm. Additionally, b = 0 . determine the critical value of the centrifugal parameter, ω c , as the value of ω where I ( t f ) = 1 .
01, denoted in figure 12 by • (blue). For each value of Q considered, we areable to approximate ω c such that the interface is circular at the end of the simulationwhen ω > ω c . Of course, when ω < ω c , we expect the interface to become circular if thesimulations were run for a longer period of time. It is interesting to note that for thetapered plate configuration discussed in § α can be chosen such that I ( t f ) will be lessthan the corresponding parallel plate configuration (see figure 8 for example). However,our results suggest that there does not exist an analogous critical taper angle for every Q such that the interface will be completely stabilised over the duration of a simulation.Returning to time-dependent injection rates, as was discussed in § § et al. (2012) and the methodologypresented in § s satisfies a linear second-order differential equation withconstant coefficients, which reduces to : s = 0 when ω = 0. We can easily solve thisequation exactly and compute the flow rate via Q = 2 πbs s to give Q ( t ) = b πω µ ´ C e b ω µ t − C e − b ω µ t ¯ ´ C e − b ω µ t + C e b ω µ t ¯ , (4.17)where C = e b tf ω µ ˆ R e b tf ω µ − R f ˙ ˆ e b ωtf µ − ˙ − , C = R − C . (4.18)To illustrate the effect of implementing this injection rate, we compare fully nonlinearsimulations of our Hele-Shaw problem using (4.17) (second row of figure 10) with thecorresponding constant injection case (first row). Imposing (4.17) has the effect of bothreducing the number and length of viscous fingers that form, and the interface appearsto be completely stabilised for ω = 30 and 40 g/(s · mL). Regarding the isoperimetricratio (figure 11( b )), we find that, similar to the parallel and tapered plate configurationsdiscussed in § I compared to the corresponding constant injection rate (figure 11( a )), in exchange for afaster growth rate for times leading up to t f . We also note that I is now a monotonically0 L. C. Morrow, T. J. Moroney, and S. W. McCue
Figure 13.
Row 1 is experimental results from Zheng et al. (2015) where the plates areseparated according to b = b t / , each produced for a different value of the control parameterˆ J , reproduced with permission from the American Physical Society. Row 2 is the correspondingnumerical simulations with parameters µ = 0 .
95 Pa s, σ = 2 . , b = 8 . × − cm/s / ,and injection rate (left to right) Q = 0 .
09, 0 .
12, 0 .
17, and 0 .
18 mL/s. The initial condition forall simulations is (5.7). The scale bar denotes a length of 5 cm. increasing function in time, and the turning point observed when Q is constant is absent.As expected, implementing (4.17) reduces I ( t f ) compared to the corresponding constantinjection case, and in particular, for ω = 30 and 40 g/(s · mL), the interface is stabilisedover the entire duration of the simulation. Finally, we are able to determine the criticalcentrifugal parameter, ω c , (denoted as ◦ in figure 12) and, as expected, implementing(4.17) results in a reduction of ω c compared to the corresponding constant injection case(denoted with • ).
5. Controlling the number of fingers In §
4, we investigated our first objective for controlling the development of viscousfingers, which involved reducing the fingering pattern when injecting a prescribed amountof viscous fluid over a finite period of time. We now turn our attention to the secondobjective, which is to control the number of non-splitting fingers that develop in a Hele-Shaw cell. Numerous theoretical and experimental investigations have been performed todetermine strategies for controlling the number of viscous fingers. Using linear stabilityanalysis, Zheng et al. (2015) proposed that if the gap thickness and injection rates areof the form b = b t α b , Q = Q t α Q , where 7 α b − α Q = 1 , (5.1)then the expected number of fingers is N ≈ d J , (5.2)where ˆ J = 6 µQ / π / σ ( α Q + 1) / b / , (5.3) ontrolling interfacial instabilities in Hele-Shaw flow J = 107 J = 146 J = 191 J = 242 α b = / , α Q = α b = , α Q = − / α b = / , α Q = − / Figure 14.
Numerical simulations with distance between plates and injection rate of the form(5.1) for different values of α Q and α b . For row 1, Q = 1 . b = 0 . t / cm · s − / ,and (left to right) σ = 0 . . For row 2, Q = 1 . t − / mL/s / , b = 0 .
12 cm and (left to right) σ = 2 . . For row 3, Q = 1 . t − / mL/s / , b = 0 . t / cm · s − / , and (left to right) σ = 0 . . Forall simulations, the initial condition is (5.7). The scale bar represents a length of 5 cm. is a dimensionless control parameter. Li et al. (2009) showed using numerical simulationsthat when α Q = − / α b = 0, the interface will tend to N -fold symmetric shapesas time increases. Furthermore, the case α Q = 0 and α b = 1 / et al. (2015), who were able to produce interfaces with differentnumbers of non-splitting fingers. To date, however, configurations where both α Q (cid:54) = 0 and α b (cid:54) = 0 have not been considered either experimentally or numerically. In this section, weperform nonlinear numerical simulations to gain insight into the feasibility of controllingthe number of non-splitting fingers when imposing a time-dependent gap thickness and/orinjection rate according to (5.1). Note for all simulations in this section, ω = 0 g/(s · mL)and µ = 1 /
12 g/(cm · s).Performing linear stability analysis on the circular solution to (2.2)-(2.7), we find that γ n γ n = n − s ˆ Q πb − n ( n + 1) b σ µs ˙ − ( n + 1)2 b d b d t , (5.4)2 L. C. Morrow, T. J. Moroney, and S. W. McCue and the most unstable mode of perturbation is n max = d ˆ µQs b πσ − µs b σ d b d t ˙ . (5.5)By setting d b/ d t = 0, (5.4) and (5.5) reduce to (3.3) and (3.4). We can infer from (5.4)that increasing the gap between the plates, d b/ d t >
0, contributes a negative offsetto γ n /γ n , resulting in a stabilising effect. By considering s = b (cid:82) Q/πb d t such that s ( t ) (cid:29) s (0) and choosing Q and b of the form (5.1), it follows that n max will beindependent of time and equal to n max = c J J = (1 + α Q − α b )6 µQ / π / σ ( α Q + 1) / b / = 1 + α Q − α b α Q ˆ J (5.6)We note the discrepancy between J and ˆ J is due to Zheng et al. (2015) possibly ignoringthe non-homogeneous term in (2.4) such that (5.4) and (5.5) reduce to (3.3) and (3.4).However, we can see from (5.6) that its inclusion is significant when α b (cid:54) = 0.In addition to linear stability analysis, Zheng et al. (2015) also performed a series ofexperiments where the gap between the plates satisfies b ∝ t / for different constantinjection rates. In figure 13, we compare these experiments with the correspondingnumerical solution to (2.2)-(2.7). Simulations are performed with initial condition s ( θ,
0) = 1 + 5 × − (cid:88) n =2 cos( n ( θ − πθ n )) , (5.7)where θ n is a uniformly generated random number between 0 and 1. Simulations areperformed on the domain 0 ≤ r ≤
15 with 0 ≤ θ < π using 800 ×
355 equally spacednodes. The interface is evolved until the mean radius of the interface is 10 cm, which isapproximately the point at which the experiments by Zheng et al. (2015) are concluded.In these experiments, it was observed that the interface develops non-splitting fingers,and our numerical simulations reproduce this morphology. In addition to this, we alsoperform simulations, shown in figure 14, for different choices of α Q and α b that satisfy(5.1) with parameters chosen such that n max varies from 6 to 9. The injection rate andplate gap width are chosen of the form Q = Q ( t + t ) α Q and b = b ( t + t ) α b where t = 0 . Q = ∞ and b = 0 at t = 0. We see that for each configurationchosen, we are able to generate interfaces whose number of fingers compare well with thenumber predicted by linear theory. Finally, we emphasise again that we are deliberatelyrunning our simulations for roughly the same time-scales as Zheng et al. (2015) does intheir experiments; for much longer scales, obviously the Hele-Shaw model would breakdown for α b > J , and predictednumber of fingers from (5.6), we perform a series of numerical simulations with differentchoices of α Q and α b . For each combination of parameters, 10 simulations are performed,and the number of fingers at t = t f are averaged and illustrated in figure 15. Alsoshown is the closest integer to the most unstable mode (black line). We see that acrossthe values of J considered, the average number of fingers that develop is consistent foreach combination of α Q and α b . In comparison to linear stability analysis, this figureindicates agreement between (5.6) and numerical simulations for parameters that giverise to n max = 1, 5, 6, 7, and 8. However, for larger values of J , the number of fingersobserved from the numerical simulations is slightly less than the number predicted bylinear stability analysis. ontrolling interfacial instabilities in Hele-Shaw flow Figure 15.
Average number of fingers that develop as a function of the control parameter J .Injection rate and plate gap thickness are of the form (5.1) with: α Q = − / α b = 0 ( IJ , blue); α Q = 0 and α b = 1 / İ , red); and α Q = − / α b = 1 /
14 ( § , yellow). The parametersused are same as figure 14. The initial condition for all simulations is (5.7). The horizontal(black) lines denote the most unstable mode of perturbation approximated from linear stabilityanalysis by (5.2) rounded to the nearest integer. Inserts are examples of numerical simulationsfor different values of J . To further investigate this apparent discrepancy, we examine the behaviour of thesolution to (2.2)-(2.7) with parameters chosen such that n max = 10, shown in the firstrow of figure 16. For small time, we see 13 fingers developing (second column) and, astime increases, several of these fingers retract resulting in 10 fingers (third column),which is the same as that predicted by linear stability analysis. However, for later times,two of the fingers (denoted by a ∗ in the fourth column) do not appear to grow as fastas their neighbours. As a result, these fingers are ‘blocked off’ and retract, resulting inthe interface developing eight fingers by the end of the simulation (fifth column). Bycomparison, the second row of figure 16 shows the numerical solution with parameterschosen such that n max = 7. Again, fingers begin to grow at various rates; however, theinteraction between the fingers appears to be not as severe as in the first row and thusthe interface maintains seven fingers. From these observations, we infer that for largervalues of J , there are more fingers that compete with each other as they grow and, inturn, this competition can result in fewer fingers than that predicted by linear stabilityanalysis.
6. Discussion
We have conducted a numerical investigation into determining how manipulating thegeometry of the classic Hele-Shaw cell experiment can be used to control viscous fingeringpatterns. By utilising a numerical scheme based on the level set method, we have beenable to compute nonlinear numerical solutions both when the gap between the plates isspatially- and time-dependent as well in the case in which the plates are rotating. Asa preliminary test of our scheme, we have shown that our numerical solutions of (2.2)-(2.7) compare well with a variety of experimental results for different injection rates andplate configurations. Subsequently, we have been able to determine new relationshipsbetween these various manipulations and their influence on interfacial instabilities whichextends well beyond the limitations of linear stability analysis and previously performedexperiments. We summarise our findings below.In § L. C. Morrow, T. J. Moroney, and S. W. McCue n max = 10 n max = 7 Figure 16.
Evolution of numerical solution to (2.2)-(2.7) for different values of n max . For row1, b = 0 .
12 cm, Q = 1 . t − / mL/s / , and σ = 0 .
808 g/s . The ∗ refer to fingers that retractdue to competition with neighbouring fingers. For row 2, b = 0 .
12 cm, Q = 1 . t − / mL/s / ,and σ = 1 .
654 g/s . Solutions for both rows are shown at times (left to right) t = 0, 2.3, 6.518.7 and 48.4 s, and initial condition is (5.7). We note that solutions here have been scaled suchthat the average radius of the interface is 1 cm. amount inviscid fluid is injected over a finite time interval. In particular, we investigatedhow imposing a time-dependent injection rate and/or tapering the plate gap (in eitherthe converging or diverging configurations) influences the morphology of the interface. Byperforming a series of numerical simulations and applying standard metrics for measuringhow round the interface is, we have shown that each of these configurations is able toproduce a less unstable interface at t = t f (for certain taper angles) than the standardconfiguration with parallel plates and constant injection rate. In other words, we canreduce viscous fingering by either tapering plates in the converging ( α >
0) or diverging( α <
0) directions and suppress the fingering further by imposing a time-dependentinjection rate that is chosen to minimise the growth rate of the most unstable mode ofperturbation (from linear stability analysis). Of all of these strategies, our results indicatethat injecting at a wisely chosen time-dependent rate when α > t ∼ t f for the converging configuration( α > § Q increases, the isoperimetric ratio I increases at a rate that is higherthan the other configurations, while the increase in circularity ratio C is not so dramatic.Thus for higher flow rates we are observing interfaces that are highly complex but whosefjords are not as deep as in the standard pattern.Another example relates to the converging case with optimal injection (fourth row offigure 4). Here, the morphology is significantly different from the standard case, and infact this interfacial pattern resembles the short flat-tipped “stubby” fingers observed bothexperimentally and numerically by Pihler-Puzovi´c et al. (2012, 2013), who considered aHele-Shaw cell where the top plate is replaced by an elastic membrane. Other closely ontrolling interfacial instabilities in Hele-Shaw flow et al. (2013); Duclou´e et al. (2017); Juel et al. (2018); Lister et al. (2013); McCue (2018); Pihler-Puzovi´c et al. (2014, 2018). These observations suggest there is a one-parameter family of solutions ( α )joining the standard Hele-Shaw problem to one where the pattern formation is similarto that produced by a deformable boundary. This is perhaps not surprising as the elasticmembrane acts like a tapered upper boundary near the interface.The morphological features for diverging plates are also interesting. Here, we noticethe interface appears significantly different from both the parallel and converging cases,with fewer and longer fingers forming over the duration of each simulation. In addition,implementing our ‘optimal’ (time-dependent) injection rate has little observable effecton this morphology.In § ω >
0. This is somewhat of an unusualresult, since for all the configurations considered in § §
5, we performed nonlinear simulations to determine whether imposing a time-dependent injection rate and/or plate gap can be used to control the number of fingerswhich develop. In particular, following the suggestion by Zheng et al. (2015), we allowedthe injection rate and gap thickness to vary according to power laws in time, withexponents α Q and α b , respectively. For a range of combinations of α Q , α b and the controlparameter J , we tested the hypothesis that (after an initial period in which various modesof perturbation grow or decay) a fixed number of non-splitting fingers emerge, n max ,which is equal to the most unstable mode. Our results support this hypothesis for thecases which predict between 5 and 8 fingers, as well as the less interesting case in whichthere are no fingers (a stable interface). For larger values of n max , the average number offingers observed in our simulations is slightly less than that predicted by linear stabilitytheory. In this parameter regime, our explanation for observing fewer than n max fingersis that, on the time-scale of our simulations, it appears that when n max is sufficientlylarge, nonlinear interactions between closely packed fingers can cause a small numberof them to retract (see Figure 15). These numerical simulations are consistent with theexperimental results of Leshchiner et al. (2010), who test a time-dependent injection rate( α Q = − / α b = 0) for a single control parameter, and Zheng et al. (2015) whotreat the lifting-plate case ( α Q = 0 and α b = 1 /
7) in some detail. Further, our findingsfor the lifting-plate case ( α Q = 0 and α b = 1 /
7) are in agreement with the very recentnumerical study of Vaquero-Stainer et al. (2019), who run their simulations for a largerrange of the control parameter J , namely up to J = 1000.From a theoretical perspective, beyond the lifetime of a normal experiment, there is aquestion about the ultimate long-time behaviour of our mathematical solutions. It wasshown numerically by Li et al. (2009) that, for the special case α Q = − / α b = 0(stationary plates with Q ∼ t − / ), the interface develops N -fold symmetry independent6 L. C. Morrow, T. J. Moroney, and S. W. McCue of the initial condition over an extremely long time period, at least for values of J whichpredict up to nine non-splitting fingers. The likely reason for this N -fold symmetriclong-time attractor is that the problem with α Q = − / α b = 0 has self-similarsolutions of the form p = t − / P ( X, Y ), where (
X, Y ) = ( x/t / , y/t / ) (Ben Amar et al. et al. (2019) alsoobserves time-dependent solutions approaching a self-similar form for the case α Q = 0and α b = 1 / N -fold symmetric. In general, it seems there are self-similar solutions of theform p = − σ/b + t − (2 α b +1) / P ( X, Y ), where (
X, Y ) = ( x/t (2 α b +1) / , y/t (2 α b +1) / ) and7 α b − α Q = 1, although these have not been computed before. For this combinationof parameters, at this stage it is not clear whether the bubble evolves to an perfect N -fold symmetric shape with N fingers predicted by (5.5) for randomly chosen modes ofperturbation in the initial condition (except for the special case α Q = − / α b = 0).Either way, these questions are worthy of further enquiry. Acknowledgements
The authors acknowledge the support of the Australian Research Council via theDiscovery Project DP140100933, as well as the computational resources provided by theHigh Performance Computing and Research Support Group at QUT. We thank AnneJuel and Draga Pihler-Puzovi´c for helpful discussions. Finally, we are grateful to theanonymous referees for their detailed reviews and insight.
REFERENCESAl-Housseiny, T. T., Christov, I. C. & Stone, H. A.
Phys. Rev. Lett. , 034502.
Al-Housseiny, T. T. & Stone, H. A.
Phys. Fluids , 092102. Al-Housseiny, T. T., Tsai, P. A. & Stone, H. A.
Nat. Phys. , 747. Alvarez-Lacalle, E, Ortın, J & Casademunt, J
Phys. Fluids , 908–924. Anjos, P H A, Alvarez, V M M, Dias, E O & Miranda, J A
Phys. Rev. Fluids , 124003. Anjos, P H A, Dias, E O & Miranda, J A
Phys. Rev. Fluids ,124004. Ben Amar, M., Hakim, V., Mashaal, M. & Couder, Y.
Phys. Fluids A , 1687–1690. Ben-Jacob, E. & Garik, P.
Nature , 523.
Ben-Jacob, E., Shmueli, H., Shochet, O. & Tenenbaum, A.
Physica A , 378–424.
Bongrand, G. & Tsai, P. A.
Phys. Rev. E , 061101. Brener, E.A., Kessler, D.A., Levine, H. & Rappei, W.J. ◦ geometry. Euro. Phys. Lett. , 161. Carrillo, L, Magdaleno, F X, Casademunt, J & Ort´ın, J
Phys. Rev. E , 6260. Carrillo, L, Soriano, J & Ortın, J
Phys. Fluids , 778–785. Chen, C. Y., Chen, C. H. & Miranda, J. A.
Phys. Rev. E , 056304. ontrolling interfacial instabilities in Hele-Shaw flow Chen, J.-D.
Exp. Fluids , 363–371. Chen, S., Merriman, B., Osher, S. & Smereka, P.
J. Comput. Phys. , 8–29.
Combescot, R. & Ben Amar, M.
Phys. Rev. Lett. , 453. Dai, W.-S. & Shelley, M. J.
Phys. Fluids A , 2131–2146. Dallaston, M. C. & McCue, S. W.
Nonlinearity , 1639–1665. DeGregoria, A. J. & Schwartz, L. W.
J. Fluid Mech. , 383–400.
Dias, E. O., Alvarez-Lacalle, E., Carvalho, M. S. & Miranda, J. A.
Phys. Rev. Lett. , 144502.
Dias, E. O. & Miranda, J.
Phys. Rev. E , 016312. Dias, E. O. & Miranda, J. A.
Phys. Rev. E , 053015. Dias, E. O., Parisio, F. & Miranda, J. A.
Phys. Rev. E , 067301. Duclou´e, L., Hazel, A. L., Pihler-Puzovi´c, D. & Juel, A.
J. Fluid Mech. . Enright, D., Fedkiw, R., Ferziger, J. & Mitchell, I.
J. Comput. Phys. , 83–116.
Fast, P. & Shelley, M. J.
J. Comput. Phys. , 117–142.
Gadˆelha, H & Miranda, J
Phys. Rev. E , 066308. Homsy, G. M.
Ann. Rev. Fluid Mech. , 271–311. Hou, T. Y., Li, Z., Osher, S. & Zhao, H.
J. Comput. Phys. , 236–252.
Jackson, S. J., Power, H., Giddings, D. & Stevens, D.
Comput. MethodsAppl. Mech. Eng. , 606–632.
Juel, A., Pihler-Puzovi´c, D. & Heil, M.
Ann. Rev. FluidMech. , 691–714. Leshchiner, A., Thrasher, M., Mineev-Weinstein, M. B. & Swinney, H. L.
Phys. Rev. E , 016206. Li, S., Lowengrub, J.S., Leo, P. H. & Cristini, V.
J. Cryst. Growth , 703–713.
Li, S., Lowengrub, J. S., Fontana, J. & Palffy-Muhoray, P.
Phys. Rev. Lett. , 174501.
Liang, S.
Phys. Rev. A , 2663. Lindner, A., Derks, D. & Shelley, M. J.
Phys. Fluids , 072107. Lins, T. F. & Azaiez, J
J. Fluid Mech. , 713–729.
Lister, J. R., Peng, G. G. & Neufeld, J. A.
Phys. Rev. Lett. , 154501.
Lu, D., Municchi, F. & Christov, I. C. arXiv preprint 1811.06960 . McCue, S. W.
J. Fluid Mech. , 1–4.
McLean, J W & Saffman, P G
J. Fluid Mech. , 455–469.
Mirzadeh, M. & Bazant, M. Z.
Phys. Rev.Lett. , 174501. L. C. Morrow, T. J. Moroney, and S. W. McCue
Moroney, T. J., Lusmore, D. R., McCue, S. W. & McElwain, S.
J. Comput. Phys. , 170–185.
Mullins, W. W. & Sekerka, R.F.
Dynamics of Curved Fronts , pp. 345–352. Elsevier.
Nase, J., Derks, D. & Lindner, A.
Phys. Fluids , 123101. Osher, S. & Sethian, J. A.
J. Comput. Phys. , 12–49. Paterson, L.
J. Fluid Mech. , 513–529.
Pihler-Puzovi´c, D., Illien, P., Heil, M. & Juel, A.
Phys. Rev.Lett. , 074502.
Pihler-Puzovi´c, D, Juel, A. & Heil, M.
Phys. Fluids , 022102. Pihler-Puzovi´c, D., Peng, G. G., Lister, J. R., Heil, M. & Juel, A.
J. Fluid Mech. , 163–191.
Pihler-Puzovi´c, D., P´erillat, R., Russell, M., Juel, A. & Heil, M.
J. Fluid Mech. ,162–183.
Rabbani, H. S., Or, D., Liu, Y., Lai, C.-Y., Lu, N. B., Datta, S. S., Stone, H. A. &Shokri, N.
Proc. Natl.Acad. Sci. p. 201800729.
Saffman, P. G. & Taylor, G. I.
Proc. R. Soc. Lond. A , 312–329.
Shelley, M. J., Tian, F. & Wlodarski, K.
Nonlinearity , 1471. Stone, H. A.
Phys. Rev.Fluids , 100507. Vaquero-Stainer, C, Heil, M, Juel, A & Pihler-Puzovi´c, D
Preprint: arXiv:1903.00903 . Witten, T. A. & Sander, L. M.
Phys. Rev. B , 5686. Zheng, Z., Kim, H. & Stone, H. A.
Phys. Rev. Lett. , 174501.
Appendix A. Numerical scheme
A.1.
The level set method
To implement the level set method, a level set function, φ ( x , t ), is constructed asa signed distance function whose zero level set describes the location of the interfacebetween the viscous and inviscid fluid regions, and φ > x ∈ R \ Ω( t ) , (A 1) φ < x ∈ Ω( t ) . (A 2)If the interface has a normal speed V n , then we wish to construct a function, F , suchthat F = V n on the interface and is continuous over the entire computational domain.Thus φ satisfies the level set equation ∂φ∂t + F |∇ φ | = 0 . (A 3)We approximate the spatial derivatives in (A 3) using a second order essentially non-oscillatory scheme, and integrate in time using second order Runge-Kutta. We choose atime step size of ∆ t = 0 . × ∆ x/ max | F | to ensure stability. Furthermore, to maintain ontrolling interfacial instabilities in Hele-Shaw flow φ as a signed distance function, re-initialisation is periodically performed by solving ∂φ∂τ + S ( φ )( |∇ φ |−
1) = 0 , (A 4)to steady state where S ( φ ) = φ a φ + ∆ x , (A 5)and τ is a pseudo time variable.A major limitation of the level set method is that solutions can suffer from volumeloss (or gain). To alleviate this problem, we implement the particle level set method,which combines the Eulerian level set method with a marker particle based Lagrangianapproach. The particle level set method, first proposed by Enright et al. (2002), extendsthe traditional level set method by placing massless marker particles around the interface.These particles are advected using the same velocity field as the level set function. Asthe particles do not suffer from mass loss, the level set function can be corrected if theparticles are found to incorrectly cross the interface. We refer the reader to Enright et al. (2002) for details on how to implement the particle level set method, as well as examplesillustrating its effectiveness. A.2. Solving for F Defining n = ∇ φ/ |∇ φ | as the outward facing normal of the interface, we have theexpression F = − b µ ∇ p · ∇ φ |∇ φ | x / ∈ Ω( t ) . (A 6)This satisfies (2.6) on the interface, and provides a continuous expression for F in theviscous fluid region. However, to solve (A 3) we require an expression for F over the entirecomputational domain. Moroney et al. (2017) proposed that the speed function can beextended into the inviscid fluid region by solving the biharmonic equation ∇ F = 0 x ∈ Ω( t ) . (A 7)By solving (A 7), this ensures that F = V n on the interface and gives a continuousexpression for F away from the interface. The sign of φ is used to determine nodes insidethe interface that need to be included in the biharmonic stencil. As such, the location ofthe interface does not need to be known explicitly, similar to the level set method itself.This velocity extension process is a variant of a thin plate spline in two dimensions. Werefer the reader to Moroney et al. (2017) for further details.A.3. Solving for pressure
To evaluate the speed function F , we must first compute the pressure field. We consider(2.2)-(2.7) in polar coordinates with p = p ( r, θ, t ) and the location of the interface is givenby r = s ( θ, t ). Thus (2.4) becomes1 r ∂∂r ˆ rb µ ∂p∂r ˙ + 1 r ∂∂θ ˆ b µ ∂p∂θ ˙ = ∂b∂t r > s ( θ, t ) . (A 8)In order to solve for the pressure at nodes that are not adjacent to the interface, wediscretise (A 8) using a standard central finite difference scheme. Denoting β = rb / µ ,0 L. C. Morrow, T. J. Moroney, and S. W. McCue the r -derivatives in (A 8) are approximated via1 r ∂∂r ˆ β ∂p∂r ˙ → r i,j ∆ r ˆ β i +1 / ,j p i +1 ,j − p i,j ∆ r − β i − / ,j p i,j − p i − ,j ∆ r ˙ , (A 9)where β i +1 / ,j = ( β i +1 ,j + β i,j ) / β i − / ,j = ( β i − ,j + β i,j ) /
2. The derivatives in the θ -direction are discretised in a similar fashion.Special care must be taken when solving for nodes adjacent to the interface. Supposethat the interface is located at r = r I where r i − ,j < r I < r i,j where the nodes r i − ,j and r i,j are in the inviscid and viscous fluid regions respectively. When discretising (A 8),we can no longer incorporate p i − ,j into our finite difference stencil as it is not in thedomain x ∈ R \ Ω( t ). Instead, we define a ghost node at r I whose value is p I . By notingthat φ is a signed distance function, the distance between r i,j and r I is computed via h = ∆ r ˇˇˇˇ φ i,j φ i − ,j − φ i,j ˇˇˇˇ . (A 10)As per Chen et al. (1997), our finite difference stencil becomes1 r ∂∂r ˆ β ∂p∂r ˙ → r i,j (∆ r + h ) ˆ β i +1 / ,j p i +1 ,j − p i,j ∆ r − ˆ β i − / ,j p i,j h ˙ + non-homog. hkkkkkkkkkkkkkkikkkkkkkkkkkkkkj r i,j h (∆ r + h ) ˆ β i − / ,j p I . (A 11)Here ˆ β i − / ,j = ( β i,j + β I ) / β I is the value of β on the interface, and is computedvia linear interpolation using β i,j and β i − ,j . When the node and interface are sufficientlyclose together such that h < ∆ r , we set p i,j = p I . A similar procedure is applied if theinterface lies between r i,j < r I < r i +1 ,j and in the azimuthal direction. The value of p I is computed from the dynamic boundary condition (2.5), where the curvature of theinterface is κ = ∇ · n .A.3.1. Far-field boundary condition
To incorporate the far-field boundary condition (2.7) into our finite difference stencil,we utilise a Dirichlet to Neumann map. This is implemented by imposing an artificialcircular boundary at r = R such that R > s ( θ, t ). By only considering the region indomain r ≥ R , we seek a solution to (2.4) of the formˆ p ( r, θ, t ) = A − Q π log r + r ∂b∂t + ∞ (cid:88) n =1 r − n p A n cos nθ + B n sin nθ q , (A 12)where A , A n , and B n are unknown, and ˆ p = ( b / µ ) p . The expansion (A 12) assumesthat b is spatially uniform in r ≥ R , and the choice of linearly tapered plate gap (4.6) isconsistent with this. Considering the value of pressure on the artificial boundary, supposethat ˆ p ( R, θ, t ) can be represented as a Fourier seriesˆ p ( R, θ, t ) = a + ∞ (cid:88) n =1 a n cos nθ + b n sin nθ, (A 13) ontrolling interfacial instabilities in Hele-Shaw flow a = 12 π (cid:90) π ˆ p ( R, θ, t ) d θ, (A 14) a n = 1 π (cid:90) π ˆ p ( R, θ, t ) cos nθ d θ, (A 15) b n = 1 π (cid:90) π ˆ p ( R, θ, t ) sin nθ d θ. (A 16)By equating (A 13) with (A 12) evaluated at r = R , we find that A = a +( Q/ π ) log R − bR / A n = R n a n and B n = R n b n .We differentiate our expression for ˆ p with respect to r and evaluate it at r = R to give ∂∂r ˆ p ( R, θ j ) = − Q πR + R ∂b∂t − ∞ (cid:88) n =1 nR p a n cos nθ j + b n sin nθ j q , (A 17) ≈ − Q πR + R ∂b∂t − ∆ θRπ m (cid:88) k =1 w jk ˆ p ( R, θ k ) , (A 18)where w jk = ∞ (cid:88) n =1 n cos( n ( θ j − θ k )) . (A 19)Defining I as the outermost index at which r = R , then our expression for ∂p/∂r isincorporated into our finite difference stencil, (A 20)1 r ∂∂r (cid:18) β ∂p∂r (cid:19) → R ∆ r (cid:40) − β I − / ,j p I,j − p I − ,j ∆ r + R (cid:34) − Q πR + R ∂b∂t − b µ ∆ θRπ m (cid:88) k =1 w jk p I,k (cid:35)(cid:41) , recalling β = rb / µ . The finite difference stencil for the derivatives in the θ -direction isnot changed. Furthermore, a similar procedure to the one presented here could be usedto model the Dirichlet boundary condition p ∼ p ∞ as r → r where p ∞ is prescribed.A.4. Numerical validation
To establish that our numerical scheme converges as the grid is refined, we considerthe initial condition (in cm) s ( θ,
0) = 0 . . θ, (A 21)where 0 ≤ θ ≤ π . Simulations are performed by employing an increasingly refined meshwith a set of parameter values that are in the range of those used elsewhere in this study.These simulations, shown in figure 17, indicate that when the grid is sufficiency refined,both the size and shape of the fingers that develop are unchanged, and convergenceappears to be achieved using 750 ×
628 equally spaced nodes. Furthermore, the bubbleappears to maintain three-fold symmetry over the duration of the simulation.We also wish to determine that our numerical scheme is able to accurately describe thebehaviour of the interface for the different plate configurations considered in this article.To do so, we perform simulations where the interface is initially a circle of radius 0.5 cm,2
L. C. Morrow, T. J. Moroney, and S. W. McCue
Figure 17.
Convergence test of numerical scheme for the evolution of a bubble with initialcondition (A 21) and Q = 1 .
25 mL/s, b = 0 . µ = 1 g/(cm · s), σ = 1 . , and t f = 12 . ≤ r ≤ . ≤ θ < π . Figure 18.
Comparison of the numerical solution to (2.2)-(2.7) (solid blue) with solutionto (A 22) (dashed red), where the gap between the plates is ( a ) b = 0 . b ) b = 0 . .
04 sin(4 πt/t f ) cm, and ( c ) of the form (4.6) with b = 0 .
005 cm, r = 7 cm,and α = 0 . t f = 9 .
81 s and the injection rate is Q = 1 . ≤ r ≤ . ≤ θ < π with µ = 1g/(cm · s), σ = 2 g/s and 750 ×
628 equally spaced nodes. and compare the evolution of the radius with the solution to s = Q πb ( s ) s − s b ( s ) ∂b∂t . (A 22)We consider three configurations. The first is the classic configuration in which the platesare parallel and stationary. The second involves parallel plates with the distance betweenthe two plates evolving according to (in cm) b ( t ) = 0 . .
04 sin ˆ πtt f ˙ ..