Self-similar and disordered front propagation in a radial Hele-Shaw channel with time-varying cell depth

TL;DR

This study investigates how time-varying cell depth in a Hele-Shaw cell affects viscous fingering, revealing delayed instability onset and the emergence of stable, self-similar interface patterns under specific perturbations.

Contribution

It provides a comprehensive linear stability analysis and numerical simulations showing the effects of dynamic cell depth on pattern formation and stability in viscous fingering.

Findings

Plate separation delays instability onset.

Self-similar interface evolution occurs under unimodal perturbations.

Complex finger patterns develop with non-unimodal perturbations.

Abstract

The displacement of a viscous fluid by an air bubble in the narrow gap between two parallel plates can readily drive complex interfacial pattern formation known as viscous fingering. We focus on a modified system suggested recently by [1], in which the onset of the fingering instability is delayed by introducing a time-dependent (power-law) plate separation. We perform a complete linear stability analysis of a depth-averaged theoretical model to show that the plate separation delays the onset of non-axisymmetric instabilities, in qualitative agreement with the predictions obtained from a simplified analysis by [1]. We then employ direct numerical simulations to show that in the parameter regime where the axisymmetrically expanding air bubble is unstable to nonaxisymmetric perturbations, the interface can evolve in a self-similar fashion such that the interface shape at a given time is simply a rescaled version of the shape at an earlier time. These novel, self-similar solutions are linearly stable but they only develop if the initially circular interface is subjected to unimodal perturbations. Conversely, the application of non-unimodal perturbations (e.g. via the superposition of multiple linearly unstable modes) leads to the development of complex, constantly evolving finger patterns similar to those that are typically observed in constant-width Hele-Shaw cells. [1] Z. Zheng, H. Kim, and H. A. Stone, Controlling viscous fingering using time-dependent strategies, Phys. Rev. Lett. 115, 174501 (2015).

Equations

$${\mathbf{u}}^{\ast }=-\frac{{b^{\ast }}^2}{12\mu }\nabla p^{\ast },\frac{\text{d}{b}^{\ast }}{\text{d}{t}^{\ast }}=-\nabla \cdot (b^{\ast }{\mathbf{u}}^{\ast })=\frac{{b^{\ast }}^3}{12\mu }{\nabla }^2{p}^{\ast }\text{ in }{r}^{\ast }>R^{\ast }(\theta ,t^{\ast }),$$

$$\frac{\partial {\mathbf{R}}^{\ast }}{\partial t^{\ast }}\cdot \mathbf{n}$$

$$p^{\ast }(R^{\ast },\theta ,t^{\ast })$$

$$Q{t}^{\ast }=\frac{b^{\ast }(t^{\ast })}{2}{\int }_0^{2\pi }{R^{\ast }}^2(\theta ,t^{\ast })d\theta .$$

$$\bar{R}(t)={\left(\frac{2t}{b(t)}\right)}^{1/2}$$

$$\bar{p}(r,t)=\frac{1}{b^3}\frac{\text{d}b}{\text{d}t}\frac{r^2-1}{4}-\frac{1}{b^3}\log r,$$

$$p_g(t)=\frac{1}{b^3}\frac{\text{d}b}{\text{d}t}\frac{{\bar{R}}^2-1}{4}-\frac{1}{b^3}\log \bar{R}+\frac{\mathcal{A}}{12Ca}\left(\frac{\mathcal{A}}{\bar{R}}+\frac{2}{b}\right),$$

$${\lambda }_k=\frac{1}{{\hat{R}}_k}\frac{\text{d}{\hat{R}}_k}{\text{d}t}=\frac{1+{\bar{R}}^{2k}}{1-{\bar{R}}^{2k}}k\left(\frac{{\mathcal{A}}^2{b}^2(1-k^2)}{12Ca{\bar{R}}^3}-\frac{1}{b}\left(\frac{1}{2}\frac{\text{d}b}{\text{d}t}-\frac{1}{{\bar{R}}^2}\right)\right)-\frac{1}{b}\left(\frac{1}{2}\frac{\text{d}b}{\text{d}t}+\frac{1}{{\bar{R}}^2}\right).$$

$${\lambda }_k=\frac{k-1}{{\bar{R}}^2}\left(\frac{1}{b}-\frac{{\mathcal{A}}^2{b}^{2}k(k+1)}{12Ca\bar{R}}-\frac{{\bar{R}}^2}{2b}\frac{\text{d}b}{\text{d}t}\frac{k+1}{k-1}\right).$$

$$\frac{\text{d}b}{\text{d}t}-\frac{(k-1)}{(k+1)}\frac{b}{t}+\frac{{\mathcal{A}}^{2}k(k-1)}{12\sqrt{2}Ca}\frac{b^{9/2}}{t^{3/2}}>0.$$

$$S_k=\frac{1}{7}-\frac{k-1}{k+1}+\frac{k(k-1)}{J}>0,$$

$$J=\frac{12\sqrt{2}Ca}{{\mathcal{A}}^2{b}_1^{7/2}}=\frac{6\mu Q^{3/2}}{{\pi }^{3/2}\gamma b_1^{\ast 7/2}}$$

$$J<J_{\mathrm{cr}}(k)=\frac{7k(k^2-1)}{6k-8},$$

$$J=\frac{7}{6}\left(3k_{\max }^2-1\right).$$

$$R(\theta ,t=t_0)=R_{\mathrm{init}}+ϵ\cos (k_{\max }\theta ),$$

$$R(\theta ,t=t_0)=R_{\mathrm{init}}+ϵ\sum _{\mathrm{unstable}k}{\hat{R}}_k\cos (k\theta +{\varphi }_k),$$

$$c(k)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }(R-⟨R⟩)\exp (-ik\theta )\text{d}\theta .$$

Full text

I Introduction

An expanding air bubble that displaces a viscous fluid within the narrow gap between two parallel plates (a Hele-Shaw cell) is unstable to non-axisymmetric perturbations beyond a critical value of the capillary number $Ca$, the ratio of viscous to surface tension forces [2, 3, 4, 5]. In this radial geometry, the unstable interface typically deforms into a set of growing fingers which evolve continuously in time through a sequence of tip-splitting events followed by competition between the newly-formed fingers. This non-linear evolution is an archetype for front-propagating, pattern forming phenomena [6].

Long-standing interest in viscous fingering stems from its similarities with a range of other front propagation phenomena, such as the growth of bacterial colonies [7] and the solidification instabilities during crystal growth [8]. In hydraulic fracture used for oil recovery, viscous fingering is promoted to ensure heterogeneous placement of particles (proppants) into fractures, thereby increasing their hydraulic conductivity, but it is also actively suppressed to avoid early breakthrough of water into adjacent production wellbores [9, 10]. This has stimulated a recent resurgence in research effort to delay the onset of viscous fingering, e.g. by manipulating the geometry of the cell, either actively [11, 12], or passively by using compliant cells [13, 14, 15]; by controlling the injection rate [16]; or by tuning the viscosity ratio of miscible fluid pairs [17].

In this paper, we consider fingering in a radial Hele-Shaw cell in which the width of the gap between the parallel bounding plates increases as a function of time. This effect is stabilising since it reduces the rate at which the axisymmetric bubble expands; in fact, if done sufficiently rapidly, the plate separation can bring the expansion of the bubble to a halt or even cause the bubble to contract. In a recent study, Zheng et al. [1] employed a simplified analysis, based on Paterson’s classical results for the system with constant gap width [3], to examine the case when the gap width increases according to the power law $b^{*}(t^{*})=b_{1}^{*}\ {t^{*}}^{1/7}$, where $b^{*}$ is the distance between the plates and $t^{*}$ is time. They not only confirmed that viscous fingering can be suppressed if the plates are separated sufficiently rapidly (i.e. for sufficiently large values of the constant $b_{1}^{*}$), but also predicted that the wavenumber $k_{\rm max}$ of the most rapidly-growing small-amplitude perturbation to the axisymmetric bubble is independent of the bubble’s radius. Furthermore, they performed experiments which showed that in the parameter range for which the axisymmetric bubble is unstable, a small number of finite-amplitude non-splitting fingers tended to develop and that the number of these fingers was remarkably close to $k_{\rm max}$.

The observation of non-splitting fingers is unusual because it contrasts with the typical sequence of tip-splitting and finger competition observed in experiments without plate lifting [4, 5]. The findings of Zheng et al. [1] suggest that, following some initial transients, the finger growth may become self-similar in the sense that the finger shape at a given time is simply a rescaled version of the shape at an earlier time. This is analogous to proportionate growth observed in biological systems, e.g. during growth of a mammal whose body parts grow at nearly the same rate and thus in direct proportion to each other [18]. In radial Hele-Shaw cells without lifting, proportionate growth of fingering (or self-similar fingering) has been observed in some experiments with miscible fluids [17]. Self-similar growth of the interface between two immiscible fluids has previously been observed [19, 20] and predicted [21, 22, 23] in a related setup where the cell consists of a disk sector. This geometry promotes the formation of a single finger which is symmetric about the sector centreline and which does not split for sufficiently small values of $Ca$ or sector angles. Other, more complex but well-defined, reproducible finger shapes were also obtained by imposing selected initial perturbations. Moreover, Li et al. [24] have computed self-similarly evolving fingering patterns in a radial geometry by appropriately varying the air injection rate, although the number of fingers observed in that study differed from the most unstable wavenumber predicted by linear stability analysis.

In this paper we revisit the scenario studied by Zheng et al. [1] and perform a full linear stability analysis to confirm that $k_{\rm max}$ is indeed independent of the bubble radius. We then perform numerical simulations of the system’s nonlinear evolution following the onset of the linear instability. We demonstrate that self-similarly evolving fingering patterns can be realised with the lifting law introduced by Zheng et al. [1], but only when the initially circular interface is perturbed with a single, linearly unstable mode. Unimodal perturbations with wavenumber $k_{\rm max}$ lead to the development of $k_{\rm max}$ distinct non-splitting, self-similar finite-amplitude fingers. If the interface is subjected to unimodal perturbations with other wavenumbers the perturbed interface typically undergoes a transient phase during which fingers split before they approach a self-similar regime. Therefore, the number of self-similar fingers that emerge from the instability does not necessarily coincide with the most unstable wavenumber predicted by the linear stability analysis. We also find that if non-unimodal, random initial perturbations are introduced, such as those occurring in a typical experiment, the interface does not evolve towards a self-similar solution, in contrast with the scenario suggested by the experiments of Zheng et al. [1]. Instead, the fingering pattern evolves continuously as the interface advances across the cell through a succession of tip-splitting and finger competition events similarly to the patterns in a radial cell with fixed parallel plates.

II Mathematical model

We describe the motion of the wetting viscous fluid occupying the narrow gap between the two bounding plates using the lubrication approximation. Using horizontal polar coordinates ($r^{*},\theta$), centred at the injection point (see figure 1), the equations for the vertically averaged velocity $\boldsymbol{u}^{*}(r^{*},\theta,t^{*})$ and the fluid pressure $p^{*}$ are then

[TABLE]

where $\mu$ is the liquid viscosity, $R^{*}(\theta,t^{*})$ is the bubble radius, and $b^{*}(t^{*})$ is the (spatially constant) gap width. Throughout this paper we employ asterisks to distinguish dimensional quantities from their non-dimensional equivalents. In the region $r^{*}<R^{*}(\theta,t^{*})$, the gas bubble has a uniform pressure $p_{g}^{*}(t^{*})$. We neglect the effects of the thin films of liquid that are left behind the advancing bubble tip, and ignore viscous normal stresses at the interface. The kinematic and dynamic boundary conditions at $r^{*}=R^{*}(\theta,t^{*})$ then become

[TABLE]

where $\boldsymbol{n}$ is the unit normal to the interface and $\gamma$ is the surface tension. The mean curvature of the interface is approximated as the sum of the in-plane and transverse curvatures, $\kappa^{*}_{\parallel}$ and $\kappa^{*}_{\perp}$, respectively. At the outer boundary of the Hele-Shaw cell we impose $p^{*}(R^{*}_{\rm outer},\theta,t^{*})=0$. Finally, for a constant injection flux $Q$, the volume of gas in the bubble is given by

[TABLE]

We study the system’s evolution starting from $t^{*}=t_{0}^{*}$ when the cell walls are separated by $b^{*}_{0}$ and the bubble has an initial radius $R^{*}=R^{*}_{\rm init}$.

In the following analysis, we non-dimensionalize in-plane lengths with $R^{*}_{\rm outer}$, the gap width with $b^{*}_{0}$, time with $2\pi R^{*2}_{\rm outer}b^{*}_{0}/Q$ and pressures with $6\mu Q/\pi b^{*3}_{0}$. Then the problem is governed by three non-dimensional parameters: the capillary number $Ca=\mu Q/(2\pi\gamma R^{*}_{\rm outer}b^{*}_{0})$, the cell aspect ratio $\mathcal{A}=b^{*}_{0}/R^{*}_{\rm outer}$ and the initial radius of the bubble $R_{\rm init}=R^{*}_{\rm init}/R^{*}_{\rm outer}$. At the initial time, when $t=t_{0}=R_{\rm init}^{2}/2$, we then have $b(t_{0})=1$. When imposing Zheng et al.’s [1] power-law behaviour for the plate separation we have $b(t)=b_{1}\ t^{1/7}$, where $b_{1}=b_{1}^{*}\ (2\pi R^{*2}_{\rm outer}/(b_{0}^{*6}Q))^{1/7}=t_{0}^{-1/7}=(R^{2}_{\rm init}/2)^{-1/7}$.

III Axisymmetric solutions and linear stability analysis

Equations (1)-(4) have an axisymmetric solution for which the bubble radius is given by

[TABLE]

while the pressures in the viscous fluid and the gas bubble are

[TABLE]

and

[TABLE]

respectively. To assess the stability of this solution to non-axisymmetric perturbations that change the shape of the interface, we assume $R(\theta,t)=\bar{R}(t)+\varepsilon\hat{R}_{k}(t)\sin(k\theta)$, where $\varepsilon\ll 1$ is the amplitude and $k>1$ is the integer wavenumber of the perturbation. A straightforward linear stability analysis then shows that the instantaneous growth rate $\lambda_{k}$ of the small-amplitude perturbation with wavenumber $k$ is given by

[TABLE]

For bubbles that are much smaller than the outer radius of the cell, we have $(1+\bar{R}^{2k})/(1-\bar{R}^{2k})\approx 1$, allowing us to approximate the growth-rate as

[TABLE]

In the absence of lifting, i.e. $b=1$ and $\text{d}b/\text{d}t=0$, this reduces to Paterson’s classical expression for the growth rate in an infinitely large Hele-Shaw cell with constant gap width [3]. A positive rate of plate separation, $\text{d}b/\text{d}t>0$, can be seen to reduce the growth rate of the perturbations and therefore stabilises the axisymmetric state.

We note that, in general, the instantaneous growth rate $\lambda_{k}$ of the small-amplitude perturbations varies with the evolving mean radius of the bubble. As a result, the range of unstable wavenumbers (i.e. wavenumbers for which $\lambda_{k}>0$) and the most unstable wavenumber (for which $d\lambda_{k}/dk=0$) generally change throughout the system’s evolution. One of the key observations made by Zheng et al. [1] is that if the gap width has a power-law behaviour of the form $b(t)\sim t^{1/7}$ then the mode with the largest positive growth rate predicted by Paterson’s analysis – which is used somewhat inconsistently because its derivation assumes that $\text{d}b/\text{d}t=0$ – remains constant.

To assess to what extent Zheng et al.’s conclusions are affected by this inconsistency we note that for an arbitrary time-varying gap width $b(t)$, equations (5) and (9) show that perturbations with wavenumber $k$ decay if

[TABLE]

For $b(t)=b_{1}\ t^{1/7}$ this condition becomes

[TABLE]

where

[TABLE]

is precisely the control parameter appearing in Zheng et al.’s analysis. The inclusion of the $\text{d}b/\text{d}t$ term in our analysis results in the appearance of an additional term in the stability criterion (11) (the constant $1/7$ on the left-hand-side of this equation) but its presence does not re-introduce a time-dependence into this condition. Thus Zheng et al.’s key observation remains unchanged. Furthermore, the condition for a perturbation with wavenumber $k$ to have a negative instantaneous growth rate (implying linear stability with respect to such perturbations) can be expressed in terms of $J$ as

[TABLE]

while, for a given, sufficiently large value of $J$, the wavenumber $k_{\rm max}$ of the perturbation with the largest positive growth rate satisfies

[TABLE]

These results again only differ slightly from those obtained by Zheng et al.’s approach which yields $J_{cr}(k)=k(k+1)$ and $J=3k_{\rm max}^{2}-1$.

Figure 2 compares the two sets of results in a plot of the stability criterion (11) for three different values of $J$. Recall that the axisymmetric state is stable when $S_{k}>0$, therefore wavenumbers for which $S_{k}<0$ represent unstable perturbations. Zheng et al.’s approach can be seen to consistently over-estimate the growth-rate of the non-axisymmetric perturbations because their analysis omits the stabilising effect of the $\text{d}b/\text{d}t$ term. However, both analyses show that the range of unstable wavenumbers and the wavenumber with the fastest growth rate increase with $J$.

IV Direct numerical simulations

The results presented above confirm that for a lifting law of the form proposed by Zheng et al. [1], the wavenumber of the most rapidly growing small-amplitude perturbation remains constant throughout the system’s evolution. To assess if this behaviour is responsible for the development of the non-splitting fingers observed in Zheng et al.’s experiments, we now conduct numerical simulations of the system’s nonlinear evolution following the onset of the linear instability. For this purpose we used an oomph-lib-based [25] finite-element discretization of the governing equations, details of which can be found in [26].

All simulations were performed for $Ca=0.3893$ and $R_{\rm init}=0.05$ as in our previous work [26, 27]. Suitable temporal and spatial convergence tests were performed to ascertain that the results presented below are fully converged.

IV.1 Unimodal perturbations and self-similar fingering

We start by perturbing the initial circular bubble of radius $R_{\rm init}=0.05$ with the single, most rapidly growing mode so that

[TABLE]

where we set the amplitude of the perturbation to $\epsilon=5\times 10^{-4}$, i.e. 1% of the initial radius. Figure 3(a) shows the time-evolution of the interface for $J=600$ (for which $k_{\rm max}=13$), using snapshots of the interface plotted at regular time intervals. Thirteen identical finite-amplitude fingers grow from the initial perturbation and, interestingly, the fingers show no signs of splitting. Instead, they remain symmetric about their radial centreline and appear to approach a self-similar shape. This is confirmed in figure  3(b), where we rescaled the radial coordinate by the radius of the finger tip, $R_{\rm tip}$, for each snapshot. This shows that, after an initial transient evolution that occurs within the first five contours shown in figure 3(b), the successive rescaled interfaces become virtually indistinguishable. Equivalent behaviour was observed in simulations for 14 different values of $J$ in the range $25\leq J\leq 1000$, where in each case we imposed the most rapidly growing single-mode perturbation, with $k_{\rm max}$ given by equation (14).

In figure 4 we plot the tip radius, $R_{\rm tip}$, as a function of time on a log-log scale for a range of $J$ values. As the tip radius increases from its initial value it rapidly approaches a power law behaviour, $R_{\rm tip}\sim t^{3/7}$, which is identical to that of the axisymmetrically growing bubble (see equation (5) and recall that $b(t)\sim t^{1/7}$). This is again consistent with the observed self-similar evolution of the interface shape – if the bubble radius $R(t,\theta)$ approaches a self-similar behaviour such that $R(t,\theta)=f(t)F_{k}(\theta)$ volume conservation requires that $f(t)\sim t^{3/7}$. We note that the curves representing $R_{\rm tip}(t)$ for different values of $J$ almost overlap, implying that the pre-factor (which is always larger than that for the axisymmetrically growing bubble) is approximately independent of $J$, despite significant variations in $k_{\rm max}$ which ranges from 8 to 17.

Next, we explore whether the interface shape evolves towards a self-similar solution when the initial single-mode perturbation has a wavenumber other than $k_{\rm max}$. For this purpose we choose a larger value of $J$, $J=1000$, in order to operate in a regime with a wide range of linearly unstable modes. We then perturb the initial circular bubble as in (15) but replace $k_{\rm max}$ with a wavenumber from the unstable range defined by equation (13). The resulting evolution of the interface is shown in figures 5(a-b) for different values of $k$, starting from a perturbation of amplitude $\epsilon=5\times 10^{-4}$. In each case, the system undergoes a transient initial evolution but ultimately reaches a regime in which a number of identical, self-similar fingers have emerged from the initial perturbation. However, the number of self-similar fingers that develop depends on the wavenumber of the initial perturbation. For the perturbation with $k=k_{\rm max}=17$ (figure 5(a)) we retain the wavenumber of the initial perturbation. If we impose an initial perturbation with $k=7$ (figure 5(b)) the finite-amplitude fingers that develop from the initial perturbation split once (see the zoomed-in region shown in the inset) before approaching a self-similar pattern with 14 fingers. This implies that the number of self-similar fingers that can emerge from the initial perturbation is not directly related to the value of $k_{\rm max}$ from the linear stability analysis.

The tip splitting events that may occur during the early transients that precede the self-similar evolution of the interface are not restricted to tip doubling. This is illustrated in figure 5(c) where we imposed an initial perturbation with $k=5$ and also increased the amplitude of the perturbation to 25% of the initial radius. The early-time growth of the initial perturbation is now followed by a tip-tripling event (see the zoomed-in region shown in the inset) before the bubble approaches a configuration with 15 self-similar fingers.

For all cases shown in figure 5 the imposed unimodal initial perturbation is symmetric about the radial centreline of a sector spanned by an angle $2\pi/k$, shown by the dashed lines. Even when the evolution of the interface involved tip-splitting into two or three fingers, the newly-formed fingers did not compete with each other to break this initial symmetry. The fact that the computations were performed with unstructured meshes which (despite being sufficiently fine to ensure mesh-independent solutions) inevitably introduce small non-axisymmetric perturbations, suggests that the self-similar fingers are linearly stable.

IV.2 Non-unimodal perturbations and disordered fingering

Unimodal initial perturbations, such as those imposed in the simulations presented in the previous section, are, of course, difficult to realise in actual experiments. The initial perturbation to a real bubble are likely to comprise a wide range of modes. To analyse this scenario we now perturb the interface with the full range of linearly-unstable wavenumbers. We introduce these perturbations by adding all of the unstable modes to the initial shape of the interface, each with a randomly generated (and uniformly distributed) relative amplitude $\hat{R}_{k}\in[0,1]$ and phase $\phi_{k}\in[0,2\pi)$,

[TABLE]

where the overall amplitude $\epsilon$ was chosen to be 0.5% of the initial radius. We performed 50 simulations for 14 values of $J$ within the range $25\leq J\leq 1000$.

A typical evolving interface is shown in figure 6 for $J=600$ for which $k_{\rm max}=13$. The snapshots of the interface are shown for different values of the average radius $\langle R\rangle$ and are plotted in different colours. When the mean radius has grown to $\langle R\rangle=0.1$ (red) the random initial perturbation has led to the development of 13 distinct fingers of varying widths and lengths. Differences in the length of adjacent fingers are amplified as the fingers grow (finger competition) and when $\langle R\rangle=0.17$ (cyan), the total number of fingers has decreased from 13 to 10. At $\langle R\rangle=0.17$ (cyan) the four largest fingers have broadened considerably so that their tips are approximately flat. When $\langle R\rangle=0.31$ (dark green) these fingers have undergone a tip-doubling followed by competition between the newly-formed fingers, while another finger has broadened sufficiently for another tip-doubling to occur. Because the fingers are not symmetric about their radial centreline, the tip-doubling events lead to pairs of fingers of different shapes and radial tip positions, so that one of the fingers screens the other and adopts a strongly asymmetric shape; see, for example, the finger identified by the arrow in figure 6. The final pattern shown in figure 6 (blue) at $\langle R\rangle=0.6$ comprises 11 distinct fingers of different shapes, which continue to evolve and interact as they grow.

We characterise the patterns in figure 6 by following [1] and compute the Fourier coefficients $c(k)$ of $R-\langle R\rangle$,

[TABLE]

To compare coefficients computed for different values of $\langle R\rangle$, we normalise the amplitudes $c(k)$ by the area under the corresponding spectral curves and denote the normalised amplitudes by $||c(k)||$. This allows us to focus on the relative amplitudes of the modes rather than their absolute value which grows with increasing radius.

Figure 7 shows that the evolution of the spectral curves is closely correlated with the evolution of the interfacial pattern. The interface shape can be seen to comprise a mixture of modes, with the number of distinct fingers observed at the interface approximately corresponding to the mode with the maximum normalised amplitude, which we label $K_{\rm max}$.

In all the simulations performed with random initial perturbations the value of $K_{\rm max}$ continued to fluctuate throughout the system’s evolution. Typically, $K_{\rm max}$ decreased significantly from its initial value at early times when many modes had comparable amplitudes (see figure 7(a)). At later times, $K_{\rm max}$ was more likely to increase due to tip-splitting events. The overall rate of change of $K_{\rm max}$ decreased with increasing time, but $K_{\rm max}$ continued to fluctuate until fingers reached the edge of the computational domain. This scenario was observed for all the values of $J$ investigated ($J\in[25,1000]$), despite the fact that the growth rate of the imposed small-amplitude initial perturbation is smaller for smaller values of $J$.

Figure 8 shows how variations in $J$ affect $K_{\rm max}$ which we evaluated when $\langle R\rangle=0.6$ – the largest value reached in all simulations (which necessarily terminate when one of the fingers reaches the outer edge of the computational domain). The error bars on the symbols quantify the standard deviation of $K_{\rm max}$ across 50 runs for each value of $J$, and highlight the disordered nature of the interface evolution. The two lines show the most unstable wavenumber $k_{\rm max}$ predicted by our linear stability analysis (solid) and by Zheng et al.’s approach (dashed). For modest values of $J$ the wavenumber $K_{\rm max}$ observed in the numerical simulations is remarkably close to $k_{\rm max}$ from the linear stability analysis. However, for larger values of $J$, $K_{\rm max}$ remains significantly below $k_{\rm max}$. Furthermore, none of the simulations that were started with random initial perturbations showed any signs of settling on a self-similar state – the fingers continued to split and then compete with each other over the entire range of the simulations.

Our results show that the system’s behaviour is fundamentally different when the initial, axisymmetrically expanding bubble is subjected to unimodal or non-unimodal perturbations to its shape. This raises the question in what way the small-amplitude non-unimodal perturbations that we deliberately introduced in the computations shown in figure 6 differ from those introduced unintentionally by the use of unstructured meshes in the simulations in figure 5, say. Clearly, the deliberately introduced unimodal perturbation imposed by equation (15) must have a larger initial amplitude than those caused by the unstructured mesh. This was ensured by performing the computations on sufficiently fine meshes. However, we also have to make sure that the growth rate of the deliberately introduced unimodal perturbation is sufficiently large so that it deforms the interface into a non-linear regime before the other modes have grown to comparable amplitude. This is, of course, easiest to achieve if the deliberately introduced perturbation is the one with the maximum growth rate, as in figure 5(a) where we applied a perturbation with $k=k_{\rm max}=17$. The initial amplitude of $\epsilon=5\times 10^{-4}$ used for the perturbation with $k=7<k_{\rm max}$ in figure 5(b) sufficed to create seven finite-amplitude fingers that subsequently underwent a single tip-splitting event before approaching the self-similar regime. However, when using that amplitude to impose a perturbation with $k=5$ (which has a much smaller but still positive growth rate) the system evolved in a disordered manner, with unequal fingers that kept splitting and competing with each other. The five-fold overall symmetry of the initial perturbation could only be retained by increasing its amplitude, as in figure 5(c) where the initial fingers undergo a tip-tripling before approaching a self-similar regime with 15 non-splitting fingers.

V Summary and discussion

We studied the dynamics of an expanding gas bubble that displaces a viscous liquid in a radial Hele-Shaw cell with a time-varying gap width, focusing on the case of a $t^{1/7}$ power-law for the plate separation. Our full linear stability analysis confirmed Zheng et al.’s [1] prediction that the wavenumber $k_{\rm max}$ of the most unstable small-amplitude perturbation to the axisymmetrically expanding bubble is independent of the bubble’s radius. We then employed numerical simulations to follow the growth of the instability into the finite-amplitude regime. This showed that for unimodal perturbations with a wavenumber that matches the most unstable wavenumber from the linear stability analysis, self-similar fingers formed on the interface. This is in stark contrast to the typical behaviour observed in Hele-Shaw cells with constant gap widths where fingers tend to split and compete with each other, resulting in complex, disordered and continuously evolving interface shapes.

Interestingly, self-similar fingers also emerged from unimodal perturbations with wavenumbers other than $k_{\rm max}$. In this case, the finite-amplitude fingers tended to pass through an initial transient period during which they split (via tip doubling or tripling) before ultimately approaching a final regime in which their shapes evolved in a self-similar fashion without any further splitting. The number of identical finite-amplitude, self-similar fingers emerging from this process depended on the initial condition. The self-similar fingers appear to be linearly stable in the sense that they persist despite inevitable small perturbations due to the use of unstructured meshes in the simulations. The fact that the self-similar solutions can only be realised from unimodal perturbations with wavenumber $k$, so that the interface was symmetric about the radial centreline of a sector spanned by an angle $2\pi/k$, suggests that they are unlikely to be observable in actual experiments (unless an experiment happens to be performed in the narrow regime where there is only a single unstable mode).

When starting the simulations from non-unimodal, random perturbations to the initial bubble shape – a scenario more representative of the situation in actual physical experiments – multiple, competing and continuously splitting fingers developed, similar to the behaviour typically found in the constant-gap system. For modest values of the control parameter $J$, the wavenumber $K_{\rm max}$ of the dominant mode contained in the Fourier spectrum of the finite-amplitude fingers tended to be close to $k_{\rm max}$ (as observed in Zheng et al.’s experiments) but for larger values of $J$ we found $K_{\rm max}$ to be significantly smaller than $k_{\rm max}$, suggesting that there is no straightforward relationship between the two quantities. In fact, it is unclear how the fact that the wavenumber of the most-unstable small-amplitude perturbation to the axisymmetric bubble remains independent of its radius could affect the behaviour of the finite-amplitude fingers that emerge from the linear instability: tip-splitting of finite-amplitude fingers occurs in the course of their non-linear evolution (and does not necessarily involve the occurrence of a secondary instability), rather than because the most unstable wavenumber of perturbations to some (far away) axisymmetric state changes. We therefore suggest to interpret the non-splitting fingers observed in this system as self-similar solutions of the governing equations. They are analogues to the steadily propagating fingers that develop in rectangular Hele-Shaw channels [2, 22, 23, 28]. We observed the systematic transient evolution of the interface towards these self-similar states via tip-doubling and tripling. By analogy with rectangular Hele-Shaw channels, this suggests that in radial Hele-Shaw cells other unstable self-similar solutions may exist which are reminiscent of the unstable Romero–Vanden-Broeck multiple-tip solutions [28].

Acknowledgements.

The authors would like to thank Andrew Hazel, Scott McCue and Liam Morrow for helpful discussions. The work of A.J. was funded through EPSRC grant EP/P026044/1. The research data supporting this publication can be found in the supplementary material.