The influence of invariant solutions on the transient behaviour of an   air bubble in a Hele-Shaw channel

J. S. Keeler; A. B. Thompson; G. Lemoult; A. Juel; A. L. Hazel

arXiv:1907.12932·physics.flu-dyn·July 1, 2020

The influence of invariant solutions on the transient behaviour of an air bubble in a Hele-Shaw channel

J. S. Keeler, A. B. Thompson, G. Lemoult, A. Juel, A. L. Hazel

PDF

TL;DR

This study explores how invariant solutions influence the transient dynamics of an air bubble in a Hele-Shaw channel, revealing complex shape transitions and bifurcations through a mathematical model aligned with experimental observations.

Contribution

It demonstrates the application of dynamical systems concepts, including invariant solutions and bifurcations, to model and analyze bubble shape transitions in fluid mechanics.

Findings

01

Existence of stable asymmetric bubble solutions across flow rates.

02

Development of a second stable solution branch above a critical flow rate.

03

Identification of unstable periodic orbits as edge states influencing transient behaviour.

Abstract

We hypothesize that dynamical systems concepts used to study the transition to turbulence in shear flows are applicable to other transition phenomena in fluid mechanics. In this paper, we consider a finite air bubble that propagates within a Hele-Shaw channel containing a depth-perturbation. Recent experiments revealed that the bubble shape becomes more complex, quantified by an increasing number of transient bubble tips, with increasing flow rate. Eventually the bubble changes topology, breaking into two or more distinct entities with non-trivial dynamics. We demonstrate that qualitatively similar behaviour to the experiments is exhibited by a previously established, depth-averaged mathematical model; a consequence of the model's intricate solution structure. For the bubble volumes studied, a stable asymmetric bubble exists for all flow rates of interest, whilst a second stable…

Equations63

\nabla \cdot (b^{3} (y) \nabla p) = 0 \mbox in

\nabla \cdot (b^{3} (y) \nabla p) = 0 \mbox in

p_{b} - p = \frac{1}{3 α Q} (\frac{1}{b ( y )} + \frac{κ}{α}) \mbox o n

n \cdot R_{t} + n \cdot e_{x} U_{b} (t) + b^{2} (y) n \cdot \nabla p = 0 \mbox o n

\frac{\partial p}{\partial y} = 0 \mbox o n

p = 0 \mbox o n

\frac{\partial p}{\partial x} = - 1 \mbox o n

\int b (y) R \cdot n \mbox d Γ = V .

b (y) = 1 - \frac{1}{2} h [tanh (s (y + w)) - tanh (s (y - w))] .

b (y) = 1 - \frac{1}{2} h [tanh (s (y + w)) - tanh (s (y - w))] .

R (\overset{u}{˙}, u, β) = 0,

R (\overset{u}{˙}, u, β) = 0,

R (\overset{u}{˙}, u, β) \equiv M [u] \overset{u}{˙} + β F [u] + G [u] = 0,

R (\overset{u}{˙}, u, β) \equiv M [u] \overset{u}{˙} + β F [u] + G [u] = 0,

u =

u =

\frac{\partial}{\partial t} =

β =

β_{c} F [u_{s}] + G [u_{s}] = 0.

β_{c} F [u_{s}] + G [u_{s}] = 0.

L (\frac{\partial u _{1}}{\partial t _{0}}, u_{1}) \equiv M [u_{s}] \frac{\partial u _{1}}{\partial t _{0}} + J [u_{s}] u_{1} = 0,

L (\frac{\partial u _{1}}{\partial t _{0}}, u_{1}) \equiv M [u_{s}] \frac{\partial u _{1}}{\partial t _{0}} + J [u_{s}] u_{1} = 0,

J [u_{s}] v_{1} \equiv D_{v_{1}} K [u_{s}] \equiv k \to 0 lim \frac{K [ u _{s} + k v _{1} ] - K [ u _{s} ]}{k},

J [u_{s}] v_{1} \equiv D_{v_{1}} K [u_{s}] \equiv k \to 0 lim \frac{K [ u _{s} + k v _{1} ] - K [ u _{s} ]}{k},

u_{1} = A (t_{1}) g (x) \mbox e^{i ω_{c} t_{0}} + \mbox c . c .,

u_{1} = A (t_{1}) g (x) \mbox e^{i ω_{c} t_{0}} + \mbox c . c .,

L (\frac{\partial u _{2}}{\partial t _{0}}, u_{2}) = - M [u_{1}] \frac{\partial u _{1}}{\partial t _{0}} - δ β_{c} F [u_{s}] - \frac{1}{2} H [u_{s}] (u_{1}, u_{1}),

L (\frac{\partial u _{2}}{\partial t _{0}}, u_{2}) = - M [u_{1}] \frac{\partial u _{1}}{\partial t _{0}} - δ β_{c} F [u_{s}] - \frac{1}{2} H [u_{s}] (u_{1}, u_{1}),

H [u_{s}] (v_{1}, v_{2}) \equiv D_{v_{2}} (J [u_{s}] v_{1}) \equiv k \to 0 lim \frac{J [ u _{s} + k v _{2} ] v _{1} - J [ u _{s} ] v _{1}}{k} .

H [u_{s}] (v_{1}, v_{2}) \equiv D_{v_{2}} (J [u_{s}] v_{1}) \equiv k \to 0 lim \frac{J [ u _{s} + k v _{2} ] v _{1} - J [ u _{s} ] v _{1}}{k} .

u_{2} = A^{2} φ_{0} \mbox e^{2 i ω_{c} t_{0}} + \mbox c . c . + ∣ A ∣^{2} φ_{1} + φ_{2},

u_{2} = A^{2} φ_{0} \mbox e^{2 i ω_{c} t_{0}} + \mbox c . c . + ∣ A ∣^{2} φ_{1} + φ_{2},

(2 i ω_{c} M [u_{s}] + J [u_{s}]) φ_{0} = - i ω_{c} M [g] g - \frac{1}{2} H [u_{s}] (g, g),

(2 i ω_{c} M [u_{s}] + J [u_{s}]) φ_{0} = - i ω_{c} M [g] g - \frac{1}{2} H [u_{s}] (g, g),

J [u_{s}] φ_{1} = - \frac{1}{2} H [u_{s}] (g, g^{*}) - \frac{1}{2} H [u_{s}] (g^{*}, g) - i ω_{c} M [g^{*}] g + i ω_{c} M [g] g^{*},

J [u_{s}] φ_{2} = - δ β_{c} F [u_{s}],

L (\frac{\partial u _{3}}{\partial t _{0}}, u_{3}) = - M [u_{2}] \frac{\partial u _{1}}{\partial t _{0}} - M [u_{1}] \frac{\partial u _{2}}{\partial t _{0}} - M [u_{s}] \frac{\partial u _{1}}{\partial t _{1}} - δ β_{c} J_{F} [u_{s}] u_{1} - \frac{1}{2} H [u_{s}] (u_{1}, u_{2}) - \frac{1}{2} H [u_{s}] (u_{2}, u_{1}) - \frac{1}{6} T [u_{s}] (u_{1}, u_{1}, u_{1}),

L (\frac{\partial u _{3}}{\partial t _{0}}, u_{3}) = - M [u_{2}] \frac{\partial u _{1}}{\partial t _{0}} - M [u_{1}] \frac{\partial u _{2}}{\partial t _{0}} - M [u_{s}] \frac{\partial u _{1}}{\partial t _{1}} - δ β_{c} J_{F} [u_{s}] u_{1} - \frac{1}{2} H [u_{s}] (u_{1}, u_{2}) - \frac{1}{2} H [u_{s}] (u_{2}, u_{1}) - \frac{1}{6} T [u_{s}] (u_{1}, u_{1}, u_{1}),

T [u] (v_{1}, v_{2}, v_{3}) \equiv D_{v_{3}} (H [u_{s}] (v_{1}, v_{2})) \equiv k \to 0 lim \frac{H [ u _{s} + k v _{3} ] ( v _{1} , v _{2} ) - H [ u _{s} ] ( v _{1} , v _{2} )}{k} .

T [u] (v_{1}, v_{2}, v_{3}) \equiv D_{v_{3}} (H [u_{s}] (v_{1}, v_{2})) \equiv k \to 0 lim \frac{H [ u _{s} + k v _{3} ] ( v _{1} , v _{2} ) - H [ u _{s} ] ( v _{1} , v _{2} )}{k} .

⟨ L u, v ⟩ = ⟨ u, L^{†} v ⟩, ⟨ u, v ⟩ = \int_{0}^{T} \int_{x \in Ω} u v^{*} \mbox d x \mbox d t_{0}

⟨ L u, v ⟩ = ⟨ u, L^{†} v ⟩, ⟨ u, v ⟩ = \int_{0}^{T} \int_{x \in Ω} u v^{*} \mbox d x \mbox d t_{0}

\overset{ν}{^} \frac{\partial A}{\partial t _{1}} + \hat{λ} A + \overset{μ}{^} A ∣ A ∣^{2} = 0,

\overset{ν}{^} \frac{\partial A}{\partial t _{1}} + \hat{λ} A + \overset{μ}{^} A ∣ A ∣^{2} = 0,

\overset{ν}{^} = ⟨ M [u_{s}] g, g^{†} ⟩, \hat{λ} = k = 0 \sum 3 ⟨ Γ_{k}, g^{†} ⟩, \overset{μ}{^} = k = 0 \sum 7 ⟨ Λ_{k}, g^{†} ⟩,

\overset{ν}{^} = ⟨ M [u_{s}] g, g^{†} ⟩, \hat{λ} = k = 0 \sum 3 ⟨ Γ_{k}, g^{†} ⟩, \overset{μ}{^} = k = 0 \sum 7 ⟨ Λ_{k}, g^{†} ⟩,

Γ_{0} = Λ_{0} = Λ_{4} = Λ_{7} = δ β_{c} J_{F} [u_{s}] g, Γ_{1} = i ω_{c} M [φ_{2}] g, Γ_{2} = \frac{1}{2} H [u_{s}] (g, φ_{2}), Γ_{3} = \frac{1}{2} H [u_{s}] (φ_{2}, g), - i ω_{c} M [φ_{0}] g^{*}, Λ_{1} = i ω_{c} M [φ_{1}] g, Λ_{2} = 2 i ω_{c} M [g^{*}] φ_{0}, Λ_{3} = \frac{1}{2} H [u_{s}] (g, φ_{1}), \frac{1}{2} H [u_{s}] (g^{*}, φ_{0}), Λ_{5} = \frac{1}{2} H [u_{s}] (φ_{1}, g), Λ_{6} = \frac{1}{2} H [u_{s}] (φ_{0}, g^{*}), \frac{1}{6} T [u_{s}] (g^{*}, g, g) + \frac{1}{6} T [u_{s}] (g, g^{*}, g) + \frac{1}{6} T [u_{s}] (g, g, g^{*}) .

Γ_{0} = Λ_{0} = Λ_{4} = Λ_{7} = δ β_{c} J_{F} [u_{s}] g, Γ_{1} = i ω_{c} M [φ_{2}] g, Γ_{2} = \frac{1}{2} H [u_{s}] (g, φ_{2}), Γ_{3} = \frac{1}{2} H [u_{s}] (φ_{2}, g), - i ω_{c} M [φ_{0}] g^{*}, Λ_{1} = i ω_{c} M [φ_{1}] g, Λ_{2} = 2 i ω_{c} M [g^{*}] φ_{0}, Λ_{3} = \frac{1}{2} H [u_{s}] (g, φ_{1}), \frac{1}{2} H [u_{s}] (g^{*}, φ_{0}), Λ_{5} = \frac{1}{2} H [u_{s}] (φ_{1}, g), Λ_{6} = \frac{1}{2} H [u_{s}] (φ_{0}, g^{*}), \frac{1}{6} T [u_{s}] (g^{*}, g, g) + \frac{1}{6} T [u_{s}] (g, g^{*}, g) + \frac{1}{6} T [u_{s}] (g, g, g^{*}) .

\frac{\partial A}{\partial t _{1}} = λ A + μ A ∣ A ∣^{2},

\frac{\partial A}{\partial t _{1}} = λ A + μ A ∣ A ∣^{2},

H [u_{s}] (v_{1}, v_{2}) \approx \frac{J [ u _{s} + k _{1} v _{2} ] v _{1} - J [ u _{s} - k _{1} v _{2} ] v _{1}}{2 k _{1}},

H [u_{s}] (v_{1}, v_{2}) \approx \frac{J [ u _{s} + k _{1} v _{2} ] v _{1} - J [ u _{s} - k _{1} v _{2} ] v _{1}}{2 k _{1}},

T [u_{0}] (v_{1}, v_{2}, v_{3}) \approx \frac{J [ u _{s} + k _{2} v _{2} + k _{2} v _{3} ] v _{1} - J [ u _{s} + k _{2} v _{2} ] v _{1} - J [ u _{s} - k _{2} v _{2} + v _{3} ] v _{1} + J [ u _{s} - k _{2} v _{2} ] v _{1}}{2 k _{2}^{2}} .

T [u_{0}] (v_{1}, v_{2}, v_{3}) \approx \frac{J [ u _{s} + k _{2} v _{2} + k _{2} v _{3} ] v _{1} - J [ u _{s} + k _{2} v _{2} ] v _{1} - J [ u _{s} - k _{2} v _{2} + v _{3} ] v _{1} + J [ u _{s} - k _{2} v _{2} ] v _{1}}{2 k _{2}^{2}} .

\mbox H 1 : Q_{c} = Q_{1} = 0.0301, ω_{c} = 3.58,

\mbox H 1 : Q_{c} = Q_{1} = 0.0301, ω_{c} = 3.58,

\mbox H 2 : Q_{c} = Q_{2} = 0.0896, ω_{c} = 10.8,

\frac{\mbox d r ^}{\mbox d t _{1}} = a \overset{r}{^} + b \overset{r}{^}^{3}, \frac{\mbox d θ}{\mbox d t _{1}} = c + d \overset{r}{^}^{2},

\frac{\mbox d r ^}{\mbox d t _{1}} = a \overset{r}{^} + b \overset{r}{^}^{3}, \frac{\mbox d θ}{\mbox d t _{1}} = c + d \overset{r}{^}^{2},

u = u_{s} + ε^{2} φ_{2} + O (ε^{4}) .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

The influence of invariant solutions on the transient behaviour of an air bubble in a Hele-Shaw channel.

J. S. Keeler

[email protected]

School of Mathematics and Manchester Centre for Nonlinear Dynamics (MCND), University of Manchester, Oxford Road, Manchester, M13 9PL, UK

A. B. Thompson

[email protected]

School of Mathematics and Manchester Centre for Nonlinear Dynamics (MCND), University of Manchester, Oxford Road, Manchester, M13 9PL, UK

G. Lemoult

[email protected]

Normandie Université, UniHavre, CNRS, UMR 6294, Laboratoire Onde et Milieux Complexes (LOMC) 53, rue de Prony, 76058 Le Havre Cedex, France.

A. Juel

[email protected]

School of Mathematics and Manchester Centre for Nonlinear Dynamics (MCND), University of Manchester, Oxford Road, Manchester, M13 9PL, UK. School of Physics and Astronomy and MCND, University of Manchester, Oxford Road, Manchester, M13 9PL, UK

A. L. Hazel

[email protected]

School of Mathematics and Manchester Centre for Nonlinear Dynamics (MCND), University of Manchester, Oxford Road, Manchester, M13 9PL, UK

[email protected]

Abstract

We hypothesize that dynamical systems concepts used to study the transition to turbulence in shear flows are applicable to other transition phenomena in fluid mechanics. In this paper, we consider a finite air bubble that propagates within a Hele-Shaw channel containing a depth-perturbation. Recent experiments revealed that the bubble shape becomes more complex, quantified by an increasing number of transient bubble tips, with increasing flow rate. Eventually the bubble changes topology, breaking into two or more distinct entities with non-trivial dynamics. We demonstrate that qualitatively similar behaviour to the experiments is exhibited by a previously established, depth-averaged mathematical model; a consequence of the model’s intricate solution structure. For the bubble volumes studied, a stable asymmetric bubble exists for all flow rates of interest, whilst a second stable solution branch develops above a critical flow rate and transitions between symmetric and asymmetric shapes. The region of bistability is bounded by two Hopf bifurcations on the second branch. By developing a method for a numerical weakly nonlinear stability analysis we show that unstable periodic orbits emanate from the Hopf bifurcation at the lower flow rate and, moreover, that these orbits are edge states that influence the transient behaviour of the system.

I Introduction

A Hele-Shaw channel consists of two parallel glass plates, separated by a distance much smaller than the width of the channel. If a trapped viscous fluid is extracted at a constant flux from one end of the channel and an air bubble is placed at the other end, then the bubble will propagate and change shape as it does so. For sufficiently large bubbles the only stable solution is for the bubble to propagate symmetrically along the centreline of the channel; a solution analogous to the symmetric semi-infinite air finger that develops when one end of the channel is left open to the atmosphere saffman1958penetration . For higher flow rates in large aspect ratio channels, propagating air fingers develop complex patterns via multiple tip splitting events as well as side-branching Tabeling1987 ; MooreJuel2002 . The onset of this complex interfacial dynamics appears to be a subcritical transition; a feature that it shares with the transition to turbulence in shear flows. Specifically, the steady symmetrically propagating solution is linearly stable for all values of the flow rate at which it has been computed, meaning that finite perturbations are required to initiate the complex dynamics. Moreover, the value of the critical dimensionless flow rate for the onset of patterns cannot be precisely determined and is very sensitive to the level of perturbation in the system: the transition occurs at lower dimensionless flow rates as the roughness of the channel walls is increasedTabeling1987 .

In this paper, we concentrate on bubble propagation in a geometrically-perturbed Hele-Shaw cell: a rectangular prism is added to the base of the channel, as sketched in figure 1. By studying bubbles, rather than air fingers, and working in a co-moving frame, we can follow long-time evolution of the system in short computational domains. This relatively simple system still exhibits a wide variety of nonlinear dynamical phenomena de2009tube ; pailha2012oscillatory ; hazel2013multiple ; jisiou2014geometry ; thompson2014multiple ; franco2016sensitivity ; franco2017propagation ; franco2017bubble and the recent experimental results of Franco-Gómez et alfranco2017propagation have shown that the system contains regions of bistability in which finite perturbations can provoke a range of possibilities for the time-evolution of an initially centred bubble. For small volume fluxes the bubble will eventually settle towards a stable asymmetric state on one side (or the other) of the depth perturbation, but if the volume flux exceeds a critical threshold the bubble shape becomes increasingly deformed before eventually breaking up into two or more distinct parts, as sketched in figure 2(b), despite theoretical predictions that a stable steady state exists for all flow rates; a feature preserved from the unperturbed Hele-Shaw channel. Thus, this model system allows us to explore whether the dynamical systems concepts applied in the study of transition to turbulence in shear flows also apply to the transient behaviour of other canonical problems in fluid mechanics.

We study the behaviour theoretically by finding invariant solutions, namely steady states and periodic orbits, of the system. We pursue the idea, hypothesised by Franco-Gómez et al franco2017propagation , that when the flux is large enough the complex time-dependent behaviour of a single bubble can be interpreted as a transient exploration of the stable manifolds of weakly unstable edge states. We consider only the dynamics before the changes in topology when the bubble breaks up into two or more separate bubbles: systems with their own dynamics that we do not pursue here.

In the context of fluid mechanics an edge state is an invariant solution of the governing equations that has only a small number of unstable eigenvalues and whose stable manifold forms the ‘boundary’ between two qualitatively different dynamical outcomes. Edge tracking techniques have been utilised in a number of different scenarios including shear flow and pipe flow kerswell2005recent ; schneider2007turbulence ; schneider2008laminar ; eckhardt2008dynamical and more recently droplet breakup gallino2018edge . In these studies the edge state was found by direct numerical simulation of the governing equations and interval bisection of the initial conditions. This methodology is easy to implement numerically but can be computationally expensive and does not reveal whether the edge state corresponds to a steady state, periodic orbit or other invariant solution.

One advantage of the Hele-Shaw system is that the only nonlinearities in the system arise due to boundary conditions on the bubble because the reduced Reynolds number is extremely small. Thus the shape of the interface gives an immediate visual representation of nonlinear behaviour. Another advantage is that by assuming that the width of the channel is large compared to its height, the behaviour of the system can be described by a depth-averaged set of equations that is more amenable to analysis than the full Navier–Stokes equationssaffman1958penetration ; taylor1959note . These reduced equations (stated in § II), often called Darcy or Hele-Shaw equations, have been used in most theoretical studies of this system, and lead to predictions for the bubble shape as viewed from above (sketched in figure 2(b)). In the unperturbed system at zero surface tension, exact steady solutions can be found by conformal mapping techniques, and for each bubble volume, there is a two-parameter family of solutions described by the centroid offset and bubble speed tanveer1987surfacetension ; tanveer1987stability ; tanveer1987newsolutions . The introduction of surface tension selects both a main, symmetric branch of linearly stable steady solutions which persists for all fluxes, along with a countably infinite sequence of unstable ‘exotic’ bubble shapes; these exist in both channels bounded by parallel side walls tanveer1987stability and in unbounded channels green2017effect . Figure 2(a) shows the bubble shape and speed for the first three solution branches in channels with a rectangular cross-section. Introducing a depth-perturbation to the system (see sketch in figure 1) allows the solution branches to interact, which results in the diverse range of observed stable steady states and time-dependent behaviour de2009tube ; pailha2012oscillatory ; hazel2013multiple ; jisiou2014geometry ; thompson2014multiple ; franco2016sensitivity ; franco2017propagation ; franco2017bubble .

We will show that the bifurcation structure of the depth-averaged equations does indeed feature bistability between steady states, and that previously unknown unstable periodic orbits influence whether the bubble breaks up or returns to a stable configuration. We have sufficient knowledge of the solution structure that, rather than edge tracking, we can develop a general numerical procedure for weakly nonlinear analysis that allows us to determine an approximation to the periodic orbits and steady states. We emphasise that the method described in this paper, although applied to our specific set of equations, can readily be adapted for other problems and we provide a numerical recipe for this procedure. The results of the weakly nonlinear analysis, along with time-dependent calculations, show that the previously unknown unstable periodic orbits are indeed edge states of the system that influence the eventual fate of the bubble.

The paper is set out as follows. Initially, in § II, we introduce the details of the system and describe the depth-averaged governing equations, non-dimensional parameters and numerical methods. In § III, we start our investigation with a range of initial-value calculations and compare these to the previous experiments franco2017propagation . In § IV, we present an extended bifurcation diagram for the invariant steady states and perform a linear stability analysis. In § V we derive our general method for determining the weakly nonlinear approximation to the invariant periodic orbits of the system by employing the method of multiple scales near the bifurcation points. Finally, in § VI, we summarise our results and briefly discuss the importance of the periodic orbits to the transient behaviour of the system and their interpretation as edge states.

II Governing Equations

The physical situation is shown in figure 1. The channel has outer dimensional width $W^{*}$ and height $H^{*}$ . The channel is filled with a viscous fluid of density $\rho$ and dynamic viscosity $\mu$ , containing a single air bubble with known volume $V^{*}=(W^{*})^{2}H^{*}V/4$ , where $V$ is the non-dimensional volume. Fluid flow is driven by withdrawing fluid at a constant flux $Q^{*}$ far ahead of the bubble. The system is non-dimensionalised based on a length scale $W^{*}/2$ in the $x$ and $y$ directions, and on the velocity scale $U_{0}^{*}=Q^{*}/(W^{*}H^{*})$ . We define the non-dimensional flow rate $Q=\mu U_{0}^{*}/\gamma,$ where $\gamma$ is the coefficient of surface tension at the air/fluid interface. Additionally we define the aspect ratio of the channel $\alpha=W^{*}/H^{*}$ . Details of the derivation of the depth-averaged equations from the Navier-Stokes equations can be found in thompson2014multiple ; franco2017propagation . The critical assumptions are that the channel aspect ratio $\alpha$ is large; the reduced Reynolds number $U_{0}^{*}W^{*}/(\rho\mu\alpha^{2})$ and Bond number $\mbox{Bo}=\rho gL^{*2}/\gamma$ are small; the bubble occupies the full height of the channel; and the component of curvature as $z$ varies corresponds to a semi-circle filling the channel height with the fluid perfectly wetting the upper and lower walls. We neglect the thin film corrections proposed by Homsy homsy1987viscous and Reinelt reinelt1987interface , because they do not change the qualitative comparison with the experiment, and it is not obvious how they should be modified due to the presence of the depth-perturbation and multiple solutions thompson2014multiple ; franco2016sensitivity ; franco2017propagation .

The nondimensional domain is shown in figure 3. The depth-averaged velocity at any point within the fluid is given by $\mathbf{\hat{u}}=-b^{2}(y)\nabla p$ , where $p=p(x,y)$ is the fluid pressure and $b(y)$ is the height of the channel (see figure 3(b)). We work in a frame moving at dimensionless speed $U_{b}(t)$ chosen at each time step so that the $x$ component of the bubble centroid remains fixed in this moving frame. The equations to be solved are

[TABLE]

Here $\Omega$ denotes the two-dimensional fluid domain (the shaded region in figure 3(a)) and $\Gamma$ the bubble boundary, $p_{b}(t)$ is the unknown pressure inside the bubble, $\textbf{R}=(x_{b},y_{b})$ is the position of the bubble interface, $\kappa$ is the curvature of this 2D interface, and $\mathbf{n}$ is a unit normal vector directed out of the bubble. Equation (1a) is the equation for the pressure in the fluid domain whilst equations (1b) and (1c) are the dynamic and kinematic conditions on the bubble surface. Our depth-averaged equations are second order in space, and hence we can apply only no-penetration conditions, (1d); likewise the condition (1b) derives from the normal component of the dynamic boundary condition at the bubble interface, and we neglect the tangential stress balance. Equation (1e) imposes the pressure to be zero far behind the bubble at $x=-L$ and equation (1f) achieves the constant volume flux at $x=L$ by imposing a constant pressure gradient. Additionally, the unknown value of $p_{b}$ is determined by ensuring the volume constraint, (1g), is satisfied. The computational domain is truncated at a fixed distance $L$ ahead of and behind the bubble. Following previous papers, we model the depth-perturbation by a smoothed profile:

[TABLE]

A typical example of this smoothed profile is shown in figure 3(b). The dimensionless depth-perturbation height $h$ can be viewed as a controlled, axially uniform perturbation of the rectangular channel. The topography profile $b(y)$ enters the equations through the bulk equation (which is no longer amenable to conformal mapping), the kinematic boundary conditions and the variable transverse curvature of the interface. The unknowns of the problem are $[p,\textbf{R},p_{b},U_{b}]$ and the control parameters are $[Q,V,h,w,s,\alpha]$ . An additional solution measure is the centroid of the bubble, denoted $\overline{y}$ , which is useful in determining the symmetry of the system. (1) are solved by a finite-element discretisation, utilising the open-source software oomph-lib heil2006oomph . Details of the implementation can be found in thompson2014multiple for the case of a propagating air-finger and franco2017propagation for a finite air bubble.

It is convenient to keep $h,w$ and $\alpha$ constant and analyse the system mainly through variations in $Q$ and $V$ , which are easily manipulated in a physical, experimental set-up. Unless otherwise specified the results presented in this paper are compared to the values used by Franco-Gómez et al franco2017propagation and hence a bubble of fixed volume $V$ is chosen so that the projected area is $A=\pi r^{2}$ with $r=0.46$ , together with a perturbation height $h=0.024$ , width $w=0.25$ , sharpness $s=10$ and aspect ratio $\alpha=40$ .

III Initial-value calculations

Franco-Gómez et al franco2017propagation performed a range of experiments which showed that an initially circular bubble would settle to a stable state below a critical flow rate whilst above this threshold the bubble becomes increasingly deformed, resulting in a variety of transient outcomes, see figure 0, including topological break up, despite evidence to suggest that an asymmetric stable state existed for all flow rates.

In order to investigate the potential range of evolution in our model we perform a number of numerical simulations designed to mimic these experiments. The nonlinear dynamics can be examined by recording the bubble shape as it propagates along the channel. We do this by solving the equations in (1) with the initial condition of a circular bubble with radius $r$ centered at $(x_{c},y_{c})=(0.0,0.01)$ and released from rest. Our simulations reveal a number of qualitatively distinct modes of propagation. The eventual propagation mode selected has a dependence on the flow rate $Q$ and the radius $r$ which sets the bubble volume $V$ . A summary of the qualitative behaviour is shown in figure 5, with regions of parameter space classified according to the characteristic shape first adopted by the bubble. As $Q$ and $V$ are varied, we find that the first emergent shape may have between one and four tips. The start of each different coloured region as $Q$ increases indicates the place in parameter space where the bubble first develops two, three or four tips for a given bubble volume. Single tipped bubbles evolve smoothly to an asymmetric steadily-propagating configuration (see inset labelled i in figure 5), while multiple-tipped bubbles may self-intersect in the simulations and thus break up into one or more fragments, in which case the simulations are terminated (see insets labelled ii,iii,iv in figure 5) . In some regimes, the bubble may exhibit oscillatory behaviour before decaying to either a steady state (see inset v in figure 5) or stable periodic oscillations (see inset vi). The nature of the oscillatory modes that tend to a stable state are dependent on the flow rate in that for lower flow rates ( $Q\approx 0.05$ ) the bubble appears to oscillate around a symmetric steady state whilst for larger flow rates ( $Q\approx 0.1$ ) the oscillations manifest as ‘waves’ above one edge of the depth-perturbation with the opposite side of the asymmetric bubble above the other edge of the depth-perturbation as seen in the iii inset in figure 5.

We can compare these numerical results to the experimental data of Franco-Gómez et al franco2017propagation shown in figure 0. Although the boundaries between regions are quantitatively different in figures 0 and 5, the phenomenological behaviour is in agreement, and the transition between different modes of propagation occur over similar ranges of $Q$ and $V$ . For nearly all values of $Q$ explored in figure 5, varying the bubble volume at fixed $Q$ allows us to access a range of qualitatively different modes of propagation.

As demonstrated by these numerical time simulations, there is a large range of possible outcomes from an initially circular bubble including oscillatory transients and multiple stable steady states. To understand this further we now perform a detailed analysis of the steady state bifurcation structure.

IV Bifurcation Structure for Steady States

In this section we analyse the steady solution space. We use the flow rate, $Q$ , as a continuation parameter for a fixed bubble volume and record the bubble speed, $U_{b}$ , and the centroid, $\overline{y}$ as convenient solution measures. We also analyse the linear stability of the numerical solutions as the solution branches are traced out. The solution space is shown in figure 6 in the $(Q,U_{b})$ and $(Q,\overline{y})$ projections respectively, with stable (unstable) solution branches denoted by solid (dashed) lines. We note that the branches in this steady solution space are connected through several bifurcations, in contrast to the $h=0$ case (see figure 2(a)), where the solution branches are disjoint.

For small $Q$ there are three distinct branches as shown in figure 6. The ‘upper’ branch (with largest $U_{b}$ ), AS1, is stable, characterised by an asymmetric bubble shape and persists for all values of $Q$ that were sampled. The other two branches, denoted as S1a and S1b, correspond to symmetric double-tipped solutions. As the value of $h$ is decreased to zero the S1a branch approaches the $m=0$ solution sketched in figure 2(b) whilst the lower part becomes the double-tipped solution $m=1$ solution. The S1 branch only exists for a small range of the parameter $Q$ , in contrast to the case with no depth-perturbation where the solution persists for all $Q$ . Both parts of the S1 branch are unstable; the segment with larger $U_{b}$ has one unstable eigenmode, while the slower segment is doubly unstable. Note that although it appears from figure 6(a) that the S1 and AS1 branches intersect at small $Q$ , these branches are actually disjoint as can be seen by the projection of the solution onto the $(Q,\overline{y})$ plane in figure 6(c).

For $Q\gtrapprox 0.01$ , we are able to compute a three-tipped symmetric solution, S2a, which is unstable. As $Q$ increases this branch undergoes two fold bifurcations. The first fold, F1, (shown in figure 6(b)) separates the S2a and S2b branches and the second fold, F2, separates the S2b and S2c branches. By decreasing $h$ to zero, the S2a branch converges to a triple-tipped solution branch, the S2b branch becomes the double-tipped solution branch and the S2c branch becomes the stable solution branch (this transition sequence is likely to involve interaction between the AS1, S1 and S2 branches). All of the S2 branches are unstable but for this bubble volume and channel geometry, the solution finally stabilises via a Hopf bifurcation H1 just beyond the second fold point. The stable branch of steady solutions that results from this Hopf bifurcation is labelled S3.

As $Q$ increases further on the S3 branch, a supercritical pitchfork bifurcation (PF) occurs where the S3 solution breaks symmetry and two stable asymmetric branches, denoted AS2, emanate from the bifurcation, see figure 6. Beyond the pitchfork, the unstable symmetric branch experiences multiple bifurcations as $Q$ is varied and becomes increasingly unstable. This collection of unstable branches is denoted by S4. The most extreme shapes on the S4 branch, with $U_{b}$ close to 1, correspond to the eight-tipped solution as seen in figure 6. Note the loops on the S2 branch around $Q\approx 0.11$ where an asymmetric branch connects to the main S4 branch via two pitchfork bifurcations. This is reminiscent of the so called ‘snakes and ladders’ bifurcations seen in other nonlinear systems, see, for example, burke2007snakes ; schneider2010snakes . This branch has at least two unstable eigenmodes and as we do not believe it affects the transient behaviour of the system we do not pursue the details of this branch here. For the largest $Q$ in our calculations ( $Q=0.15$ ), two distinct asymmetric states persist in the $(Q,U_{b})$ projection, (four in the $(Q,\overline{y})$ projection). The AS1 branch is always stable but the stable AS2 branch experiences a Hopf bifurcation, denoted by H2, and becomes unstable, the resulting branch labelled AS3.

In contrast with the $h=0$ solution space, where only one solution is ever stable, in our case there is a finite-width region of bistability (with respect to steady states) between the H1 and H2 bifurcations. The AS1 solution is always stable but the nature of the second stable solution changes from symmetric (the S2 branch) to asymmetric (the AS2 branch), via the pitchfork bifurcation. The symmetric stable solution (S2) has been observed before in this system, see Franco-Gómez et al franco2017propagation , but its transition to an asymmetric stable solution as $Q$ increases, and indeed the presence of the two Hopf bifurcations, is a new observation.

The results presented here are for a fixed volume, $V$ , and fixed depth-perturbation height, $h$ but we can use bifurcation tracking calculations to investigate the robustness of the location and order of the bifurcations F2, H1, PF and H2 that bound the bistable region. The locations of these bifurcations as functions of $V$ and $h$ are shown in figure 7. We find that the ordering of the F2, H1, PF and H2 bifurcations are very robust to changes in volume and, with the exception of H2, their position is fixed, see figure 7(a). As $V$ decreases, H2 migrates to larger values of $Q$ and hence the width of the bistable region increases. All the bifurcations are more sensitive to $h$ than $V$ , and when $h<h_{c}\approx 0.012$ , the change in stability occurs at the fold, F2, instead of the Hopf, H1, see figure 7(b). A fold-Hopf bifurcation occurs when $h=h_{c}$ . The dynamics of the system near this type of co-dimension two bifurcation are very complex and can lead to the appearance of invariant tori and heteroclinic orbits, see, for example, kuznetsov2013elements .

The time-dependent behaviour explored in § III can now be interpreted in terms of the bifurcation structure. For all flow rates, there is an asymmetric stable mode of propagation, which corresponds to the AS1 branch. The region of bistability appears to coincide with the solid symbols in figure 5 which denote that the bubble is evolving towards a steady solution. The multi-tipped modes of propagation may relate to the presence of the unstable S1, S2 and S4 branches, while the oscillating modes of propagation may be expected to arise due to the Hopf bifurcations. However, the linear stability analysis does not reveal the centre, nor the stability and size of the periodic orbits emanating from the Hopf bifurcations and thus the criticality of H1 and H2 remains unknown.

There are at least three possible approaches to determine the periodic orbits. If the periodic orbit is stable, we expect that for suitable initial conditions, the periodic orbit is attracting and hence initial-value simulations will eventually converge towards the periodic orbit. These simulations may be expensive if the convergence rate is slow, and in any case will not capture unstable periodic orbits. The second approach is to calculate the periodic orbits directly by solving an extended system of equations as proposed by, for example, net2015continuation . This direct solution will capture unstable periodic orbits, but leads to a significantly larger system of equations and requires a discretization in both space and time. The third approach is to perform a numerical weakly nonlinear stability analysis near the Hopf points, leading to an analytic normal form of the perturbation equations near the Hopf points with numerically computed coefficients, providing an approximate expression for the periodic orbits in terms of eigenfunctions. This latter approach is the least computationally expensive and has the advantage of providing semi-analytic approximations for the amplitude, period and location of the periodic orbits and we pursue this analysis in the next section.

V Weakly Nonlinear Stability

The aim of this section is to perform a local analysis near the Hopf bifurcation points to obtain approximations for the location, stability and size of the periodic orbits. In linear stability theory the growth of the perturbation occurs on a single time-scale and is determined by the eigenvalues alone and not the amplitude of the perturbation. In weakly nonlinear stability theory near a marginally stable solution (i.e. bifurcation point) the growth/decay of a perturbation is assumed to happen over two time scales, the eigenvalue providing the growth/decay on the fast scale $t$ , whilst the amplitude varies on a slower time-scale. The main objective of a weakly nonlinear analysis is to obtain an evolution equation for this unknown amplitude function. This method has been applied to a number of specific examples in physics, including shear flows, shallow water waves, thermoacoustics and magnetohydrodynamics, see, for example fujimura1989equivalence ; bera_poiseuille ; orchiniweaklynonlinear ; sanchez2006amplitude ; ghidaouiwaterwavenonlinear ; laroze2010amplitude ; laroze2009amplitude . In these examples the nature of the equations analysed means that the algebra often becomes very complicated. Our approach (inspired by sanchez2006amplitude ) is to perform the analysis for a general set of equations that can be applied to a wide range of systems so that an analytic approximation for the periodic orbits and steady states can be obtained. The analysis is based on a continuous set of equations that are independent of the nature of a discretisation procedure. We then describe a versatile numerical procedure that can be implemented to obtain semi-analytic expressions for the periodic orbits and steady states. We examine a set of equations of the form

[TABLE]

where $\mathcal{R}$ is a nonlinear function that depends on the state variables $u\in\mathcal{U}$ , where $\mathcal{U}$ is an appropriately defined Hilbert space, time derivatives $\dot{u}$ and a parameter $\beta\in\mathbb{R}$ . The state variable, $u$ , will in general depend continuously on spatial coordinates, $x$ and temporal coordinates $t$ . To proceed further it is assumed that the set of equations can be separated so that

[TABLE]

where $\mathcal{M}$ is a linear mass operator and $\mathcal{F}$ , $\mathcal{G}$ are nonlinear operators on the state variables independent of the parameter $\beta$ . The form of the equations in (4), whilst not completely general, are representative of a large number of physical systems where time derivatives appear in linear combinations, including the equations that form the subject of this paper, (1).

We are interested in the nonlinear evolution of a perturbation near a Hopf point. A standard linearisation procedure about a steady solution, $u_{s}$ , will yield a generalised eigenproblem that can be solved to find the eigenmodes, denoted $g$ , and corresponding eigenvalues, denoted $s$ . A Hopf point occurs at a particular value of the parameter, $\beta_{c}$ say, such that a single pair of eigenvalues are located on the imaginary axis, i.e. of the form $s=\pm\mathrm{i}\omega_{c}$ . Near this point of marginal stability the oscillations will occur on a fast time scale and will be modulated by a time-dependent amplitude function that operates on a slower time scale. Following sanchez2006amplitude ; orchiniweaklynonlinear ; drazin2004hydrodynamic we define a fast time scale, $t_{0}=t$ and a slow time scale, $t_{1}=\varepsilon^{2}t$ , where $\varepsilon\ll 1$ is a small unfolding parameter, and employ the method of multiple scales. The solution, $u$ , is now a function of two different time scales as well as the spatial variables denoted by $x$ . The solution, time derivative and parameter $\beta$ are therefore expanded in ascending powers of $\varepsilon$ :

[TABLE]

The parameter $\delta=\pm 1$ signifies which ‘side’ of the Hopf bifurcation is being analysed. Substituting the expressions in (5) into (4) yields a sequence of linear problems that can be solved at each order of $\varepsilon$ .

Examining terms of $\mathcal{O}(\varepsilon^{0})$ in the expansion gives the following equation to solve for the steady base solution,

[TABLE]

With $u_{s}$ known, we continue further and at $\mathcal{O}(\varepsilon)$ we have

[TABLE]

where the operator $\mathcal{J}[u_{s}]$ is the Jacobian operator defined by its action on $v_{1}\in\mathcal{U}$ by the Fréchet derivative,

[TABLE]

where $\mathcal{K}[u]\equiv\beta_{c}\mathcal{F}[u]+\mathcal{G}[u]$ . We separate $u_{1}$ into spatial and temporal parts so that $u_{1}(t_{0},t_{1},x)=A(t_{1})g(x)\mbox{e}^{st_{0}}$ and equation (7) becomes a generalised eigenproblem. Its solution determines the linear stability of the steady solutions solved in (6). There are an infinite number of solutions to (7) but as we are only interested in periodic solutions arising from a Hopf bifurcation, we choose $u_{1}$ to be

[TABLE]

where $A(t_{1})$ is the undetermined amplitude of the perturbation depending on the slow time scale $t_{1}$ . The function ${g}(x)$ is a complex eigenfunction of the generalised eigenvalue problem, $\mathrm{i}\omega_{c}\mathcal{M}{g}=\mathcal{J}{g}$ . Proceeding to $\mathcal{O}(\varepsilon^{2})$ the problem to be solved is

[TABLE]

where $\mathcal{H}[u_{s}]$ represents the bilinear Hessian operator defined its action on $v_{1},v_{2}\in\mathcal{U}$ and using (8) by

[TABLE]

The terms on the right hand side of (10) are either formed of products of $u_{1}$ or have no time dependence. It can be shown that the products of $u_{1}$ , defined in (9), introduce time dependent terms proportional to $\mbox{exp}(2\mathrm{i}\omega_{c}t_{0})$ in (10) and these do not resonate with the homogeneous solutions of the operator $\mathcal{L}$ . By examining the form of ${u}_{1}$ in (9) a particular solution to (10) is sought in the form

[TABLE]

where the functions $\varphi_{i}$ are to be determined. Substituting (12) into (10) and then equating coefficients of $A^{2}$ , $|A|^{2}$ and constant terms leads to three linear equations that can be solved to find the undetermined functions ${\varphi}_{i}$ in (12). These equations can be stated as

[TABLE]

where the star denotes complex conjugation. The amplitude function $A(t_{1})$ remains undetermined at this order. Therefore the analysis continues to the next order, $\mathcal{O}(\varepsilon^{3})$ . The corresponding equation to be solved is

[TABLE]

where $\mathcal{J}_{\mathcal{F}}[u_{s}]u_{1}\equiv D_{u_{1}}\mathcal{F}[u_{s}]$ and the operator $\mathcal{T}[u]$ is the trilinear third-order differential operator defined by its action on $v_{1},v_{2},v_{3}\in\mathcal{U}$ and using (11) by

[TABLE]

At this order the products formed on the right-hand side of (14) will result in time-dependent terms proportional to $\mbox{exp}(\mathrm{i}\omega_{c}t_{0})$ and these are resonant with homogeneous solutions of the operator $\mathcal{L}$ . Therefore a solvability condition is invoked. In linear operator theory, by the Fredholm alternative (see, for example, boyce1969elementary ), a linear system of the form $\mathcal{L}u=f$ will either have (a) a unique solution $u$ , or (b) a non-trivial solution to $\mathcal{L}^{\dagger}v=0$ . For condition (a) to hold $\langle v^{\dagger},f\rangle=0$ . The dagger superscripts indicate the adjoint problem. The adjoint problem and inner product on the Hilbert space, $U$ , are defined by

[TABLE]

for any functions $u$ and $v$ belonging to the solution function space, where $\Omega$ is the domain of definition of the spatial variable $x$ , and $T$ is the period of oscillation in the fast time scale $t_{0}$ . Applying the Fredholm alternative to (14) yields a solvability condition which provides a constraint on the amplitude function $A(t_{1})$ . This can be written as

[TABLE]

which is the weakly nonlinear Landau equation. The coefficients, $\hat{\nu}$ , $\hat{\lambda}$ and $\hat{\mu}$ are defined by

[TABLE]

where $g^{\dagger}$ is the eigenfunction corresponding to the adjoint problem, $\mathcal{L}^{\dagger}g^{\dagger}=0$ . The values of $\Lambda_{k}$ and $\Gamma_{k}$ are defined as

[TABLE]

The values of $\hat{\nu},\hat{\lambda},\hat{\mu}$ can be obtained once the functions $\varphi_{i}$ and the eigenmodes $g$ and $g^{\dagger}$ have been calculated. For convenience the eigenmodes $g$ and $g^{\dagger}$ are normalised so that $|\langle\mathcal{M}g,g^{\dagger}\rangle|=1$ . Therefore, $|\nu|=1$ . Furthermore it is more convenient to write the Landau equation as

[TABLE]

where $\lambda=-\hat{\lambda}/\hat{\nu}$ and $\mu=-\hat{\mu}/\hat{\nu}$ .

The analysis up to this point has been based on a spatially and temporally continuous set of equations. We now describe a procedure, independent of the choice of discretisation, appropriate for a spatially discretised version of the equations in (3). In the description of the algorithm that follows, J and M are the matrix representations of $\mathcal{J}$ and $\mathcal{M}$ respectively and bold symbols are the discretised vectors of their continuous counterpart. It is assumed that efficient numerical linear solvers (for example, SuperLu li2005overview ), generalised eigensolvers (for example, Trilinos herouxtrilnos ). and continuation algorithms are available. We note that for high-dimensional systems, these calculations are computationally demanding, but they need to be performed only once for each bifurcation point. The procedure is as follows:

Find a Hopf bifurcation as $\beta$ is varied by a continuation method as described by, for example, kuznetsov2013elements . 2. 2.

Calculate the solution of (6), $\textbf{u}_{s}$ , using Newton’s method at $\beta=\beta_{c}$ . 3. 3.

Calculate the eigenfunction g and eigenvalue $s_{c}=\mathrm{i}\omega_{c}$ using an eigensolver to solve $\textbf{J}[\textbf{u}_{s}]\textbf{g}=s\textbf{M}[\textbf{u}_{s}]\textbf{g}$ . 4. 4.

Calculate the Hessian products on the right hand side of (13). If an exact expression for the Jacobian matrix is available, the Hessian products can be numerically computed using the central-difference formula

[TABLE]

where $k_{1}$ is a small finite difference parameter. 5. 5.

Calculate the functions $\bm{\varphi}_{0},\bm{\varphi}_{1}$ and $\bm{\varphi}_{2}$ in (13) using a linear solver. 6. 6.

Calculate the adjoint eigenfunction, $g^{\dagger}$ by solving $\textbf{J}^{T}[\textbf{u}_{s}]\textbf{g}^{\dagger}=-s\textbf{M}^{T}[\textbf{u}_{s}]\textbf{g}^{\dagger}$ . 7. 7.

Calculate the Hessian and third-order products in (18) using (20) by approximating the third-order operator (again assuming an exact Jacobian) as

[TABLE] 8. 8.

Finally combine the quantities in (18) to calculate the values of $\hat{\nu},\hat{\lambda}$ and $\hat{\mu}$ given in (17), defining the discrete inner product as $\langle\textbf{f},\textbf{g}\rangle=\textbf{f}\textbf{g}^{T}$ .

We now return our focus to the specific equations in (1). These are discretised by the finite-element method, utilising the open-source software oomph-lib heil2006oomph in which the functions required for the above procedure have all been implemented. It is convenient to choose the independent parameter $\beta=Q^{-1}$ and, to ensure the finite difference formulas stated above converge, $k_{1}=10^{-4}$ and $k_{2}=10^{-3}$ . When $V=\pi r^{2}$ , with $r=0.46$ , $\alpha=40$ , $w=0.25$ and $h=0.024$ , this procedure was implemented and the critical values of $Q_{c}$ , $\omega_{c}$ and the Landau coefficients were found to 3 significant figures to be

[TABLE]

We have checked that the quadrant in the complex plane of which $\lambda$ and $\mu$ lie (setting $\delta$ constant) does not change as parameters are varied thus not affecting the criticality of the bifurcation.

At this point, we wish to analyse solutions to (19). We note that a decomposition into modulus argument form, i.e. $A=\hat{r}(t)\exp(\mathrm{i}\theta)$ converts the equation into a pair of ODEs:

[TABLE]

where $a=\delta\Re(\lambda)$ , $b=\Re(\mu)$ , $c=\delta\Im(\lambda)$ and $d=\Im(\mu)$ . The invariant solutions occur when $\mbox{d}{\hat{r}}/\mbox{d}{t_{1}}=0$ . The solutions corresponding to $\hat{r}=0$ are the steady states of the fully nonlinear steady equations, while the solution corresponding to $\hat{r}=\sqrt{-a/b}$ is an invariant periodic orbit of the fully nonlinear system. By examining the signs of $a$ and $b$ at each Hopf point and performing a stability analysis on the equations in (23) it is straightforward to show that H1 is subcritical; unstable periodic orbits exist when $Q>Q_{1}$ for $Q$ near to $Q_{1}$ . Similarly, H2 is supercritical so stable periodic orbits exist when $Q>Q_{2}$ .

Near the H1 and H2 points the steady state solution of (3), corresponding to $\hat{r}=0$ in (23), can be written in terms of the state variables as

[TABLE]

This weakly nonlinear approximation can be compared to the fully nonlinear steady state bifurcation structure of § IV. The dotted lines in the inset panels of figure 8 are the approximations given by (24). The approximation near the H2 bifurcation point is excellent in this projection; the fully nonlinear and weakly nonlinear curves are visually indistinguishable. For the H1 bifurcation point the weakly nonlinear steady state approximation deviates significantly away from the solution branch as we move away from H1 in the $(Q,U_{b})$ projection. This difference is likely due to the existence of the fold, F2, in the immediate vicinity of H1. As $h$ is varied, to approach the fold-Hopf point, a weakly nonlinear analysis should include the eigenmode associated with the zero eigenvalue and hence the form of the perturbation equation in (19) would differ.

The periodic orbits, corresponding to $\hat{r}=\sqrt{-a/b}$ , can be stated as

[TABLE]

where $\theta=\theta_{0}$ at $t=0$ and $\Delta=\theta_{0}+\varepsilon^{2}(c-ad/b)t_{0}$ represents the $\mathcal{O}(\varepsilon^{2})$ correction to the phase of the oscillations. The location of the centre of oscillations (ignoring time-dependent terms in (25)) of these periodic orbits are shown as dotted lines in the main panel of figure 8. The circular markers in this diagram indicate the amplitude of the periodic orbits. The initial phase $\theta_{0}$ will not alter the amplitude or period of the periodic orbit.

An example of the weakly nonlinear bubble shapes of the periodic orbits near the H1 point is shown in figure 9. The oscillations are characterised by an oscillation about the symmetric steady solution and the oscillations are visible on both sides of the bubble. It can be seen that after half a period the bubble has ‘flipped’ from its initial shape. The weakly nonlinear equations in (23) predict that if the initial amplitude, $r_{0}$ , of a perturbation is smaller than a critical amplitude $\hat{r}$ , then the system will evolve towards $r=0$ , i.e. the S3 steady state. In contrast if $r_{0}>\hat{r}$ then the trajectory will diverge to $r\to\infty$ . The H1 periodic orbit can therefore be considered an edge state of the reduced system as it forms the boundary between stable and unstable behaviour.

The orbits near H1 are unstable so as noted in § IV we are unable to precisely capture the fully nonlinear periodic orbit by time-integration alone. However by choosing the weakly nonlinear solution in (25) as an initial condition to the fully nonlinear time-dependent problem a qualitative comparison can be made. A characteristic feature of the unstable periodic orbits is that the perturbation oscillates all the way around the edge of the bubble (see figure 9). These oscillations will cause oscillations in the centroid of the bubble. Therefore a useful way to compare the fully nonlinear simulations and weakly nonlinear approximations is to measure the time signal of $\overline{y}$ so that the amplitude and period can be compared. A result of one such simulation when $\varepsilon=0.01$ is shown in figure 10 in the fully nonlinear simulations (a), and weakly nonlinear simulations, (b). The period and amplitude of oscillations of both regimes are nearly identical in both regimes. At this value of $\epsilon$ , we are very close to the H1 point and hence the oscillations in the fully nonlinear regime will take a long time to decay. For larger values of $\varepsilon$ we find that depending on the initial phase of the initial condition, the bubble will either break up, evolve to the AS1 solution or continue to oscillate towards the S3 solution, see figure 11. This is further evidence that the H1 orbit is an edge state in the fully nonlinear regime. In contrast to the weakly nonlinear regime there are three possible outcomes with the fully nonlinear H1 periodic edge state forming the dividing line between stable symmetric evolution and either breaking up or being attracted towards the stable asymmetric states.

An example of the weakly nonlinear bubble shapes of the periodic orbits near the H2 point is shown in figure 12; now the underlying steady state is asymmetric. The stable periodic orbit near the H2 point is characterised by oscillations appearing on the single side of the bubble that is over the edge of the depth-perturbation. These periodic orbits are linearly stable and hence time-integration of the fully nonlinear system should give a good indication of the period and amplitude of the orbit as $t$ progresses. As the centroid of the bubble shape does not vary significantly, more appropriate time-signals are $p_{b}$ and $U_{b}$ , shown in figure 13. The period and amplitude of the weakly nonlinear periodic orbits demonstrate excellent agreement with the fully nonlinear time-dependent simulations.

VI Conclusion

We have investigated the transient dynamics and invariant solutions of a finite air bubble propagating in a perturbed Hele-Shaw channel, a model system in which to explore complex transition scenarios in fluid mechanics. In § III, time-dependent simulations confirmed that a number of different modes of propagation are possible, including multi-tipped solutions and oscillatory solutions as the flow-rate increases. These simulations and an analysis of the steady solution space revealed a finite width region of bistability (with respect to steady states). The demarcation of this bistable region by two Hopf bifurcations explains the oscillatory solutions observed from the time simulations. The first Hopf bifurcation, H1, is subcritical, resulting in the emergence of unstable periodic orbits and the second, H2, is supercritical, resulting in stable periodic orbits. An interesting feature is that the nature of the bistability changes as $Q$ increases. Initially ( $Q$ immediately larger than $Q_{1}$ ) a slower symmetric steady state, S3, co-exists with an ever-present asymmetric steady state, AS1. A transition occurs at the pitchfork bifurcation as the stable AS2 asymmetric branch emerges. Finally due to the supercritical Hopf bifurcation, H2, a stable periodic orbit appears as an alternative stable invariant solution to the asymmetric AS1 branch.

The criticality of the Hopf bifurcations was determined and periodic orbits were approximated using a weakly nonlinear stability analysis in § V. The method described here applies to a range of governing equations and is independent of the choice of discretisation. It has been implemented in a general form in the open-source library oomph-lib. The approximate periodic orbits were in good agreement with the fully nonlinear transient simulations. The oscillations of the unstable orbit near H1 manifest themselves as oscillations all the way around the bubble that cause an oscillatory change in the centroid, $\overline{y}$ , of the bubble. For flow rates just above the H2 bifurcation point, stable periodic orbits exist. The oscillations now appear as ‘waves’ on the side of the bubble above the depth-perturbation only, with the lower edge of the bubble ‘resting’ on the lower limits of the depth-perturbation. In this case the offset of the bubble centroid remains non-zero and constant and the underlying steady state is asymmetric. These stable, asymmetric oscillations have similar characteristics to the oscillations found for air-fingers by de2009tube ; pailha2012oscillatory ; hazel2013multiple ; thompson2014multiple and can be viewed as the finite-air bubble analogue. The oscillations are driven by interactions between the edge of the bubble and the edge of the constriction, which is typically only possible for asymmetric bubble shapes. In contrast the unstable periodic orbits near H1 have not been seen before and appear to be a new phenomena in this system. In this case the oscillations may be driven by underdamping in the mechanism that restores a symmetric bubble to the centre of the channel after perturbation; the generic mechanism requires an interaction between bubble deformation and surface tension and is described in franco2017propagation .

The periodic orbits near the H1 bifurcation point, although unstable, are important in understanding the full time-dependent behaviour of the system. Figure 11(a) is intriguing evidence that these unstable periodic orbits are edge states of the system (see, for example kerswell2005recent ; duguet2008transition ). Slight changes in the initial condition can result in either the bubble returning to the steady S2 state or breaking up/evolving to the AS1 solution. The evidence presented here shows that the unstable periodic orbit is on the boundary between these two forms of transient behaviour. We have quantified the edge state by approximating the unstable periodic orbits by a weakly nonlinear approximation but others have calculated these by using edge-tracking techniques, including shear flow kerswell2005recent ; duguet2008transition , and droplet pinch-off, gallino2018edge which rely on two distinct types of transient behaviour being present. An added complication in this case is that are three possible eventual outcomes; the bubble returning to the S2 steady state, breaking up, or evolving to the AS1 symmetric state. We would expect to see evidence of the periodic orbits and the transitions in the bistable region in any experimental results. The periodic orbit near H1, despite being unstable, is therefore highly influential in the underlying dynamics of a bubble and understanding the basin boundary and a full comparison of the sensitivity of bubble break up to initial conditions, especially close to bifurcation points, is the subject of a combined experimental and theoretical study currently underway.

VII Acknowledgments

We acknowledge funding from the Engineering and Physical Sciences Research Council through grant number EP/P026044/1.

Bibliography40

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) Saffman PG, Taylor G. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Phil. Trans. R. Soc. Lond. 245 , 312–329.
2(2) Tabeling P, Zocchi G, Libchaber A. 1987 An experimental study of the Saffman-Taylor instability.. J. Fluid Mech. 177 , 67–82.
3(3) Moore MG, Juel A, Burgess JM, Mc Cormick WD, Swinney HL. 2002 Fluctuations in viscous fingering.. Phys. Rev. E 65 , 030601(R).
4(4) de Lózar A, Heap A, Box F, Hazel AL, Juel A. 2009 Tube geometry can force switchlike transitions in the behavior of propagating bubbles. Phys. Fluids 21 , 101702.
5(5) Pailha M, Hazel AL, Glendinning PA, Juel A. 2012 Oscillatory bubbles induced by geometrical constraint. Phys. Fluids 24 , 021702.
6(6) Hazel AL, Pailha M, Cox SJ, Juel A. 2013 Multiple states of finger propagation in partially occluded tubes. Phys. Fluids 25 , 021702.
7(7) Jisiou M, Dawson G, Thompson A, Mohr S, Fielden P, Hazel A, Juel A. 2014 Geometry-induced oscillations of finite bubbles in microchannels. Procedia IUTAM 11 , 81–88.
8(8) Thompson AB, Juel A, Hazel AL. 2014 Multiple finger propagation modes in Hele-Shaw channels of variable depth. J. Fluid Mech. 746 , 123–164.