1. Introduction
Interface dynamics in a Hele-Shaw cell – where the motion of the viscous fluids is confined to the narrow gap between two closely spaced parallel glass plates – is a problem of considerable interest both from a theoretical standpoint as a moving free boundary problem and from a practical perspective in view of its connection to other important physical systems, such as flows in porous media, dendritic crystal growth and directional solidification (Pelcé Reference Pelcé1988). The Hele-Shaw system is particularly interesting when one fluid is much less viscous than the other and surface tension effects are neglected, for in this case the problem becomes quite tractable mathematically and many steady and time-dependent solutions have been found since the pioneering work by Saffman & Taylor (Reference Saffman and Taylor1958), Saffman (Reference Saffman1959) and Taylor & Saffman (Reference Taylor and Saffman1959). Recently, deep mathematical connections have been discovered between Hele-Shaw flows and other problems in mathematical physics, such as integrable systems, random matrix theory and quantum gravity (Mineev-Weinstein, Wiegmann & Zabrodin Reference Mineev-Weinstein, Wiegmann and Zabrodin2000). There is also a close relation between interface dynamics in a Hele-Shaw cell and an important growth model known as Loewner evolution (Hohlov et al. Reference Hohlov, Howison, Huntingford, Ockendon and Lacey1994; Gubiec & Szymczak Reference Gubiec and Szymczak2008; Zabrodin Reference Zabrodin2009; Durán & Vasconcelos Reference Durán and Vasconcelos2011). Motivated by these findings, interest in Hele-Shaw flows (or Laplacian growth, as it is also known) has grown well beyond its original hydrodynamics context and there is now an extensive literature on the subject and its mathematical ramifications; for an overview, see e.g. the recent monograph by Gustafsson, Teodorescu & Vasil’ev (Reference Gustafsson, Teodorescu and Vasil’ev2014).
Yet, despite all these developments, the Hele-Shaw system continues to surprise and reveal more mathematical structures underlying its dynamics. Recent investigations by Vasconcelos, Marshall & Crowdy (Reference Vasconcelos, Marshall and Crowdy2014) – motivated in part by the problem of interface dynamics in a Hele-Shaw cell – led to the discovery of a new class of special functions (called the secondary Schottky–Klein prime functions) associated with planar multiply connected domains. These functions are particularly useful in tackling potential-theory problems involving multiply connected domains with mixed boundary conditions. One such problem – and the main theme of the present paper – is the motion of bubbles in a Hele-Shaw channel, in which case the velocity potential satisfies Dirichlet boundary conditions on the bubble interfaces and Neumann conditions on the channel walls.
In this paper, the formalism of the secondary prime functions is used to construct exact solutions for the problem of multiple bubbles steadily translating in a Hele-Shaw channel, both for a finite assembly of bubbles and for a periodic stream of bubbles with an arbitrary number of bubbles per period cell. The problem of multiple fingers penetrating into the channel and moving together with an assembly of bubbles is also analysed as a particular case of the multi-bubble solutions (when some of the bubbles become infinitely elongated). In all cases, the solutions are given in terms of a conformal mapping
$z({\it\zeta})$
from a multiply connected circular domain in an auxiliary complex
${\it\zeta}$
-plane to the flow region exterior to the bubbles in the complex
$z$
-plane. This mapping is written as the sum of two analytic functions – corresponding to the complex potentials in the laboratory and co-moving frames – that map the circular domain onto slit domains. Analytical formulae for these slit maps are obtained in terms of the secondary Schottky–Klein (S–K) prime functions, which then allows us to obtain an explicit solution for the desired mapping
$z({\it\zeta})$
.
In the case of a finite assembly of bubbles, a generalised method of images is used at first to construct the relevant complex potentials and then the resulting expressions (containing an infinite product of terms) are recast in terms of the secondary prime functions. This function-theoretic formulation is more advantageous in that not only does it have a firmer mathematical basis (the theory of compact Riemann surfaces and their associated prime functions), but it can also handle more general cases, such as periodic solutions, that are not easily tackled by the heuristic method of images.
Solutions for multiple steady bubbles in a Hele-Shaw channel were first obtained by the author (Vasconcelos Reference Vasconcelos1994, Reference Vasconcelos2001) for the cases when the bubbles either are symmetrical about the channel centreline or have fore-and-aft symmetry. In such cases, the fluid region can be reduced by virtue of symmetry to a simply connected domain, whereby the desired mappings can be constructed via the Schwarz–Christoffel formula. Solutions for an arbitrary number of steady bubbles in an unbounded cell were obtained by Crowdy (Reference Crowdy2009a ) in terms of the (primary) Schottky–Klein prime function. Crowdy (Reference Crowdy2009b ) also considered the case of a finite assembly of steady bubbles in a channel with no assumed symmetry, but it was subsequently found that the second family of prime functions used in his approach was not correctly defined (D. G. Crowdy, 2011, personal communication). Exact solutions for this problem were later obtained by Green & Vasconcelos (Reference Green and Vasconcelos2014) using an alternative method based on the generalised Schwarz–Christoffel mapping for multiply connected domains. Their solution is expressed in the form of an indefinite integral whose integrand consists of products of (primary) S–K prime functions and which contains several accessory parameters that need to be determined numerically. The solutions reported here for multiple Hele-Shaw bubbles in a channel are based on an entirely different approach and have the advantage that they are given by an explicit analytical formula in terms of the secondary prime functions, with no accessory parameters whatsoever. Furthermore, they can be used to generate new solutions for multiple fingers moving together with an assembly of bubbles, as will be seen later.
In the case of a periodic array of bubbles, solutions were first obtained by Burgess & Tanveer (Reference Burgess and Tanveer1991) for a single stream of symmetrical bubbles. This class of solutions was later extended to include the case of multiple bubbles per period cell under certain symmetry assumptions (Vasconcelos Reference Vasconcelos1994; Silva & Vasconcelos Reference Silva and Vasconcelos2011, Reference Silva and Vasconcelos2013). The new family of periodic solutions reported here is much more general in that it describes an arbitrary stream of groups of bubbles, with no symmetry restriction on the geometrical arrangement of the bubbles within a period cell. Here again the solutions are given in analytical form in terms of the secondary prime functions, making the computation of the bubble shapes a rather simple task, once the preimage domain in the
${\it\zeta}$
-plane is specified.
In light of the existence of this large class of exact solutions for multiple bubbles, it can be argued that the variety of forms observed by Maxworthy (Reference Maxworthy1986), in his experiments on bubbles rising in an inclined Hele-Shaw cell, is in part a manifestation of this multitude of solutions. Further studies would of course be required for a more direct comparison between theory and experiments, but it is worth pointing out that a good agreement was already obtained for the case of a small bubble at the nose of a larger bubble (Ikeda & Maxworthy Reference Ikeda and Maxworthy1990; Vasconcelos Reference Vasconcelos2000). It is to be noted, however, that not all exact solutions reported here are expected to have experimental counterparts, since they correspond to an idealised model where surface tension and three-dimensional thin-film effects are neglected.
It is also important to emphasise that obtaining analytical solutions for steady Hele-Shaw flows naturally paves the way to finding time-dependent solutions, as the form of the steady solutions often suggests an ansatz for the time-dependent ones (see e.g. Bensimon et al. Reference Bensimon, Kadanoff, Liang, Shraiman and Tang1986). In this context, the problem of finding exact solutions for steady Hele-Shaw flows is of physical and mathematical interest not only on its own merit but also because it serves as a starting point to study more general interfacial problems. (The possibility of extending the steady solutions reported herein to the time-dependent case will be briefly discussed towards the end of the paper.)
The analysis presented here for Hele-Shaw bubbles might also find applications in other related problems, such as hollow vortices and streamer discharges in a strong electric field. A hollow vortex is a vortex whose fluid in the interior is at rest (and hence at constant pressure), and so it can be viewed as a bubble with non-zero circulation. The formalism of the primary S–K prime functions has been used to find solutions for a pair of translating hollow vortices (Crowdy, Llewellyn Smith & Freilich Reference Crowdy, Llewellyn Smith and Freilich2013) as well as for a von Kármán street of hollow vortices (Crowdy & Green Reference Crowdy and Green2011). It is thus possible that the present method of analysis, involving multiple bubbles, may be adapted to study more general configurations of hollow vortices. In the case of steady streamers in strong electric fields (Meulenbroek, Rocco & Ebert Reference Meulenbroek, Rocco and Ebert2004), the governing equations are identical to those for Hele-Shaw flows, with a streamer corresponding to a bubble or finger, and so the solutions presented here are likely to be relevant for the problem of multiple interacting streamers (Khalid, Brau & Ebert Reference Luque, Brau and Ebert2008).
The paper is organised as follows. In § 2, the problem of an assembly of a finite number of steady bubbles in a Hele-Shaw channel is formulated in terms of a conformal mapping
$z({\it\zeta})$
from a circular domain in an auxiliary complex
${\it\zeta}$
-plane to the fluid region in the physical
$z$
-plane. The Schottky groups associated with this circular domain and their corresponding Schottky–Klein prime functions are discussed in § 3. The formalism of the secondary prime functions is then used, in § 4, to obtain an analytical solution for the mapping
$z({\it\zeta})$
. For ease of presentation, here the problem is first solved by the method of images and then the results are recast in terms of the prime functions. Configurations with multiple fingers moving together with a group of bubbles are discussed in § 5, as a particular case of the general solution for an assembly of bubbles. In § 6, the case of a periodic array of steady bubbles in a Hele-Shaw channel is considered. Here a fully fledged function-theoretic approach is used to construct an explicit solution for the corresponding mapping
$z({\it\zeta})$
in terms of the secondary prime functions. We conclude the paper by briefly discussing, in § 7, the main features of our results as well as other possible applications of the analysis presented herein.
2. Finite assembly of bubbles: problem formulation
Here we consider the problem of an assembly of a number
$M$
of bubbles of a fluid of negligible viscosity translating uniformly with speed
$U$
parallel to the
$x$
axis in a horizontal Hele-Shaw channel filled with a viscous fluid; see figure 1(a) for a schematic for the case
$M=2$
. The fluid velocity at infinity in front of and behind the bubbles is denoted by
$V$
, where
$U>V$
. To avoid a proliferation of factors of
${\rm\pi}$
in our expressions, it is assumed that the channel has width equal to
${\rm\pi}$
; and without loss of generality the far-field velocity is set to unity, i.e.
$V=1$
. It is also assumed that the pressure inside each bubble is constant and surface tension effects are neglected, so that the viscous fluid pressure on each bubble boundary is taken to be constant (i.e. equal to the pressure inside the bubble).
The flow domains for a finite assembly of bubbles in: (a) the physical
$z$
-plane; (b) the plane of the complex potential in the laboratory frame; (c) the plane of the complex potential in the moving frame; and (d) the auxiliary complex
${\it\zeta}$
-plane.
2.1. The complex potentials
As is well known (Howison Reference Howison1992), the motion of a viscous fluid in a Hele-Shaw cell can be described by a complex potential
$$\begin{eqnarray}F(z,t)={\it\phi}(x,y,t)+\text{i}{\it\psi}(x,y,t)\end{eqnarray}$$
where
$z=x+\text{i}y$
and the velocity potential
${\it\phi}(x,y,t)$
is given by Darcy’s law:
$$\begin{eqnarray}{\it\phi}(x,y,t)=-\frac{b^{2}}{12{\it\mu}}p(x,y,t).\end{eqnarray}$$
Here
$b$
is the cell gap,
${\it\mu}$
is the fluid viscosity,
$p(x,y,t)$
is the pressure and
${\it\psi}(x,y,t)$
is the stream function conjugated to
${\it\phi}(x,y,t)$
.
Since it is assumed that the bubbles all move with a constant speed
$U$
, the complex potential
$F(z,t)$
has a trivial dependence on the time
$t$
, namely
$$\begin{eqnarray}F(z,t)=w(z-Ut),\end{eqnarray}$$
for some function
$w(z)$
to be determined. Note that
$w(z)$
represents the complex potential in the laboratory frame expressed in terms of the coordinates of a frame of reference co-travelling with the bubbles. It is therefore convenient to introduce a second complex potential,
${\it\tau}(z)$
, defined by
$$\begin{eqnarray}{\it\tau}(z)=w(z)-Uz,\end{eqnarray}$$
which describes the flow in the frame of reference co-moving with the bubbles. Henceforth, it will be understood that the variable
$z=x+\text{i}y$
labels points in the co-moving frame.
Now let
$D_{z}$
denote the fluid region in the
$z$
-plane exterior to the bubbles, and denote by
$\partial D_{j}$
, for
$j=1,\ldots ,M$
, the boundary of the
$j$
th bubble; see figure 1(a). The complex potential
$w(z)$
must be analytic in
$D_{z}$
and satisfy the following boundary conditions:
$$\begin{eqnarray}\displaystyle & \text{Im}[w(z)]=0\quad \text{on }y=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{Im}[w(z)]={\rm\pi}\quad \text{on }y={\rm\pi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{Re}[w(z)]=c_{j}\quad \text{for }z\in \partial D_{j}, & \displaystyle\end{eqnarray}$$
$\{c_{j}\mid j=1,\ldots ,M\}$
is a set of real constants. Conditions (2.5) and (2.6) simply state that the channel walls at
$y=0$
and
$y={\rm\pi}$
are streamlines of the flow, whereas (2.7) follows from the fact that the pressure is constant on each bubble boundary (i.e. the constants
$c_{j}$
are related to the values of the pressure in each bubble). Furthermore, as the far-field flow is uniform with unity velocity, it follows that
$w(z)$
has the following asymptotic behaviour at infinity: 
$$\begin{eqnarray}w(z)\sim z\quad \text{for }|x|\rightarrow \infty .\end{eqnarray}$$
From (2.5)–(2.8) one then concludes that the flow domain,
$D_{w}$
, in the
$w$
-plane is a horizontal strip of width
${\rm\pi}$
containing
$M$
vertical slits in its interior, where each slit corresponds to a bubble in the
$z$
-plane; see figure 1(b).
Let us now consider the boundary conditions satisfied by
${\it\tau}(z)$
. As the bubble boundaries
$\partial D_{j}$
are streamlines of the flow in the co-moving frame, one has that
$$\begin{eqnarray}\text{Im}[{\it\tau}(z)]=d_{j}\quad \text{for }z\in \partial D_{j},\end{eqnarray}$$
where
$\{d_{j}\mid j=1,\ldots ,M\}$
is a set of real constants (corresponding to the values of the streamfunction in the moving frame on each bubble boundary). Furthermore, conditions (2.5), (2.6) and (2.8) imply, in view of (2.4), that
$$\begin{eqnarray}\displaystyle & \text{Im}[{\it\tau}(z)]=0\quad \text{on }y=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{Im}[{\it\tau}(z)]=(1-U){\rm\pi}\quad \text{on }y={\rm\pi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}{\it\tau}(z)\sim (1-U)z\quad \text{for }|x|\rightarrow \infty .\end{eqnarray}$$
From (2.9)–(2.12) one readily sees that the flow domain,
$D_{{\it\tau}}$
, in the
${\it\tau}$
-plane is a strip of width
$(U-1){\rm\pi}$
, with
$M$
horizontal slits in its interior; see figure 1(c).
2.2. Conformal mapping
We shall seek a solution for the free boundary problem defined in § 2.1 in terms of a conformal mapping
$z({\it\zeta})$
from a bounded circular domain
$D_{{\it\zeta}}$
in an auxiliary complex
${\it\zeta}$
-plane to the fluid region
$D_{z}$
. To be specific, let
$D_{{\it\zeta}}$
be the domain obtained from the unit disk by excising
$M$
non-overlapping smaller disks. A schematic of
$D_{{\it\zeta}}$
is shown in figure 1(d) for the triply connected case
$M=2$
.
Label the unit circle by
$C_{0}$
and the
$M$
inner circular boundaries by
$C_{1},\ldots ,C_{M}$
; and let
${\it\delta}_{j}$
and
$q_{j}$
denote respectively the centres and radii of the circles
$C_{j}$
,
$j=0,\ldots ,M$
. (Note that
${\it\delta}_{0}=0$
and
$q_{0}=1$
.) The mapping
$z({\it\zeta})$
is chosen such that the unit circle
$C_{0}$
maps to the channel walls, whilst the inner circles
$C_{1},\ldots ,C_{M}$
map to the bubble boundaries
$\partial D_{1},\ldots ,\partial D_{M}$
. This implies that
$z({\it\zeta})$
will necessarily have two logarithmic singularities, denoted by
${\it\zeta}_{1}$
and
${\it\zeta}_{2}$
, on the unit circle, which map in the
$z$
-plane to the two end points of the channel,
$x=\mp \infty$
, respectively. By the degrees of freedom afforded by the Riemann–Koebe mapping theorem (Goluzin Reference Goluzin1969), we can set
${\it\zeta}_{1}=-1$
and
${\it\zeta}_{2}=1$
. With this choice, the upper unit semicircle,
$C_{0}^{+}$
, maps to the upper channel wall (
$y={\rm\pi}$
) and the lower unit semicircle,
$C_{0}^{-}$
, maps to the lower wall (
$y=0$
); see figure 1.
The extended flow domain
$F_{0}$
in the auxiliary complex
${\it\zeta}$
-plane.
It is expedient to augment the flow region by reflecting the original channel in its lower wall, thus generating an extended channel defined by
$D_{z}\cup \overline{D_{z}}$
, where a bar denotes complex conjugation. (Note that this extended channel contains
$2M$
bubbles where each bubble in the lower half-channel is the mirror image of a corresponding bubble in the upper half-channel.) Accordingly, the extended flow domain, denoted by
$F_{0}$
, in the auxiliary
${\it\zeta}$
-plane is obtained by adding to
$D_{{\it\zeta}}$
its reflection in
$C_{0}$
:
$$\begin{eqnarray}F_{0}=D_{{\it\zeta}}\cup {\it\varphi}_{0}(D_{{\it\zeta}}),\end{eqnarray}$$
where
${\it\varphi}_{0}({\it\zeta})=1/\overline{{\it\zeta}}$
defines reflection in
$C_{0}$
. In addition, a branch cut must be inserted along
$C_{0}^{+}$
so that the lower and upper sides of this cut map to the upper and lower walls of the extended channel,
$y=\pm {\rm\pi}$
, respectively. A schematic of
$F_{0}$
is shown in figure 2 for the case
$M=2$
.
Let us denote by
$C_{-j}$
the reflection of the circle
$C_{j}$
in
$C_{0}$
, i.e.
$C_{-j}={\it\varphi}_{0}(C_{j})$
. The region
$F_{0}$
defined in (2.13) thus corresponds to the exterior of the circles
$C_{j}$
and
$C_{-j}$
, for
$j=1,\ldots ,M$
. For convenience of notation, define the following set of labels
$$\begin{eqnarray}K_{M}\equiv \{-M,\ldots ,-1,1,\ldots ,M\},\end{eqnarray}$$
and let
$\mathscr{C}_{M}$
denote the set of circles bounding
$F_{0}$
, that is,
$$\begin{eqnarray}\mathscr{C}_{M}\equiv \{C_{k}\mid k\in K_{M}\}.\end{eqnarray}$$
Next, introduce the functions
$W({\it\zeta})$
and
$T({\it\zeta})$
through the following compositions:
$$\begin{eqnarray}\displaystyle & W({\it\zeta})\equiv w(z({\it\zeta})), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & T({\it\zeta})\equiv {\it\tau}(z({\it\zeta})). & \displaystyle\end{eqnarray}$$
$w=W({\it\zeta})$
conformally maps
$F_{0}$
onto a slit strip domain in the
$w$
-plane defined by
$D_{w}\cup \overline{D_{w}}$
. Similarly, the function
${\it\tau}=T({\it\zeta})$
maps
$F_{0}$
onto the slit strip domain
$D_{{\it\tau}}\cup \overline{D_{{\it\tau}}}$
in the
${\it\tau}$
-plane. Note, in particular, that both
$W({\it\zeta})$
and
$T({\it\zeta})$
must have logarithmic singularities at
${\it\zeta}=\pm 1$
. These singularities act as a sink and a source, respectively, for the flows generated by these complex potentials – a fact that will be exploited in § 4 to compute explicit formulae for
$W({\it\zeta})$
and
$T({\it\zeta})$
via the method of images. Once these functions are known, the desired mapping function
$z({\it\zeta})$
that describes the bubble shapes is then given by 
$$\begin{eqnarray}z({\it\zeta})=\frac{1}{U}[W({\it\zeta})-T({\it\zeta})],\end{eqnarray}$$
as follows from (2.4).
3. Schottky groups and the Schottky–Klein prime function
For any circular region
$F_{0}$
as defined in (2.13), one can define a so-called Schottky group generated by the Möbius transformations that ‘pair’ the circles
$C_{-j}$
and
$C_{j}$
. Associated with this Schottky group and its subgroups there can be defined special transcendental functions, called primary and secondary Schottky–Klein prime functions. These functions naturally appear in the context of Hele-Shaw flows with multiples bubbles, as will be seen in §§ 4–6, and so it was thought desirable to present here a brief introduction to the S–K prime functions.
3.1. The primary S–K prime function
Consider the circular domain
$F_{0}$
defined in (2.13). For
$k=-M,\ldots ,M$
, denote by
${\it\varphi}_{k}({\it\zeta})$
the reflection map in the circle
$C_{k}$
which is defined by
$$\begin{eqnarray}{\it\varphi}_{k}({\it\zeta})={\it\delta}_{k}+\frac{q_{k}^{2}}{\overline{{\it\zeta}}-\overline{{\it\delta}}_{k}}.\end{eqnarray}$$
Now, introduce the following Möbius maps:
$$\begin{eqnarray}{\it\theta}_{k}({\it\zeta})\equiv {\it\varphi}_{k}(1/\overline{{\it\zeta}})={\it\delta}_{k}+\frac{q_{k}^{2}{\it\zeta}}{1-\overline{{\it\delta}}_{k}{\it\zeta}}.\end{eqnarray}$$
Note that
${\it\theta}_{k}({\it\zeta})$
consists of the composition of a reflection in
$C_{0}$
followed by a reflection in
$C_{k}$
, i.e.
${\it\theta}_{k}={\it\varphi}_{k}{\it\varphi}_{0}$
. Alternatively, reflection in
$C_{k}$
can be expressed in terms of
${\it\theta}_{k}$
as
${\it\varphi}_{k}={\it\theta}_{k}{\it\varphi}_{0}$
, or more explicitly:
$$\begin{eqnarray}{\it\varphi}_{k}({\it\zeta})={\it\theta}_{k}(1/\overline{{\it\zeta}}).\end{eqnarray}$$
One important property of the maps
${\it\theta}_{k}({\it\zeta})$
is the following relation:
$$\begin{eqnarray}{\it\theta}_{-k}({\it\zeta})={\it\theta}_{k}^{-1}({\it\zeta})=\frac{1}{\overline{{\it\theta}_{k}(1/\bar{{\it\zeta}})}},\end{eqnarray}$$
where
${\it\theta}_{k}^{-1}({\it\zeta})$
denotes the inverse of
${\it\theta}_{k}({\it\zeta})$
. This relation can be derived using geometrical arguments, or it can be verified directly.
Now it may be verified that
${\it\theta}_{j}({\it\zeta})$
, for
$j=1,\ldots ,M$
, maps the interior of
$C_{-j}$
onto the exterior of
$C_{j}$
. Conversely, the map
${\it\theta}_{-j}({\it\zeta})$
maps the exterior of
$C_{j}$
onto the interior of
$C_{-j}$
. The set
${\it\Theta}_{0}$
consisting of all compositions of the maps
${\it\theta}_{k}({\it\zeta})$
,
$k\in K_{M}$
, defines what is called a classical Schottky group. The region
$F_{0}$
(which we recall consists of the exterior of the circles
$\mathscr{C}_{M}$
) is called a fundamental region of the group
${\it\Theta}_{0}$
and the maps
$\{{\it\theta}_{j}\mid j=1,\ldots ,M\}$
are called the fundamental generators of the group. For any given Schottky group
${\it\Theta}_{0}$
and fundamental region
$F_{0}$
, the S–K prime function,
${\it\omega}({\it\zeta},{\it\alpha})$
, can be defined for any two points
${\it\zeta},{\it\alpha}\in F_{0}$
.
The S–K prime function admits the following infinite product representation (Baker Reference Baker1897):
$$\begin{eqnarray}{\it\omega}({\it\zeta},{\it\alpha})=({\it\zeta}-{\it\alpha})\mathop{\prod }_{{\it\theta}\in {\it\Theta}_{0}^{\prime \prime }}\frac{({\it\zeta}-{\it\theta}({\it\alpha}))({\it\alpha}-{\it\theta}({\it\zeta}))}{({\it\zeta}-{\it\theta}({\it\zeta}))({\it\alpha}-{\it\theta}({\it\alpha}))}.\end{eqnarray}$$
where
${\it\Theta}_{0}^{\prime \prime }\subset {\it\Theta}_{0}$
is the set such that for all
${\it\theta}\in {\it\Theta}_{0}$
, excluding the identity, either
${\it\theta}$
or
${\it\theta}^{-1}$
(but not both) is contained in
${\it\Theta}_{0}^{\prime \prime }$
. For example, if
${\it\theta}_{1}{\it\theta}_{-2}$
is included in
${\it\Theta}_{0}^{\prime \prime }$
, then
${\it\theta}_{2}{\it\theta}_{-1}$
must be excluded.
For later use, it is convenient to quote here the following relation:
$$\begin{eqnarray}\frac{{\it\omega}({\it\zeta},{\it\alpha})}{{\it\omega}({\it\zeta},{\it\gamma})}=C({\it\alpha},{\it\gamma})\mathop{\prod }_{{\it\theta}\in {\it\Theta}_{0}}\frac{{\it\zeta}-{\it\theta}({\it\alpha})}{{\it\zeta}-{\it\theta}({\it\gamma})},\end{eqnarray}$$
where the prefactor
$C({\it\alpha},{\it\gamma})$
depends on the values of
${\it\alpha}$
and
${\it\gamma}$
(but not on the point
${\it\zeta}$
). Note in particular that the product in (3.6) is over the entire group
${\it\Theta}_{0}$
. For a derivation of this formula, see, e.g. DeLillo & Kropf (Reference DeLillo and Kropf2010); it is also implied by an alternative representation of
${\it\omega}({\it\zeta},{\it\alpha})$
given by Marshall (Reference Marshall2005).
The S–K prime function is intimately connected with the theory of compact Riemann surfaces (Fay Reference Fay2008), but for the present purposes it suffices to think of it as a special computable function. In this context, a few remarks on the numerical computation of the S–K prime functions are in order. The infinite products (3.5) and (3.6) are known to converge if the circular boundaries of
$D_{{\it\zeta}}$
are either centred on the real axis or are sufficiently well separated (Baker Reference Baker1897; DeLillo & Kropf Reference DeLillo and Kropf2010). But even when this product formula converges it may be inefficient for numerical computation. An alternative numerical scheme that relies on a representation of the prime function in terms of much more rapidly convergent Fourier–Laurent sums has been developed by Crowdy & Marshall (Reference Crowdy and Marshall2007). This scheme (when applicable) will be used in our numerical computations; see below.
3.2. Secondary S–K prime functions
Given a Schottky group
${\it\Theta}_{0}$
as defined above, a family of Schottky subgroups
${\it\Theta}_{N}\subset {\it\Theta}_{0}$
, for
$N=1,\ldots ,M$
, can be defined, and prime functions can naturally be associated to them. These so-called secondary prime functions were introduced by Vasconcelos et al. (Reference Vasconcelos, Marshall and Crowdy2014), as building blocks for constructing conformal mappings for mixed slit domains, and are briefly reviewed here.
For a given
$N$
, with
$1\leqslant N\leqslant M$
, define
${\it\Theta}_{N}$
to be the set of all elements
${\it\theta}\in {\it\Theta}_{0}$
such that
${\it\theta}$
contains only combinations of an even number of the maps
${\it\theta}_{k}$
,
$k=-N,\ldots ,N$
,
$k\neq 0$
, but may contain any number (even or odd) of the other maps
${\it\theta}_{k}$
, i.e. for
$N+1\leqslant |k|\leqslant M$
. For example, for the case
$M=3$
and
$N=2$
, one can show that the group
${\it\Theta}_{N}$
is generated by the following maps:
$\{{\it\theta}_{3},{\it\theta}_{1}^{2},{\it\theta}_{1}{\it\theta}_{2},{\it\theta}_{1}{\it\theta}_{-2},{\it\theta}_{1}{\it\theta}_{3}{\it\theta}_{-1}\}$
; it can indeed be verified that these maps and their inverses generate only (and all) combinations of an even number of the maps
${\it\theta}_{\pm 1}$
and
${\it\theta}_{\pm 2}$
, but any number of the maps
${\it\theta}_{\pm 3}$
can appear.
It may be verified that the set
${\it\Theta}_{N}$
defined above is itself a Schottky group, which is obviously a subgroup of the original group
${\it\Theta}_{0}$
. Associated with the group
${\it\Theta}_{N}$
one can define a corresponding prime function. To avoid confusion with the primary S–K prime function
${\it\omega}({\it\zeta},{\it\alpha})$
introduced in § 4.1 (and associated with the original group
${\it\Theta}_{0}$
), this secondary S–K prime function associated with
${\it\Theta}_{N}$
is denoted by
${\it\Omega}_{N}({\it\zeta},{\it\alpha})$
. This function admits a product representation as in (3.5), with the only difference that the product is over the set
${\it\Theta}_{N}^{\prime \prime }$
whose definition mirrors that of the set
${\it\Theta}_{0}^{\prime \prime }$
.
Now fix an integer
$l$
such that
$1\leqslant l\leqslant N$
. From the definition of
${\it\Theta}_{N}$
, one can verify that any element
${\it\theta}$
of the original group
${\it\Theta}_{0}$
is either an element of
${\it\Theta}_{N}$
or a composite map of the form
${\it\psi}\circ {\it\theta}_{l}$
, for some
${\it\psi}\in {\it\Theta}_{N}$
. For example, for the case
$M=3$
,
$N=2$
and
$l=1$
, we have that
${\it\theta}_{2}\notin {\it\Theta}_{N}$
, but
${\it\theta}_{2}={\it\theta}_{2}({\it\theta}_{-1}{\it\theta}_{1})=({\it\theta}_{2}{\it\theta}_{-1}){\it\theta}_{1}={\it\psi}{\it\theta}_{1}$
, where
${\it\psi}={\it\theta}_{2}{\it\theta}_{-1}\in {\it\Theta}_{N}$
. Using this decomposition property of the group
${\it\Theta}_{0}$
, together with identity (3.6), one can establish the following relation between the primary and secondary prime functions:
$$\begin{eqnarray}\frac{{\it\omega}({\it\zeta},{\it\alpha})}{{\it\omega}({\it\zeta},{\it\gamma})}=\tilde{C}({\it\alpha},{\it\gamma})\frac{{\it\Omega}_{N}({\it\zeta},{\it\alpha}){\it\Omega}_{N}({\it\zeta},{\it\theta}_{l}({\it\alpha}))}{{\it\Omega}_{N}({\it\zeta},{\it\gamma}){\it\Omega}_{N}({\it\zeta},{\it\theta}_{l}({\it\gamma}))},\end{eqnarray}$$
where the prefactor
$\tilde{C}({\it\alpha},{\it\gamma})$
depends only on
${\it\alpha}$
and
${\it\gamma}$
. (This relation holds irrespective of the choice of
$l$
, for
$1\leqslant l\leqslant N$
.)
Note, in particular, that for
$N=M$
the subgroup
${\it\Theta}_{M}$
consists of all even combinations of the maps
${\it\theta}_{k}$
,
$k\in K_{M}$
; here the generators of the group are the Möbius maps
$\{{\it\theta}_{1}{\it\theta}_{k}\mid k\in K_{M}\}$
(excluding the identity) which map the interior of the circles
$C_{-k}$
onto the exterior of their images in
$C_{1}$
. This subgroup and its associated S–K prime function,
${\it\Omega}_{M}({\it\zeta},{\it\alpha})$
, play a crucial role in constructing solutions for a finite assembly of bubbles in a Hele-Shaw channel, as discussed in § 4. The function
${\it\Omega}_{M-1}({\it\zeta},{\it\alpha})$
, on the other hand, appears in the problem of periodic arrays of bubbles to be discussed in § 6.
4. Solutions for a finite assembly of bubbles
In this section the formalism of the S–K prime functions is used to construct explicit formulae for the complex potentials
$W({\it\zeta})$
and
$T({\it\zeta})$
introduced in § 2.2. A crucial step in this task is the computation of the infinite sets of images (associated with the sink and source at
${\it\zeta}=\pm 1$
) that are necessary to enforce the appropriate boundary conditions on the circles
$\mathscr{C}_{M}$
bounding the flow region
$F_{0}$
in the
${\it\zeta}$
-plane. The location of these images can be expressed in terms of the action of one of the groups
${\it\Theta}_{N}$
on the positions of the original source and sink. The specific group required (i.e. the value of
$N$
) depends on whether the flow is described in the laboratory frame or in the co-moving frame. We begin by considering the complex potential
$T({\it\zeta})$
in the co-moving frame.
4.1. The function
$T({\it\zeta})$
Recall that in the co-moving frame each bubble is a streamline of the complex potential
${\it\tau}(z)$
, thus implying that the circles
$\mathscr{C}_{M}$
are streamlines of the flow generated by
$T({\it\zeta})$
in the
${\it\zeta}$
-plane. These boundary conditions can be satisfied with a judicious choice of images, as discussed below.
Consider first the source at
${\it\zeta}=-1$
. Reflection of this source in each of the
$2M$
circles in
$\mathscr{C}_{M}$
yields a set of
$2M$
image sources at the positions
${\it\varphi}_{k}={\it\theta}_{k}(-1)$
,
$k\in K_{M}$
; see (3.3). Here we used the fact that the image of a point source with respect to a streamline circle is a source of the same strength and located at the corresponding reflection point.
Now, for any given image source,
${\it\theta}_{k}(-1)$
, for some
$k\in K_{M}$
, its subsequent reflection in a circle
$C_{k^{\prime }}$
,
$k^{\prime }\neq k$
, yields another source at the point
${\it\varphi}_{k^{\prime }}({\it\theta}_{k}(-1))={\it\theta}_{k^{\prime }}(1/\overline{{\it\theta}_{k}(-1)})={\it\theta}_{k^{\prime }}({\it\theta}_{-k}(-1))$
, where we used (3.3) and (3.4). Generalising this argument, one can show that reflection in the circles
$\mathscr{C}_{M}$
of the first set of images
$\{{\it\theta}_{k}(-1)\mid k\in K_{M}\}$
generates second-level image sources at the points
$\{{\it\theta}_{k_{1}}\circ {\it\theta}_{k_{2}}(-1)\mid k_{1},k_{2}\in K_{M},~k_{1}+k_{2}\neq 0\}$
. More generally, it may be verified that after
$m$
reflections of the original source in
$\mathscr{C}_{M}$
, one obtains a set of sources whose locations are given by
$\{{\it\theta}_{k_{1}}\circ \cdots \circ {\it\theta}_{k_{m}}(-1)\mid k_{j}\in K_{M},~k_{j}+k_{j+1}\neq 0,~j=1,\ldots ,M\}$
. Continuing this procedure ad infinitum, one obtains an infinite set of image sources located at the following points:
$\{{\it\theta}(-1)\mid {\it\theta}\in {\it\Theta}_{0}\}$
. A similar procedure for the sink at
${\it\zeta}=1$
yields a system of image sinks at the points
$\{{\it\theta}(1)\mid {\it\theta}\in {\it\Theta}_{0}\}$
.
The velocity potential
$T({\it\zeta})$
produced by the set of sources and sinks computed above is given by
$$\begin{eqnarray}T({\it\zeta})=(1-U)\log \left(\mathop{\prod }_{{\it\theta}\in {\it\Theta}_{0}}\frac{{\it\zeta}-{\it\theta}(-1)}{{\it\zeta}-{\it\theta}(1)}\right),\end{eqnarray}$$
where the prefactor was determined from the requirement that the logarithmic singularities of
$T({\it\zeta})$
at
${\it\zeta}=\pm 1$
have the appropriate strength. In other words, when going around either one of these singularities (from one side of the cut to the other) the jump in
$T({\it\zeta})$
must equal
$\text{i}2{\rm\pi}(U-1)$
, which corresponds to the width of the extended strip domain in the
${\it\tau}$
-plane.
Using (3.6), one can rewrite (4.1) in terms of the S–K prime functions:
$$\begin{eqnarray}T({\it\zeta})=(1-U)\log \left(\frac{{\it\omega}({\it\zeta},-1)}{{\it\omega}({\it\zeta},1)}\right)+b,\end{eqnarray}$$
where
$b$
is an immaterial constant. Alternative derivations of this formula directly from the properties of the S–K prime function were given by Crowdy (Reference Crowdy2009b
) and Green & Vasconcelos (Reference Green and Vasconcelos2014). The derivation presented above is arguably more intuitive in that it is based on the well-known method of images.
For later use, it is convenient to make use of (3.7) and rewrite (4.2) in terms of the secondary prime functions
${\it\Omega}_{M}({\it\zeta},{\it\alpha})$
:
$$\begin{eqnarray}T({\it\zeta})=(1-U)\log \left(\frac{{\it\Omega}_{M}({\it\zeta},-1){\it\Omega}_{M}({\it\zeta},{\it\theta}_{l}(-1))}{{\it\Omega}_{M}({\it\zeta},1){\it\Omega}_{M}({\it\zeta},{\it\theta}_{l}(1))}\right),\end{eqnarray}$$
where
$1\leqslant l\leqslant M$
and an overall additive constant was omitted since its value is not relevant for the flow.
4.2. The function
$W({\it\zeta})$
Here we start by recalling that the bubbles’ boundaries are equipotentials of the complex potential
$w(z)$
in the laboratory frame, and so the circles
$\mathscr{C}_{M}$
must be equipotentials of the flow described by
$W({\it\zeta})$
. Using a similar approach as in § 4.1, one can readily compute the system of images required to satisfy these boundary conditions. The only difference to bear in mind is that the image of a source with respect to an equipotential circle is a sink, and vice versa.
Consider the source at
${\it\zeta}=-1$
. Its first set of images with respect to reflections in the circles
$\mathscr{C}_{M}$
correspond to
$2M$
sinks at the positions
$\{{\it\theta}_{k}(-1)\mid k\in K_{M}\}$
. A second round of reflections of these sinks in
$\mathscr{C}_{M}$
then yields sources at the points
$\{{\it\theta}_{k_{1}}\circ {\it\theta}_{k_{2}}(-1)\mid k_{1},k_{2}\in K_{M},k_{1}+k_{2}\neq 0\}$
. Upon further reflections of these sinks in
$\mathscr{C}_{M}$
one gets a third-level set of sources, and so on, where at each successive level of reflection sources generate sinks and vice versa. In other words, after a sequence of an even number of reflections of the original source one gets back a source, whereas an odd number of such reflections produces a sink. The system of images associated with the primary source at
${\it\zeta}=-1$
thus consists of the following two infinite sets:
-
(i) sources at the points
$\{{\it\theta}(-1)\mid {\it\theta}\in {\it\Theta}_{M}\}$
; -
(ii) sinks at the points
$\{{\it\theta}\circ {\it\theta}_{l}(-1)\mid {\it\theta}\in {\it\Theta}_{M}\}$
;
where we recall that
${\it\Theta}_{M}$
is the set of all even combinations of the maps
${\it\theta}_{k}$
,
$k\in K_{M}$
, and
$1\leqslant l\leqslant M$
is an arbitrary integer.
Similarly, associated with the sink at
${\it\zeta}=1$
one finds the following system of images:
-
(iii) sinks at the points
$\{{\it\theta}(1)\mid {\it\theta}\in {\it\Theta}_{M}\}$
; -
(iv) sources at the points
$\{{\it\theta}\circ {\it\theta}_{l}(1)\mid {\it\theta}\in {\it\Theta}_{M}\}$
.
Note that in writing down the locations of the images in sets (ii) and (iv) above, use was made of the fact that that any combination of an odd number of the maps
${\it\theta}_{k}$
,
$k\in K_{M}$
, can be written as
${\it\psi}\circ {\it\theta}_{l}$
for some
${\it\psi}\in {\it\Theta}_{M}$
, as discussed in § 3.2.
Given the sets of sources and sinks in (i)–(iv) above, it then follows that the resulting complex potential
$W({\it\zeta})$
is
$$\begin{eqnarray}W({\it\zeta})=\log \left(\mathop{\prod }_{{\it\theta}\in {\it\Theta}_{M}}\frac{({\it\zeta}-{\it\theta}(-1))({\it\zeta}-{\it\theta}({\it\theta}_{l}(1)))}{({\it\zeta}-{\it\theta}(1))({\it\zeta}-{\it\theta}({\it\theta}_{l}(-1)))}\right).\end{eqnarray}$$
where the prefactor (unity) was chosen so that the width of the extended channel in the
$w$
-plane is equal to
$2{\rm\pi}$
. Now using (3.6), one can rewrite (4.4) in terms of the secondary S–K prime functions
${\it\Omega}_{M}({\it\zeta},{\it\alpha})$
:
$$\begin{eqnarray}W({\it\zeta})=\log \left(\frac{{\it\Omega}_{M}({\it\zeta},-1){\it\Omega}_{M}({\it\zeta},{\it\theta}_{l}(1))}{{\it\Omega}_{M}({\it\zeta},1){\it\Omega}_{M}({\it\zeta},{\it\theta}_{l}(-1))}\right),\end{eqnarray}$$
where again an irrelevant additive constant was omitted. Once the complex potentials
$W({\it\zeta})$
and
$T({\it\zeta})$
have been obtained, the mapping function
$z({\it\zeta})$
immediately follows from (2.18), as discussed next.
4.3. The conformal mapping
$z({\it\zeta})$
After substituting (4.3) and (4.5) into (2.18) and performing some simplification, one finds
$$\begin{eqnarray}z({\it\zeta})=\log \left(\frac{{\it\Omega}_{M}({\it\zeta},-1)}{{\it\Omega}_{M}({\it\zeta},1)}\right)+\left(1-\frac{2}{U}\right)\log \left(\frac{{\it\Omega}_{M}({\it\zeta},{\it\theta}_{l}(-1))}{{\it\Omega}_{M}({\it\zeta},{\it\theta}_{l}(1))}\right).\end{eqnarray}$$
Recall that
$1\leqslant l\leqslant M$
is an integer that can be chosen arbitrarily; in specific computations it is convenient to set
$l=1$
.
The coordinates (
$x_{j},y_{j}$
) of each bubble interface
$\partial D_{j}$
,
$j=1,\ldots ,M$
, are thus given in parametric form by
$$\begin{eqnarray}x_{j}(s)+\text{i}y_{j}(s)=z({\it\delta}_{j}+q_{j}\exp (\text{i}s)),\quad 0\leqslant s<2{\rm\pi},\end{eqnarray}$$
with
$z({\it\zeta})$
as in (4.6). Note that all the geometrical information about the bubble configuration described by the solution above is encapsulated in the prescription of the preimage domain
$D_{{\it\zeta}}$
. This domain is characterised by its
$3M$
conformal moduli
$\{{\it\delta}_{j},q_{j}\mid j=1,\ldots ,M\}$
, which correspond physically to the area and centroids of each of the
$M$
bubbles. Thus, once the conformal moduli of
$D_{{\it\zeta}}$
are prescribed, a solution for a specific assembly of bubbles is obtained.
Let us recall here that solutions with
$U=2$
are special in the sense that all other bubble assemblies with different values of
$U$
can be obtained from the
$U=2$
solutions by a simple rescaling of coordinates (possibly followed by a rigid translation). This fact was first noticed by Millar (Reference Millar1992) in the context of the Taylor–Saffman solution for a single bubble and later shown to be also valid for multiple symmetrical bubbles (Vasconcelos Reference Vasconcelos2001). In fact, this property holds in general, i.e. for an assembly of steady bubbles with no assumed symmetry, as can be shown by the following argument, which parallels that presented by Vasconcelos (Reference Vasconcelos2001), albeit with a different notation.
First note that the complex potential
$W({\it\zeta})$
given in (4.5) does not depend explicitly on the velocity
$U$
, whereas the complex potential
$T({\it\zeta})$
in the moving frame depends on
$U$
only via a prefactor; see (4.3). For a given domain
$D_{{\it\zeta}}$
, one can then write
$$\begin{eqnarray}\displaystyle & W_{U}({\it\zeta})=W_{2}({\it\zeta}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & T_{U}({\it\zeta})=(U-1)T_{2}({\it\zeta}), & \displaystyle\end{eqnarray}$$
$U$
. Inserting these equations into (2.18), one can show after some manipulation that the coordinates
$(x_{j}^{(U)},y_{j}^{(U)})$
for the
$j$
th bubble in an assembly with speed
$U$
can be obtained from the coordinates
$(x_{j}^{(2)},y_{j}^{(2)})$
of the corresponding bubble in the assembly with
$U=2$
by the following relations: 
$$\begin{eqnarray}\displaystyle & x_{j}^{(U)}(s)=(1+{\it\rho})x_{j}^{(2)}(s)-{\it\rho}c_{j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & y_{j}^{(U)}(s)=(1-{\it\rho})x_{j}^{(2)}(s)-{\it\rho}d_{j}, & \displaystyle\end{eqnarray}$$
${\it\rho}=1-(2/U)$
. Here
$c_{j}=\text{Re}[W_{2}({\it\zeta})]$
and
$d_{j}=\text{Im}[T_{2}({\it\zeta})]$
for
${\it\zeta}\in C_{j}$
; see (2.7) and (2.9). One then sees that a bubble with
$U>2$
(i.e.
$0<{\it\rho}<1$
) is obtained from the corresponding shape with
$U=2$
by a stretching along the
$x$
direction and a contraction along the
$y$
direction (followed by a rigid translation); whereas for
$1<U<2$
(i.e.
$-1<{\it\rho}<0$
) the opposite occurs: contraction along the
$x$
-axis and stretching along the
$y$
-axis. For this reason we shall restrict ourselves to the case
$U=2$
when computing specific examples of bubble configurations, as discussed next.
4.4. Examples
As already mentioned in the Introduction, solutions for multiple steady bubbles in a Hele-Shaw channel were recently found by Green & Vasconcelos (Reference Green and Vasconcelos2014) in terms of an indefinite integral whose integrand consists of a product of primary S–K prime functions and which contains several accessory parameters that need to be determined numerically. The solution (4.6), in contradistinction, is expressed as an explicit analytical formula in terms of the secondary prime functions
${\it\Omega}_{M}({\it\zeta},{\it\alpha})$
, with no accessory parameters. Given a domain
$D_{{\it\zeta}}$
, the corresponding bubble shapes can be readily obtained upon computation of the relevant prime functions.
Examples of assemblies of three symmetrical bubbles (
$M=3$
). In (a) the bubbles are symmetric about the centreline, whereas in (b) the bubbles have fore-and-aft symmetry. The conformal moduli of
$D_{{\it\zeta}}$
are as follows:
${\it\delta}_{1}=0.4$
,
${\it\delta}_{2}=0.0$
,
${\it\delta}_{3}=-0.6$
,
$q_{1}=0.1$
,
$q_{2}=0.1$
,
$q_{3}=0.3$
in (a); and
${\it\delta}_{1}=\text{i}0.7$
,
${\it\delta}_{2}=\text{i}0.1$
,
${\it\delta}_{3}=-\text{i}0.5$
,
$q_{1}=q_{2}=0.3$
,
$q_{3}=0.3$
in (b).
In this context, it is important to point out that the numerical scheme developed by Crowdy & Marshall (Reference Crowdy and Marshall2007) for the computation of the primary S–K function – which avoids the infinite product and relies on a more rapidly convergent Laurent series, as discussed in § 3.2 – can be easily adapted for the evaluation of the secondary prime functions
${\it\Omega}_{M}({\it\zeta},{\it\alpha})$
; see Vasconcelos et al. (Reference Vasconcelos, Marshall and Crowdy2014) for details. The numerical computation of
${\it\Omega}_{M}({\it\zeta},{\it\alpha})$
can thus be performed in an efficient manner for domains of arbitrary connectivity. Using this method, we have reproduced at considerably less computational cost the specific solutions reported by Green & Vasconcelos (Reference Green and Vasconcelos2014). Other examples of multi-bubble configurations are discussed below, where it is assumed that
$U=2$
.
In the particular case that the bubbles either are symmetrical about the centreline or have fore-and-aft symmetry, solutions can be obtained by reducing the flow region to a simply connected domain and then applying the standard Schwarz–Christoffel formula (Vasconcelos Reference Vasconcelos2001). These symmetrical solutions can be recovered from our formula by simply prescribing a domain
$D_{{\it\zeta}}$
with the appropriate symmetry. More precisely, centreline symmetry is enforced by choosing the centres of all circles
$C_{j}$
on the real axis, whereas bubbles with fore-and-aft symmetry are obtained by placing the centres of
$C_{j}$
on the imaginary axis. Two examples of assemblies of symmetrical bubbles are shown in figure 3. (In this and in the following figures, we have departed from our original convention and chosen to normalise the channel width to unity for convenience of presentation.)
More generally, if
$D_{{\it\zeta}}$
is reflectionally symmetric about the real (imaginary) axis, then the resulting bubble configuration has centreline (fore-and-aft) symmetry – but not all bubbles will necessarily have the symmetry of the overall solution (this happens only in the two cases just mentioned). For instance, in figure 4 we show examples of three-bubble assemblies in which the configuration as a whole has either centreline or fore-and-aft symmetry, but where only one of the bubbles (the largest one in each case) possesses the overall symmetry of the solution, whilst the other two bubbles are totally asymmetric.
Examples of assemblies of three bubbles in which the configuration as a whole has either centreline symmetry (a) or fore-and-aft symmetry (b), but where only one of the bubbles (the largest one in each case) has the overall symmetry of the solution. The conformal moduli of
$D_{{\it\zeta}}$
were chosen as follows:
${\it\delta}_{1}=-0.2$
,
${\it\delta}_{2}=0.25+\text{i}0.12$
,
${\it\delta}_{3}=0.25-\text{i}0.12$
,
$q_{1}=0.3$
,
$q_{2}=q_{3}=0.1$
in (a); and
${\it\delta}_{1}=-\text{i}0.2$
,
${\it\delta}_{2}=0.2+\text{i}0.25$
,
${\it\delta}_{3}=-0.2+\text{i}0.25$
,
$q_{1}=0.3$
,
$q_{2}=q_{3}=0.1$
in (b).
An example of an assembly of three bubbles in a general asymmetric configuration. Here the conformal moduli of
$D_{{\it\zeta}}$
are
${\it\delta}_{1}=0.2+\text{i}0.2$
,
${\it\delta}_{2}=-0.4+\text{i}0.2$
,
${\it\delta}_{3}=-\text{i}0.3$
,
$q_{1}=q_{2}=q_{3}=0.25$
.
A domain
$D_{{\it\zeta}}$
that entails no symmetry yields, of course, a completely asymmetric bubble configuration. An example of this case is given in figure 5 for an assembly of three asymmetric bubbles. Solutions for a higher number
$M$
of bubbles can be handled in a similar manner.
5. Solutions for multiple fingers and bubbles
In the instance that some of the bubbles within a multi-bubble solution become infinitely elongated, whilst the other bubbles remain of finite area, one obtains a situation where multiple fingers penetrate into the channel with an assembly of bubbles moving ahead of the fingers. To be specific, consider a situation with
$p$
fingers and
$M$
bubbles. Let us number the fingers from the bottom up and denote by
${\rm\pi}{\it\lambda}_{j}$
the width of the
$j$
th finger. Similarly, let us denote by
${\rm\pi}{\it\Delta}_{j}$
, for
$j=1,\ldots ,p+1$
, the widths of the fluid gaps separating adjacent fingers or between a finger and a channel wall, with the gaps also numbered from the bottom up. A schematic of the flow domain
$D_{z}$
in the
$z$
-plane is shown in figure 7 for the case
$p=2$
and
$M=1$
.
As before, it is convenient to work in an extended channel,
$D_{z}\cup \overline{D_{z}}$
, containing
$2p$
fingers and
$2M$
bubbles, where each additional interface is the reflection in the real axis of an interface in the original channel. The solution to this multifinger problem can be obtained from the multi-bubble solution given in § 4 by starting with an assembly of
$2(M+p)$
bubbles in the extended channel and then taking the limit in which
$2p$
bubbles become infinitely elongated so as to yield the desired fingers. This limit can be easily obtained in the
${\it\zeta}$
-plane, as follows. The
$p$
pairs of circles
$C_{j}$
and
$C_{-j}$
, for
$j=M+1,\ldots ,M+p$
, corresponding to the
$2p$
bubbles that will become fingers should coalesce into a single circle that orthogonally intersects the unit circle
$C_{0}$
and that encloses the point
${\it\zeta}=-1$
, which was the preimage of
$x=-\infty$
in the multi-bubble solution. The other circles
$C_{k}$
,
$k\in M$
, remain as they are (and so they will map to
$2M$
bubbles of finite area). The resulting flow domain, denoted by
$\tilde{F}_{0}$
, is shown in figure 6 for the case
$M=1$
and
$p=2$
.
The flow region
$\tilde{F}_{0}$
in the auxiliary complex
${\it\zeta}$
-plane for the case of multiple fingers and bubbles. Here is represented the case of
$p=2$
fingers and
$M=1$
bubbles in the original channel.
Before proceeding further, let us establish some notation. Label by
$C_{M+1}$
the circle orthogonal to
$C_{0}$
, and let
${\it\delta}_{M+1}$
and
$q_{M+1}$
denote its centre and radius, respectively. From the orthogonality condition one has
$$\begin{eqnarray}|{\it\delta}_{M+1}|^{2}=1+q_{M+1}^{2}.\end{eqnarray}$$
The Möbius map
${\it\theta}_{M+1}({\it\zeta})$
associated with reflection in
$C_{M+1}$
, see (3.2), is given by
$$\begin{eqnarray}{\it\theta}_{M+1}({\it\zeta})=\frac{{\it\delta}_{M+1}-{\it\zeta}}{1-\overline{{\it\delta}}_{M+1}{\it\zeta}}.\end{eqnarray}$$
It can be verified that this map is of order two:
${\it\theta}_{M+1}^{2}=1$
. (Here 1 denotes the identity map.) One may also verify that
$C_{M+1}$
is invariant under
${\it\theta}_{M+1}$
in the following sense:
$$\begin{eqnarray}C_{M+1}^{(in)}={\it\theta}_{M+1}(C_{M+1}^{(out)}),\end{eqnarray}$$
where
$C_{M+1}^{(in)}$
and
$C_{M+1}^{(out)}$
denote the segments of
$C_{M+1}$
that are inside and outside
$C_{0}$
, respectively. In fact, one can show that
${\it\theta}_{M+1}$
maps the interior of
$C_{M+1}$
onto its exterior. In analogy with the group
${\it\Theta}_{M}$
introduced in § 3.2, we define the Schottky group
$\tilde{{\it\Theta}}_{M}$
as the set consisting of all even combinations of the maps
$\{{\it\theta}_{j}\mid j=1,\ldots ,M+1\}$
. (The generators of this group are the maps
$\{{\it\theta}_{M+1}{\it\theta}_{k}~|~k\in K_{M}\}$
, which send the interior of
$C_{-k}$
onto the exterior of their images in
$C_{M+1}$
.)
The flow domain
$D_{z}$
for multiple fingers and bubbles in a Hele-Shaw channel.
As discussed above, the orthogonal circle
$C_{M+1}$
is to be mapped to the
$2p$
fingers, whereas the circles
$C_{k}$
,
$k\in K_{M}$
, map to the
$2M$
bubbles. This implies, in particular, that the function
$z({\it\zeta})$
must have
$2p$
logarithmic branch points on
$C_{M+1}$
, corresponding to the left end points of the fingers (i.e.
$x=-\infty$
). Let us denote by
${\it\zeta}_{j}$
the singularities that lie on
$C_{M+1}^{(in)}$
and by
${\it\zeta}_{j}^{\ast }$
those lying on
$C_{M+1}^{(out)}$
, where
$$\begin{eqnarray}{\it\zeta}_{j}^{\ast }\equiv {\it\theta}_{M+1}({\it\zeta}_{j}),\quad j=1,\ldots ,p+1.\end{eqnarray}$$
Note, in particular, that
${\it\zeta}_{1}$
and
${\it\zeta}_{p+1}$
are the points of intersection between
$C_{0}$
and
$C_{M+1}$
, and so
${\it\zeta}_{1}^{\ast }={\it\zeta}_{1}$
and
${\it\zeta}_{p+1}^{\ast }={\it\zeta}_{p+1}$
, as these are the two fixed points of the map
${\it\theta}_{M+1}({\it\zeta})$
. Recall also that by construction the circle
$C_{M+1}$
encloses the point
${\it\zeta}=-1$
, which was originally mapped to
$x=-\infty$
. This condition can be fulfilled without loss of generality by setting
$q_{M+1}=1$
and
${\it\delta}_{M+1}=-\sqrt{2}$
; this choice will be implied in the remainder of this section.
From the preceding discussion, it then follows that solutions with multiple fingers and bubbles can be obtained from the multi-bubble solutions (4.6) by replacing the original logarithmic singularity at
${\it\zeta}=-1$
with
$2p$
logarithmic singularities at the points
${\it\zeta}_{j},{\it\zeta}_{j}^{\ast }\in C_{M+1}$
, for
$j=1,\ldots ,p+1$
. In other words, the solution (for
$U=2$
) given in (4.6) becomes
$$\begin{eqnarray}z({\it\zeta})=\log \left(\frac{\displaystyle \mathop{\prod }_{j=1}^{p+1}[\tilde{{\it\Omega}}_{M}({\it\zeta},{\it\zeta}_{j})\tilde{{\it\Omega}}_{M}({\it\zeta},{\it\zeta}_{j}^{\ast })]^{{\it\alpha}_{j}/2}}{\tilde{{\it\Omega}}_{M}({\it\zeta},1)}\right),\end{eqnarray}$$
where
$\tilde{{\it\Omega}}_{M}({\it\zeta},{\it\alpha})$
is the S–K prime function defined over the group
$\tilde{{\it\Theta}}_{M}$
and the parameters
$0\leqslant {\it\alpha}_{j}\leqslant 1$
must satisfy the condition
$$\begin{eqnarray}\mathop{\sum }_{j=1}^{p+1}{\it\alpha}_{j}=1,\end{eqnarray}$$
to ensure single-valuedness of
$z({\it\zeta})$
in
$\tilde{F}_{0}$
. Note in particular that the widths of the fluid gaps (relative to the channel width) are given by
${\it\Delta}_{j}={\it\alpha}_{j}/2$
, for
$j=1,\ldots ,p+1$
, whereas the parameters
${\it\zeta}_{j}$
,
$j=2,\ldots ,p$
, control the relative widths of the individual fingers. From (5.6) it follows that
$\sum _{j=1}^{p+1}{\it\Delta}_{i}=1/2$
, which implies that the combined relative width of the fingers is
$\sum _{j=1}^{p}{\it\lambda}_{j}=1/2$
, as required by fluid mass conservation (recall that we have set
$U=2$
).
Once the parameters
$\{{\it\alpha}_{j}\mid j=1,\ldots ,p\}$
and
$\{{\it\zeta}_{j}\mid j=2,\ldots ,p\}$
are specified and the conformal moduli
$\{{\it\delta}_{j},q_{j}\mid j=1,\ldots ,M\}$
of the domain
$\tilde{D}_{{\it\zeta}}$
are prescribed, a specific solution for
$p$
fingers moving together with an assembly of
$M$
bubbles in a Hele-Shaw channel is obtained. The shapes of the different interfaces correspond to the images under the mapping (5.5) of the circles
$C_{j}$
,
$j=1,\ldots ,M+1$
. For example, the coordinates of the
$j$
th finger are given by
$$\begin{eqnarray}x_{j}^{(f)}(s)+y_{j}^{(f)}(s)=z(-\sqrt{2}+\text{e}^{\text{i}s}),\quad j=1,\ldots ,p,\end{eqnarray}$$
where
$\arg ({\it\zeta}_{j}+\sqrt{2})<s<\arg ({\it\zeta}_{j+1}+\sqrt{2})$
. A similar expression to that in (4.7) is obtained for the bubble coordinates.
A related solution for one symmetrical finger with an assembly of symmetrical bubbles in front of it was obtained before (Vasconcelos Reference Vasconcelos1999) using Schwarz–Christoffel methods. The present solutions are much more general in that they describe any number of fingers and bubbles with no symmetry assumption and, furthermore, have the advantage of being given by an explicit mapping function from which the shapes of the interfaces can be easily computed. An example with one asymmetric finger and one asymmetric bubble (
$p=M=1$
) is shown in figure 8, while figure 9 shows a configuration with two fingers and one bubble (
$p=2$
,
$M=1$
).
An asymmetric finger with a small asymmetric bubble in front of it. Here
$p=M=1$
and the parameters are
${\it\delta}_{1}=0.2$
,
$q_{1}=0.2$
and
${\it\alpha}_{1}=0.75$
.
A configuration with two fingers and one bubble (
$p=2$
and
$M=1$
) corresponding to the parameters
${\it\delta}_{1}=\text{i}0.2$
,
$q_{1}=0.2$
,
${\it\alpha}_{1}=0.25$
and
${\it\alpha}_{2}=0.5$
.
Before leaving this section, it is perhaps worth mentioning that in the case when no bubble is present (
$M=0$
), the S–K prime function becomes a monomial, i.e.
$\tilde{{\it\Omega}}({\it\zeta},{\it\alpha})={\it\zeta}-{\it\alpha}$
, yielding a solution (for
$p$
fingers only) of the form
$$\begin{eqnarray}z({\it\zeta})=\log \left(\frac{\displaystyle \mathop{\prod }_{j=1}^{p+1}[({\it\zeta}-{\it\zeta}_{j})({\it\zeta}-{\it\zeta}_{j}^{\ast })]^{{\it\alpha}_{j}/2}}{{\it\zeta}-1}\right).\end{eqnarray}$$
This solution describes (in a different representation) the multifinger solutions obtained by Vasconcelos (Reference Vasconcelos1998). Of particular interest is the case of a single finger (
$p=1$
), for which the expression above reduces to
$$\begin{eqnarray}z({\it\zeta})=-\log ({\it\zeta}-1)+\frac{{\it\alpha}}{2}\log ({\it\zeta}-{\it\zeta}_{1})+\left(1-\frac{{\it\alpha}}{2}\right)\log ({\it\zeta}-{\it\zeta}_{2}),\end{eqnarray}$$
where
$0<{\it\alpha}<1$
and
${\it\zeta}_{2}=\overline{{\it\zeta}_{1}}=(-1+\text{i})/\sqrt{2}$
. This recovers (in a different representation) the asymmetric finger obtained by Taylor & Saffman (Reference Taylor and Saffman1959) as a generalisation of their previous solution (Saffman & Taylor Reference Saffman and Taylor1958) for a symmetrical finger (
${\it\alpha}=1$
).
6. Periodic array of bubbles
In this section we consider the case of a periodic array of bubbles steadily moving in a Hele-Shaw channel. The problem is formulated in a general set-up where it is supposed that there is an arbitrary number of bubbles per period cell and no assumption is made as to the geometrical arrangement of the bubbles within a period cell; see figure 10 for a schematic. This contrasts with all previous periodic solutions (Burgess & Tanveer Reference Burgess and Tanveer1991; Vasconcelos Reference Vasconcelos1994; Silva & Vasconcelos Reference Silva and Vasconcelos2011, Reference Silva and Vasconcelos2013) where symmetry requirements are imposed a priori. As in the case of a finite assembly of bubbles discussed in § 4, exact solutions for the present problem are obtained in terms of a conformal mapping from a circular domain to the flow region exterior to the bubbles (within a period cell). Here, however, a fully fledged function-theoretic approach is used, whereby the desired mapping functions are obtained by directly exploiting the properties of the secondary prime functions.
6.1. Problem formulation
Consider a periodic assembly of bubbles moving with speed
$U$
in a Hele-Shaw channel. Here there are no factors of
${\rm\pi}$
to worry about and so we set the width of the channel to unity. The average fluid velocity,
$V$
, across the channel in the
$x$
-direction is also normalised to unity, i.e.
$V=1$
. The streamwise period is denoted by
$L$
, and it is assumed that there are
$M-1$
bubbles per period cell. Because of the flow periodicity, the problem can be restricted to the fluid region,
$D_{z}$
, within one period cell. A schematic of
$D_{z}$
is shown in figure 11(a) for the case
$M=3$
.
Periodic array of bubbles in a Hele-Shaw channel.
Now let
$D_{{\it\zeta}}$
be our usual circular domain with
$M$
inner circles; see figure 11(b). We shall seek a conformal mapping
$z({\it\zeta})$
from
$D_{{\it\zeta}}$
to
$D_{z}$
, such that the unit circle
$C_{0}$
maps to the upper edge of the period cell (
$y=1$
), the inner circle
$C_{M}$
maps to the lower edge (
$y=0$
), and the other inner circles
$C_{j}$
,
$j=1,\ldots ,M-1$
, map to the bubble boundaries
$\partial D_{j}$
. In addition, we insert a branch cut from a point on
$C_{M}$
to a point on
$C_{0}$
such that the two sides of this cut map to the left and right edges of the period cell, respectively; see figure 11(b). The specific ‘path’ of the branch cut is not relevant as it only affects the choice of period cell in the
$z$
-plane. (In the examples discussed below, the branch cut is placed on the positive real axis for convenience.)
6.2. General solutions
As before, let
$W({\it\zeta})$
and
$T({\it\zeta})$
denote the complex potentials in the laboratory and co-moving frames, respectively. Carrying out an analysis analogous to that presented in § 2.2, one can show that
$w=W({\it\zeta})$
maps
$D_{{\it\zeta}}$
onto a ‘rectangular’ region,
$D_{w}$
, in the
$w$
-plane which is bounded by two horizontal edges located at
$\text{Im}\,w=0$
and
$\text{Im}\,w=1$
(the images of
$C_{0}$
and
$C_{M}$
) and by two ‘curved’ lateral edges (the image of the branch cut from
$C_{0}$
to
$C_{M}$
), and which contains
$M-1$
slits in its interior (the images of
$C_{j}$
,
$j=1,\ldots ,M-1$
); see figure 11(c). In the same vein, one may verify that
${\it\tau}=T({\it\zeta})$
maps
$D_{{\it\zeta}}$
onto a ‘rectangular’ domain,
$D_{{\it\tau}}$
, in the
${\it\tau}$
-plane with
$M-1$
horizontal slits; see figure 11(d).
The flow domain within a period cell in: (a) the
$z$
-plane, (b) the
${\it\zeta}$
-plane, (c) the
$w$
-plane and (d) the
${\it\tau}$
-plane.
To obtain
$W({\it\zeta})$
, let us first introduce the following transformation:
$$\begin{eqnarray}S=\exp (-\text{i}{\it\lambda}w),\quad {\it\lambda}\in \mathbb{R}_{{>}0},\end{eqnarray}$$
which maps
$D_{w}$
onto a domain in a subsidiary
$S$
-plane that consists of a concentric annulus with
$M-1$
radial slits; see figure 12(a). The real parameter
${\it\lambda}$
allows us the freedom to vary the modulus of the annulus in the
$S$
-plane.
Now, it was shown by Vasconcelos et al. (Reference Vasconcelos, Marshall and Crowdy2014) that the function
$$\begin{eqnarray}S({\it\zeta})=\frac{{\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{M}({\it\alpha})){\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\alpha}))}{{\it\Omega}_{M-1}({\it\zeta},{\it\alpha}){\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\theta}_{M}({\it\alpha})))},\end{eqnarray}$$
where
$1\leqslant l\leqslant M-1$
and
${\it\alpha}\in {\it\varphi}_{0}(D_{{\it\zeta}})$
, maps the circular domain
$D_{{\it\zeta}}$
conformally onto a concentric annulus with
$M-1$
radial slits. Here
$C_{0}$
is mapped to the outer circumference of the annulus and
$C_{M}$
maps to the inner circumference, whereas
$C_{j}$
,
$j=1,\ldots ,M-1$
, map to the slits.
Using (6.1) and (6.2), one finds that the map
$w=W({\it\zeta})$
is of the form
$$\begin{eqnarray}W({\it\zeta})=\text{i}K\log \left[\frac{{\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{M}({\it\alpha})){\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\alpha}))}{{\it\Omega}_{M-1}({\it\zeta},{\it\alpha}){\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\theta}_{M}({\it\alpha})))}\right],\end{eqnarray}$$
where
$K$
is a real constant. The value of
$K$
is determined from the condition that the ‘rectangular’ domain
$D_{w}$
has a height equal to
$1$
, that is,
$$\begin{eqnarray}\text{Im}[W(A)-W(D)]=1,\end{eqnarray}$$
where the points
$A$
and
$D$
correspond to the intersections of the branch cut with the circles
$C_{0}$
and
$C_{M}$
, respectively; see figure 11(b) and 11(c).
The flow domain in the subsidiary
$S$
-plane (a) and
${\it\eta}$
-plane (b).
Similarly, to obtain the mapping
${\it\tau}=T({\it\zeta})$
one first applies an exponential transformation
$$\begin{eqnarray}{\it\eta}=\exp (\text{i}{\it\lambda}{\it\tau}),\quad {\it\lambda}\in \mathbb{R}_{{>}0},\end{eqnarray}$$
which maps
$D_{{\it\tau}}$
onto a domain
$D_{{\it\eta}}$
consisting of a concentric annulus with
$M-1$
concentric circular-arc slits; see figure 12(b). Next, we recall that, as shown by Crowdy & Marshall (Reference Crowdy and Marshall2007), the function
$$\begin{eqnarray}{\it\eta}({\it\zeta})=\frac{{\it\omega}({\it\zeta},{\it\theta}_{M}({\it\alpha}))}{{\it\omega}({\it\zeta},{\it\alpha})},\end{eqnarray}$$
where
${\it\alpha}\in {\it\varphi}_{0}(D_{{\it\zeta}})$
, maps
$D_{{\it\zeta}}$
onto the desired annular slit domain
$D_{{\it\eta}}$
. (Here
$C_{0}$
maps to the outer circumference of the annulus,
$C_{M}$
maps to the inner circumference, and
$C_{j}$
,
$j=1,\ldots ,M-1$
, map to the circular slits.)
Using (6.5) and (6.6), one then finds that
$$\begin{eqnarray}T({\it\zeta})=\text{i}K^{\prime }\log \left[\frac{{\it\omega}({\it\zeta},{\it\theta}_{M}({\it\alpha}))}{{\it\omega}({\it\zeta},{\it\alpha})}\right],\end{eqnarray}$$
where the prefactor
$K^{\prime }$
is determined from the requirement that the height of the domain
$D_{{\it\tau}}$
is equal to
$U-1$
:
$$\begin{eqnarray}\text{Im}[T(D)-T(A)]=U-1.\end{eqnarray}$$
In view of (3.7), one can rewrite (6.7) in terms of the secondary prime functions
${\it\Omega}_{M-1}({\it\zeta},{\it\alpha})$
:
$$\begin{eqnarray}T({\it\zeta})=\text{i}K^{\prime }\log \left[\frac{{\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{M}({\it\alpha})){\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\theta}_{M}({\it\alpha})))}{{\it\Omega}_{M-1}({\it\zeta},{\it\alpha}){\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\alpha}))}\right],\end{eqnarray}$$
where an additive constant was omitted.
Inserting (6.3) and (6.9) into (2.18), and performing some straightforward rearrangements, then yields the desired mapping
$z({\it\zeta})$
:
$$\begin{eqnarray}z({\it\zeta})=\frac{\text{i}K_{-}}{U}\log \left[\frac{{\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{M}({\it\alpha}))}{{\it\Omega}_{M-1}({\it\zeta},{\it\alpha})}\right]+\frac{\text{i}K_{+}}{U}\log \left[\frac{{\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\alpha}))}{{\it\Omega}_{M-1}({\it\zeta},{\it\theta}_{l}({\it\theta}_{M}({\it\alpha})))}\right],\end{eqnarray}$$
where
$$\begin{eqnarray}K_{\pm }=K\pm K^{\prime }.\end{eqnarray}$$
A periodic stream of bubbles with one bubble per period cell (
$M=2$
). The conformal moduli of
$D_{{\it\zeta}}$
are as follows:
${\it\delta}_{1}=-0.6$
,
${\it\delta}_{2}=0.0$
,
$q_{1}=0.3$
,
$q_{2}=0.15$
; the corresponding period is
$L=3.49$
. The lateral edges of the period cell,
$y=\pm L/2$
(dashed lines), are equipotentials of the flow.
Using the degrees of freedom allowed by the Riemann–Koebe mapping theorem (Goluzin Reference Goluzin1969), we can place the centre of the circle
$C_{M}$
at the origin, i.e.,
${\it\delta}_{M}=0$
. Its radius,
$q_{M}$
, is then a free parameter that essentially controls the period
$L$
. The remaining
$3M-3$
parameters, corresponding to the conformal moduli
$\{{\it\delta}_{j},q_{j}\mid j=1,\ldots ,M-1\}$
of
$D_{{\it\zeta}}$
, determine the centroid and area of the
$M-1$
bubbles in a period cell. Thus, once the domain
$D_{{\it\zeta}}$
is prescribed, a specific solution for a periodic assembly of bubbles can be readily computed from (6.10).
6.3. Examples
As seen above, obtaining specific solutions for a periodic array of Hele-Shaw bubbles requires computation of the secondary prime functions
${\it\Omega}_{N}({\it\zeta},{\it\alpha})$
for
$N=M-1$
. This can be done by using the infinite product (3.5) defined over the appropriate group
${\it\Theta}_{M-1}$
. (Alternative numerical schemes to compute
${\it\Omega}_{N}({\it\zeta},{\it\alpha})$
are known only for the cases
$N=M$
and
$N=1$
; see Vasconcelos et al. (Reference Vasconcelos, Marshall and Crowdy2014).) For the present purposes, it suffices to truncate the infinite product at the fourth-level maps (i.e. maps involving up to four generators of the group). Including higher order maps makes no discernible difference within the scale of the figures.
An example of a staggered two-file stream of bubbles. Here
$M=3$
and the conformal moduli of
$D_{{\it\zeta}}$
are
${\it\delta}_{1}=-0.1+\text{i}0.5$
,
${\it\delta}_{2}=-0.2-\text{i}0.3$
,
${\it\delta}_{3}=0.0$
,
$q_{1}=q_{2}=0.2$
,
$q_{3}=0.1$
. The period is
$L=2.82$
and two period cells are shown.
An example of a periodic array of bubbles with only one bubble per period cell is shown in figure 13. It follows from symmetry considerations that in this case the bubble shape has fore-and-aft symmetry, so that the vertical lines at
$x=\pm L/2$
are equipotentials of the flow. In particular, it is worth noting that when the bubble is also symmetrical about the centreline, our solution recovers the stream of symmetrical bubbles obtained by Burgess & Tanveer (Reference Burgess and Tanveer1991). As discussed before, periodic solutions with more general symmetrical arrangements, such as those found by Silva & Vasconcelos (Reference Silva and Vasconcelos2011, Reference Silva and Vasconcelos2013), can easily be reproduced in our formalism by simply choosing
$D_{{\it\zeta}}$
with the appropriate symmetry. Furthermore, the analytical solution given in (6.10) can handle asymmetric configurations with equal ease, as illustrated below.
An example of a staggered two-file array of unequal bubbles is shown in figure 14. For the choice of parameters used in this case (see caption of figure) the period is
$L=2.82$
, and two period cells are shown in the figure. Bubble configurations with a higher number of bubbles per period cell can be computed in a similar manner. An example with three asymmetric bubbles per unit cell is given in figure 15, where again two period cells are shown in the figure.
7. Discussion
The motion of an assembly of bubbles in a Hele-Shaw channel is a free boundary problem, which is made more difficult by the fact that the relevant field (the velocity potential) is defined over a multiply connected domain and on whose boundaries it satisfies mixed boundary conditions. This means that in the complex potential plane the flow region is a slit strip domain of mixed type: the bubbles are described by vertical slits, whilst the channel walls correspond to horizontal lines. Recently, a formalism based on the secondary S–K prime function was developed to construct a large class of such mixed slit maps (Vasconcelos et al. Reference Vasconcelos, Marshall and Crowdy2014).
An example of a periodic array of bubbles with three asymmetric bubbles per period cell (
$M=4$
) corresponding to the following choice of parameters:
${\it\delta}_{1}=0.2+\text{i}0.5$
,
${\it\delta}_{2}=-0.2-\text{i}0.5$
,
${\it\delta}_{3}=-0.38$
,
${\it\delta}_{4}=0.0$
,
$q_{1}=q_{2}=q_{3}=0.2$
,
$q_{4}=0.1$
. Here the period is
$L=2.76$
and two period cells are shown.
Here the formalism of the secondary prime functions was used to compute exact solutions for multiple bubbles steadily translating in a Hele-Shaw channel in various configurations: (i) a finite assembly of bubbles; (ii) multiple fingers moving together with an assembly of bubbles; and (iii) a periodic array of bubbles. In all cases considered, analytical formulae in terms of the secondary prime functions were obtained for the conformal mapping from a circular domain to the corresponding flow region in the physical plane. Several examples of specific solutions for these distinct arrangements were given. Taken together, the results reported here represent the complete set of solutions for multiple steady bubbles and fingers in a horizontal Hele-Shaw channel when surface tension is neglected.
A variant of the Hele-Shaw problem that has received less attention is the case when the cell is rotated about the centreline away from the horizontal, which introduces a gravitational potential transverse to the centreline. To the best of our knowledge, the only known solution for this situation (in the channel geometry) was obtained by Brener, Levine & Tu (Reference Brener, Levine and Tu1991) for the case of an asymmetric finger. It would be interesting to investigate whether this solution can be extended for the case of multiple bubbles in a rotated channel. The additional complication here is that the flow region in the complex potential plane consists of a strip with slanted slits (rather than vertical ones), and explicit conformal mappings to this type of slit domain are not yet known. Injection flows in an unbounded Hele-Shaw cell in the presence of a uniform gravitational field have also been considered using Schwarz function-type methods (Entov, Etingof & Kleinbock Reference Entov, Etingof and Kleinbock1995; McDonald Reference McDonald2011). It remains to be seen whether these alternative approaches could be adapted to the channel geometry and generalised to the multiply connected case.
Another possible extension of the present research would be to consider time-evolving bubbles in a Hele-Shaw cell. Time-dependent solutions for Hele-Shaw flows date back to the early work by Saffman (Reference Saffman1959), who found a solution for a finger growing from a flat interface. Since then, more general solutions have been found both for the growth of fingers in a channel (Dawson & Mineev-Weinstein Reference Dawson and Mineev-Weinstein1994) and for an expanding bubble in an unbounded cell (Shraiman & Bensimon Reference Shraiman and Bensimon1984; Bensimon & Pelcé Reference Bensimon and Pelcé1986; Howison Reference Howison1986). Several exact solutions for doubly connected time-dependent Hele-Shaw flows have also been obtained both in an unbounded cell (Richardson Reference Richardson1994, Reference Richardson1996a
; Crowdy Reference Crowdy2002; Dallaston & McCue Reference Dallaston and McCue2012) and in the channel geometry (Richardson Reference Richardson1996b
; Crowdy & Tanveer Reference Crowdy and Tanveer2004). Richardson (Reference Richardson2001) has obtained time-dependent solutions for injection flows in an unbounded cell involving a multiply connected fluid region, where the solution is given by a conformal mapping written in terms of a ratio of Poincaré theta series associated with a Schottky group (the group
${\it\Theta}_{0}$
in our notation). Crowdy & Marshall (Reference Crowdy and Marshall2004) used the formalism of the Schottky–Klein prime functions to construct conformal maps to multiply connected quadrature domains, of which the Hele-Shaw flow domains considered by Richardson (Reference Richardson2001) are a subclass. It would be interesting to investigate whether this approach could be generalised to other geometries. Recently, time-dependent solutions have also been obtained for injection flows in an unbounded Hele-Shaw cell in the presence of both a flat plate (Marshall Reference Marshall2015a
) and a circular cylinder (Marshall Reference Marshall2015b
), in which case the flow domain becomes doubly connected after the fluid has fully engulfed the obstacle. Extension of these results to the case of multiple obstacles – thus implying multiply connected domains – should in principle be possible with the help of the secondary prime functions.
Particularly relevant to the work presented herein is the class of time-dependent solutions for a single bubble in a Hele-Shaw channel recently obtained by Vasconcelos & Mineev-Weinstein (Reference Vasconcelos and Mineev-Weinstein2014). Their solution is a direct extension to the time-dependent case of the steady solution for a single bubble in a channel (Taylor & Saffman Reference Taylor and Saffman1959; Tanveer Reference Tanveer1987), and so it is natural to expect that the steady solutions reported here should also lead to a generalisation for time-dependent flows with multiple bubbles. Constructing time-dependent solutions for moving bubbles in a Hele-Shaw cell is of special importance in connection with the so-called selection problem, which concerns the question of why a bubble or a finger advances with a unique velocity – twice that of the background fluid velocity – despite the fact that there is a continuum of steady solutions for zero surface tension. It was generally accepted that surface tension was necessary for selecting the pattern (Combescot et al.
Reference Combescot, Dombre, Hakim, Pomeau and Pumir1986; Hong & Langer Reference Hong and Langer1986; Shraiman Reference Shraiman1986; Tanveer Reference Tanveer1987). However, recent studies (Mineev-Weinstein Reference Mineev-Weinstein1998; Vasconcelos & Mineev-Weinstein Reference Vasconcelos and Mineev-Weinstein2014; Khalid, McDonald & Vanden-Broeck Reference Khalid, McDonald and Vanden-Broeck2015) have demonstrated, by using time-dependent exact solutions without surface tension, that the steady solution with
$U=2$
is the only attractor for the non-singular unsteady solutions, thus showing that velocity selection in a Hele-Shaw cell is inherently determined by the zero-surface-tension dynamics. It has been conjectured (Vasconcelos & Mineev-Weinstein Reference Vasconcelos and Mineev-Weinstein2014) that a similar scenario holds in domains of arbitrary connectivity. The steady solutions presented here – in terms of the secondary S–K prime functions – should thus serve as a starting point to address the general selection problem for multiple bubbles. Work in this direction is currently underway. It is also hoped that the methods employed in the present work can be adapted to study other related physical systems, such as equilibrium configurations of multiple hollow vortices and the formation of multiple streamers in strong electric fields.
Acknowledgements
The author is appreciative of the hospitality of the Department of Mathematics at Imperial College London (ICL) where this work was carried out. He acknowledges financial support from a scholarship from the Conselho Nacional de Desenvolvimento Cientifico e Tecnologico (CNPq, Brazil) under the Science Without Borders programme for a sabbatical stay at ICL. Helpful discussions with D. Crowdy are gratefully acknowledged. The author would also like to thank J. Marshall for a critical reading of the manuscript and for his help in generating the data used in figures 14 and 15.
