2.1 Generalised Hele-Shaw model

We consider a generalised one-phase model of Hele-Shaw flow where the gap between the plates is either spatially or temporally dependent such that $b \to b(\boldsymbol {x}, t)$ and the Hele-Shaw plates can rotate with angular velocity $\bar {\omega }$ . We suppose an inviscid bubble is injected into the viscous fluid at rate $Q(t)$ , and denote the domain occupied by the inviscid fluid by $\Omega (t)$ . The interface separating the inviscid bubble and the viscous fluid is denoted by $\partial \Omega (t)$ . As is commonplace with Hele-Shaw flows, we consider a depth-averaged model that comes about from averaging Stokes flow over the gap between the plates, which itself is assumed to be small.

Denoting the pressure, viscosity and density of the viscous fluid by P, $\mu $ and $\rho $ respectively, and denoting the angular velocity of the Hele-Shaw cell by $\bar {\omega }$ , the governing equations for the velocity of the viscous fluid are

(2.1) $$ \begin{align} &\boldsymbol{v} = -\frac{b^2}{12 \mu} (\boldsymbol{\nabla} P - \bar{\omega}^2 \rho r \boldsymbol{e}_r), \quad\boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t), \end{align} $$

(2.2) $$ \begin{align} &\boldsymbol{\nabla} \cdot (b \boldsymbol{v}) = - \frac{\partial b}{\partial t}, \quad\boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t), \end{align} $$

where $r = |\boldsymbol {x}|$ . Equation (2.1) is Darcy’s law modified to include the rotational effects of the Hele-Shaw cell, while (2.2) ensures that the mass of the fluid is conserved. Defining a reduced pressure according to $p = P - \bar {\omega } \rho r^2 / 2$ and then substituting (2.1) into (2.2) generates the governing equation for pressure,

(2.3) $$ \begin{align} \boldsymbol{\nabla} \cdot \bigg( \frac{b^3}{12 \mu} \boldsymbol{\nabla} p \bigg) = \frac{\partial b}{\partial t},\quad \boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t). \end{align} $$

When the gap between the plates is both spatially and temporally uniform, (2.3) reduces to Laplace’s equation $\nabla ^2 p = 0$ . We have two boundary conditions on the fluid–fluid interface given by

(2.4) $$ \begin{align} p &= -\gamma \bigg( \kappa + \frac{2}{b} \bigg) - \frac{\rho \bar{\omega}^2 r^2}{2}, \quad\boldsymbol{x} \in \partial \Omega(t) , \end{align} $$

(2.5) $$ \begin{align} v_n &= -\frac{b^2}{12 \mu} \frac{\partial p}{\partial n}, \quad\boldsymbol{x} \in \partial \Omega(t) , \end{align} $$

where $\gamma $ is the surface tension, $\kappa $ is the signed curvature of $\partial \Omega $ , $\rho $ is the density of the viscous fluid, and $v_n$ is the normal speed of the interface. The dynamic boundary condition (2.4) incorporates both the effects of surface tension and the rotation of the Hele-Shaw plates. The kinematic boundary condition (2.5) relates the velocity of the viscous fluid with the normal velocity of the interface. We also have the far-field boundary condition

(2.6) $$ \begin{align} \frac{b^3}{12 \mu}\frac{\partial p}{\partial r} &\sim -\frac{Q}{2 \pi r} + \frac{1}{2} r \frac{\partial b}{\partial t}, \quad r \to \infty, \end{align} $$

which acts as a source/sink term at infinity. For $Q> 0$ ( $Q < 0$ ), this condition corresponds to the bubble area expanding (contracting) at rate Q. The inclusion of the $\partial b / \partial t$ in (2.6) comes about from the nonhomogeneous term in (2.3), and ensures the change of rate of volume of the bubble is Q.

To nondimensionalise (2.3)-(2.6), we introduce $r_0$ as the average initial radius of the bubble and $Q_0$ as the average injection rate over the duration of a simulation, and $b_0 = b(0,0)$ . Then space, time, gap width, pressure and velocity are scaled according to

(2.7) $$ \begin{align} \hat{\boldsymbol{x}}=\frac{\boldsymbol{x}}{r_0}, \quad \hat{t}=\frac{t}{T}, \quad \hat{b}=\frac{b}{b_0}, \quad \hat{p}=\frac{b_0^2T}{12\mu r_0^2}p, \quad \hat{\boldsymbol{v}}=\frac{T}{r_0}\boldsymbol{v}, \end{align} $$

respectively, where T is a representative time-scale. Dropping the hats and retaining our original variable names, (2.3)–(2.6) become

(2.8) $$ \begin{align} &\nabla \cdot ( b^3 \nabla p ) = \frac{\partial b}{\partial t},\, \qquad\qquad\qquad\boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t),\end{align} $$

(2.9) $$ \begin{align} &p = -\sigma \bigg( \kappa + \frac{2 R_0}{b} \bigg) - \omega^2 r^2, \qquad\boldsymbol{x} \in \partial \Omega(t), \end{align} $$

(2.10) $$ \begin{align} &v_n = -b^2 \frac{\partial p}{\partial n}, \,\quad\qquad\qquad\qquad\boldsymbol{x} \in \partial \Omega(t), \end{align} $$

(2.11) $$ \begin{align} &b^3 \frac{\partial p}{\partial r} \sim -\frac{Q}{2 \pi r} + \frac{1}{2} r \frac{\partial b}{\partial t} \ \qquad r \to \infty, \end{align} $$

where $\sigma = b_0^2 T\gamma / 12 \mu r_0^3$ , $R_0 = r_0/b_0$ , and $\omega ^2 = \rho b_0^2 T\bar {\omega }^2 / 24 \mu $ . For this configuration, an appropriate time-scale could be $T=r_0^2b_0/Q_0$ , in which case the dimensionless average injection rate would become $Q_0=1$ .

In addition to this expanding bubble problem, we shall be concerned with two other scenarios, namely the blob geometry, where viscous fluid occupies $\Omega (t)$ and is withdrawn from a point or the cell rotates around a perpendicular axis, and the channel geometry, where viscous fluid occupies a long rectangular channel and is displaced by the inviscid fluid that is injected at one end. For these two scenarios, modifications to (2.8)–(2.11) will be made as appropriate.

2.2 Complex variable formulation

Before outlining our numerical scheme in Section 3, we take some time to illustrate some of the mathematical properties of the Hele-Shaw problem, especially in the special case of zero surface tension. This mathematical exploration, which relies heavily on complex variable theory and conformal mapping, is based on many previous studies in this spirit [Reference Crowdy and Tanveer29–Reference Cummings, Howison and King31, Reference Dallaston and McCue34, Reference Hohlov and Howison61, Reference Hohlov, Howison, Huntingford, Ockendon and Lacey62, Reference Howison69–Reference Howison71, Reference Mineev-Weinstein92, Reference Mineev-Weinstein, Wiegmann and Zabrodin93]. To keep the discussion contained and to connect with numerical simulations described later in this paper, we restrict ourselves to examples of three geometries (the expanding bubble problem, the blob problem and the channel geometry).

In the standard set-up in which $b=1$ (parallel, stationary plates) and $\omega =0$ (no rotation), equations 2.8–2.11 respectively reduce to

(2.12) $$ \begin{align} &\nabla^2 p = 0, \quad \ \, \qquad\qquad\boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t), \end{align} $$

(2.13) $$ \begin{align}&p = -\sigma\kappa, \, \qquad\qquad\boldsymbol{x} \in \partial \Omega(t), \end{align} $$

(2.14) $$ \begin{align} &v_n = - \frac{\partial p}{\partial n}, \ \ \ \quad\qquad\boldsymbol{x} \in \partial \Omega(t), \end{align} $$

(2.15) $$ \begin{align} &\frac{\partial p}{\partial r} \sim -\frac{Q}{2 \pi r}, \qquad r \to \infty. \end{align} $$

It is instructive to reformulate this problem using complex variable methods as follows. Given the fluid pressure p satisfies Laplace’s equation (2.12), it can be interpreted as the negative real part of an analytic function $W(z,t)=-p(x,y,t)+\mathrm {i}\psi (x,y,t)$ of the complex variable $z=x+\mathrm {i}y$ . Here, W is acting as a complex potential, while $\psi $ is a streamfunction.

Further, there exists a time-dependent conformal map $z=f(\zeta ,t)$ from the unit disc in the plane of an auxiliary variable $\zeta $ to the fluid region in the z-plane (that is, $\mathbb {R}^2 \backslash \Omega (t)$ ) and the unit circle $|\zeta |=1$ to the fluid interface $\partial \Omega (t)$ , as depicted schematically in Figure 1. The map will be univalent (that is, one-to-one) and analytic in the unit disc except for at a single point, which we choose to be $\zeta =0$ , that represents $z\to \infty $ . In the limit $\zeta \rightarrow 0$ , we have $f\sim a(t)/\zeta $ . By fixing a rotational degree of freedom we force $a(t)$ to be real. Now the complex potential $W(z,t)$ is also an analytic function of $\zeta $ and so we write $w(\zeta ,t)=W(f(\zeta ,t),t)$ . Given the far-field condition (2.15), which implies $W\sim (Q/2\pi )\log z$ as $|z|\rightarrow \infty $ , we now have the local behaviour $w\sim -(Q/2\pi )\log \zeta $ as $\zeta \rightarrow 0$ .

Figure 1 A schematic diagram indicating the time-dependent complex mapping $f(\zeta ,t)$ from the auxiliary $\zeta $ plane to the physical $z (=x+\mathrm i y)$ plane. The interface $\delta \Omega (t)$ is the image of the unit circle $|\zeta |=1$ , and the complex representation of the unit normal vector can be expressed in terms of the derivative of f.

To formulate the kinematic condition (2.14) in terms of the map f, it is useful to introduce some complex variable equivalents of standard concepts from vector algebra. Firstly, a complex number can be used to represent a vector (with components given by the real and imaginary parts), such as the normal to an interface. The unit normal to the unit circle $|\zeta |=1$ is $n^{(\zeta )} = \zeta $ , and given the interface $\partial \Omega $ is the image of the unit circle under $z = f(\zeta ,t)$ , the normal $n^{(z)}$ in the z-plane is found by rotating $n^{(\zeta )}$ by the argument of $f_\zeta $ , thus $n^{(z)} = \zeta f_\zeta /|\zeta f_\zeta |$ (see Figure 1). Secondly, the equivalent of the dot product between two complex numbers a and b is $\Re \{a\overline b\}$ . The time derivative of the map $f_t$ at a point on the unit disc gives a velocity vector of a point on the interface $\partial \Omega $ . Therefore the normal velocity $v_n$ of the interface $\delta \Omega $ as a function of $\zeta $ is given by $\Re \{f_t\overline {\zeta f_\zeta }\}/|\zeta f_\zeta |$ , while the normal derivative $\partial p/\partial n$ is given by $-\Re \{\zeta W_z f_\zeta \}/|f_\zeta | = \Re \{\zeta w_\zeta \}/|f_\zeta |$ . This calculation allows the kinematic condition to be represented as

$$ \begin{align*}\Re\{f_t\overline{\zeta f_\zeta}\}=\Re\{\zeta w_\zeta\}, \quad |\zeta|=1. \end{align*} $$

Now to reformulate the dynamic condition (2.13) we note the curvature on the fluid boundary can be written as $\Re \{\zeta (\zeta f_\zeta )_\zeta \overline {\zeta f_\zeta }\}/|\zeta f_\zeta |^3$ on the unit circle. Given (2.13) and the logarithmic behaviour of w as $\zeta \rightarrow 0$ , we can write $w=-(Q/2\pi )\log \zeta -\sigma \mathcal {K}(\zeta ,t)$ , where $\mathcal {K}(\zeta ,t)$ is an analytic function of $\zeta $ in the unit disc whose real part on the unit circle $|\zeta |=1$ is given by the curvature $\kappa $ . Combining these ideas, we arrive at the single governing equation

(2.16) $$ \begin{align} \Re\{f_t\overline{\zeta f_\zeta}\}=-\frac{Q}{2\pi}-\sigma\Re\{\zeta \mathcal{K}_\zeta\}, \quad |\zeta|=1. \end{align} $$

This equation is often referred to as the Polubarinova–Galin equation, especially when surface tension is ignored [Reference Gustafsson, Vasiĺev, Galdi, Heywood and Rannacher60].

We shall briefly outline five illustrative examples, chosen to demonstrate the key features for the special case in which surface tension is ignored. We shall later refer back to these examples when we implement our numerical scheme with surface tension included. First, we shall consider the mapping $f=a(t)/\zeta +b(t)\zeta ^N$ , $N\geq 2$ [Reference Howison69]. By substituting into (2.16) with $\sigma =0$ , we find that a and b must satisfy the coupled system of ordinary differential equations (ODEs) $a\dot {a}-Nb\dot {b}=Q/2\pi $ , $N\dot {a}b-a\dot {b}=0$ . Say, for definiteness, $a(0)=1$ and $b(0)=\epsilon $ , then the second of these equations gives $b=\epsilon a^N$ , while the first equation integrates to give the time-dependence

$$ \begin{align*}t=\frac{\pi}{Q}(a^2-1-\epsilon^2N(a^{2N}-1)). \end{align*} $$

Thus we have an exact solution, as shown by the red dashed curves in the top panel of Figure 2(a). The innermost curve is the initial bubble boundary. Here we have chosen $N=5$ , $Q=2\pi $ and $\epsilon =0.01$ , so this inner curve is a circle with a sixfold perturbation. As time increases, the bubble expands and starts to develop six small fingers.

Figure 2 Solutions with and without surface tension. The red dashed curves in (a)–(e) are zero-surface-tension solutions described by the five examples in Section 2.2. The solid blue curves are numerical solutions, including surface tension, for the same initial conditions. For the examples in (a), (c) and (d), the zero-surface-tension solutions involve a form of finite-time blow-up characterised by cusps forming on the interface; the inclusion of surface 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 of a 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 $\zeta $ -plane with the unit circle and critical points and logarithmic singularities indicated by solid red dots and black diamonds, respectively. Note that in (b) we do not plot the critical points at $t=0$ as they are outside the field of view here. For (a)–(e), as time increases, the critical points and logarithmic singularities move towards the unit circle (colour available online).

It is of interest to track the critical points, $\zeta =\zeta ^*$ , which are the points at which $f_\zeta =0$ . For this example, $\zeta ^*=(1/\epsilon N a^{N-1})^{1/(N+1)}$ . Clearly there are $N+1$ critical points that are equally spaced along a circle that is outside the unit circle in the $\zeta $ -plane. As time evolves, each of these critical points moves in a straight line towards the origin and intersects the unit circle when $|\zeta ^*|=1$ , that is, $a=1/(\epsilon N)^{1/(N-1)}$ . We can compute the exact time when this occurs, namely,

$$ \begin{align*}t^*=\frac{\pi}{Q}\bigg(\frac{N-1}{N(\epsilon N)^{2/(N-1)}}-1+\epsilon^2N \bigg)\,. \end{align*} $$

For the case in Figure 2(a), for which $N=5$ , we can see the six critical points (red dots) in the bottom panel moving towards the unit circle. At $t=t^*$ , there is finite-time blow-up, characterised by six cusps of order $3/2$ along the bubble boundary, which we can see in the top panel of Figure 2(a) (a cusp of order $3/2$ is characterised by a curvature singularity that appears locally like two branches of $y^2=x^3$ meeting at a cusp, suitably scaled and 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 the map $f=a(t)/\zeta +\log (\zeta -d(t))-\log (\zeta +d(t))$ [Reference Howison69]. Again, the geometry is an expanding bubble, but this time the behaviour is qualitatively different. By substituting this map for f into (2.16) with $\sigma =0$ , we find that a and d satisfy the coupled system

$$ \begin{align*}\dot{a}=\frac{Q}{2\pi}\frac{ad^4-a+4d}{a^2d^4-(2d-a)^2}, \quad \dot{d}=-\frac{Q}{2\pi}\frac{d(d^4-1)}{a^2d^4-(2d-a)^2}. \end{align*} $$

For initial conditions, we choose both $a(0)>d(0)>1$ so the denominators are initially positive, which means that $\dot {a}>0$ and $\dot {d}<0$ for small time. To determine the precise time-dependent behaviour, one option is to integrate this system numerically, which demonstrates that $a(t)$ continues to increase while $d(t)$ decreases towards $d=1$ as t increases. To make progress analytically, by dividing one equation by the other, we can also derive a first-order ODE with exact solution

$$ \begin{align*}a=\frac{1}{d}\bigg[\ln\bigg(\frac{(d(0)^2-1)(d^2+1)}{(d(0)^2+1)(d^2-1)}\bigg)+a(0)d(0)\bigg], \end{align*} $$

which shows that $a\sim -\ln (d-1)$ as $d\rightarrow 1^+$ . Clearly the map has logarithmic singularities at $\zeta _s=\pm d$ and critical points where $\zeta ^*=\pm d\sqrt {a}/\sqrt {a-2d}$ . For sufficiently large time, we have $a-2d>0$ and so both $\zeta _s$ and $\zeta ^*$ lie on the real $\zeta $ -axis and move towards the unit circle as $t\rightarrow \infty $ . Since $|\zeta ^*|\sim d+d^2/a$ as $a\rightarrow \infty $ , we see the critical points $\zeta ^*$ are further away from the unit circle than the logarithmic singularities $\zeta _s$ , and therefore finite-time blow-up is avoided (each logarithmic singularity asymptotes to the unit circle, but does not cross it; see Howison et al. [Reference Howison, Ockendon and Lacey72], and the discussion on the channel problem at the end of this section).

This second example is illustrated in Figure 2(b) using $a(0)=4$ , $d(0)=10/3$ and $Q=2\pi $ . In the top panel, the exact 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 on the positive and negative x-axes. In the bottom panel, both the critical points (red dots) and logarithmic singularities (black dots) are indicated in the $\zeta $ -plane. As just described, even though the critical points move towards the unit circle, they are further away than the logarithmic singularities; therefore, the logarithmic singularities have the effect of shielding the critical points and preventing blow-up. An interesting observation is that each of the two fjords appears to take the shape of a classical Saffman–Taylor finger [Reference Saffman and Taylor118]. To see this, note that on the unit circle near $\zeta =1$ we can derive a local analysis by setting $\zeta =1+\mathrm {i}\eta $ , so that $f\sim a+\log (1-d+\mathrm {i}\eta )$ . Taking real and imaginary parts and then eliminating $\eta $ gives $x=\mathrm {constant}-\ln (\cos y)$ , which is the famous Saffman–Taylor finger shape with width $\pi $ .

The third example is probably the best-known example of an exact solution in Hele-Shaw flow. Here suppose the geometry is such that there is a blob of fluid in the Hele-Shaw cell, occupying a region $\Omega (t)$ , surrounded by inviscid fluid. If the fluid is withdrawn from a point in space (the origin, say), then the blob boundary contracts and we have the less viscous fluid displacing the more viscous fluid. The governing equations (2.12)–(2.14) apply in $\Omega (t)$ , while (2.12) is replaced by $\partial p/\partial r\sim Q/2\pi r$ as $r\rightarrow 0$ . The map we have in mind here for this example is the quadratic map $f=a(t)\zeta +b(t)\zeta ^2$ , where the initial conditions $a(0)=1$ , $b(0)=\epsilon \ll 1$ correspond to an initial blob boundary that is a perturbed circle [Reference Polubarinova-Kochina113]. Substituting this map into the Polubarinova–Galin equation with $\sigma =0$ leads to a coupled system of integrable odes with the exact solution

$$ \begin{align*}a^2b=\epsilon,\quad t=\frac{\pi}{Q}(1-a^2+2(\epsilon^2-b^2)). \end{align*} $$

There is one critical point $\zeta ^*=-a^3/2\epsilon $ , which is initially located in the $\zeta $ -plane at $-1/2\epsilon $ and moves towards the origin as time evolves and a decreases. At $t=t^*$ this critical point hits the unit circle, where

$$ \begin{align*}t^*=\frac{\pi}{Q}\bigg(1+2\epsilon^2-\frac{9}{8}(2\epsilon)^{2/3} \bigg), \end{align*} $$

causing a cusp of order $3/2$ to form on the blob boundary. This behaviour is illustrated in Figure 2(c) for $\epsilon =0.1$ and $Q=2\pi $ . in the left panel, the blob boundary is represented by the dashed red curves. Initially, this curve appears circular, as the perturbation is very small. For intermediate times, the left portion of the boundary begins contracting faster than the remainder of the boundary until the cusp forms, corresponding to finite-time blow-up. On the right panel of Figure 2(c) the location of the fixed point in the $\zeta $ -plane is represented (red dots) for the same three times that the boundary is drawn for in the left panel. Here we see that the critical point touches the unit circle at the precise time that finite-time blow-up occurs. Note that for polynomial maps like the one in this example, it has been proven that a cusp will always form before the interface reaches the sink [Reference Hohlov and Howison61]; other explicit solutions exist whose boundary evolves to the location of the sink before or at the same time as cusp formation (see [Reference Richardson115], for example). The time reversibility of the system (2.12)–(2.15) in the absence of surface tension ( $\sigma =0$ ) implies that the only initial condition that will lead to the removal of all fluid is a disc centred on the sink, that is, $f(\zeta ,t) = a(t)\zeta $ .

For completeness we include two more examples, providing only the key details. These examples are for the geometry of flow in a Hele-Shaw channel, which we fix to be $2\pi $ units wide. For the fourth example, the map takes $f=-\log \zeta +a(t)+b(t)\zeta $ , with initial conditions $a(0)=0$ , $b(0)=\epsilon \ll 1$ , corresponding to a slightly perturbed flat interface [Reference Howison, Ockendon and Lacey72]. 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

$$ \begin{align*}a-\ln b=-\ln\epsilon, \quad t=\frac{\pi}{Q}(2a-b^2+\epsilon^2). \end{align*} $$

There are critical points at $\zeta ^*=1/b$ that move towards the origin as b increases and ultimately intersect the unit circle at $t^*=\pi (2\ln (1/\epsilon )-1+\epsilon ^2)/Q$ , at which time a $3/2$ cusp forms on the interface. This example is illustrated in Figure 2(d), where the interface profiles are shown in the top panel as red dashed curves. In the bottom panel the critical points are indicated (red dots). Here $\epsilon =0.1$ and so the critical point is initially at $\zeta ^*=10$ and ultimately hits the unit circle at $t=t^*$ .

Finally, the fifth example is for $f=-\log \zeta +a(t)+\alpha \log (\zeta +d(t))$ with initial conditions $a(0)=0$ , $d(0)=1/\epsilon \gg 1$ [Reference Howison, Ockendon and Lacey72], which is analogous to the second example above. There is a critical point at $\zeta ^*=d/(\alpha -1)$ and a logarithmic singularity at $\zeta _s=-d$ . We shall not include the details here, but it is possible to derive a coupled system of odes for $a(t)$ and $d(t)$ that can be solved numerically or reduced further analytically by dividing one by the other. For $0<\alpha <2$ the singularity $\zeta _s$ is always smaller in magnitude, and therefore closer to the unit circle, than $\zeta ^*$ . As a result, it turns out that d is a decreasing function and that $\dot {d}\rightarrow 0^+$ as $d\rightarrow 1^+$ . Therefore, $\zeta ^*$ and $\zeta _s$ do not intersect the unit circle but in fact approach it asymptotically as $t\rightarrow \infty $ . As such, there is no cusp formation, and instead the proximity of $\zeta _s$ to the unit circle results in the interface forming a long finger, whose width is $(2-\alpha )\pi $ . In Figure 2(e) we present an example with $\epsilon =0.2$ , $\alpha =1.2$ . The interface in the top panel, given by the dashed red curves, clearly approaches a finger in shape. The bottom panel shows the critical point (red dots) moving towards the unit circle, but always further away than the logarithmic singularity (black dots).

In each of Figure 2(a) (top panel), (b) (top panel), (c) (left panel), (d) (top panel) and (e) (top panel), we have included numerical solutions drawn as solid blue curves that are computed using the same initial conditions but with a nonzero value of surface tension. While we discuss these (more physically realistic) solutions at various points later in the paper, it is worth repeating here that nonzero surface tension is required in the model in order to relate to the real physics of a Hele-Shaw experiment. The historical interest in complex singularities and finite-time blow-up is therefore mostly of a mathematical nature.

3.1 The level-set method

To numerically solve (2.8)–(2.11), following the methodology of Osher and Sethian [Reference Osher and Sethian104], we construct a level-set function $\phi $ such that the fluid–fluid interface $\partial \Omega $ is the zero level set of $\phi $ or

$$ \begin{align*} \partial \Omega(t) = \lbrace \boldsymbol{x} \,| \, \phi(\boldsymbol{x}, t) = 0 \rbrace. \end{align*} $$

If the interface has the normal speed $v_n$ , then we wish to construct a speed function, F, such that $v_n = F$ on $\boldsymbol {x} \in \partial \Omega (t)$ and is continuous over the entire computational domain. Thus $\phi $ satisfies the level-set equation

(3.1) $$ \begin{align} \frac{\partial \phi}{\partial t} + F | \boldsymbol{\nabla} \phi | = 0. \end{align} $$

We discuss how F is computed in Section 3.2. To solve (3.1), we approximate the spatial derivatives using a second-order essentially nonoscillatory scheme (see the papers by Osher and Fedkiw [Reference Osher and Fedkiw102, Ch. 3] and Sethian [Reference Sethian120, Ch. 6] for details), and integrate in time using second-order total variation diminishing Runge–Kutta, which is performed by taking two forward Euler steps

$$ \begin{align*} \tilde{\phi}^{(n + 1)} &= \phi^{(n)} - \Delta t F^{(n)} | \boldsymbol{\nabla} \phi^{(n)} |,\\ \phi^{(n + 2)} &= \tilde{\phi}^{(n+1)} - \Delta t F^{(n+1)} | \boldsymbol{\nabla} \tilde{\phi}^{(n+1)} |, \end{align*} $$

and then taking an averaging step $\phi ^{(n+1)} = ( \phi ^{(n)} + \phi ^{(n + 2)} ) / 2$ . We note that the inclusion of the second-order curvature term, $\kappa $ , in the dynamic boundary condition (2.4) would typically require $\Delta t \sim \Delta x^2$ . However, we find that for the results presented in this work, the surface tension parameter is sufficiently small such that we can maintain numerical stability by choosing $\Delta t = \Delta x / (4 \max |F|)$ .

The level-set function $\phi $ is initialised as a signed distance function satisfying

(3.2) $$ \begin{align} \phi = \begin{cases} d &\mbox{if } \boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t), \\ -d & \mbox{if } \boldsymbol{x} \in \Omega(t), \end{cases} \end{align} $$

where d is the minimum distance between $\boldsymbol {x}$ and $\partial \Omega $ , via the method of crossing times [Reference Osher and Fedkiw102, Ch. 7]. That is, we advect $\phi $ in the normal direction to the interface by solving (3.1) with $F = 1$ and determine the point in time where each value of $\phi $ crosses from positive to negative. This process is repeated for $F = -1$ . To reduce numerical error, which can result when the gradient of $\phi $ becomes excessively small or large, we periodically perform reinitialisation in order to keep $\phi $ approximately equal to a signed distance function, which satisfies $|\boldsymbol {\nabla } \phi | = 1$ , over the duration of a simulation. Reinitialisation is performed by solving

(3.3) $$ \begin{align} \frac{\partial \phi}{\partial \tau} + S(\phi) (|\boldsymbol{\nabla} \phi| - 1) = 0, \end{align} $$

where

$$ \begin{align*} S(\phi) = \frac{\phi}{\sqrt{\phi^2 + \Delta x^2}}, \end{align*} $$

to steady state. Here $\tau $ is a pseudo-time variable where $\Delta \tau = \Delta x / 5$ . We find that performing reinitialisation every five time steps is sufficiently frequent.

While the level-set method has successfully been used as a framework for studying a variety of moving boundary problems, a limitation of the method is that it can suffer from volume loss or gain in regions where the mesh is underresolved. In an effort to alleviate this problem, Enright et al. [Reference Enright, Fedkiw, Ferziger and Mitchel42] proposed the particle level-set method, 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. [Reference Enright, Fedkiw, Ferziger and Mitchel42] for a more comprehensive description, as well as examples illustrating the effectiveness of the particle level-set method.

The method works by placing massless particles in the regions where $\phi>0$ and $\phi <0$ , that is, on both sides of the interface, which are referred to positive and negative particles, respectively. We denote $r_p$ as the minimum distance between the interface and the particles’ location. The marker particles are advected according to

(3.4) $$ \begin{align} \frac{{d}\boldsymbol{x}_p}{{d} t} = F \boldsymbol{n} , \end{align} $$

where $\boldsymbol {x}_p$ is the location of the particle and $\boldsymbol {n} = \boldsymbol {\nabla } \phi / |\boldsymbol {\nabla } \phi |$ is a unit vector that reduces to the outward-facing normal 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 four adjacent nodes to the particle defined as $\phi _p(\boldsymbol {x}) = s_p(r_p - |\boldsymbol {x} - \boldsymbol {x}_p|)$ , where $s_p = 1$ and $-1$ for positive and negative particles, respectively. Using $\phi _p$ , $\phi $ is corrected using

$$ \begin{align*}\phi = \begin{cases} \phi^+ &\mbox{if } |\phi^+| \le |\phi^-| ,\\ \phi^- &\mbox{if } |\phi^+|> |\phi^-|, \end{cases}\end{align*} $$

where $\phi ^+ = \max (\phi _p, \phi ^+)$ and $\phi ^- = \min (\phi _p, \phi ^-)$ . This procedure is performed both when the level-set equation (3.1) is solved and when reinitialisation is performed.

3.2 Solving for F

From the kinematic boundary condition (2.10), we have the expression for the speed function

(3.5) $$ \begin{align} F = -b^2 \boldsymbol{\nabla} p \cdot \boldsymbol{n}, \quad \boldsymbol{x} \in \mathbb{R}^2 \backslash \Omega(t), \end{align} $$

recalling that F is required to solve the level-set equation (3.1) and $\boldsymbol {n} = \boldsymbol {\nabla } \phi / |\boldsymbol {\nabla } \phi |$ reduces to the outward-facing normal on the interface. Thus by solving (3.5), we find that $F = v_n$ on $\boldsymbol {x} \in \partial \Omega (t)$ , recalling that $v_n$ is the normal speed of the interface. Further, (3.5) provides a continuous expression for F in the viscous fluid region. The derivatives in (3.5) are evaluated using central differencing. However, to solve (3.1) we require an expression for F over the entire computational domain. It was proposed by Moroney et al. [Reference Moroney, Lusmore, McCue and McElwain96] that the speed function be extended into the inviscid fluid region by solving the biharmonic equation

(3.6) $$ \begin{align} \nabla^4 F = 0, \quad \boldsymbol{x} \in \Omega(t). \end{align} $$

This ensures that $F = v_n$ is satisfied on the interface and gives a continuous expression for F away from the interface. For the purpose of discretisation, the sign of $\phi $ is used to determine nodes inside the interface that need to be included in the biharmonic stencil. This discretisation results in a symmetric system of linear equations that is solved using LU decomposition. As such, the location of the interface does not need to be known explicitly, similar to the level-set method itself. This velocity extension process is a variant of a thin plate spline in two dimensions [Reference Bookstein15]. To illustrate the velocity extension process, we consider an example F defined in the region $\Omega (t)$ , shown in Figure 3(a). The red line represents the fluid–fluid interface $\partial \Omega (t)$ . Figure 3(b) shows F after the biharmonic equation (3.6) is solved. We see that this velocity extension process gives a differentiable expression for F over the entire computational domain, which can be used to solve (3.1).

Figure 3 An illustration of the velocity extension process used to extend F to be defined over the entire computational domain. (a) shows an example function that is undefined where $\boldsymbol {x} \in \Omega $ . (b) shows this regions being “filled in’ by solving the biharmonic equation (3.6). This gives us a differentiable extension of F over the entire domain (colour available online).

3.3 Solving for pressure

To evaluate the speed function F, we must first compute the pressure field. We consider (2.8)–(2.11) in polar coordinates with $p = p(r, \theta , t)$ , and the location of the interface is given by $r = s(\theta , t)$ . Thus (2.8) becomes

(3.7) $$ \begin{align}\frac{1}{r} \frac{\partial}{\partial r} \bigg( r b^3 \frac{\partial p}{\partial r} \bigg) + \frac{1}{r^2} \frac{\partial}{\partial \theta} \bigg( b^3 \frac{\partial p}{\partial \theta} \bigg) = \frac{\partial b}{\partial t}, \quad r > s(\theta, t). \end{align} $$

For nodes away from the interface, the derivatives in (3.7) are discretised using a standard five-point stencil, illustrated in Figure 4(a). Denoting $\beta = r b^3$ , the r-derivatives in (3.7) are approximated via

(3.8) $$ \begin{align} \frac{1}{r} \frac{\partial}{\partial r} \bigg( \beta \frac{\partial p}{\partial r} \bigg) \to \frac{1}{r_{i,j} \Delta r} \bigg( \beta_{i+1/2,j} \frac{p_{i+1,j} - p_{i,j}}{\Delta r} - \beta_{i-1/2,j} \frac{p_{i,j} - p_{i-1,j}}{\Delta r} \bigg), \end{align} $$

where $\beta _{i+1/2,j} = (\beta _{i+1,j} + \beta _{i,j})/2$ and $\beta _{i-1/2,j} = (\beta _{i-1,j} + \beta _{i,j})/2$ . The derivatives in the $\theta $ -direction are discretised in a similar fashion.

Figure 4 An illustration of how the pressure of the viscous fluid is solved for using the finite-difference method. (a) For nodes away from the interface, we solve for pressure (3.7) using a standard five-point finite-difference stencil (3.8). (b) However, this stencil cannot be used for nodes adjacent to the interface, as in this case the southern node is not in the domain $\boldsymbol {x} \in \mathbb {R}^2 \backslash \Omega (t)$ , and thus cannot be used in our stencil. Instead we impose a ghost node, denoted by the red dots, on the interface whose location corresponds to the point where $\phi =0$ . The value at this ghost node is determined from the dynamic boundary condition (2.9). This leads to the nonstandard finite-difference stencil (3.9) (colour available online).

As illustrated in Figure 4(b), special care must be taken when solving for nodes adjacent to the interface. Suppose that the interface is located at $r = r_I$ with $r_{i-1,j} < r_I < r_{i,j}$ , where the nodes $r_{i-1,j}$ and $r_{i,j}$ are in the inviscid and viscous fluid regions, respectively. When discretising (3.7), we can no longer incorporate $p_{i-1,j}$ into our finite-difference stencil as it is not in the domain $\boldsymbol {x} \in \mathbb {R}^2 \backslash \Omega (t)$ . Instead, we define a ghost node at $r_I$ (denoted by the red dots in Figure 4(b)) whose value is $p_I$ . By noting that $\phi $ is a signed distance function, the distance between $r_{i,j}$ and $r_I$ is computed via

$$ \begin{align*} h = \Delta r \bigg| \frac{\phi_{i,j}}{\phi_{i-1,j} - \phi_{i,j}} \bigg|. \end{align*} $$

As per Chen et al. [Reference Chen, Merriman, Osher and Smereka25], our finite-difference stencil becomes

(3.9) $$ \begin{align} \frac{1}{r} \frac{\partial}{\partial r} \bigg( \beta \frac{\partial p}{\partial r} \bigg) & \approx \frac{2}{r_{i,j}(\Delta r + h)} \bigg( \beta_{i+1/2,j} \frac{p_{i+1,j} - p_{i,j}}{\Delta r} -\hat{\beta}_{i-1/2,j} \frac{p_{i,j}}{h} \bigg) \nonumber \\ &\quad+ \frac{2}{r_{i,j} h (\Delta r + h)} \hat{\beta}_{i-1/2,j} p_I. \end{align} $$

Here $\hat {\beta }_{i-1/2,j} = (\beta _{i,j} + \beta _I)/2$ where $\beta _I$ is the value of $\beta $ on the interface. When the node and interface are sufficiently close together such that $h < \Delta r^2$ , we set $p_{i,j} = p_I$ . A similar procedure is applied if the interface lies between $r_{i,j} < r_I < r_{i+1,j}$ and in the azimuthal direction.

The value of $p_I$ is computed from the dynamic boundary condition (2.9). To determine the value of $p_I$ , we first compute the curvature of $\phi $ via $\kappa = \boldsymbol {\nabla } \cdot \boldsymbol {n}$ over the entire computational domain, recalling that $\boldsymbol {n} = \boldsymbol {\nabla } \phi / |\boldsymbol {\nabla } \phi |$ . The value of $\kappa $ on the interface is computed via linear interpolation using $\kappa _{i,j}$ and $\kappa _{i-1,j}$ , leading to

$$ \begin{align*} p_I = -\sigma \bigg( \kappa_{i,j} - \frac{h(\kappa_{i,j}-\kappa_{i-1,j})}{\Delta r} + \frac{2}{b(r_I, \theta_j)} \bigg) - \omega^2 r_I^2. \end{align*} $$

For more information about solving both elliptic and parabolic problems in irregular domains using the finite-difference method in conjunction with level-set functions, we refer the reader to Coco and Russo [Reference Coco and Russo26] and Gibou et al. [Reference Gibou, Fedkiw, Cheng and Kang53] (and the references therein). Once the finite-difference stencil has been formed, the resulting system of linear equations is solved using the LU decomposition.

3.3.1 Far-field boundary condition

To incorporate the far-field boundary condition (2.11) 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(\theta , t)$ . By only considering the region $s(\theta , t) \le r \le R$ , we seek a solution to (2.8) of the form

(3.10) $$ \begin{align} \hat{p}(r, \theta, t) = A_0 - \frac{Q}{2 \pi} \log r + \frac{r^2}{4} \frac{\partial b}{\partial t} + \sum_{n=1}^{\infty} r^{-n} ( A_n \cos n \theta + B_n \sin n \theta ), \end{align} $$

where $A_0$ , $A_n$ , and $B_n$ are unknown, and $\hat {p} = b^3 p$ . The expansion (3.10) assumes that b is spatially uniform in $r \ge R$ so that $\hat {p}$ satisfies the appropriate Poisson equation. Considering the value of pressure on the artificial boundary, suppose that $\hat {p}(R, \theta , t)$ can be represented as a Fourier series

(3.11) $$ \begin{align} \hat{p}(R, \theta, t) = a_0 + \sum_{n = 1}^{\infty} a_n \cos n \theta + b_n \sin n \theta, \end{align} $$

where

$$ \begin{align*} a_{0} &= {\frac{1}{2\pi} }\int_{0}^{2\pi} \hat{p}(R, \theta,t) \, d \theta, \quad a_n = \frac{1}{\pi} \int_{0}^{2\pi} \hat{p}(R, \theta,t) \cos n \theta \, d\theta,\\[5pt] b_n &= \frac{1}{\pi} \int_{0}^{2\pi} \hat{p}(R, \theta,t) \sin n \theta \, d\theta. \end{align*} $$

By equating (3.11) with (3.10) evaluated at $r = R$ , we find that

$$ \begin{align*}A_0 = a_0 + (Q/2 \pi) \log R - \dot{b}R^2 / 4, \quad A_n = R^n a_n\quad \text{and} \quad B_n = R^n b_n.\end{align*} $$

To incorporate our expression for $\hat {p}$ into our finite-difference stencil, we differentiate (3.11) with respect to r and evaluate it at $r = R$ to give

$$ \begin{align*} \frac{\partial}{\partial r} \hat{p}(R, \theta_j) = -\frac{Q}{2 \pi R} + \frac{R}{2} \frac{\partial b}{\partial t} - \sum_{n = 1}^{\infty} \frac{n}{R} ( a_n \cos n \theta_j + b_n \sin n \theta_j ). \end{align*} $$

By approximating the integrals in our expressions for $a_n$ and $b_n$ as

$$ \begin{align*} a_n \approx \frac{\Delta \theta}{\pi} \sum_{k=1}^{m} \hat{p}(R, \theta_k) \cos n \theta_k \quad \textrm{and} \quad b_n \approx \frac{\Delta \theta}{\pi} \sum_{k=1}^{m} \hat{p}(R, \theta_k) \sin n \theta_k, \end{align*} $$

we have

$$ \begin{align*} \frac{\partial}{\partial r} \hat{p}(R, \theta_j)\approx -\frac{Q}{2 \pi R} + \frac{R}{2}\frac{\partial b}{\partial t} - \frac{\Delta \theta}{R \pi} \sum_{k = 1}^{m} v_{jk} \hat{p}(R, \theta_k), \end{align*} $$

where

$$ \begin{align*} v_{jk} = \sum_{n = 1}^{\infty} n \cos (n (\theta_k - \theta_j)). \end{align*} $$

Defining I as the outermost index at which $r = R$ , then our expression for $\partial p / \partial r$ is incorporated into our finite-difference stencil,

$$ \begin{align*} \frac{1}{r} \frac{\partial}{\partial r} \bigg( \beta \frac{\partial p}{\partial r} \bigg) \to \frac{2}{R \Delta r} \bigg[ -\beta_{I-1/2,j} \frac{ p_{I,j} - p_{I-1,j}}{\Delta r} + R \bigg\lbrace -\frac{Q}{2 \pi R} + \frac{R}{2} \frac{\partial b}{\partial t} - b^3 \frac{\Delta \theta}{R \pi} \sum_{k = 1}^{m} w_{jk} p_{I,k} \bigg \rbrace \bigg] , \end{align*} $$

recalling that $\beta = r b^3$ .

As an aside, we note this procedure can be adapted to model a Dirichlet boundary condition of the form $p \sim p_{\infty }$ as $r \to \infty $ where $p_\infty $ is prescribed. This type of boundary condition would be more appropriate for the model of a Stefan problem [Reference Morrow, King, Moroney and McCue98], where now p would represent temperature that is prescribed in the far field. Assuming that b is constant, then (3.10) becomes

$$ \begin{align*} p = p_\infty + A_0\ln r + \sum_{n=1}^{\infty} r^{-n} ( A_n \cos n \theta + B_n \sin n \theta ). \end{align*} $$

By following the same steps outlined above, we have

$$ \begin{align*} \frac{\partial}{\partial r} p(R, \theta_j) &= \frac{a_0 - p_\infty}{R \ln R} - \sum_{n = 1}^{\infty} \frac{n}{R} ( a_n \cos n \theta_j + b_n \sin n \theta_j ). \end{align*} $$

This expression is then incorporated into our finite-difference stencil as per usual.

3.4 General algorithm

We summarise our numerical algorithm as follows.

1. Step 1 For a given initial interface $s(\theta , 0)$ , initialise $\phi $ as a signed distanced function such that it satisfies (3.2) using the method of crossing times.

2. Step 2 Place marker particles around the interface, noting which side of the interface they are on, as well as their minimum distance from the interface.

3. Step 3 Solve for pressure in the domain $\boldsymbol {x} \in \mathbb {R}^2 \backslash \Omega $ using the procedure described in Section 3.3.

4. Step 4 Compute F according to (3.5), and then extend F into the region $\boldsymbol {x} \in \Omega $ by solving the biharmonic equation (3.6).

5. Step 5 Using F, update both $\phi $ and the marker particles by solving (3.1) and (3.4), respectively.

6. Step 6 Correct $\phi $ (if necessary) using the marker particles.

7. Step 7 Reinitialise the level-set function by solving (3.3), and then correct $\phi $ using the marker particles.

8. Step 8 Repeat steps 2–7 until $t = t_f$ .

4.1 Expanding bubble problem

We first consider the standard Hele-Shaw problem in which the inviscid bubble is injected into the viscous fluid, while the plates are parallel and stationary such that $b = 1$ and $\omega = 0$ .

4.1.1 Numerical validation

As a preliminary test for our scheme, we demonstrate that numerical solutions converge for a sufficiently refined grid. To do so, we perform simulations with the initial condition

(4.1) $$ \begin{align} s(\theta, 0) = 1 + \varepsilon \cos m \theta, \end{align} $$

where $\varepsilon = 0.1$ , $m = 6$ , on the computational domain $0 \le r \le 7.5$ and $0 \le \theta < 2 \pi $ , performed using an increasingly refined mesh. These simulations, shown in Figure 5, indicate that the interfacial profiles are converging as the mesh is refined, and that grid independence is achieved, at this scale, using $750 \times 628$ equally spaced nodes. Further, our solutions appear to maintain sixfold symmetry over the duration of the simulation.

Figure 5 Convergence test of numerical scheme for the evolution of a bubble with initial condition (4.1), where solutions are computed using (a) $350 \times 293$ , (b) $550 \times 461$ , (c) $750 \times 628$ , and (d) $850 \times 712$ equally spaced nodes. Additionally, $Q = 1$ , $\sigma = 5 \times 10^{-4}$ , $\omega = 0$ and $t_f = 100$ . Simulations are performed on the domain $0 \le r \le 7.5$ and $0 \le \theta < 2\pi $ . Solutions are plotted in time intervals of $t = 10$ .

We also demonstrate that our use of the Dirichlet-to-Neumann map, described in Section 3.3.1, results in the bubble’s volume changing at rate Q. To do so, we consider three different injection rates. The first is the constant injection rate $Q = 1$ , the second is the sinusoidal injection rate,

$$ \begin{align*} Q = 1 + 0.2 \sin (4\pi t/t_f), \end{align*} $$

and the third is the piecewise rate

$$ \begin{align*} Q = \begin{cases} 0.8 & \quad t \leq t_f/2, \\ 1.2 & \quad t> t_f/2. \end{cases} \end{align*} $$

We find that the rate of change of volume computed from the numerical simulations compares well with the corresponding exact rate of expansion, shown in Figure 6(a). The relative error, shown in Figure 6(b), suggests that over the duration of a simulation, we experience mass loss of only approximately $0.1\%$ . This result suggests that the Dirichlet-to-Neumann map correctly ensures that the bubble expands at the correct rate for both constant and time-dependent Q.

Figure 6 (a) The rate of change of volume, $\dot {V}$ , of numerical solution with constant (red), periodic (blue) and piecewise (yellow) injection rates. Dotted black lines denote the corresponding exact rate of change of volume. (b) Corresponding relative error. Simulations are performed on the domain $0 \le r \le 7.5$ and $0 \le \theta < 2 \pi $ using $750 \times 628$ equally spaces nodes with the initial condition (4.1). The surface tension parameter is $\sigma = 5 \times 10^{-4}$ and the final time of simulations is $t = 100$ (colour available online).

4.1.2 Effect of surface tension

Surface tension, modelled via the inclusion of the signed curvature term in the dynamic boundary condition (2.9), 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) despite initial conditions that are arbitrarily close; hence the zero-surface-tension problem is ill-posed. Adding surface tension removes the possibility of cusp formation; in Figure 2, we have included solutions for small but nonzero surface tension, computed using the method described in Section 3, for each of the cases described in Section 2.2. With nonzero surface tension, the interface remains smooth and solutions exist for all time, or, in the case of the finite blob (Figure 2)(c), until the interface intersects the point at which fluid is being withdrawn.

Linear stability analysis [Reference Paterson107] indicates that increasing the surface tension parameter $\sigma $ acts to make the interface less unstable, and nonlinear numerical simulations and experiments [Reference Chen24, Reference Dai and Shelley33, Reference Hou, Lowengrub and Shelley66] show that increasing the injection rate, which is mathematically equivalent to decreasing the dimensionless surface tension $\sigma $ , results in an increase in the number of fingers that develop. We perform numerical simulations for $Q=1$ with values of surface tension varying over several orders of magnitude with the initial condition

(4.2) $$ \begin{align} s(\theta, 0) = 1 + \varepsilon ( \cos m \theta + \sin n \theta ), \end{align} $$

where $\varepsilon = 0.1$ , $m = 3$ and $n = 2$ . These numerical simulations, shown in Figure 7, are able to reproduce the key morphological features of the Saffman–Taylor instability. For each value of $\sigma $ considered, the interface is unstable, and the sinusoidal perturbations in the initial condition (4.2) grow and evolve into viscous fingers. As $\sigma $ is decreased, the interface becomes more unstable and the number of fingers that develop over the duration of a simulation increases due to tip-splitting occurring (see $(1)$ in Figure 7(b), for example). Additionally, our simulations are able to produce so called “shielding” behaviour, where neighbouring fingers can block one another off, which in turn results in fingers retracting (denoted by $(2)$ in Figure 7(c)). This behaviour is known to occur experimentally (see [Reference Paterson107] 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 (denoted by $(3)$ in Figure 7(d)). Again this behaviour has been observed in the experiments (see [Reference Chen24]).

Figure 7 The numerical solution to (2.8)–(2.11) with $Q = 1$ for values of surface tension parameter $\sigma $ equal to (a) $1.5 \times 10^{-3}$ , (b) $5 \times 10^{-4}$ , (c) $1.75 \times 10^{-4}$ and (d) $6.25 \times 10^{-5}$ with initial condition (4.2). Our numerical scheme 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 $0 \le r \le 10$ and $0 \le \theta < 2 \pi $ using $1000 \times 628$ equally spaced nodes. Solutions are plotted in time intervals of $t = 5$ up to $t_f = 55$ .

4.2 Tapered plates

One of the more popular modifications to the Hele-Shaw cell (particularly in recent years) is to consider that the plates of the Hele-Shaw cell are no longer parallel but instead linearly tapered (either converging or diverging) in the direction of the flow. The concept was first introduced by Zhao et al. [Reference Zhao, Casademunt, Yeung and Maher140], who considered tapered plates in a channel geometry, discussed further in Section 4.4. Imposing tapered plates in radial geometry has been studied using linear stability analysis [Reference Al-Housseiny, Christov and Stone3], weakly-nonlinear stability analysis [Reference Anjos, Alvarez, Dias and Miranda7], experimentally [Reference Al-Housseiny, Christov and Stone3, Reference Bongrand and Tsai14] and by numerical simulations [Reference Jackson, Stevens, Power and Giddings73]. Results indicate that tapering the gap between the plates such that they either converge or diverge can have a stabilising or destabilising effect on the interface, depending on the injection rate. The numerical scheme presented in Section 3 was used by us in Morrow et al. [Reference Morrow, Moroney and McCue99] to study how tapering the plates while injecting at either a constant or time-dependent injection rate can be used to reduce the development of viscous fingering patterns.

Here we perform numerical simulations where the gap between the plates is linearly tapered according to

(4.3) $$ \begin{align} b(r) = \begin{cases} b_\infty + \alpha (r - r_B) & \text{if } r \le r_B, \\ b_\infty & \text{if } r> r_B, \end{cases} \end{align} $$

where $\alpha $ is the gradient of the taper ( $\alpha =0$ being the unmodified Hele-Shaw cell). Experiments were recently performed by Bongrand and Tsai [Reference Bongrand and Tsai14], who considered a Hele-Shaw cell, where the plate gap is (4.3) with $\alpha> 0$ 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. In our nondimensionalisation, for which the time-scale is set by the injection rate, a smaller (dimensional) injection rate corresponds to a larger (dimensionless) surface tension value $\sigma $ . To confirm that our numerical scheme produces solutions consistent with these experiments, we perform simulations with two different values of $\sigma $ , shown in Figure 8, to model two different (dimensional) injection rates. The initial condition of these simulations is (4.2), where $0 \le \theta < 2 \pi $ , $\varepsilon = 5 \times 10^{-3}$ , $m = 5$ and $n = 4$ . For a larger $\sigma $ (Figure 8(a)), we indeed find that the interface is stabilised, and remains circular over the duration of the simulation. For a smaller $\sigma $ (Figure 8(b)), we find that the interface is unstable; in particular, these fingers appear distinct from traditional Saffman–Taylor fingers (see Figure 7, for example). This behaviour is consistent with the results of Bongrand and Tsai [Reference Bongrand and Tsai14], who described the interface as becoming “wavy” as it expanded. These results suggest that our numerical simulations are producing the correct behaviour when the plates are of the form of (4.3).

Figure 8 Numerical simulations where the gap between the plates linearly tapered according to (4.3) with $Q = 1$ , $b_\infty = 1/15$ , $R_0 = 8/3$ , $r_B = 7$ and $\alpha = 2/15$ , with surface tension $\sigma $ equal to (a) 6 and (b) 1. Simulations are performed on the domain $0 \le r \le 7.5$ and $0 \le \theta < 2 \pi $ using $750 \times 628$ equally spaced nodes. Solutions are plotted in time intervals of 5.6. up to $t_f = 44.8$ .

4.3 Blob problem

In this subsection we now consider the complementary problem to (2.8)–(2.11), for which a blob of viscous fluid now occupies the region $\Omega (t)$ and the inviscid fluid is in $\mathbb {R}^2 \backslash \Omega (t)$ . This problem is typically studied by considering the withdrawal of the viscous fluid, which in turn causes fingers to develop inward (as in Figure 2(c), for example). However, popular modifications, including considering the gap between the plates to be time-dependent or rotating the entire Hele-Shaw cell, have also received interest. In this subsection we explain how the numerical scheme presented in Section 3 is modified to solve for these variations, and show that our simulations are consistent with previous experimental and numerical results.

For the case in which the inviscid bubble is injected into the viscous fluid (Sections 4.1–4.2, say), the velocity of the fluid is driven by a sink term in the far field given by (2.11). For the blob problem considered here, the withdrawal of the viscous fluid is incorporated into the model via a sink at the origin. To include this sink into our numerical model, we follow Hou et al. [Reference Hou, Lowengrub and Shelley66] and introduce a smoothed Dirac delta function in (2.8) such that

(4.4) $$ \begin{align} \boldsymbol{\nabla} \cdot ( b^3 \boldsymbol{\nabla} p ) = \frac{\partial b}{\partial t} + S, \end{align} $$

where

(4.5) $$ \begin{align} S = \begin{cases} \dfrac{Q}{b \bar{r}_0^2} \bigg( 1 + \cos \dfrac{\pi r}{\bar{r}_0} \bigg) & \text{if } r \leq \bar{r}_0, \\ 0 & \text{if } r> \bar{r}_0. \end{cases} \end{align} $$

As shown by Hou et al. [Reference Hou, Lowengrub and Shelley66], this choice of source term ensures that the rate of change of volume of the viscous fluid is Q. For all results in this subsection, we use $\bar {r}_0 = 0.05$ . We note that it is straightforward to extend (4.5) to include multiple sink/source points.

When solving for the governing equation for pressure (2.8) for the case where an inviscid fluid is injected into a viscous one (Sections 4.1-4.2), we consider the model in polar coordinates as it is simpler to incorporate the far-field boundary condition (2.11) via the Dirichlet-to-Neuman map (Section 3.3.1) on a circle. However, when the inviscid region surrounds the viscous fluid, we no longer have a far-field boundary condition and, as such, it is more convenient to solve for pressure in Cartesian coordinates, although of course either coordinate system could be used. The discretisation of the spatial derivatives (4.4) is performed as was described in Section 3.3. That is, we use central finite differences for nodes away from the interfaces, and incorporating a ghost value of p for nodes adjacent to the interface. We refer to the papers by Chen et al. [Reference Chen, Merriman, Osher and Smereka25] and Gibou et al. [Reference Gibou, Fedkiw, Cheng and Kang53, Reference Gibou, Fedkiw and Osher54] for more details on implementing this stencil in Cartesian coordinates.

In a similar way to that described in Section 3.2, once we have computed F via (3.5), we extend it into the region $\boldsymbol {x} \in \mathbb R^2 \backslash \Omega (t)$ by solving the biharmonic equation (3.6). When solving (2.8)–(2.11), the speed function was known outside the interface and was 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 forming our biharmonic stencil, we apply homogeneous Neumann boundary conditions on each of the boundaries. We illustrate this process in Figure 9, which shows an example F before (Figure 9(a)) and after (Figure 9(b)) the velocity extension process. Again we see that by solving the biharmonic equation we obtain a smooth expression for F over the entire computational domain.

Figure 9 An illustration of the velocity extension process where the viscous fluid is surrounded by the inviscid bubble. As discussed in Section 3.2, this is done by solving the biharmonic equation in the region where $\boldsymbol {x} \in \Omega (t)$ , where now we apply homogeneous Neumann boundary conditions on the edge of the computational domain (colour available online).

4.3.1 Withdrawal of viscous fluid

We consider the case where the viscous fluid is withdrawn such that as the interface contracts, viscous fingers form inward. As described in Section 2.2, exact solutions in the absence of surface tension ( $\sigma = 0$ ) are ill-posed, 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 effects of surface tension preventing these cusps from forming [Reference Ceniceros, Hou and Si20, Reference Kelly and Hinch75]. Experimental results [Reference Paterson107, Reference Thomé, Rabaud, Hakim and Couder130] show that the fingers that form exhibit morphological features distinct from traditional Saffman–Taylor fingers, in that these fingers do not appear to undergo tip-splitting but instead appear to be in competition to “race” towards the sink.

Figure 10 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.1, 0)$ where $\sigma = 8 \times 10^{-4}$ and $t_f = 1.88$ . For (b), initial condition is (4.6) where $\sigma = 1.6 \times 10^{-6}$ and $t_f = 1.45$ . A black dot denotes the region where $r \le 0.05$ . Simulations are performed on the domain $-1.15 \le x \le 1.15$ and $-1.15 \le y \le 1.15$ using $400 \times 400$ equally spaced nodes.

As well as the comparison with the zero-surface-tension solution made in Figure 2(c), we perform two further different simulations for this configuration. The first is with an initially circular interface centred at $(0,-0.1)$ shown in Figure 10(a). This simulation shows that as the interface contracts, it becomes nonconvex until a single finger develops that tends towards the origin. This behaviour compares well with previous numerical simulations performed using the boundary integral method [Reference Kelly and Hinch75]. For the second simulation, shown in Figure 10(b), we consider a perturbed circle centred at the origin of the form

(4.6) $$ \begin{align} s(\theta, 0) = 1 + \varepsilon ( \cos 3 \theta + \sin 7 \theta + \cos 15 \theta + \sin 25 \theta ), \end{align} $$

where $\varepsilon = 5 \times 10^{-3}$ . We find that the interface initially develops numerous short fingers. These fingers do not appear to exhibit the same morphological features as the case in which the inviscid bubble is injected such as tip-splitting and feathering (see Figure 7, for example), but instead the number of fingers remains constant. Due to the pressure differential between the sink and the boundary of the bubble, the velocity of one of the fingers rapidly increases, and the simulation is stopped when this finger reaches the origin. We note that this behaviour compares well with experimental results (see [Reference Thomé, Rabaud, Hakim and Couder130] for example), as well as numerical simulations in [Reference Chen, Huang and Miranda23].

4.3.2 Lifting plates

A popular modification to the blob problem is to consider the case where the upper plate is uniformly lifted in time such that $b \to b(t)$ . The volume of viscous fluid remains constant ( $Q = 0$ ) so that when the plates are separated, viscous fingers developing inward. For this problem, the Hele-Shaw approximation itself can only remain valid for as long as the gap $b(t)$ is sufficiently small. For example, in dimensional terms, we must be very careful about using the model when the gap width is of the same order as the important length scales in the lateral direction.

For this lifting-plates configuration, the governing equation for pressure (4.4) becomes Poisson’s equation. However, pressure can be reduced via $P = p - \dot {b} |\boldsymbol {x}|^2/(4 b^3)$ , essentially moving the nonhomogeneous term to the dynamic boundary condition (2.9). This transformation reduces (2.8) to Laplace’s equation, meaning this configuration can also be solved numerically with the boundary integral method [Reference Zhao, Li, Ying, Belmonte, Lowengrub and Li141, Reference Zhao, Niroobakhsh, Lowengrub and Li142].

The lifting-plates problem was first considered mathematically by Shelley et al. [Reference Shelley, Tian and Wlodarski122]. In the absence of surface tension ( $\sigma = 0$ ), they [Reference Shelley, Tian and Wlodarski122] argued that the generic behaviour is that a cusp will develop in a finite time. When surface tension is included, numerical simulations suggested a relationship between the number of fingers that develop and the surface tension parameter. In particular, it was shown that the number of fingers is a monotonically decreasing function of time. This behaviour was later shown to be consistent with experimental results [Reference Lindner, Derks and Shelley83, Reference Nase, Derks and Lindner100]. As a point of comparison, for the traditional Hele-Shaw configuration, discussed in Section 4.1, the number of fingers typically increases with time due to tip-splitting. We note that the case in which the inviscid bubble is injected into the viscous fluid while the plates are separated has also received attention [Reference Morrow, Moroney and McCue99, Reference Vaquero-Stainer, Heil, Juel and Pihler-Puzović133, Reference Zheng, Kim and Stone143].

We perform a simulation using the linearly increasing gap between the plates $b = 1 + t$ for different values of $\sigma $ . Note that this time-dependence with $\dot {b}=1$ effectively chooses the appropriate time-scale T in (2.7). When $\sigma = 10^{-4}$ (top row of Figure 11), we find that the interface quickly destabilises and approximately 15 fingers develop by $t = 1$ . As time increases further, neighbouring fingers begin to merge with each other and the overall number of fingers significantly decreases, and we find only around five fingers remain by the conclusion of the simulation. For the lower value of surface tension $\sigma = 5 \times 10^{-5}$ (bottom row of Figure 11), we find the number of fingers that initially develop has increased compared to when $\sigma = 10^{-4}$ , about 25 at $t = 1$ . However as the interface contracts, again fingers begin to merge and only eight fingers remain at $t = 4$ . We perform these simulations over a longer time period (not shown), which reveals that the interface will become circular when the gap between the plates is sufficiently large. This behaviour is consistent both with previous experimental and numerical results [Reference Lindner, Derks and Shelley83, Reference Nase, Derks and Lindner100, Reference Shelley, Tian and Wlodarski122].

Figure 11 Time evolution of viscous fluid where plates are separated according to $b = 1 + t$ with initial condition (4.6). Solutions are plotted at times (left to right) $t = 0$ , 1, 1.8, 2.8, and 4. Simulations are performed on the domain $-1.1 \le x \le 1.1$ and $-1.1 \le y \le 1.1$ using $400 \times 400$ equally spaced nodes.

4.3.3 Rotating plates

While Saffman–Taylor fingers traditionally form due to the injection/withdrawal of one immiscible fluid into another, it is known that these fingers can also be triggered by body forces. The two most commonly studied body forces are gravity and centrifugal forces. For the latter, when the entire Hele-Shaw cell is rotated, this rotation results in the dense viscous fluid being propelled outward, which in turn leads to finger formation [Reference Schwartz119]. Experimental [Reference Carrillo, Magdaleno, Casademunt and Ortín17] and numerical [Reference Alvarez-Lacalle, Gadêlha and Miranda4, Reference Folch, Alvarez-Lacalle, Or´tın and Casademunt46, Reference Paiva, Lira and Andrade105] studies reveal that the interface patterns are distinct from the traditional Saffman–Taylor instability. That is, fingers appear more “stretched-out” and generally do not undergo tip-splitting. We note that our model (2.8)–(4.4) ignores the effect of Coriolis forces, however several studies have investigated its effect on interfacial dynamics [Reference Alvarez-Lacalle, Gadêlha and Miranda4, Reference Schwartz119, Reference Waters and Cummings137]. Further, we consider the case where the inviscid bubble is injected into the viscous fluid while the plates are rotated in Morrow et al. [Reference Morrow, Moroney and McCue99], where we show that the angular velocity has a stabilising effect on the interface.

The incorporation of the centrifugal term is straightforward. While body forces appear in the governing equation for the velocity of the viscous fluid (2.1), scaling pressure means that the angular velocity term can be moved to the dynamic boundary condition (2.9) (this can be done for any conservative body force $\boldsymbol {f}$ satisfying $\boldsymbol {\nabla } \times \boldsymbol {f} = \boldsymbol {0}$ ). The rotating Hele-Shaw cell has previously been studied numerically using boundary integral [Reference Schwartz119] and diffusive interface [Reference Chen, Huang and Miranda23, Reference Paiva, Lira and Andrade105] techniques. We perform simulations where the volume of viscous fluid between the plates is constant ( $Q = 0$ ) and the Hele-Shaw plates are rotated with $\omega =1$ for different values of $\sigma $ , shown in Figure 12. Note that this choice of $\omega $ effectively fixes the appropriate time-scale T in (2.7).

Figure 12 Numerical simulation of a rotating Hele-Shaw cell with $\omega = 1$ , $Q = 0$ , and $\sigma $ equal to (a) $ 10^{-2}$ , (b) $5 \times 10^{-3}$ , (c) $2.5 \times 10^{-3}$ , and (d) $10^{-3}$ . Corresponding final time of simulations is $t = 0.61$ , $0.55$ , $0.425$ , and $0.36$ . Initial condition for each simulation is (4.6). Simulations are performed on the domain $-2 \le x \le 2$ and $-2 \le y \le 2$ using $500 \times 500$ equally spaced nodes.

For each value of $\sigma $ considered in Figure 12, we find that the interface is unstable, and the fingers that develop are distinct from traditional Saffman–Taylor fingers. In particular, the fingers that develop do not appear to tip-split but instead remain constant. Additionally, the number of fingers that develop increases as the surface tension parameter is decreased such that when $\sigma = 10^{-2}$ (Figure 12(a)), seven fingers form, while for $\sigma = 10^{-3}$ (Figure 12(d)), the number of fingers is approximately 21. These results are consistent with experimental results [Reference Carrillo, Magdaleno, Casademunt and Ortín17] and numerical simulations [Reference Alvarez-Lacalle, Gadêlha and Miranda4, Reference Paiva, Lira and Andrade105].

4.4 Channel geometry

In Sections 4.1–4.3, we considered Hele-Shaw flow in radial geometry, where the bubble–fluid interface is completely immersed by an infinite amount of viscous or inviscid fluid. In this subsection we focus on another well-studied version of the Hele-Shaw cell in channel (or rectangular) geometry, where the shape of the cell is a narrow rectangle of infinite length and width L. As discussed in Section 2.2, for the zero-surface-tension case, exact solutions are known to exist, which may involve a type of blow-up in finite time with a cusp forming on the boundary (as in Figure 2(d)).

This channel problem dates back to the work of Saffman and Taylor [Reference Saffman and Taylor118], who showed that when an inviscid bubble is injected into the channel filled with a viscous fluid, typically a single finger develops that propagates through the channel (see the numerical solution in Figure 2(e)). Since the work of Saffman and Taylor, extensive research has been carried out determining how the parameters of the model influence the width (relative to L) and speed of this finger. In particular, for a fixed injection rate, as the surface tension parameter is increased, it is established that the width and speed of the finger increase and decrease, respectively. We refer to the papers by Homsy [Reference Homsy63] and Saffman [Reference Saffman117] (and references therein) for a comprehensive overview of the problem. In this subsection (and Section 2.2) we restrict ourselves to the classic configuration where b is constant, but we note that, similar to the radial problem, our scheme can easily be used to study nonstandard cases, such as those in the papers [Reference Al-Housseiny, Christov and Stone3, Reference Franco-Gómez, Thompson, Hazel and Juel49, Reference Thompson, Juel and Hazel131, Reference Zhao, Casademunt, Yeung and Maher140].

As with the blob problem discussed in Section 4.3, it is more convenient to consider this problem in Cartesian coordinates such that $p \to p(x, y, t)$ and the interface is given by $x = f(y, t)$ . Similar to the radial case (see (2.15), for example), the velocity of the fluid is driven by the sink term in the far field

(4.7) $$ \begin{align} b^3 \frac{\partial p}{\partial x} \sim \frac{Q}{L} \quad \textrm{as } x \to \pm \infty. \end{align} $$

Equation (4.7) is incorporated into our finite-difference stencil using a Dirichlet-to- Neumann map. We do not provide full details, but note that the procedure is similar to that described in Section 3.3.1, where we impose an artificial boundary at $x = X$ , and seek a solution to (2.8) of the form

$$ \begin{align*} \hat{p}(x, y, t) = \frac{Q}{L}x + \sum_{n=0}^{\infty} A_n \textrm{e}^{-\lambda_n y} \cos \lambda_n x, \end{align*} $$

where $\lambda _n = 2 \pi n / L$ and $A_n$ is to be determined. Additionally, we impose $\partial p / \partial y = 0$ on $y = \pm L/2$ .

We perform a series of simulations using the same parameters as those of DeGregoria and Schwartz [Reference DeGregoria and Schwartz37], who studied this problem using a boundary-integral approach, to demonstrate that our solutions are consistent with the expected behaviour. We choose the initial condition

$$ \begin{align*} f(y, 0) = \varepsilon \cos 2 \pi y, \end{align*} $$

where $\varepsilon = 0.05$ , and perform simulations over a range of values of $\sigma $ , shown in Figure 13. In this figure $L=1$ , which sets the length scale $r_0$ in (2.7), while $Q=1$ fixes T in these equations. For low $\sigma $ (first row of Figure 13), we find that as the bubble expands, a finger grows, which is unstable and split into two. This is consistent with the results of DeGregoria and Schwartz [Reference DeGregoria and Schwartz37], and this behaviour is also observed experimentally by Tabeling et al. [Reference Tabeling, Zocchi and Libchaber124] when the injection rate is sufficiently large. For larger values of $\sigma $ (second to fourth rows of Figure 13), a single stable finger propagates through the channel whose speed and width decrease and increase as $\sigma $ increases. Again, this behaviour is consistent with previous experimental and numerical results [Reference DeGregoria and Schwartz37, Reference Tabeling, Zocchi and Libchaber124].

Figure 13 Numerical simulations in channel geometry with values of the surface tension parameter (top to bottom) $\sigma = 2 \times 10^{-4}$ , $5 \times 10^{-4}$ , $7 \times 10^{-3}$ , and $2 \times 10^{-2}$ . The initial condition for all simulations is $f(x,0) = 0.05 + \cos (2 \pi y)$ . Simulations are performed on the domain $-0.5 \le x \le 0.5$ and $-0.5 \le y \le 6$ using $100 \times 650$ equally spaced nodes. Solutions are plotted in time intervals of $t = 0.5$ up to $t_f = 3$ .