A review of one-phase Hele-Shaw flows and a level-set method for non-standard configurations
Liam C. Morrow, Timothy J. Moroney, Michael C. Dallaston, Scott W. McCue
AA REVIEW OF ONE - PHASE H ELE –S HAW FLOWS AND ALEVEL - SET METHOD FOR NON - STANDARD CONFIGURATIONS
A P
REPRINT
Liam C. Morrow , , Timothy J. Moroney , Michael C. Dallaston and Scott W. McCue Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, United Kingdom School of Mathematical Sciences, Queensland University of Technology, Brisbane, QLD, 4001, AustraliaJanuary 20, 2021 A BSTRACT
The classical model for studying one-phase Hele–Shaw flows is based on a highly nonlinearmoving boundary problem with the fluid velocity related to pressure gradients via a Darcy-typelaw. In a standard configuration with the Hele–Shaw cell made up of two flat stationary plates,the pressure is harmonic. Therefore, conformal mapping techniques and boundary integralmethods can be readily applied to study the key interfacial dynamics, including the Saffman–Taylor instability and viscous fingering patterns. As well as providing a brief review of thesekey issues, we present a flexible numerical scheme for studying both standard and non-standardHele–Shaw flows. Our method consists of using a modified finite difference stencil in con-junction with the level set method to solve the governing equation for pressure on complicateddomains and track the location of the moving boundary. Simulations show that our methodis capable reproducing the distinctive morphological features of the Saffman–Taylor instabil-ity on a uniform computational grid. By making straightforward adjustments, we show howour scheme can easily be adapted to solve for a wide variety of configurations, including caseswhere the gap between the plates is linearly tapered, the plates are separated in time, and theentire Hele–Shaw cell is rotated at a given angular velocity.
Keywords:
Hele–Shaw flow; Saffman–Taylor instability; viscous fingering patterns; movingboundary problem; conformal mapping; level-set method.
Viscous fingering patterns that form in a Hele–Shaw flow are very well studied in fluid dynamics. These patterns,which arise due to the Saffman–Taylor instability [112], occur when a viscous fluid that fills a gap between twonarrowly separated parallel plates is displaced by a less viscous fluid, which is injected into the cell. Providedthese two fluids are immiscible, an interface forms which is usually unstable and develops visually striking pat-terns characterised by their branching morphology. As the governing equation for the velocity of the viscous a r X i v : . [ phy s i c s . f l u - dyn ] J a n review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT fluid is the same as Darcy’s law, Hele–Shaw flow can be thought of as a model for flow through a homoge-neous porous medium. Further, the Hele–Shaw framework has also been used to model interfacial instabilitiesappearing in other scenarios including bacterial colony growth [12], crystal solidification [74], random walks,[77], and the flow of electrolytes [47, 89]. We refer the reader to [17, 58, 82, 130] for a historical summary andcomprehensive overview of the broad applicability of the Hele–Shaw model.If we assume that the viscosity of the less viscous fluid (air, say) can be neglected entirely, then the classicalmodel for flow in the more viscous fluid, Ω( t ) , is the one-phase moving boundary problem v = − b µ ∇ p, ∇ · v = 0 , x ∈ Ω( t ) , (1)where v is the fluid velocity (averaged across the distance b between the parallel Hele–Shaw plates), p is thefluid pressure and µ is the constant viscosity, together with the boundary conditions p = − γκ + constant , v n = − b µ ∂p∂n , x ∈ ∂ Ω( t ) , (2)where γ is the surface tension, κ is the signed curvature of ∂ Ω , defined to be positive if the interface is convexfrom the fluid side, and v n is the normal speed of the interface. Typically the flow is driven by injection orwithdrawal of fluid through a point or at infinity. This original model for two immiscible fluids is described bySaffman & Taylor [112] in 1958, except that for the most part they kept track of both fluids.The one-phase Hele–Shaw model that we are concerned with has been applied to three main configurations,namely an expanding bubble of air into an infinite body of fluid, a contracting finite blob of fluid, and thedisplacement of viscous fluid in a Hele–Shaw channel. In each of these three scenarios, the fluid boundary isunstable (the Saffman–Taylor instability), and a typical outcome involves portions of the interface propagatingincreasingly faster than other portions and in some cases leading to a striking fingering pattern at the boundary.For the special zero-surface tension case (also known as Laplacian growth), a host of mathematical studies basedmostly on conformal mapping, conserved moments and the Baiocchi transform have highlighted the possiblescenarios for this ill-posed model, including exact solutions and finite-time blow-up [28–30, 33, 56, 57, 64–66, 72, 86, 87]. For the more physically realistic nonzero surface tension case (which is well-posed), the broaderstrategies to study this problem include stability analysis [88, 101], small-surface-tension asymptotics [18, 122],employing harmonic moments and conserved quantities [73] and fully numerical methods mostly with boundaryintegral methods [19, 32, 61, 63, 70, 95] but also the level set formulation [62]. While we shall devote much ofour attention in this article to certain non-standard variations of (1)-(2), we do not provide any commentary onhow the boundary conditions (2) may be altered by considering additional physical effects on the boundary apartfrom surface tension, including the effect of a dynamic contact angle and the related issue of kinetic undercooling[5, 6, 21, 34, 35, 40, 100, 106, 111, 132]. Similarly, we do not review non-Newtonian flows, which themselvesare well-studied [9, 43, 45, 71, 84, 110]. Finally, our focus is on time-dependent problems and so we are notintending to review the extensive literature on travelling wave problems involving a steadily propagating finger[20, 26, 48, 49, 60, 85, 111, 117, 120, 126] or bubble [54, 59, 81, 119, 121, 123, 128, 129].In recent years, there has been increased interest in studying how variations to the classic Hele–Shaw modelinfluence the development of viscous fingering patterns. Many of these studies consider the effect of imposinga time-dependent injection rate, specifically to control or reduce the growth of fingers [10, 11, 27, 37–39, 76].Further, much attention has been devoted to manipulating the geometry of the Hele–Shaw cell. One of the earliestexamples of this approach is by Zhao et al. [133], who considered the classic Saffman–Taylor experiment [112]and linearly tapered the gap between the plates in the direction of the fluid flow. Since this experiment, numerous2 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT studies have been performed to give further insight into how the taper angle influences viscous fingering [1, 2,13, 69, 93]. Other popular physical alterations to the Hele–Shaw cell include uniformly separating the plates intime [42, 78, 94, 116, 127, 135], rotating the entire Hele–Shaw cell at a given angular velocity [7, 16, 113], orreplacing one of the plates with an elastic membrane [3, 31, 80, 83, 102–105]. All of these configurations havebeen shown to produce patterns distinct from traditional Saffman–Taylor fingering.One of the most commonly used analytical tools for studying both standard and non-standard Hele–Shaw flow islinear stability analysis. For the standard configuration, Paterson [101] showed that modes of perturbation to thecircular solution becomes successively unstable as the bubble expands, predicting the most unstable wave numberfor a given bubble radius. Further, linear stability analysis has also been used to derive injection rates to control[75] or minimise [39] the development of viscous fingering. For non-standard Hele–Shaw flow, linear stabilityanalysis gives insight into how manipulating the geometry of the cell influences the development of viscousfingers, including when the plates are tapered [1, 2], rotating [16], or are being separated [116]. While linearstability analysis is a flexible tool that leads to analytic predictions, it only provides an accurate description ofsolutions for small time. As such, this restriction increases the need for flexible and accurate numerical methodsthat can be used to understand the full nonlinear behaviour of these problems.Computing numerical solutions to Hele–Shaw flow (and related moving boundary problems) can be a challengingtask, as interfacial patterns develop which then requires solving the governing equations in complicated movingdomains. Such approaches can be classified as either front tracking, where the interface is solved for explicitly, orfront capturing, where the interface is represented implicitly. For the classic Hele–Shaw problem, as the pressureof the viscous fluid satisfies Laplace’s equation, the most popular choice is the boundary integral method (alsoknown as the boundary element method), which is classified as a front tracking method. In particular, theboundary integral method has been used to solve the classic one-phase Hele–Shaw problem [32, 75, 76], as wellas two-phase flow [68, 108], time-dependent plate gap [116, 134], and Hele–Shaw flow in channel geometry[36]. However, for non-standard Hele–Shaw configurations, the pressure may no longer be harmonic and theboundary integral method becomes a less suitable tool. Another disadvantage of front tracking methods is thatthe mesh may need to be regenerated as the interface evolves, in which case care must be taken to avoid meshdistortion effects.A popular alternative to the boundary integral method is the level set method, which represents the interfaceimplicitly as the zero level set of a higher dimensional hypersurface [98]. A commonly cited advantage of thelevel set method is that it can easily handle complex interfacial behaviour such as the merging and splittingof interfaces. Another, more pertinent, advantage of the level set method is that it can describe the formationof complex interfacial patterns (such as those that occur in Hele–Shaw flow) on a uniform grid, eliminatingthe need to re-mesh as the interface evolves. One of the most significant drawbacks of the level set method isthat it can suffer from mass loss/gain in regions where the mesh is under resolved. However, this issue can bemitigated by using the particle level set method [41], which uses massless marker particles to correct the locationof the interface when mass is lost/gained. The level set method is a popular tool for studying moving boundaryproblems in fluid dynamics, and has been used to investigate interfacial instabilities that occur in Stefan problems[24, 51] and Hele–Shaw flow [62, 79]. We refer the reader to [52, 97, 115] for more information about the levelset method, including details regarding implementation and applications.While the level set method can be used to implicitly represent the location of the interface, to numericallysimulate Hele–Shaw flow we are also required to determine the pressure within the viscous fluid, which involvessolving a partial differential equation in a complicated domain that changes in time. When applying the boundary3 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT integral method for the classic Hele–Shaw problem, the solution to Laplace’s equation can be expressed interms of Green’s functions. As such, the problem is reformulated as an integro-differential equation, and nodesneed only be placed on the fluid-fluid interface. An alternative choice is to solve for the pressure using thefinite difference method, which can be modified to solve problems on complex domains when coupled withlevel set functions that describe the location of the interface [25, 50]. An advantage of this approach is thatthe finite difference method can be easily adapted to wide classes of partial differential equations. Further,while the boundary integral method can easily handle non-trivial far-field boundary conditions, their inclusioninto the finite difference stencil is not so straightforward. One solution to overcome this is to use a very largecomputational domain, but this in turn results in significantly longer computational times. Another, more elegant,solution is to use a Dirichlet-to-Neumann map [53], which has been shown to accurately capture the far-fieldboundary condition even when the interface is relatively close to the curve on which the Dirichlet-to-Neumannmap is applied [91–93].In this work, we provide a brief review of the one-phase Hele–Shaw model, touching on the mathematical con-sequences of including or excluding surface tension in the model. We focus on the three well-studied scenarios,namely an expanding bubble, a contracting blob and displacement of fluid in a linear channel. Our approachis to write down a generalised model that allows for a number variations of the standard approach, including atime-dependent flow rate, a spatially and/or time-dependent gap between the plate, or rotating plates. We thenpresent a flexible numerical framework for solving this generalised model. Our scheme is based on the work ofChen [23] and Hou et al. [62], and uses a level set based approach to track the location of the liquid-air interface.There are several novel aspects of our numerical framework. The first is that our scheme overcomes the limi-tations of the boundary integral method in that it can easily solve Hele–Shaw flow in non-homogeneous media,i.e. where the plates are not parallel. Second, by representing the interface implicitly by a higher dimensionallevel set function, we are able represent the complicated interfacial patterns easily on a uniform mesh. By per-forming a series of simulations, we show that our numerical solutions are able to reproduce the morphologicalfeatures of viscous fingering in a Hele–Shaw cell. Further, by making straightforward adjustments, we show thatour scheme can easily be modified for a wide range of non-standard Hele–Shaw configurations, including wherethe plates are linearly tapered, uniformly separated in time, or rotated. For all the configurations considered, ournumerical solutions are shown to compare well with previous simulations and experiments.
We consider a generalised model of Hele–Shaw flow where the gap between the plates is either spatially ortemporally dependent such that b → b ( x , t ) and the Hele–Shaw plates can rotate with angular velocity ˆ ω .Suppose an inviscid bubble is injected into the viscous fluid at rate Q ( t ) , and we denote the domain occupiedby the inviscid fluid to be Ω( t ) , and the interface separating the two immiscible fluids ∂ Ω( t ) . We consider adepth averaged model which comes about from averaging Stokes flow over the gap between the plates, which isassumed to be small. 4 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Denoting P , µ and ρ as the pressure, viscosity and density of the viscous fluid, respectively, and denoting theangular velocity of the Hele–Shaw cell by ˆ ω , the governing equations for the velocity of the viscous fluid are v = − b µ ( ∇ P − ˆ ω ρr e r ) , x ∈ R \ Ω( t ) , (3) ∇ · ( b v ) = − ∂b∂t , x ∈ R \ Ω( t ) , (4)where r = | x | . Equation (3) is Darcy’s law modified to include the rotational effects of the Hele–Shaw cell, while(4) ensures that the mass of the fluid is conserved. Defining a reduced pressure according to p = P − ˆ ωρr / and then substituting (3) into (4) generates the governing equation for pressure, ∇ · (cid:18) b µ ∇ p (cid:19) = ∂b∂t , x ∈ R \ Ω( t ) . (5)When the gap between the plates is both spatially and temporally uniform, (5) reduces to Laplace’s equation ∇ p = 0 . We have two boundary conditions on the fluid-fluid interface given by p = − γ (cid:18) κ + 2 b (cid:19) − ρ ˆ ω r , x ∈ ∂ Ω( t ) , (6) v n = − b µ ∂p∂n , x ∈ ∂ Ω( t ) , (7)where γ is the surface tension, κ is the signed curvature of ∂ Ω , ρ is the density of the viscous fluid, and v n isthe normal speed of the interface. The dynamic boundary condition (6) incorporates both the effects of surfacetension as well as the rotation of the Hele–Shaw plates. The kinematic boundary condition (7) relates the velocityof the viscous fluid with the normal velocity of the interface. We also have the far-field boundary condition b µ ∂p∂r ∼ − Q πr + 12 r ∂b∂t , r → ∞ , (8)which acts as a source/sink term at infinity. For Q > and Q < , this corresponds to the bubble area expandingand contracting at rate Q , respectively.To nondimensionalise (5)-(8), we introduce r as the average initial radius of the bubble and Q as the averageinjection rate over the duration of a simulation, and b = b (0 , , then space, time, pressure, and velocity arescaled according to x = r ˆ x , t = r b Q ˆ t, b = b ˆ b, p = 12 µQ b ˆ p, v = Q r b ˆ x . Retaining our original variable names, (5)-(8) become ∇ · (cid:0) b ∇ p (cid:1) = ∂b∂t , x ∈ R \ Ω( t ) , (9) p = − σ (cid:18) κ + 2 R b (cid:19) − ωr , x ∈ ∂ Ω( t ) , (10) v n = − b ∂p∂n , x ∈ ∂ Ω( t ) , (11) b ∂p∂r ∼ − Q πr + 12 r ∂b∂t r → ∞ , (12)where σ = b γ/ µQ r , R = r /b , and ω = ρ ˆ ω r b / µQ .In addition to this expanding bubble problem, we shall also be concerned with two other scenarios, namely theblob geometry, where viscous fluid occupies Ω( t ) and is withdrawn from a point or the cell rotates around aperpendicular axis, and the channel geometry, where viscous fluid occupies a long rectangular channel and isdisplaced by the inviscid fluid which is injected at one end. For these two scenarios, modifications to (9)-(12)will be made as appropriate. 5 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
In the standard case in which b = 1 (parallel, stationary plates) and ω = 0 (no rotation), equations (9)-(12)reduce to ∇ p = 0 , x ∈ R \ Ω( t ) , (13) p = − σκ, x ∈ ∂ Ω( t ) , (14) v n = − ∂p∂n , x ∈ ∂ Ω( t ) , (15) ∂p∂r ∼ − Q πr , r → ∞ . (16)It is instructive to reformulate this problem using complex variable methods as follows. Given the fluid pressure p satisfies Laplace’s equation (13), it can be interpreted as the negative real part of an analytic function W ( z, t ) = − p ( x, y, t ) + i ψ ( x, y, t ) of the complex variable z = x + i y . Here, W is acting as a complex potential, while ψ is a streamfunction.Further, there exists a conformal map z = f ( ζ, t ) from the unit disc of the ζ -plane to the fluid region in the z -plane (i.e., R \ Ω( t ) ). The map will be univalent (that is, one-to-one) and analytic in the unit disc except for ata single point, which we choose to be ζ = 0 , that represents z = ∞ . In the limit ζ → , we have f ∼ a ( t ) /ζ .By fixing a rotational degree of freedom we force a ( t ) to be real. Now the complex potential W ( z, t ) is alsoan analytic function of ζ and so we write w ( ζ, t ) = W ( f ( ζ, t ) , t ) . Given the far-field condition (16), whichimplies W ∼ ( Q/ π ) log z as | z |→ ∞ , we now have the local behaviour w ∼ − ( Q/ π ) log ζ as ζ → .To rewrite the kinematic condition (15), we note that a unit normal vector can be represented in complex numbersby ζf ζ / | ζf ζ | . Therefore the normal velocity v n on the unit circle given by (cid:60){ f t ζf ζ } / | ζf ζ | . This calculationallows the kinematic condition to be represented as (cid:60){ f t ζf ζ } = (cid:60){ ζw ζ } , | ζ | = 1 . Now to reformulate the dynamic condition (14) we note the curvature on the fluid boundary can be written as (cid:60){ ζ ( ζf ζ ) ζ ζf ζ } / | ζf ζ | on the unit circle. Given (14) and the logarithmic behaviour of w as ζ → , we canwrite w = − ( Q/ π ) log ζ − σ K ( ζ, t ) , where K ( ζ, t ) is an analytic function of ζ in the unit disc whose real parton the unit circle | ζ | = 1 is given by the curvature κ . Combining these ideas, we arrive at the single governingequation (cid:60){ f t ζf ζ } = − Q π − σ (cid:60){ ζ K ζ } | ζ | = 1 . (17)This equation is often referred to as the Polubarinova–Galin equation, especially when surface tension is ignored[55].We shall briefly outline five illustrative examples, chosen to demonstrate the key features for the special case inwhich surface tension is ignored. We shall later refer back to these examples when we implement our numericalscheme with surface tension included. First, we shall consider the mapping f = a ( t ) /ζ + b ( t ) ζ N , N ≥ [64]. By substituting into (17) with σ = 0 , we find that a and b must satisfy the coupled system of ODEs a ˙ a − N b ˙ b = Q/ π , N ˙ ab − a ˙ b = 0 . Say, for definiteness, a (0) = 1 and b (0) = (cid:15) , then the second of theseequations gives b = (cid:15)a N , while the first equation integrates to give the time-dependence t = πQ (cid:0) a − − (cid:15) N ( a N − (cid:1) . Thus we have an exact solution, as shown by the red dashed curves in the top panel of Figure 1(a). The innermostcurve is the initial bubble boundary. Here we have chosen N = 5 , Q = 2 π and (cid:15) = 0 . , so this inner curve is acircle with a six-fold perturbation. As time increases, the bubble expands and starts to develop six small fingers.6 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 1:
Solutions with and without surface tension. The red dashed curves in (a)-(e) are zero-surface-tensionsolutions described by the five examples in Section 2.2. The solid blue curves are numerical solutions, includingsurface tension, for the same initial conditions. For the examples in (a), (c) and (d), the zero-surface-tensionsolutions involve a form of finite-time blow-up characterised by cusps forming on the interface; the inclusion ofsurface tension regularises these singularities, allowing the full solution to continue past these blow-up times.For the examples in (b) and (e), blow-up in the zero-surface-tension solution is prevented by the presence ofa logarithmic singularity; here the numerical solution with small surface tension remains close to the zero-surface-tension solution for small time and then deviates away so that the long-time behaviour is different.Each case includes a sketch of the ζ -plane with the unit circle and critical points and logarithmic singularitiesindicated by solid red dots and black diamonds, respectively. Note in (b) we do not plot the critical points at t = 0 as they are outside the field of view here. A PREPRINT
It is of interest to track the critical points, ζ = ζ ∗ , which are the points at which f ζ = 0 . For this example, ζ ∗ = (1 /(cid:15)N a N − ) / ( N +1) . Clearly there are N + 1 critical points that are equally spaced along a circle whichis outside the unit circle in the ζ -plane. As time evolves, each of these critical points moves in a straight linetowards the origin and intersects the unit circle when | ζ ∗ | = 1 , ie a = 1 / ( (cid:15)N ) / ( N − . We can compute theexact time at which this occurs, namely t ∗ = πQ (cid:18) N − N ( (cid:15)N ) / ( N − − (cid:15) N (cid:19) . For the case in Figure 1(a), for which N = 5 , we can see the six critical points (red dots) in the bottom panelmoving towards the unit circle. At t = t ∗ , there is finite-time blow-up, characterised by six cusps of order / along the bubble boundary, which we can see in the top panel of Figure 1(a) (a cusp of order / is characterisedby a curvature singularity that appears locally like two branches of y = x meeting at a cusp, suitably scaledand rotated). The solution cannot be continued past t = t ∗ as the conformal mapping ceases to be univalent.The second zero-surface-tension example we consider is for f = a ( t ) /ζ + log( ζ − d ( t )) − log( ζ + d ( t )) [64]. Again, the geometry is an expanding bubble, however this time the behaviour is qualitatively different. Bysubstituting this map for f into (17) with σ = 0 , we find that a and d satisfy the coupled system: ˙ a = Q π ad − a + 4 da d − (2 d − a ) , ˙ d = − Q π d ( d − a d − (2 d − a ) . For initial conditions, we choose both a (0) > d (0) > so the denominators are initially positive, whichmeans that ˙ a > and ˙ d < for small time. To determine the precise time-dependent behaviour, one optionis to integrate this system numerically, which demonstrates that a ( t ) continues to increase while d ( t ) decreasestowards d = 1 as t increases. To make progress analytically, by dividing one equation by the other, we can alsoderive a first-order ode with exact solution a = 1 d (cid:20) ln (cid:18) ( d (0) − d + 1)( d (0) + 1)( d − (cid:19) + a (0) d (0) (cid:21) , which shows that a ∼ − ln( d − as d → + . Clearly the map has logarithmic singularities at ζ s = ± d andcritical points where ζ ∗ = ± d √ a/ √ a − d . For sufficiently large time, we have a − d > and so both ζ s and ζ ∗ lie on the real ζ -axis and move towards the unit circle as t → ∞ . Since | ζ ∗ |∼ d + d /a as a → ∞ , wesee the critical points ζ ∗ are further away from the unit circle than the logarithmic singularities ζ s , and thereforefinite-time blow-up is avoided.This second example is illustrated in Figure 1(b) using a (0) = 4 , d (0) = 10 / , Q = 2 π . In the top panel, theexact solution is denoted by the dashed red curves. Initially the bubble boundary looks oval in shape on this scale.As time increases the interface expands, leaving two fjords behind, centred both the positive and negative x -axes.In the bottom panel, both the critical points (red dots) and logarithmic singularities (black dots) are indicated inthe ζ -plane. As just described, even though the critical points move towards the unit circle, they are further awaythan the logarithmic singularities; therefore, the logarithmic singularities have the effect of shielding the criticalpoints and preventing blow-up. An interesting observation is that each of the two fjords appears to take the shapeof a classical Saffman-Taylor finger. To see this, note that on the unit circle near ζ = 1 we can derive a localanalysis by setting ζ = 1 + i η , so that f ∼ a + log(1 − d + i η ) . Taking real and imaginary parts and theneliminating η gives x = constant − log(cos y ) , which is the famous finger shape with width π .The third example is probably the most well-known example of an exact solution in Hele–Shaw flow. Heresuppose the geometry is such that there is a blob of fluid in the Hele–Shaw cell, occupying a region Ω( t ) ,surrounded by inviscid fluid. If the fluid is withdrawn from a point in space (the origin, say), then the blob8 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT boundary contracts and we have the less viscous fluid displacing the more viscous fluid. The governing equations(13)-(15) apply in Ω( t ) , while (13) is replaced by ∂p/∂r ∼ Q/ πr as r → . The map we have in mind herefor this example is the quadratic map f = a ( t ) ζ + b ( t ) ζ , where the initial conditions a (0) = 1 , b (0) = (cid:15) (cid:28) correspond to an initial blob boundary that is a perturbed circle [107]. Substituting this map into thePolubarinova–Galin equation with σ = 0 leads to a coupled system of integrable odes with the exact solution a b = (cid:15), t = πQ (cid:0) − a + 2( (cid:15) − b ) (cid:1) . There is one critical point ζ ∗ = − a / (cid:15) , which is initially located in the ζ -plane at − / (cid:15) and moves towardsthe origin as time evolves and a decreases. At t = t ∗ this critical point hits the unit circle, where t ∗ = πQ (cid:18) (cid:15) −
98 (2 (cid:15) ) / (cid:19) , causing a cusp of order / to form on the blob boundary. This behaviour is illustrated in Figure 1(c) for (cid:15) = 0 . and Q = 2 π . On the left panel, the blob boundary is represented by the dashed red curves. Initially, this curveappears circular, as the perturbation is very small. For intermediate times, the left portion of the boundary beginscontracting faster than the remainder of the boundary until the cusp forms, corresponding to finite-time blow-up.On the right panel of Figure 1(c) the location of the fixed point in the ζ -plane is represented (red dots) for thesame three times that the boundary is drawn for in the left panel. Here we see that the critical point touches theunit circle at the precise time that finite-time blow-up occurs. Note that for polynomial maps like the one in thisexample, it has been proven that a cusp will always form before the interface reaches the sink [56]; other explicitsolutions exist whose boundary evolves to the location of the sink before or at the same time as cusp formation(see [109], for example).For completeness we include two more examples, providing only the key details. These examples are for thegeometry of flow in a Hele–Shaw channel, which we fix to be π units wide. For the fourth example, the maptakes for the f = − log ζ + a ( t ) + b ( t ) ζ , with initial conditions a (0) = 0 , b (0) = (cid:15) (cid:28) , corresponding to aslightly perturbed flat interface [67]. This case is analogous to the first and third examples above. The functions a and b satisfy the coupled system of odes with an exact solution a − ln b = − ln (cid:15), t = πQ (cid:0) a − b + (cid:15) (cid:1) . There are critical points at ζ ∗ = 1 /b which move towards the origin as b increases and ultimately intersect theunit circle at t ∗ = π (2 ln(1 /(cid:15) ) − (cid:15) ) /Q , at which time a / cusp forms on the interface. This exampleis illustrated in Figure 1(d), where the interface profiles are shown in the top panel as red dashed curves. In thebottom panel the critical points are indicated (red dots). Here (cid:15) = 0 . and so the critical point is initially at ζ ∗ = 10 and ultimately hits the unit circle at t = t ∗ .Finally, the fifth example is for f = − log ζ + a ( t ) + α log( ζ + d ( t )) with initial conditions a (0) = 0 , d (0) =1 /(cid:15) (cid:29) [67], which is analogous to the second example above. There is a critical point at ζ ∗ = d/ ( α − anda logarithmic singularity at ζ s = − d . We shall not include the details here, but it is possible to derive a coupledsystem of odes for a ( t ) and d ( t ) which can be solved numerically or reduced further analytically by diving oneby the other. For < α < the singularity ζ s is always smaller in magnitude, and therefore closer to the unitcircle, than ζ ∗ . As a result, it turns out that d is a decreasing function and that ˙ d → + as d → + . Therefore,neither ζ ∗ nor ζ s intersect the unit circle but in fact approach it asymptotically as t → ∞ . As such, there is nocusp formation, and instead the proximity of ζ s to the unit circle results in the interface forming a long finger,whose width is (2 − α ) π . In Figure 1(e) we present an example with (cid:15) = 0 . , α = 1 . . The interface in the top9 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT panel, given by the dashed red curves, clearly approaches a finger in shape. The bottom panel shows the criticalpoint (red dots) moving towards the unit circle, but always further away than the logarithmic singularity (blackdots).
To numerically solve (9)-(12), we construct a level set function φ such that the fluid-fluid interface ∂ Ω is thezero level set of φ or ∂ Ω( t ) = { x | φ ( x , t ) = 0 } . If the interface has the normal speed V n , then we wish to construct a speed function, F , such that V n = F on x ∈ ∂ Ω( t ) , and is continuous over the entire computational domain. Thus φ satisfies the level set equation ∂φ∂t + F | ∇ φ | = 0 . (18)We discuss how F is computed in Section 3.2. To solve (18), we approximate the spatial derivatives usinga second order essentially non-oscillatory scheme (see Osher and Fedkiw [96, chapter 3] and Sethian [114,chapter 6] for details), and integrate in time using second order total variation diminishing (TVD) Runge-Kutta(RK), which is performed by taking two forward Euler steps ˜ φ ( n +1) = φ ( n ) − ∆ tF ( n ) | ∇ φ ( n ) | ,φ ( n +2) = ˜ φ ( n +1) − ∆ tF ( n +1) | ∇ ˜ φ ( n +1) | , and then take an averaging step φ ( n +1) = (cid:0) φ ( n ) + φ ( n +2) (cid:1) / . To maintain accuracy and numerical stability, wechoose ∆ t = ∆ x/ (4 max | F | ) .The level set function φ is initialised as a signed distance function satisfying φ = d if x ∈ R \ Ω( t ) − d if x ∈ Ω( t ) , (19)where d is the minimum distance between x and ∂ Ω , via the method of crossing times [96, chapter 7]. Thatis, we advect φ in the normal direction to the interface by solving (18) with F = 1 and determine the point intime where each value of φ crosses from positive to negative. This process is repeated for F = − . To reducenumerical error which can result when the gradient of φ becomes excessively small or large, we maintain φ as a signed distance function, which satisfies | ∇ φ | = 1 , over the duration of a simulation. Re-initialisation isperformed every five time steps by solving ∂φ∂τ + S ( φ )( | ∇ φ |−
1) = 0 , (20)where S ( φ ) = φ √ φ + ∆ x , to steady state. Here τ is a pseudo time variable where ∆ τ = ∆ x/ .While the level set method has successfully been used as a framework for studying a variety of moving bound-ary problems, a limitation of the method is that it can suffer from volume loss or gain in regions where the10 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT mesh is underresolved. In an effort to alleviate this problem, Enright et al. [41] proposed the particle level setmethod, which combines the Eulerian based level set method with a marker particle based Lagrangian approach.We briefly describe the algorithm here, and refer the reader to Enright et al. [41] for a more comprehensivedescription, as well as examples illustrating the effectiveness of the particle level set method.The method works by placing massless particles on both sides of the interface, where positive and negativeparticles are assigned in the regions φ > and φ < , respectively. We denote r p as the minimum distancebetween the interface and particle’s location The marker particles are advected according to ∂ x p ∂t = F n , where x p is the location of the particle and n = ∇ φ/ | ∇ φ | is a unit vector that reduces to the outward facingnormal on the interface. If a particle crosses the interface, this indicates that mass has been lost (or gained).We mitigate this error by locally rebuilding the interface by constructing a local level set function from the fouradjacent nodes to the particle defined as φ p ( x ) = s p ( r p − | x − x p | ) , where s p = 1 and − for positive andnegative particles, respectively. Using φ p , φ is corrected using φ = φ + if | φ + |≤ | φ − | φ − if | φ + | > | φ − | , where φ + = max( φ p , φ + ) and φ − = min( φ p , φ − ) . This procedure is performed both when the level setequation (18) is solved as well as when re-initialisation is performed. F From the kinematic boundary condition (11), we have the expression for the speed function F = − b ∇ p · n x ∈ R \ Ω( t ) , (21)recalling F is required to solve the level set equation (18) and n = ∇ φ/ | ∇ φ | reduces to the outward facingnormal on the interface. Equation (21) satisfies (11) on the interface, or V n = F on x ∈ ∂ Ω( t ) where V n isthe normal speed of the interface, and provides a continuous expression for F in the viscous fluid region. Thederivatives in (21) are evaluated using central differencing. However, to solve (18) we require an expression for F over the entire computational domain. It was proposed by Moroney et al. [90] that the speed function can beextended into the inviscid fluid region by solving the biharmonic equation ∇ F = 0 x ∈ Ω( t ) . (22)This ensures that F = V n is satisfied on the interface and gives a continuous expression for F away from theinterface. For the purpose of discretisation, the sign of φ is used to determine nodes inside the interface that needto be included in the biharmonic stencil. This discretisation results in a symmetric system of linear equationswhich is solved using LU decomposition. As such, the location of the interface does not need to be knownexplicitly, similar to the level set method itself. This velocity extension process is a variant of a thin plate splinein two dimensions [14]. To illustrate the velocity extension process, we consider an example F defined in theregion Ω( t ) , shown in Figure 2 ( a ) . The red line represented the fluid-fluid interface ∂ Ω( t ) . Figure 2 ( b ) shows F after the biharmonic equation (22) is solved. We see that this velocity extension process gives a differentiableexpression for F over the entire computational domain, which can be used to solve (18).11 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 2:
An illustration of the velocity extension process used to extend F to be defined over the entirecomputational domain. ( a ) shows an example function that is undefined where x ∈ Ω . ( b ) shows this regionsbeing ‘filled-in’ by solving the biharmonic equation (22). This gives us a differentiable extension of F over theentire domain. To evaluate the speed function F , we must first compute the pressure field. We consider (9)-(12) in polarcoordinates with p = p ( r, θ, t ) and the location of the interface is given by r = s ( θ, t ) . Thus (9) becomes r ∂∂r (cid:18) rb ∂p∂r (cid:19) + 1 r ∂∂θ (cid:18) b ∂p∂θ (cid:19) = ∂b∂t r > s ( θ, t ) . (23)For nodes away from the interface, the derivatives in (23) are discretised using a standard 5-point stencil, illus-trated in Figure 3 ( a ) . Denoting β = rb , the r -derivatives in (23) are approximated via r ∂∂r (cid:18) β ∂p∂r (cid:19) → r i,j ∆ r (cid:18) β i +1 / ,j p i +1 ,j − p i,j ∆ r − β i − / ,j p i,j − p i − ,j ∆ r (cid:19) , (24)where β i +1 / ,j = ( β i +1 ,j + β i,j ) / and β i − / ,j = ( β i − ,j + β i,j ) / . The derivatives in the θ -direction arediscretised in a similar fashion.As illustrated in Figure 3 ( b ) , 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 theinviscid and viscous fluid regions respectively, illustrated in Figure 3. When discretising (23), we can no longerincorporate p i − ,j into our finite difference stencil as it is not in the domain x ∈ R \ Ω( t ) . Instead, we define aghost node at r I (denoted by the red dots in Figure 3 ( b ) ) whose value is p I . By noting that φ is a signed distancefunction, the distance between r i,j and r I is computed via h = ∆ r (cid:12)(cid:12)(cid:12)(cid:12) φ i,j φ i − ,j − φ i,j (cid:12)(cid:12)(cid:12)(cid:12) . As per Chen et al. [24], our finite difference stencil becomes r ∂∂r (cid:18) β ∂p∂r (cid:19) ≈ r i,j (∆ r + h ) (cid:18) β i +1 / ,j p i +1 ,j − p i,j ∆ r − ˆ β i − / ,j p i,j h (cid:19) + 2 r i,j h (∆ r + h ) ˆ β i − / ,j p I . (25)12 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 3:
An illustration of how the pressure of the viscous fluid is solved for using the finite difference method,represented in Cartesian coordinates for visualisation purposes. ( a ) For nodes away from the interface, we solvefor pressure (23) using a standard 5-point finite difference stencil (24). ( b ) However, this stencil cannot be usedfor nodes adjacent to the interface, as in this case the southern node is not in the domain x ∈ R \ Ω( t ) , and thuscannot be used in our stencil. Instead we impose a ghost node, denoted by the red dots, on the interface whoselocation corresponds to the point where φ = 0 . The value at this ghost node is determined from the dynamicboundary condition (10). This leads to the non-standard finite difference stencil (25). Here ˆ β i − / ,j = ( β i,j + β I ) / where β I is the value of β on the interface. When the node and interface aresufficiently close together such that h < ∆ r , we set p i,j = p I . A similar procedure is applied if the interfacelies 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 (10). To determine the value of p I , we firstcompute the curvature of φ via κ = ∇ · n over the entire computational domain, recalling n = ∇ φ/ | ∇ φ | . Thevalue of κ on the interface is computed via linear interpolation using κ i,j and κ i − ,j leading to p I = − σ (cid:18) κ i,j − h ( κ i,j − κ i − ,j )∆ r + 2 b ( r I , θ j ) (cid:19) − ωr I . For more information about solving both elliptic and parabolic problems in irregular domains using the finitedifference method in conjunction with level set functions, we refer the reader to Coco and Russo [25], Gibouet al. [50] (and references therein). Once the finite difference stencil has been formed, the resulting system oflinear equations is solved using LU decomposition.
To incorporate the far-field boundary condition (12) into our finite difference stencil, we utilise a Dirichlet-to-Neumann map. This map is implemented by imposing an artificial circular boundary at r = R such that R > s ( θ, t ) . By only considering the region s ( θ, t ) ≤ r ≤ R , we seek a solution to (9) of the form ˆ p ( r, θ, t ) = A − Q π log r + r ∂b∂t + ∞ (cid:88) n =1 r − n ( A n cos nθ + B n sin nθ ) , (26)where A , A n , and B n are unknown, and ˆ p = b p . The expansion (26) assumes that b is spatially uniform in r ≥ R so that ˆ p satisfies the appropriate Poisson equation. Considering the value of pressure on the artificial13 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT boundary, suppose that ˆ 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θ, (27)where a = 12 π (cid:90) π ˆ p ( R, θ, t ) d θ, a n = 1 π (cid:90) π ˆ p ( R, θ, t ) cos nθ d θ, b n = 1 π (cid:90) π ˆ p ( R, θ, t ) sin nθ d θ. By equating (27) with (26) 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 .To incorporate our expression for ˆ p into our finite difference stencil, we differentiate (27) with respect to r andevaluate it at r = R to give ∂∂r ˆ p ( R, θ j ) = − Q πR + R ∂b∂t − ∞ (cid:88) n =1 nR ( a n cos nθ j + b n sin nθ j ) . By approximating the integrals in our expressions for a n and b n as a n ≈ ∆ θπ m (cid:88) k =1 ˆ p ( R, θ k ) cos nθ k and b n ≈ ∆ θπ m (cid:88) k =1 ˆ p ( R, θ k ) sin nθ k , then ∂∂r ˆ p ( R, θ j ) ≈ − Q πR + R ∂b∂t − ∆ θRπ m (cid:88) k =1 v jk ˆ p ( R, θ k ) , where v jk = ∞ (cid:88) n =1 n cos( n ( θ k − θ j )) . Defining I as the outermost index at which r = R , then our expression for ∂p/∂r is incorporated into our finitedifference stencil, (28) 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 .As an aside, we note this procedure can be adapted to model a Dirichlet boundary condition of the form p ∼ p ∞ as r → ∞ where p ∞ is prescribed. This type of boundary condition would be more appropriate for the model ofa Stefan problem [92], where now p would represent temperature which is prescribed in the far-field. Assumingthat b is constant, then (26) becomes p = p ∞ + A ln r + ∞ (cid:88) n =1 r − n ( A n cos nθ + B n sin nθ ) , By following the same steps outlined above, we have ∂∂r p ( R, θ j ) = a − p ∞ R ln R − ∞ (cid:88) n =1 nR ( a n cos nθ j + b n sin nθ j ) . This expression is then incorporated into our finite difference stencil as per usual.14 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
We summarise our numerical algorithm as follows:
Step 1
For a given initial interface s ( θ, , choose φ such that it satisfies (19), and then initialise φ as a signeddistance function using the method of crossing times. Step 2
Place marker particles around the interface, noting which side of the interface they are on, as well astheir minimum distance from the interface.
Step 3
Solve for pressure in the domain x ∈ R \ Ω using the procedure described in Section 3.3. Step 4
Compute F according to (21), and then extend F into the region x ∈ Ω by solving the biharmonicequation (22). Step 5
Using F , update both φ by solving (18) and the marker particles. Step 6
Correct φ (if necessary) using the marker particles. Step 7
Re-initialise the level set function by solving (20), and then correct φ using the marker particles. Step 8
Repeat steps 2-7 until t = t f . In this section, we present a selection of results demonstrating the capabilities of the numerical scheme presentedin Section 3. We show how our framework can be modified to solve for a wide range of different configurationsof the Hele–Shaw cell, and provide examples illustrating that simulations are capable of producing solutionsconsistent with previous experimental and numerical results.
We first consider the standard Hele–Shaw problem in which the inviscid bubble is injected into the viscous fluidwhile the plates are parallel and stationary such that b = 1 and ω = 0 . As a preliminary test for our scheme, we demonstrate that numerical solutions converge for a sufficiently refinedgrid. To do so, we perform simulations with the initial condition s ( θ,
0) = 1 + ε cos mθ, (29)where ε = 0 . , m = 6 , on the computational domain ≤ r ≤ . and ≤ θ < π , and are performed usingan increasingly refined mesh. These simulations, shown in Figure 4, indicate that the interfacial profiles areconverging as the mesh is refined, and that grid independence is achieved, at this scale, using × equallyspaced nodes. Further, our solutions appear to maintain six fold symmetry over the duration of the simulation.We also demonstrate that our use of the Dirichlet-to-Neumann map, described in Section 3.3.1, results in thebubble’s volume changing at rate Q . To do so, we consider three different injection rates. The first is theconstant injection rate Q = 1 , the second is the sinusoidal injection rate, Q = 1 + 0 . πt/t f ) ,
15 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 4:
Convergence test of numerical scheme for the evolution of a bubble with initial condition (29), wheresolutions are computed using ( a ) 350 × , ( b ) 550 × , ( c ) 750 × , and ( d ) 850 × equally spacednodes. Additionally, Q = 1 , σ = 5 × − , ω = 0 , and t f = 100 . Simulations are performed on the domain ≤ r ≤ . and ≤ θ < π . Solutions are plotted in time intervals of t = 10 . and the third is the piecewise rate Q = (cid:40) . t ≤ t f / . t > t f / . We find that the rate of change of volume computed from the numerical simulations compares well with thecorresponding exact rate of expansion, shown in Figure 5. This suggests that the Dirichlet-to-Neumann mapcorrectly ensures that the bubble expands at the correct rate for both constant and time-dependent Q . Surface tension, modelled via the inclusion of the signed curvature term in the dynamic boundary condition(10), acts to regularise the Hele–Shaw model. As shown in Section 2.2, in the absence of surface tension,exact solutions exist that exhibit very different behaviour (finite-time cusp formation or finger formation) despiteinitial conditions that are arbitrarily close; hence the zero-surface-tension-problem is ill-posed. Adding surfacetension removes the possibility of cusp formation; in Figure 1, we have included the solutions for small butnonzero surface tension, computed using the method described in Section 3, for each of the cases described in16 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 5:
The rate of change of volume, ˙ V , of numerical solution with constant (red), periodic (blue) andpiecewise (yellow) injection rates. Dotted black lines denote the corresponding exact rate of change of volume.Simulations are performed on the domain ≤ r ≤ . and ≤ θ < π using × equally spaces nodeswith the initial condition (29). The surface tension parameter is σ = 5 × − and final time of simulations is t = 100 . Section 2.2. With nonzero surface tension, the interface remains smooth and solutions exist for all time, or, inthe case of the finite blob (Figure 1(c)), until the interface intersects the point at which fluid is being withdrawn.Linear stability analysis [101] indicates that increasing the surface tension parameter σ acts to make the interfaceless unstable, and nonlinear numerical simulations and experiments [23, 32, 62] show that increasing Q , whichis mathematically equivalent to decreasing σ , results in an increase in the number of fingers that develop. Weperform numerical simulations for a fixed injection rate with values of surface tension varying over several ordersof magnitude with the initial condition s ( θ,
0) = 1 + ε (cos mθ + sin nθ ) , (30)where ε = 0 . , m = 3 , and n = 2 . Our numerical simulations, shown in Figure 6, are able to reproducethe key morphological features of the Saffman–Taylor instability. For each value of σ considered, the interfaceis unstable, and the sinusoidal perturbations in the initial condition (30) grow into fingers. As σ is decreased,the interface becomes more unstable and the number of fingers that develop over the duration of a simulationincreases due to tip-splitting occurring (see (1) in Figure 6 ( b ) for example). Additionally, our simulations areable to produce so called ‘shielding’ behaviour, where neighbouring fingers can block one another off, which inturn results in fingers retracting (denoted by (2) in Figure 6 ( c ) ). This behaviour is known to occur experimentally(see Figure 1 in [101] for example). Finally, when surface tension is sufficiently small, ‘feathering’ can occur,when a finger does not strictly tip-split, but instead develops ripples along one of its sides as it expands (denotedby (3) in Figure 6 ( d ) ). Again this behaviour has been observed in experiments (see Figure 3 in [23]). One of the more popular modifications to the Hele–Shaw cell (particularly in recent years) is to consider that theplates of the Hele–Shaw cell are no longer parallel but instead linearly tapered (either converging or diverging)17 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 6:
The numerical solution to (9)-(12) with Q = 1 for values of surface tension parameter σ ( a )1 . × − , ( b ) 5 × − , ( c ) 1 . × − , and ( d ) 6 . × − with initial condition (30). Our numericalscheme captures different morphological behaviour including (1) tip-splitting, (2) shielding, and (3) feathering,all of which have been observed experimentally. Simulations are performed on the domain ≤ r ≤ and ≤ θ < π using × equally spaced nodes. Solutions are plotted in time intervals of t = 5 up to t f = 55 .
18 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 7:
Numerical simulations where the gap between the plates linearly tapered according to (31) with Q = 1 , b ∞ = 1 / , R = 8 / , r = 7 , and α = 2 / , with surface tension σ ( a ) ( b )
1. Simulations areperformed on the domain ≤ r ≤ . and ≤ θ < π using × equally spaced nodes. Solutions areplotted in time intervals of 5.6. up to t f = 44 . . in the direction of the flow. The concept was first introduced in Zhao et al. [133], who considered tapered platesin channel geometry, discussed further in Section 4.4. Imposing tapered plates in radial geometry has beenstudied using linear stability analysis [2], weakly-nonlinear stability analysis [8], experimentally [2, 13], andby numerical simulations [69]. Results indicate that tapering the gap between the plates either such that theyconverge or diverge can have a stabilising or de-stabilising effect on the interface depending on the injection rate.We use the numerical scheme presented in Section 3 to study how tapering the plates while injecting at either aconstant or time dependent injection rate can be used to reduce the development of viscous fingers in Morrowet al. [93].We perform numerical simulations where the gap between the plates is linearly tapered according to b ( r ) = b ∞ + α ( r − r B ) if r ≤ r B ,b ∞ if r > r B , (31)where α is the gradient of the taper ( α = 0 being the unmodified Hele–Shaw cell). Experiments were recentlyperformed by Bongrand and Tsai [13], who considered a Hele–Shaw cell where the plate gap is (31) where α > for different injection rates. Bongrand and Tsai showed that if the injection rate is sufficiently small,the interface is completely stabilised over the duration of the experiment. To confirm that our numerical schemeproduces solutions consistent with these experiments, we perform simulations with several different injectionrates, shown in Figure 7. The initial condition of these simulations is (30) where ≤ θ < π , ε = 5 × − , m = 5 , and n = 4 . For a low injection rate (Figure 7 ( a ) ), we indeed find that the interface is stabilised, andremains circular over the duration of the simulation. For a faster injection rate (Figure 7 ( b ) ), we find that theinterface is unstable, and in particular, morphologically speaking these fingers appear distinct from traditionalSaffman–Taylor fingers (see Figure 6 for example). This behaviour is consistent with the results of Bongrandand Tsai [13], who described the interface as becoming ‘wavy’ as it expanded. This suggests our numericalsimulations are producing the correct behaviour when the plates are of the form of (31).19 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
In this subsection, we now consider the complementary problem to (9)-(12), for which a blob of viscous fluid nowoccupies the region Ω( t ) and the inviscid fluid is in R \ Ω( t ) . This problem is typically studied by consideringthe withdrawal of the viscous fluid, which in turn causes fingers to develop inward (as in Figure 1(c), for exam-ple). However, popular modifications including considering the gap between the plates to be time-dependent orrotating the entire Hele–Shaw cell have also received interest. In this subsection, we explain how the numericalscheme presented in Section 3 can be modified to solve for these variations, and show that our simulations areconsistent with previous experimental and numerical results.For the case where the inviscid bubble is injected into the viscous fluid, the velocity of the fluid is driven by asink term in the far-field given by (12). For the blob problem, the withdrawal of the viscous fluid is incorporatedinto the model via a sink at the origin. To include this sink into our numerical model, we follow Hou et al. [62]and introduce a smoothed Dirac delta function in (9) such that ∇ · (cid:0) b ∇ p (cid:1) = ∂b∂t + S, (32)where S = Qb ¯ r (cid:18) πr ¯ r (cid:19) if r ≤ ¯ r , if r > r . (33)As shown by Hou et al. [62], this choice of source term ensures the rate of change of volume of the viscous fluidis Q . For all results in this subsection, we use ¯ r = 0 . . We note that it is straightforward to extend (33) toinclude multiple sink/source points.When solving for the governing equation for pressure (9) for the case where in inviscid fluid is injected into aviscous one, we consider the model in polar coordinates as it is simpler to incorporate the far-field boundarycondition (12) via the Dirichlet-to-Neuman map (Section 3.3.1) on a circle. However, when the inviscid regionsurrounds the viscous fluid, we no longer have a far-field boundary condition and, as such, it is more convenientto solve for pressure in Cartesian coordinates, although of course either coordinate system could be used. Thediscretisation of the spatial derivatives (32) is performed as was described in Section 3.3. That is, we use centralfinite differences for nodes away from the interfaces, and incorporating a ghost value of p for nodes adjacent tothe interface. We refer the reader to Chen et al. [24] and Gibou et al. [50, 51] for more details on implementingthis stencil in Cartesian coordinates.Similarly to Section 3.2, once we have computed F via (21), we extend it into the region x ∈ R \ Ω( t ) by solvingthe biharmonic equation (22). When solving (9)-(12), the speed function was known outside the interface andwas extended in. We now have an expression for F inside the interface that needs to be extended outward.This means we require boundary conditions on each of the four computational boundaries. When formingour biharmonic stencil, we apply homogeneous Neumann boundary conditions on each of the boundaries. Weillustrate this process in Figure 8, which shows an example F before (Figure 8 ( a ) ) and after (Figure 8 ( b ) ) thevelocity extension process. Again we see that by solving the biharmonic equation we obtain a smooth expressionfor F over the entire computational domain. We consider the case where the viscous fluid is withdrawn such that as the interface contracts, viscous fingersform inward. As described in Section 2.2, exact solutions in the absence of surface tension ( σ = 0 ) are ill-posed,20 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 8:
An illustration of the velocity extension process where the viscous fluid is surrounded by the inviscidbubble. As discussed in Section 3.2, this is done by solving the biharmonic equation in the region where x ∈ Ω( t ) , where now we apply homogeneous Neumann boundary conditions on the edge of the computationaldomain. and unphysical cusps can form before the interface reaches the sink (the point at which liquid is withdrawn).Numerical simulations performed using the boundary integral method have investigated the regularising effectsof surface tension preventing these cusps from forming [19, 70]. Experimental results [101, 124] show that thefingers that form exhibit morphological features distinct from traditional Saffman–Taylor fingers, in that thesefingers do not appear to undergo tip-splitting but instead appear to be in competition to ‘race’ toward the sink.As well as the comparison with the zero surface tension solution made in Figure 1(c), we perform two furtherdifferent simulations for this configuration. The first is with an initially circular interface centred at (0 , − . shown in Figure 9 ( a ) . This simulation shows that as the interface contracts, it becomes non-convex until asingle finger develops that tends towards the origin. This behaviour compares well with previous numericalsimulations performed using the boundary integral method [70]. For the second simulation, shown in Figure 9 ( b ) ,we consider a perturbed circle centred at the origin of the form s ( θ,
0) = 1 + ε (cos 3 θ + sin 7 θ + cos 15 θ + sin 25 θ ) , (34)where ε = 5 × − . We find that the interface initially develops numerous short fingers. These fingers donot appear to exhibit the same morphological features as the case where the inviscid bubble is injected such astip-splitting and feathering (see Figure 6 for example), but instead the number of fingers remains constant. Dueto the pressure differential between the sink and the boundary of the bubble, the velocity of one of the fingersrapidly increases, and the simulation is stopped when this finger reaches the origin. We note that this behaviourcompares well with experimental results (see Figure 15 in [124] for example), as well as numerical simulationsin [22]. A popular modification to the blob problem is to consider the case where the upper plate is uniformly lifted intime such that b → b ( t ) . The volume of viscous fluid remains constant ( Q = 0 ) so that when the plates areseparated, this results in fingers developing inward. For this configuration, the governing equation for pressure(32) becomes Poisson’s equation. However, pressure can be reduced via P = p − ˙ b | x | / (4 b ) essentially21 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 9:
Numerical simulation where the viscous fluid is withdrawn at a point located at the origin where Q = 1 . For ( a ) , initial condition is a circle of radius unity centred at (0 . , where σ = 8 × − and t f = 1 . . For ( b ) , initial condition is (34) where σ = 1 . × − and t f = 1 . . Black dot denotes regionwhere r ≤ . . Simulations are performed on the domain − . ≤ x ≤ . and − . ≤ y ≤ . using × equally spaced nodes. moving the non-homogeneous term to the dynamic boundary condition (10). This reduces (9) to Laplace’sequation, meaning this configuration can be solved numerically with the boundary integral method [134]. Thelifting plate problem was first considered mathematically by Shelley et al. [116], who demonstrated a relationshipbetween the number of fingers that develop and the surface tension parameter. In particular, it was shown that thenumber of fingers is a monotonically decreasing function of time. As a point of comparison, for the traditionalHele–Shaw configuration, discussed in Section 4.1, the number of fingers typically increases with time due totip-splitting. This behaviour was later shown to be consistent with experimental results [78, 94]. We note that thecase where the inviscid bubble is injected into the viscous fluid while the plates are separated has also receivedattention [93, 127, 135].We perform a simulation using the linearly increasing gap between the plates b = 1 + t for different values of σ . When σ = 10 − (row one of Figure 10), we find that the interface quickly destabilises and approximately 15fingers develop by t = 1 . As time increases further, neighbouring fingers begin to merge with each other and theoverall number of fingers significantly decreases, and we find only around five fingers remain by the conclusionof the simulation. For the lower value of surface tension σ = 5 × − (row two of Figure 10), we find thenumber of fingers that initially develop has increased compared to when σ = 10 − , about 25 at t = 1 . Howeveras the interface contracts, again fingers begin to merge and only 8 fingers remain at t = 4 . We perform thesesimulations over a longer time period (not shown) which reveals that the interface will become circular when thegap between the plates is sufficiently large. This behaviour is consistent both with previous experimental andnumerical results [78, 94, 116]. While Saffman–Taylor fingers traditionally form due to the injection/withdrawal of one immiscible fluid intoanother, it is known that these fingers can also be triggered by body forces. The two most commonly studied body22 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT σ = 10 − σ = 5 × − Figure 10:
Time evolution of viscous fluid where plates are separated according to b = 1 + t with initialcondition (34). Solutions are plotted at times (left to right) t = 0 , 1, 1.8, 2.8, and 4. Simulations are performedon the domain − . ≤ x ≤ . and − . ≤ y ≤ . using × equally spaced nodes. forces are gravity and centrifugal forces. For the latter, when the entire Hele–Shaw cell is rotated, this resultsin the dense viscous fluid being propelled outward which in turn results in fingers forming [113]. Experimental[16] and numerical [4, 44, 99] studies reveal that the interface patterns are distinct from the traditional Saffman–Taylor instability. That is, fingers appear more ‘stretched-out’ and generally go not undergo tip-splitting. Wenote that our model (9)-(32) ignores the effect of Coriolis forces, however several studies have investigated itseffect on interfacial dynamics [4, 113, 131]. Further, we consider the case where the inviscid bubble is injectedinto the viscous fluid while the plates are rotated in Morrow et al. [93], where we show that the angular velocityhas a stabilising effect on the interface.The incorporation of the centrifugal term is straightforward. While body forces appear in the governing equationfor the velocity of the viscous fluid (3), scaling pressure means that the angular velocity term can be moved to thedynamic boundary condition (10) (this can be done for any conservative body force f satisfying ∇ × f = ).The rotating Hele–Shaw cell has previously been studied numerically using boundary integral [113] and diffusiveinterface [22, 99] techniques. We perform simulations where the amount of viscous fluid is constant ( Q = 0 )and the Hele–Shaw plates are rotated ( ω > ) for different values of σ , shown in Figure 11. For each valueof σ considered, we find that the interface is unstable, and the fingers that develop are distinct from traditionalSaffman–Taylor fingers. In particular, the fingers that develop do not appear to tip-split but instead remainconstant. Additionally, the number of fingers that develop increases as the surface tension parameter is decreasedsuch that when σ = 10 − (Figure 11 ( a ) ), 7 fingers form, while for σ = 10 − (Figure 11 ( d ) ), the number offingers is approximately 21. These results are consistent with experimental results [16] and numerical simulations[4, 99]. 23 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 11:
Numerical simulation of a rotating Hele–Shaw cell with ω = 1 , Q = 0 , and σ ( a ) 10 − , ( b )5 × − , ( c ) 2 . × − , and ( d ) 10 − . Corresponding final time of simulations is t = 0 . , . , . , and . . Initial condition for each simulation is (34). Simulations are performed on the domain − ≤ x ≤ and − ≤ y ≤ using × equally spaced nodes. In subsections 4.1-4.3, we considered Hele–Shaw flow in radial geometry, where the bubble-fluid interface iscompletely immersed by an infinite amount of viscous or inviscid fluid. In this subsection, we focus on anotherwell studied version of the Hele–Shaw cell is in channel (or rectangular) geometry, where the shape of the cell isa narrow rectangle of infinite length and width L . As discussion in subsection 2.2, for the zero-surface-tensioncase, exact solutions are known to exist which may involve a type of blow-up in finite time with a cusp formingon the boundary (as in Figure 1(d)).This channel problem dates back to the work of Saffman and Taylor [112], who showed that when an inviscidbubble is injected into the channel filled with a viscous fluid, typically a single finger develops that propagatesthrough the channel (see the numerical solution in Figure 1(e)). Since the work of Saffman and Taylor, extensiveresearch has been carried out determining how the parameters of the model influence the width (relative to L )24 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT and speed of this finger. In particular, for a fixed injection rate, as the surface tension parameter is increased, itis established that the width and speed of the finger increase and decrease, respectively. We refer the reader toHomsy [58] and Saffman [111] (and references therein) for a comprehensive overview of the problem. In thissubsection (and subsection 2.2), we restrict ourselves to the classic configuration where b is constant, but wenote that similar to the radial problem, our scheme can easily be used to study non-standard cases, such as thosein the papers [2, 46, 125, 133].As with the blob problem discussed in Section 4.3, it is more convenient to consider this problem in Cartesiancoordinates such that p → p ( x, y, t ) and the interface is given by x = f ( y, t ) . Similar to the radial case (see(16), for example), the velocity of the fluid is driven by the sink term in the far-field b ∂p∂x ∼ QL as x → ±∞ . (35)Equation (35) is incorporated into our finite difference stencil using a Dirichlet-to-Neumann map. We do notprovide full details, but note the procedure is similar to that described in Section 3.3.1, where we impose anartificial boundary at x = X , and seek a solution to (9) of the form ˆ p ( x, y, t ) = QL x + ∞ (cid:88) n =0 A n e − λ n y cos λ n x, where λ n = 2 πn/L and A n is to be determined. Additionally, we also impose ∂p/∂y = 0 on y = ± L/ .As as the numerical results in Figure 1, we perform a series of simulations using the same parameters as thoseof DeGregoria and Schwartz [36], who studied this problem using a boundary-integral approach, to demonstrateour solutions are consistent with the expected behaviour. We choose the initial condition f ( y,
0) = ε cos 2 πy, where ε = 0 . , and perform simulations over a range of values of σ , shown in Figure 12. For low σ (row oneof Figure 12), we find that as the bubble expands, a finger grows, which is unstable and split into two. This isconsistent with the results of DeGregoria and Schwartz [36], and this behaviour is also observed experimentallyby Tabeling et al. [118] when the injection rate is sufficiently large. For larger values of σ (rows two to four ofFigure 12), a single stable finger propagates through the channel whose speed and width decreases and increasesas σ increases. Again, this behaviour is consistent with previous experimental and numerical results [36, 118]. In this article, we have reviewed a suite of Hele-Shaw configurations with two immiscible fluids separated by asharp interface. Our focus is on one-phase models, which arise by assuming one fluid is much less viscous thanthe other (indeed we assume the less viscous fluid is inviscid and ignore its contribution). For the standard Hele-Shaw configuration with parallel stationary plates, we have summarised how complex variable and conformalmapping techniques can be applied to the zero-surface-tension model to deduce a variety of exact analyticalresults. The three geometries we have focussed on involve a bubble expanding into a body of viscous fluid, ablob of fluid withdrawn from a point or viscous fluid displaced by an inviscid fluid in a Hele-Shaw channel.Despite the drawbacks of Hele-Shaw models without surface tension in terms of physical applicability, thesecomplex variable approaches are very well studied by applied mathematicians and have motivated numerouspapers on moving boundary problems in general.We have also reviewed a series of alterations to the standard one-phase Hele-Shaw model. For these alterationsapplied in various combinations with the three geometries, we have presented a flexible numerical scheme based25 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT
Figure 12:
Numerical simulations in channel geometry with values of the surface tension parameter (top tobottom) σ = 2 × − , × − , × − , and × − . Initial condition for all simulations is f ( x,
0) =0 .
05 + cos(2 πy ) . Simulations are performed on the domain − . ≤ x ≤ . and − . ≤ y ≤ using × equally spaced nodes. Solutions are plotted in time intervals of t = 0 . up to t f = 3 . on the level set method. We have shown that our scheme is capable of reproducing the complicated interfacialpatterns that form in Hele–Shaw flow while using a uniform computational grid. By making straightforward,appropriate adjustments to the scheme, we have been able to solve for a wide range of configurations. We havepresented a selection of some of the more well-studied configurations, including the expanding bubble problem,linearly tapered plates, the withdrawal of fluid from a viscous blob, time-dependent plate gap, rotating Hele-Shawcell, and flow in a channel geometry. For all of these configurations, we have demonstrated that our simulationscompare well with previous experimental and numerical results.While we have considered a range of different Hele–Shaw configurations in this article, this is by no means an ex-haustive list. Using our numerical scheme, opportunities exist to study configurations which have not previouslybeen considered either experimentally and numerically. For example, while the linearly tapered configuration,discussed in Section 4.2, has received significant attention [1, 8, 13, 69], including our own study in Morrow et al.[93], an open question is to determine the effect of tapering the plates for the corresponding blob problem. Ourscheme could be used to gain insight the effect of the taper angle on viscous finger development when the fluidis withdrawn compared to the parallel plate case discussed in Section 4.3. Further, additional physical effects onthe interface between fluids could be easily included, such as kinetic undercooling [6, 34] or dynamic wettingeffects [15, 100]. Further adjustments could be made to apply the scheme to study controlling instabilities inHele-Shaws with an elastic membrane [102, 103, 105] or with an external electrical field [47, 89] and muchmore. Acknowledgements
SWM acknowledges the support of the Australian Research Council via the Discovery Projects DP140100933.He would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospi-26 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT tality during the programme Complex Analysis: Techniques, Applications and Computations where part of thework on this paper was undertaken. This programme was supported by the EPSRC grant EP/R014604/1. He isgrateful for the generous support of the Simons Foundation who provided further financial support for his visitto the Isaac Newton Institute via a Simons Foundation Fellowship.
References [1] T. T. Al-Housseiny and H. A. Stone. Controlling viscous fingering in tapered Hele-Shaw cells.
Physics ofFluids , :092102, 2013.[2] T. T. Al-Housseiny, P. A. Tsai, and H. A. Stone. Control of interfacial instabilities using flow geometry. Nature Physics , :747, 2012.[3] T. T. Al-Housseiny, I. C. Christov, and H. A. Stone. Two-phase fluid displacement and interfacial insta-bilities under elastic membranes. Physical Review Letters , :034502, 2013.[4] E. Alvarez-Lacalle, H. Gadˆelha, and J. A. Miranda. Coriolis effects on fingering patterns under rotation. Physical Review E , :026305, 2008.[5] P. H. A. Anjos and J. A. Miranda. Radial viscous fingering: Wetting film effects on pattern-formingmechanisms. Physical Review E , :053003, 2013.[6] P. H. A. Anjos, E. O. Dias, and J. A. Miranda. Kinetic undercooling in Hele-Shaw flows. Physical ReviewE , :043019, 2015.[7] P. H. A. Anjos, V. M. M. Alvarez, E. O. Dias, and J. A. Miranda. Rotating Hele-Shaw cell with a time-dependent angular velocity. Physical Review Fluids , :124003, 2017.[8] P. H. A. Anjos, E. O. Dias, and J. A. Miranda. Fingering instability transition in radially tapered Hele-Shaw cells: Insights at the onset of nonlinear effects. Physical Review Fluids , :124004, 2018.[9] G. Aronsson and U Janfalk. On Hele-Shaw flow of power-law fluids. European Journal of AppliedMathematics , :343–366, 1992.[10] R. Arun, S. T. M. Dawson, P. J. Schmid, A. Laskari, and B. J. McKeon. Control of instability by injectionrate oscillations in a radial Hele-Shaw cell. Physical Review Fluids , :123902, 2020.[11] T. H. Beeson-Jones and A. W. Woods. Control of viscous instability by variation of injection rate in a fluidwith time-dependent rheology. Journal of Fluid Mechanics , 829:214–235, 2017.[12] E. Ben-Jacob and P. Garik. The formation of patterns in non-equilibrium growth.
Nature , :523–530,1990.[13] G. Bongrand and P. A. Tsai. Manipulation of viscous fingering in a radially tapered cell geometry. PhysicalReview E , :061101, 2018.[14] F. L. Bookstein. Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Trans-actions on Pattern Analysis and Machine Intelligence , :567–585, 1989.[15] F. P. Bretherton. The motion of long bubbles in tubes. Journal of Fluid Mechanics , :166–188, 1961.[16] L. Carrillo, F. X. Magdaleno, J. Casademunt, and J. Ort´ın. Experiments in a rotating Hele-Shaw cell. Physical Review E , :6260, 1996.[17] J. Casademunt. Viscous fingering as a paradigm of interfacial pattern formation: Recent results and newchallenges. Chaos , :809–824, 2004. 27 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT [18] H. D. Ceniceros and T. Y. Hou. The singular perturbation of surface tension in Hele-Shaw flows.
Journalof Fluid Mechanics , :251–272, 2000.[19] H. D. Ceniceros, T. Y. Hou, and H. Si. Numerical study of Hele-Shaw flow with suction. Physics ofFluids , :2471–2486, 1999.[20] S. J. Chapman. On the role of Stokes lines in the selection of Saffman-Taylor fingers with small surfacetension. European Journal of Applied Mathematics , :513–534, 1999.[21] S. J. Chapman and J. R. King. The selection of Saffman–Taylor fingers by kinetic undercooling. Journalof Engineering Mathematics , :1–32, 2003.[22] C.-Y. Chen, Y.-S. Huang, and J. A. Miranda. Radial Hele-Shaw flow with suction: fully nonlinear patternformation. Physical Review E , :053006, 2014.[23] J.-D. Chen. Radial viscous fingering patterns in Hele-Shaw cells. Experiments in Fluids , :363–371,1987.[24] S. Chen, B. Merriman, S. Osher, and P. Smereka. A simple level set method for solving Stefan problems. Journal of Computational Physics , :8–29, 1997.[25] A. Coco and G. Russo. Finite-difference ghost-point multigrid methods on Cartesian grids for ellipticproblems in arbitrary domains. Journal of Computational Physics , :464–501, 2013.[26] R. Combescot, T. Dombre, V. Hakim, Y. Pomeau, and A. Pumir. Shape selection of Saffman-Taylorfingers. Physical Review Letters , :2036, 1986.[27] ´I. M. Coutinho and J. A. Miranda. Control of viscous fingering through variable injection rates and time-dependent viscosity fluids: Beyond the linear regime. Phys. Rev. E , :063102, 2020.[28] D. Crowdy and S. Tanveer. The effect of finiteness in the Saffman–Taylor viscous fingering problem. Journal of Statistical Physics , 114:1501–1536, 2004.[29] L. J. Cummings and J. R. King. Hele–Shaw flow with a point sink: generic solution breakdown.
EuropeanJournal of Applied Mathematics , :1–37, 2004.[30] L. M. Cummings, S. D. Howison, and J. R. King. Two-dimensional Stokes and Hele-Shaw flows with freesurfaces. European Journal of Mechanics , :635–680, 1999.[31] C. Cuttle, D. Pihler-Puzovi´c, and A. Juel. Dynamics of front propagation in a compliant channel. Journalof Fluid Mechanics , 886:A20, 2020.[32] W.-S. Dai and M. J. Shelley. A numerical study of the effect of surface tension and noise on an expandingHele–Shaw bubble.
Physics of Fluids , :2131–2146, 1993.[33] M. C. Dallaston and S. W. McCue. New exact solutions for Hele-Shaw flow in doubly connected regions. Physics of Fluids , 24:052101, 2012.[34] M. C. Dallaston and S. W. McCue. Bubble extinction in Hele-Shaw flow with surface tension and kineticundercooling regularization.
Nonlinearity , :1639–1665, 2013.[35] M. C. Dallaston and S. W. McCue. Corner and finger formation in Hele–Shaw flow with kinetic under-cooling regularisation. European Journal of Applied Mathematics , :707–727, 2014.[36] A. J. DeGregoria and L. W. Schwartz. A boundary-integral method for two-phase displacement in Hele-Shaw cells. Journal of Fluid Mechanics , :383–400, 1986.[37] E. O. Dias and J. Miranda. Control of radial fingering patterns: A weakly nonlinear approach. PhysicalReview E , :016312, 2010. 28 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT [38] E. O. Dias, F. Parisio, and J. A. Miranda. Suppression of viscous fluid fingering: A piecewise-constantinjection process.
Physical Review E , :067301, 2010.[39] E. O. Dias, E. Alvarez-Lacalle, M. S. Carvalho, and J. A. Miranda. Minimization of viscous fluid fingering:a variational scheme for optimal flow rates. Physical Review Letters , :144502, 2012.[40] U. Ebert, B. Meulenbroeck, and L. Sch¨afer. Convective stabilization of a Laplacian moving boundaryproblem with kinetic undercooling. SIAM Journal on Applied Mathematics , :292–310, 2007.[41] D. Enright, R. Fedkiw, J. Ferziger, and I. Mitchell. A hybrid particle level set method for improvedinterface capturing. Journal of Computational Physics , :83–116, 2002.[42] A. Eslami, R. Basak, and S. M. Taghavi. Multiphase viscoplastic flows in a nonuniform Hele-Shaw cell:a fluidic device to control interfacial patterns. Industrial & Engineering Chemistry Research , :4119–4133, 2020.[43] P. Fast, L. Kondic, M. J. Shelley, and P. Palffy-Muhoray. Pattern formation in non-Newtonian Hele–Shawflow. Physics of Physics , :1191, 2001.[44] R. Folch, E. Alvarez-Lacalle, J. Ort´ın, and J. Casademunt. Pattern formation and interface pinch-off inrotating Hele-Shaw flows: A phase-field approach. Physical Review E , :056305, 2009.[45] J. V. Fontana, E. O. Dias, and J. A. Miranda. Controlling and minimizing fingering instabilities in non-Newtonian fluids. Physical Review E , :013016, 2014.[46] A. Franco-G´omez, A. B. Thompson, A. L. Hazel, and A. Juel. Sensitivity of Saffman–Taylor fingers tochannel-depth perturbations. Journal of Fluid Mechanics , :343–368, 2016.[47] T. Gao, M. Mirzadeh, P. Bai, K. M. Conforti, and M. Z. Bazant. Active control of viscous fingering usingelectric fields. Nature Communications , :1–8, 2019.[48] B. P. J. Gardiner, S. W. McCue, M. C. Dallaston, and T. J. Moroney. Saffman-Taylor fingers with kineticundercooling. Physical Review E , :023016, 2015.[49] B. P. J. Gardiner, S. W. McCue, and T. J. Moroney. Discrete families of Saffman-Taylor fingers with exoticshapes. Results in Physics , :103–104, 2015.[50] F Gibou, E P Fedkiw, L.-T. Cheng, and M Kang. A second-order-accurate symmetric discretization of thePoisson equation on irregular domains. Journal of Computational Physics , :205–227, 2002.[51] F. Gibou, R. Fedkiw, R. Caflisch, and S. Osher. A level set approach for the numerical simulation ofdendritic growth. Journal of Scientific Computing , :183–199, 2003.[52] F. Gibou, R. Fedkiw, and S. Osher. A review of level-set methods and some recent applications. Journalof Computational Physics , :82–109, 2018.[53] D. Givoli. Numerical methods for problems in infinite domains , volume 33 of
Studies in Applied Mechan-ics . Elsevier, 2013.[54] C. C. Green, C. J. Lustri, and S. W. McCue. The effect of surface tension on steadily translating bub-bles in an unbounded Hele-Shaw cell.
Proceedings of the Royal Society A: Mathematical, Physical andEngineering Sciences , :20170050, 2017.[55] B. Gustafsson and A Vasil’ev. Conformal and potential analysis in Hele-Shaw cells . Birkh¨auser-Verlag,2006.[56] Y. E. Hohlov and S. D. Howison. On the classification of solutions to the zero-surface-tension model forHele–Shaw free-boundary flows.
Quarterly of Applied Mathematics , :777–789, 1993.29 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT [57] Y. E. Hohlov, S. D. Howison, C. Huntingford, J. R. Ockendon, and A. A. Lacey. A model for nonsmoothfree boundaries in Hele–Shaw flows.
The Quarterly Journal of Mechanics and Applied Mathematics , :107–128, 1994.[58] G. M. Homsy. Viscous fingering in porous media. Annual Review of Fluid Mechanics , :271–311, 1987.[59] D. C. Hong and F. Family. Bubbles in the Hele–Shaw cell — pattern selection and tip perturbations. Physical Review A , :5253–5259, 1988.[60] D. C. Hong and J. S. Langer. Analytic theory of the selection mechanism in the Saffman-Taylor problem. Physical Review Letters , :2032, 1986.[61] T. Y. Hou, J. S. Lowengrub, and M. J. Shelley. Removing the stiffness from interfacial flow with surfacetension. Journal of Computational Physics , :312–338, 1994.[62] T. Y. Hou, Z. Li, S. Osher, and H. Zhao. A hybrid method for moving interface problems with applicationto the Hele–Shaw flow. Journal of Computational Physics , :236–252, 1997.[63] T. Y. Hou, J. S. Lowengrub, and M. J. Shelley. Boundary Integral Methods for Multicomponent Fluidsand Multiphase Materials. Journal of Computational Physics , :302–362, 2001.[64] S. D. Howison. Bubble growth in porous media and Hele–Shaw cells. Proceedings of the Royal Societyof Edinburgh Section A: Mathematics , :141–148, 1986.[65] S. D. Howison. Cusp Development in Hele–Shaw Flow with a Free Surface. SIAM Journal on AppliedMathematics , :20–26, 1986.[66] S. D. Howison. Fingering in Hele-Shaw cells. Journal of Fluid Mechanics , :439–453, 1986.[67] S. D. Howison, J. R. Ockendon, and A. A. Lacey. Singularity development in moving boundary problems. The Quarterly Journal of Mechanics and Applied Mathematics , :343–354, 1985.[68] S. J. Jackson, D. Stevens, H. Power, and D. Giddings. A boundary element method for the solution offinite mobility ratio immiscible displacement in a Hele-Shaw cell. International Journal for NumericalMethods in Fluids , :521–551, 2015.[69] S. J. Jackson, H. Power, D. Giddings, and D. Stevens. The stability of immiscible viscous fingeringin Hele-Shaw cells with spatially varying permeability. Computer Methods in Applied Mechanics andEngineering , :606–632, 2017.[70] E. D. Kelly and E. J. Hinch. Numerical simulations of sink flow in the Hele-Shaw cell with small surfacetension. European Journal of Applied Mathematics , :533–550, 1997.[71] L. Kondic, M. J. Shelley, and P. Palffy-Muhoray. Non-Newtonian Hele–Shaw flow and the Saffman–Taylor instability. Physical Review Letters , :1433–1436, 1998.[72] A. A. Lacey. Moving boundary problems in the flow of liquid through porous media. ANZIAM Journal , :171–193, 1982.[73] A. Leshchiner, M. Thrasher, M. B. Mineev-Weinstein, and H. L. Swinney. Harmonic moment dynamicsin Laplacian growth. Physical Review E , :016206, 2010.[74] S. Li, J.S. Lowengrub, P. H. Leo, and V. Cristini. Nonlinear theory of self-similar crystal growth andmelting. Journal of Crystal Growth , :703–713, 2004.[75] S. Li, J. S. Lowengrub, and P. H. Leo. A rescaling scheme with application to the long-time simulation ofviscous fingering in a Hele–Shaw cell. Journal of Computational Physics , :554–567, 2007.30 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT [76] S. Li, J. S. Lowengrub, J. Fontana, and P. Palffy-Muhoray. Control of viscous fingering patterns in a radialHele-Shaw cell.
Physical Review Letters , :174501, 2009.[77] S. Liang. Random-walk simulations of flow in Hele Shaw cells. Physical Review A , :2663, 1986.[78] A. Lindner, D. Derks, and M. J. Shelley. Stretch flow of thin layers of Newtonian liquids: Fingeringpatterns and lifting forces. Physics of Fluids , :072107, 2005.[79] T. F. Lins and J. Azaiez. Resonance-like dynamics in radial cyclic injection flows of immiscible fluids inhomogeneous porous media. Journal of Fluid Mechanics , :713–729, 2017.[80] J. R. Lister, G. G. Peng, and J. A. Neufeld. Viscous control of peeling an elastic sheet by bending andpulling. Physical Review Letters , :154501, 2013.[81] C. J. Lustri, C. C. Green, and S. W. McCue. Selection of a Hele-Shaw bubble via exponential asymptotics. SIAM Journal on Applied Mathematics , :289–311, 2020.[82] K. V. McCloud and J. V. Maher. Experimental perturbations to Saffman-Taylor flow. Physics Reports , :139–185, 1995.[83] S. W. McCue. Short, flat-tipped, viscous fingers: Novel interfacial patterns in a Hele-Shaw channel withan elastic boundary. Journal of Fluid Mechanics , :1–4, 2018.[84] S. W. McCue and J. R. King. Contracting bubbles in Hele-Shaw cells with a power-law fluid. Nonlinearity , :613–641, 2011.[85] J. W. McLean and P. G. Saffman. The effect of surface tension on the shape of fingers in a Hele Shaw cell. Journal of Fluid Mechanics , :455–469, 1981.[86] M. Mineev–Weinstein. Selection of the Saffman–Taylor finger width in the absence of surface tension:an exact result. Physical Review Letters , :2113–2116, 1998.[87] M. Mineev–Weinstein, P. B. Wiegmann, and A. Zabrodin. Integrable structure of interface dynamics. Physical Review Letters , :5106–5109, 2000.[88] J. Miranda and M. Widom. Radial fingering in a Hele-Shaw cell: a weakly nonlinear analysis. Physica D , :315–328, 1998.[89] M. Mirzadeh and M. Z. Bazant. Electrokinetic control of viscous fingering. Physical Review Letters , :174501, 2017.[90] T. J. Moroney, D. R. Lusmore, S. W. McCue, and D. L. S. McElwain. Extending fields in a level setmethod by solving a biharmonic equation. Journal of Computational Physics , :170–185, 2017.[91] L. C. Morrow, M. C. Dallaston, and S. W. McCue. Interfacial dynamics and pinch-off singularities foraxially symmetric darcy flow. Physical Review E , :053109, 2019.[92] L. C. Morrow, J. R. King, T. J. Moroney, and S. W. McCue. Moving boundary problems for quasi-steadyconduction limited melting. SIAM Journal on Applied Mathematics , :2107–2131, 2019.[93] L. C. Morrow, T. J. Moroney, and S. W. McCue. Numerical investigation of controlling interfacial insta-bilities in non-standard Hele-Shaw configurations. Journal of Fluid Mechanics , :1063–1097, 2019.[94] J. Nase, D. Derks, and A. Lindner. Dynamic evolution of fingering patterns in a lifted Hele–Shaw cell. Physics of Fluids , :123101, 2011.[95] Q. Nie and F. R. Tian. Singularities in Hele–Shaw flows. SIAM Journal on Applied Mathematics , :34–54, 1998.[96] S. Osher and R. Fedkiw. Level set methods and dynamic implicit surfaces , volume 153. Springer, 2003.31 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT [97] S. Osher and R. P. Fedkiw. Level set methods: an overview and some recent results.
Journal of Compu-tational Physics , :463–502, 2001.[98] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based onHamilton-Jacobi formulations. Journal of Computational Physics , :12–49, 1988.[99] A. S. S. Paiva, S. H. A. Lira, and R. F. S. Andrade. Non-linear effects in a closed rotating radial Hele-Shawcell. AIP Advances , :025121, 2019.[100] C.-W. Park and G. M. Homsy. Two-phase displacement in Hele Shaw cells: theory. Journal of FluidMechanics , :291–308, 1984.[101] L. Paterson. Radial fingering in a Hele Shaw cell. Journal of Fluid Mechanics , :513–529, 1981.[102] D. Pihler-Puzovi´c, P. Illien, M. Heil, and A. Juel. Suppression of complex fingerlike patterns at theinterface between air and a viscous fluid by elastic membranes. Physical Review Letters , :074502,2012.[103] D. Pihler-Puzovi´c, R. P´erillat, M. Russell, A. Juel, and M. Heil. Modelling the suppression of viscousfingering in elastic-walled Hele-Shaw cells. Journal of Fluid Mechanics , :162–183, 2013.[104] D Pihler-Puzovi´c, A. Juel, and M. Heil. The interaction between viscous fingering and wrinkling inelastic-walled Hele-Shaw cells. Physics of Fluids , :022102, 2014.[105] D. Pihler-Puzovi´c, G. G. Peng, J. R. Lister, M. Heil, and A. Juel. Viscous fingering in a radial elastic-walled Hele-Shaw cell. Journal of Fluid Mechanics , :163–191, 2018.[106] N. B. Pleshchinskii and M. Reissig. Hele-shaw flows with nonlinear kinetic undercooling regularization. Nonlinear Analysis , :191–203, 2002.[107] P. Y. Polubarinova-Kochina. On motion of the contour of an oil layer. In Dokl. Akad. Nauk SSSR , volume , pages 254–257, 1945.[108] H. Power, D. Stevens, K. A. Cliffe, and A. Golin. A boundary element study of the effect of surface dis-solution on the evolution of immiscible viscous fingering within a Hele-Shaw cell. Engineering Analysiswith Boundary Elements , :1318–1330, 2013.[109] S. Richardson. Hele–Shaw flows with a free boundary produced by the injection of fluid into a narrowchannel. Journal of Fluids Mechanics , :609–618, 1972.[110] J. E. Sader, D. Y. C. Chan, and B. D. Hughes. Non-Newtonian effects on immiscible viscous fingering ina radial Hele–Shaw cell. Physical Review E , :420–432, 1994.[111] P. G. Saffman. Viscous fingering in Hele-Shaw cells. Journal of Fluid Mechanics , :73–94, 1986.[112] P. G. Saffman and G. I. Taylor. The penetration of a fluid into a porous medium or hele-shaw cell con-taining a more viscous liquid. Proceedings of the Royal Society of London. Series A. Mathematical andPhysical Sciences , :312–329, 1958.[113] L. W. Schwartz. Instability and fingering in a rotating Hele–Shaw cell or porous medium. Physics ofFluids , :167–169, 1989.[114] J. A. Sethian. Level set methods and fast marching methods: evolving interfaces in computational ge-ometry, fluid mechanics, computer vision, and materials science , volume 3. Cambridge University Press,Cambridge, UK, 1993.[115] J. A. Sethian and P. Smereka. Level set methods for fluid interfaces.
Annual Review of Fluid Mechanics , :341–372, 2003. 32 review of one-phase Hele–Shaw flows and a level-set method for non-standard configurations A PREPRINT [116] M. J. Shelley, F. Tian, and K. Wlodarski. Hele-Shaw flow and pattern formation in a time-dependent gap.
Nonlinearity , :1471, 1997.[117] B. I. Shraiman. Velocity selection and the Saffman–Taylor problem. Physical Review Letters , :2028–2031, 1986.[118] P. Tabeling, G. Zocchi, and A. Libchaber. An experimental study of the Saffman-Taylor instability. Jour-nal of Fluid Mechanics , :67–82, 1987.[119] S. Tanveer. New solutions for steady bubbles in a Hele–Shaw cell. Physics of Fluids , :651–658, 1987.[120] S. Tanveer. Analytic theory for the selection of a symmetric Saffman-Taylor finger in a Hele-Shaw cell. Physics of Fluids , :1589–1605, 1987.[121] S. Tanveer. Analytic theory for the determination of velocity and stability of bubbles in a Hele-Shaw cell. Theoretical and Computational Fluid Dynamics , :135–163, 1989.[122] S. Tanveer. Surprises in viscous fingering. Journal of Fluid Mechanics , :273–308, 2000.[123] S. Tanveer and P. G. Saffman. Stability of bubbles in a Hele–Shaw cell. Physics of Fluids , :2624–2635,1987.[124] H. Thom´e, M. Rabaud, V. Hakim, and Y. Couder. The Saffman–Taylor instability: From the linear to thecircular geometry. Physics of Fluids , :224–240, 1989.[125] A. B. Thompson, A. Juel, and A. L. Hazel. Multiple finger propagation modes in Hele-Shaw channels ofvariable depth. Journal of Fluid Mechanics , :123–164, 2014.[126] J.-M. Vanden-Broeck. Fingers in a Hele-Shaw Cell with surface tension. Physics of Fluids , 26:2033,1983.[127] C. Vaquero-Stainer, M. Heil, A. Juel, and D. Pihler-Puzovi´c. Self-similar and disordered front propagationin a radial Hele-Shaw channel with time-varying cell depth.
Physical Review Fluids , :064002, 2019.[128] G. Vasconcelos. Exact solutions for steady bubbles in a Hele-Shaw cell with rectangular geometry. Journalof Fluid Mechanics , :175–198, 2001.[129] G. L. Vasconcelos. Multiple bubbles in a Hele–Shaw cell. Physical Review E , :R3306–R3309, 1994.[130] A. Vasil’ev. From the hele-shaw experiment to integrable systems: A historical overview. Compl. Anal.Oper. Theor. , :551–585, 2009.[131] S. L. Waters and L. J. Cummings. Coriolis effects in a rotating Hele-Shaw cell. Physics of Fluids , :048101, 2005.[132] X. Xie. Rigorous results in existence and selection of Saffman–Taylor fingers by kinetic undercooling. European Journal of Applied Mathematics , :63–116, 2019.[133] H. Zhao, J. Casademunt, C. Yeung, and J. V. Maher. Perturbing Hele-Shaw flow with a small gap gradient. Physical Review A , :2455, 1992.[134] M. Zhao, X. Li, W. Ying, A. Belmonte, J. Lowengrub, and S. Li. Computation of a Shrinking Interface ina Hele-Shaw Cell. SIAM Journal on Scientific Computing , :B1206–B1228, 2018.[135] Z. Zheng, H. Kim, and H. A. Stone. Controlling viscous fingering using time-dependent strategies. Phys-ical Review Letters ,115