## Abstract

This paper presents new solutions, in analytical form, for the shapes of an assembly of steady co-travelling bubbles in a Hele-Shaw cell. The associated velocity field is also derived. The assembly can consist of any finite number of bubbles. The solutions are expressed in terms of Schottky–Klein prime functions and are derived by exploiting results on multiply connected conformal slit mappings.

## 1. Introduction

A Hele-Shaw cell is an apparatus in which a thin layer of viscous fluid is trapped between two glass plates. The evolution of the free interfaces in such a situation has constituted an important theoretical paradigm in the study of free surface flows for over a century. It has become a standard model of what have been dubbed ‘Laplacian growth processes’ owing to the fact that the governing field equation (for a scalar velocity potential) in the fluid region is Laplace's equation: the fluid interface then evolves with a velocity that is proportional to the local gradient of this scalar field.

It is a remarkable fact that, even though the free boundary problem is highly nonlinear, a wide variety of analytical solutions to it have been found both for steady and unsteady cases. S. D. Howison (www.maths.ox.ac.uk/howison/Hele-Shaw) has compiled a comprehensive list of references, up to 1998, of many important contributions to the study of Hele-Shaw flows. One of the earliest solutions for a steady bubble in a Hele-Shaw cell was found by Taylor & Saffman (1959). They studied the steady translation of a single bubble in a channel (so it is assumed that side walls are present in the cell) under the assumption that the bubble is centred on the channel centreline and is reflectionally symmetric about it. A bubble that is small compared with the channel width is found to be close to elliptical in shape. Their analysis was later generalized by Tanveer (1987) who established the functional form of the solutions for a bubble in a channel with no assumed symmetries. Combescot & Dombre (1988) have retrieved the asymmetric single-bubble solution of Tanveer (1987) using different methods based on Riemann–Hilbert problems. Other analytical solutions for steady co-travelling bubbles in a Hele-Shaw channel include those for an infinite stream of steady bubbles due to Burgess & Tanveer (1991; who made use of arguments of Schwarz–Christoffel type) and a related study by Vasconcelos (1994), who used related methods to find other solutions for infinite streams of bubbles in a channel. A further study by Vasconcelos can be found in Vasconcelos (2001).

In this paper, we study the motion of a *finite number* of bubbles steadily translating, in the absence of side walls, in an unbounded Hele-Shaw cell. Of course, side walls are usually present in an experimental apparatus, so it is useful to think of our analysis as being applicable when all the bubbles are far from any boundary walls (so that the influence of the latter can be neglected). Our methods differ from all previous studies of this problem in that they enable us to find steady solutions for any (finite) number of interacting bubbles with no assumed geometrical symmetries and no single (or double) periodicity conditions (that is, we do not need infinitely many bubbles to perform the analysis). The solution class allows the specification of the centroid and area of each bubble. In more recent work, the author has extended the analysis here to *include* the effects of side walls on the motion and shape of the bubble assembly (Crowdy 2009).

It is worth mentioning a second possible application of our results. Recently, it has been found that a convincing free boundary model (Meulenbroek *et al*. 2004) for the evolution of steady ‘streamers’ in strong electric fields (also known as Lozansky–Firsov streamers (Lozansky & Firsov 1973)) is mathematically equivalent to the Hele-Shaw free boundary problem studied here. Specifically, the ionization front in a numerical simulation of multiple steady streamer evolution appears to be remarkably well fit by a periodic array of steadily translating Saffman–Taylor fingers (Luque *et al*. 2008). It is possible, therefore, that the results herein might also be relevant to that quite separate physical problem.

## 2. Mathematical formulation

Let *D* be the region, in a two-dimensional (*x*, *y*)-plane, exterior to *M*+1 bubbles each with finite area. Let the bubbles be denoted and their boundaries be . Figure 1 shows a schematic. We will suppose that each bubble is travelling at a steady speed *U* and that the fluid velocity far away from the assembly is in the *x*-direction with speed *U*_{0}. Without loss of generality, we set *U*_{0}=1. The velocity field ** u** is derived from a velocity potential

*Φ*so that

**=∇**

*u**Φ*.

*Φ*must be single valued in the region

*D*since, within the Hele-Shaw model, it is proportional to the fluid pressure.

The mathematical problem is then to solve(2.1)everywhere in *D*. Let be a complex potential of the flow. If each bubble is assumed to be of fixed area (so that there are no sources or sinks in each bubble), then it follows that(2.2)where we have used the divergence theorem. On writing this in complex form, it follows that(2.3)where d*s* is an arc length element around the boundary and square brackets around a quantity denote its change after traversing a circuit of ∂*D*_{j}. From these conditions, and since *Φ* is known to be single valued in *D*, we deduce that is itself single valued in *D* (in addition to being analytic there).

Let us move to a frame of reference co-travelling with the bubbles at speed *U*. Let the complex potential in this frame be and . Since, in this co-travelling frame, the bubble boundaries are streamlines then(2.4)where is a set of real constants. In addition, we must have that(2.5)where is a set of real constants (related to the values of the constant pressure in each bubble). In the far field, we have to impose that(2.6)

## 3. Conformal mapping

We define a *circular domain D*_{ζ} in a complex parametric *ζ*-plane to be a domain whose boundaries are all circles. Circular domains are a canonical class of planar domains because every planar domain is conformally equivalent to some circular domain (Goluzin 1969). In particular, the fluid domain *D* is conformally equivalent to some such domain.

Take *D*_{ζ} to be the *M*+1 connected circular domain in a *ζ*-plane consisting of the unit disc with *M* smaller discs excised from its interior. The outer boundary of *D*_{ζ} is the unit circle which we label *C*_{0}. Label the *M* inner boundaries of *D*_{ζ} as *C*_{1}, …, *C*_{M}. For let the centre and radius of *C*_{k} be *δ*_{k} and *q*_{k}, respectively. Figure 1 illustrates a triply connected case. Since the domain *D* is unbounded, there will be some point *ζ*=*β*, strictly inside *D*_{ζ}, at which *z*(*ζ*) has a simple pole singularity. By the degrees of freedom of the Riemann mapping theorem, we can pick *β* as we like. That leaves a single rotational degree of freedom in the Riemann mapping theorem still to be set.

On introduction of the composed function , the boundary conditions take the form(3.1)The two functions to be found are *W*(*ζ*) and *z*(*ζ*). *W*(*ζ*) is analytic and single valued everywhere in *D*_{ζ} except for a simple pole at *ζ*=*β* (recall that *z*(*ζ*) shares all these properties). The boundary conditions (3.1) have a particularly special form which will allow us to construct the two functions we seek by reference to some prior results. Two of the constants, *c*_{0} and *d*_{0} say, can be arbitrarily specified. Then, since both *W*(*ζ*) and *z*(*ζ*) must be single valued in *D*_{ζ}, the values of the constants will be determined once *W*(*ζ*) and *z*(*ζ*) have been found.

The relevant prior results are due to Crowdy & Marshall (2006), who have shown how to derive the functional form of conformal mappings from multiply connected circular domains (such as *D*_{ζ}) to all the various types of canonical multiply connected slit domains (Nehari 1952). Among these is the class of *parallel slit domains* (Nehari 1952; Crowdy & Marshall 2006). Such mappings take circular domains *D*_{ζ} to unbounded regions in which each boundary circle *C*_{j} (for ) maps to a slit of finite length. All the slit images of the circles are parallel and their common angle of inclination can be freely specified.

The mapping from *D*_{ζ} to a collection of *M*+1 horizontal slits (parallel to the real axis) is given by Crowdy & Marshall (2006)(3.2)where(3.3)*α* is some point inside *D*_{ζ} and *ω*(.,.) is the so-called *Schottky*–*Klein prime function* (Baker 1897; Hejhal 1972) associated with *D*_{ζ}. Its definition, and how to compute this function, is described in detail in the appendix A. Its definition depends purely on the parameters characterizing the domain *D*_{ζ} and, in many cases, it can be readily evaluated by means of an infinite product formula given in (A 10).

Similarly, the mapping from *D*_{ζ} to a collection of *M*+1 vertical slits (parallel to the imaginary axis) is given by Crowdy & Marshall (2006)(3.4)Suppose that, near *ζ*=*β*,(3.5)for some parameter *a* which, using a remaining rotational degree of freedom of the Riemann mapping theorem, can be assumed real. Since, as *ζ*→*β*,(3.6)we therefore pose that(3.7)(Note that it is important to let *α* equal *β after* derivatives have been taken.) It is clear that, being a real multiple of a horizontal slit mapping, such a function satisfies the conditions(3.8)for some set of constants as well as the condition that as *ζ*→*β*. It is also a single valued function in *D*_{ζ}. On each of the circles , d*W*/d*ζ* has two simple zeros (at the points on the circles corresponding to the ends of the slits). These correspond, physically, to stagnation points of the flow which appear on each bubble boundary. We could also add any complex constant to the expression (3.7)—corresponding to modifying the choices of *c*_{0} and *d*_{0}—but this would be inconsequential because these constants do not affect the final solution for the velocity field.

In the same spirit, we pose that(3.9)where *A* is a real constant to be determined. Since the right-hand side is a real multiple of a vertical slit mapping, we clearly have(3.10)for some set of constants . The right-hand side of (3.9) also has a simple pole at *ζ*=*β* (as required) and is single valued in *D*_{ζ}. Since *W*(*ζ*) is single valued in *D*_{ζ}, it follows that *z*(*ζ*) is too. Indeed, we deduce(3.11)which can be rearranged as(3.12)Making use of (3.6), we must have(3.13)implying *A*=−*a*. Thus,(3.14)where we have simply added an arbitrary constant (which does not affect the fact that the function satisfies all the required properties). (3.14) is the conformal mapping we seek.

Let us count the undetermined parameters in (3.14). Since we are free to specify *β* and since *a* is real, there are three remaining real degrees of freedom in the choice of *a* and . There are also 3*M* real degrees of freedom in the parameters . This makes a total of 3*M*+3 parameters. However, we expect to be able to specify the centroid and area of each of the *M*+1 bubbles in the cell, which makes a total of 3*M*+3 conditions. The number of parameters in the mapping therefore equals the number of conditions, so the counting is consistent.

To summarize the key results, it has been determined that the complex potential *W*(*ζ*) and conformal mapping *z*(*ζ*) from *D*_{ζ} to the region exterior to a configuration of steadily translating bubbles are given by(3.15)where *G*_{0} is defined in (3.3) and *β* is some chosen point in *D*_{ζ}. Knowledge of these two functions determines both the free surface bubble shapes and the flow field around them.

## 4. Elliptical solutions: *M*=0

Taylor & Saffman (1959) established that a single steadily translating bubble away from any channel walls is elliptical in shape (this result is also derived, using alternative methods, in Entov *et al*. (1993)). We now demonstrate that formulae (3.15) are consistent with this result thereby providing a check on the mathematical solutions.

In the simply connected case when *M*=0, we have(4.1)so that(4.2)Substitution of (4.2) into (3.15) then yields(4.3)Let us fix the centroid of the bubble by setting(4.4)Then,(4.5)The choice of *β* is arbitrary (provided it is in the unit *ζ*-disc), so taking the limit of (4.5) as *β*→0 gives(4.6)which is well known to be the conformal mapping from the interior of the disc |*ζ*|≤1 to the exterior of an ellipse centred at *z*=0. Altering *U* changes the eccentricity of the elliptical bubble: as *U*→1, the map tends to which maps to the region exterior to a vertical slit (aligned with the *y*-axis); as *U*→∞, the map tends to which maps to the region exterior to a horizontal slit.

In this way, we see that formulae (3.15) generalize this single-bubble result to the case of an arbitrary number of bubbles.

## 5. Two-bubble solutions: *M*=1

When *M*=1 we can, without loss of generality, take *D*_{ζ} to be the concentric annulus *ρ*<|*ζ*|<1 (we then no longer have the freedom to specify *β* arbitrarily). In this case, the Schottky–Klein prime function (see the appendix A for details) is(5.1)where(5.2)and . Thus(5.3)so that(5.4)where *P*′(*ζ*, *ρ*) denotes the derivative of *P*(*ζ*, *ρ*) with respect to its first argument. It follows from (3.15) that the conformal mapping in this case is(5.5)where we define(5.6)and is some constant.

In what follows, we will focus on bubbles of equal size that are symmetrically disposed about the origin (although the solution class just described captures much more general geometrical configurations). It is possible to show, using arguments based on the Schwarz reflection principle (we omit the details), that if the choice is made then the geometrical configuration will be an identical pair of vertically aligned bubbles of equal area, which are reflections of each other in the real axis and which each have reflectional symmetry about the imaginary axis. It can also be argued that, for this choice of *β*, the point *ζ*=−*β* corresponds to the point midway between the two symmetric bubbles. We therefore fix by insisting that the two bubbles are symmetrically disposed with respect to the origin *z*=0 and, to do this, we enforce the condition . Thus, two symmetric bubbles whose centres are aligned with the imaginary axis and reflectionally symmetric about both the real and imaginary axes is given by(5.7)This depends on three real parameters: *a*, *ρ* and *U*. Specifying the area of each of the bubbles fixes *a*. Setting *ρ* can be thought of as specifying the separation of the two symmetric bubbles. Figure 2 shows an array of bubbles in a horizontal field travelling at different velocities *U*. The area of each bubble is *π* and we have chosen the fixed value *ρ*=0.4. A typical streamline distribution is shown, for *U*=2 and 4, in figures 3 and 4. Note that the bubbles become increasingly elongated in the horizontal direction as their speed (in that direction) increases.

Finally, we mention that the functions *P*(*ζ*, *ρ*) and *K*(*ζ*, *ρ*) can be related, respectively, to the Weierstrass sigma and zeta functions (Whittaker & Watson 1927) so that the above solutions can be re-expressed in terms of those special functions if desired.

## 6. Three-bubble solutions: *M*=2

We now calculate some solutions involving three interacting bubbles. On use of (3.3) we obtain(6.1)where *ω*_{2}(.,.) denotes the derivative of the Schottky–Klein prime function *ω*(.,.) with respect to its second argument. Substitution of these into (3.15) gives the conformal mapping as(6.2)where is some constant.

Figure 5 shows both two- and three-bubble solutions of similar kind to highlight the fact that the present solution scheme gives exact solutions for any finite collection of bubbles. To compute the three-bubble solution, the infinite product representation (A 10) of the Schottky–Klein prime function was employed for its evaluation. The derivative *ω*_{2}(.,.) was computed numerically using a centred finite-difference scheme.

Solutions involving additional bubbles can be computed in an exactly similar way. To incorporate larger numbers of bubbles, it is simply necessary to excise greater numbers of circles from the unit disc to form the relevant circular domain *D*_{ζ}. Then, the appropriate Schottky–Klein prime function can be constructed in the manner described in the appendix A. Formulae (3.15) give analytical expressions for the solutions in all cases.

## 7. Discussion

The solution class embodied in formulae (3.15) allows the speed *U* of the bubbles to be arbitrarily specified, even for given area and centroid of each bubble. It is a matter for future investigation to study questions, whose answers are expected to involve exponential asymptotics, concerning how regularization effects—such as surface tension on the bubble boundaries (Tanveer 1986; Combescot & Dombre 1988)—will answer the ‘selection problem’ determining the speed of a finite assembly of bubbles. The existence of exact solutions to the unregularized problem is expected to be an important component in such an analysis. Related to this matter, the linear stability of the solutions (given some regularization) is also of interest.

In addition to solving a nonlinear free boundary problem, the work has a number of points of mathematical interest. First, the new solutions are expressed in terms of a special transcendental function called the Schottky–Klein prime function (Baker 1897; Hejhal 1972). This function is a natural one to use when studying problems involving multiply connected geometries and has, for example, been used by Crowdy & Marshall (2004) in the related problem of constructing multiply connected quadrature domains. The author (Crowdy 2008) has given a recent survey of the applications of this important, but little used, special function. A second mathematical point of interest is that the solutions are derived by importing more general results on the construction of multiply connected conformal slit mappings previously found in Crowdy & Marshall (2006). Finally, it is worth mentioning that our approach has features in common with the solution of the different problem of finding the irrotational uniform potential flow past an array of circular cylindrical obstacles (Crowdy 2006).

## Acknowledgments

The author acknowledges support from a 2004 Philip Leverhulme Prize in Mathematics, an EPSRC Advanced Research Fellowship and partial support from the European Science Foundation's MISGAM and HCAA research networks. He also acknowledges discussions with S. Tanveer.

## Footnotes

- Received June 17, 2008.
- Accepted September 22, 2008.

- © 2008 The Royal Society