## Abstract

Exact solutions are presented for a doubly periodic array of steadily moving bubbles in a Hele-Shaw cell when surface tension is neglected. It is assumed that the bubbles are either symmetrical with respect to the channel centreline or have fore-and-aft symmetry, or both, so that the relevant flow domain can be reduced to a simply connected region. By using conformal-mapping techniques, a general solution with any number of bubbles per unit cell is obtained in integral form. Several examples are given, including solutions for multi-file arrays of bubbles in the channel geometry and doubly periodic solutions in an unbounded cell.

## 1. Introduction

Since the seminal papers by Saffman & Taylor (1958) and Taylor & Saffman (1959), the motion of bubbles in a Hele-Shaw cell, where the fluids are confined between two closely spaced parallel plates, has attracted a great deal of attention. Because the problem is amenable to analytic treatment, particularly in the case when surface tension effects are neglected, many exact solutions have been found for both steady and time-dependent flows. A comprehensive list of references, up to 1998, on Hele-Shaw flows was compiled by S. D. Howison (http://people.maths.ox.ac.uk/howison/Hele-Shaw). Hele-Shaw flows have also recently been recognized as having deep connections with other areas of mathematical physics (Mineev-Weinstein *et al.* 2000; Gustafsson & Vasil’ev 2006). Exact solutions for Hele-Shaw flows are thus of interest, not only for their eventual physical significance, but also on their own right, for they are solutions to a difficult free-boundary problem. Some known exact solutions for the motion of bubbles and fingers in a Hele-Shaw cell were briefly reviewed by Vasconcelos (2001). More recently, other exact solutions have been found, such as time-dependent solutions for Hele-Shaw flows around a wedge (Cummings 1999; Richardson 2001*b*,*c*; Vasconcelos 2007) and solutions for multiple steady bubbles in a Hele-Shaw cell with no assumed symmetry obtained by Crowdy (2009*a*,*b*).

Exact solutions for a periodic stream of bubbles in a Hele-Shaw cell have also been found. The first of such solutions was obtained by Burgess & Tanveer (1991) for a stream of bubbles in channel geometry, with one bubble per unit cell. In the limit that the period is taken to infinity, their solution recovers the solution originally given by Taylor & Saffman (1959) for a single bubble in a channel. Exact solutions for a stream of bubbles in an unbounded Hele-Shaw cell were subsequently found by Vasconcelos (1993). The Burgess–Tanveer solution was later generalized by Vasconcelos (1994) for a wider class of solutions, with an arbitrary number of symmetrical bubbles per unit cell. Here again, taking the period to infinity yields solutions with an arbitrary but finite number of symmetrical bubbles in a Hele-Shaw channel (Vasconcelos 2001). It is worth mentioning that a periodic solution in a Hele-Shaw channel can be extended to the entire plane by successive reflections at the channel walls, thus yielding a doubly periodic solution in an unbounded cell. In this context, exact solutions for doubly periodic flows in a Hele-Shaw cell are of interest, not only from a mathematical viewpoint, for they correspond to a more general situation, but also from a physical perspective because, in certain cases, they can be realized as a stream of bubbles in a Hele-Shaw channel and hence may (in principle) be susceptible to experimental investigation. A brief experimental study of a stream of bubbles rising in an inclined Hele-Shaw cell was reported by Maxworthy (1986).

In the present paper, we report a rather general class of exact solutions for a doubly periodic array of steadily moving bubbles in a Hele-Shaw cell. Our generic solutions are valid in an unbounded Hele-Shaw cell, but there is a large subclass of solutions that correspond to a stream of bubbles in a Hele-Shaw channel. To render the problem analytically tractable, we neglect surface tension effects and suppose that the bubbles either are symmetrical with respect to the channel centreline or have fore-and-aft symmetry (or both), so that the flow domain can be reduced to a simply connected unit cell. By using conformal-mapping techniques, a general solution for any number of bubbles per unit cell is obtained in integral form. Our solution represents the most general periodic solution for a stream of bubbles in a Hele-Shaw cell, in the sense that all previously known periodic solutions are particular cases of the solutions reported here. Furthermore, our solutions include several novel configurations that were not obtained before, such as the case of a staggered two-file array of bubbles in a channel which somewhat resembles the ‘zipper’-like flows of red cells in capillaries (Sugihara-Seki & Fu 2005). Also presented are solutions with ‘mixed symmetry’ in the channel geometry where, within a given bubble configuration, some bubbles have fore-and-aft symmetry and others are symmetrical about the channel centreline. Examples are also given of solutions that cannot be restricted to the channel geometry and hence must necessarily be considered in an unbounded Hele-Shaw cell.

## 2. Mathematical formulation

We consider the problem of a doubly periodic array of bubbles moving with a constant velocity *U* along the *x*-direction in an unbounded Hele-Shaw cell. For convenience, we shall work in a reference frame moving with the bubbles in which the flow is steady and the bubbles stationary. We assume that the solutions have streamwise and spanwise periods of 2*L* and 2*a*, respectively, and so we reduce the problem to a rectangular unit cell, as shown in figure 1. We place the origin of our system of coordinates at the centre of the unit cell, so that its upper and lower edges are at *y*=±*a* and the left and right edges at *x*=±*L*. We shall assume furthermore that the bubble configuration is symmetrical with respect to reflections about the following four axes: *x*=0, *x*=*L*, *y*=0 and *y*=*a*. With this assumption, we only need to consider the problem in a reduced unit cell corresponding to one quarter, say, the upper right quarter, of the original unit cell; see figure 2*a*.

Now let denote the domain occupied by the fluid within the reduced unit cell. For simplicity of notation, we shall denote by a generic bubble–fluid interface. As is usual for Hele-Shaw flows, we introduce the complex potential *W*(*z*)=*ϕ*+i*ψ*, where *z*=*x*+i*y*, *ϕ*(*x*,*y*) is the velocity potential and *ψ*(*x*,*y*) is the associated stream function. We recall that the velocity potential (in the moving frame) is given by *ϕ*(*x*,*y*)=−(*b*^{2}*p*/12*μ*)−*Ux*, where *b* is the gap between the cell plates, *μ* is the viscosity and *p* is the fluid pressure. The complex potential *W*(*z*) must be analytic in the fluid domain and satisfy the appropriate boundary conditions, as we now discuss.

The pressure inside a bubble is assumed to be constant, and surface tension effects are neglected. The complex potential *W* must then satisfy the following condition on a fluid–bubble interface :2.1where *C* is a complex-valued constant. (The value of the constant *C* depends on the specific bubble one considers and is related to the position of the bubble centroid.) The complex potential must also satisfy specific conditions on the boundaries of the unit cell, as follows. Given that the upper and lower edges of our unit cell are streamlines of the flow and that the left and right edges are equipotentials, we must then have
2.2
2.3
2.4and
2.5
where Im and Re denote real and imaginary parts, respectively. Here, *V* is the average velocity across the unit cell in the *x*-direction in the laboratory frame and is a positive constant whose physical meaning will become clear shortly. An explicit expression for in terms of *L* and *V* , and the total area *J* occupied by the bubbles (within a unit cell), can be obtained. We refer the reader to Burgess & Tanveer (1991) for the details and simply quote the result here,2.6The complex potential *W*(*z*) defined above can be seen as a conformal mapping from the fluid region in the *z*-plane onto the corresponding flow domain in the *W*-plane. In view of equations (2.1)–(2.5), one readily sees that the flow domain in the *W*-plane corresponds to a rectangle, , −(*U*−*V*)*a*<*ψ*<0, as shown in figure 2*b*. Bubbles placed along the left and right edges of the unit cell in the *z*-plane are mapped onto horizontal slits in the *W*-plane, whereas bubbles on the streamlines *y*=0 and *a* are mapped onto segments along the lines *ψ*=−*a*(*U*−*V*) and *ψ*=0, respectively; see figure 2*b*. Note that the horizontal slits in figure 2*b* correspond to the respective halves of the bubbles on the left and right edges that are within the reduced unit cell, and hence such slits always penetrate into the rectangular domain in the *W*-plane.

As discussed in detail elsewhere (Tian & Vasconcelos 1994; Vasconcelos 2001), steady Hele-Shaw flows have an interesting rotation invariance in the following sense: if a curve is a solution for a bubble moving with constant velocity *U* along the *x*-direction in an unbounded Hele-Shaw cell, then the curve obtained from a rotation of about the origin by an angle *α* is also a solution with the same velocity *U*. Alternatively, one can view the rotated solution as one in which the bubble itself is kept fixed while its velocity is rotated by an angle *α*. In view of the rectangular geometry of our unit cell, it is particularly relevant for us to consider here a rotation by 90^{°}, where the bubbles move with the same velocity *U* as in the original solution, but now in the *y*-direction. In this case, the region occupied by the fluid in the rotated solution is exactly the same as that for the original problem, the main difference being that the upper and lower edges of the unit cell which are streamlines in the original solution become equipotentials in the rotated problem, whereas the left and right edges which were originally equipotentials become streamlines.

As shown by Vasconcelos (2001), the complex potential for the rotated problem (in the corresponding moving frame) is given by2.7From equations (2.1)–(2.5) and equation (2.7), it then follows that the rotated complex potential satisfies the following boundary conditions:2.8
2.9
2.10
2.11and
2.12
From equations (2.12), one now sees that the parameter corresponds to the average fluid velocity (in the laboratory frame) in the *y*-direction of the rotated flow. The flow domain for the rotated problem in the -plane corresponds to a rectangle, , , with horizontal slits corresponding to the bubbles placed on *y*=0 and *a*, as shown in figure 2*c*. (For the same reasons as discussed before, these horizontal slits always penetrate into the bounded domain in the -plane.)

Consider now the conformal mapping *z*(*ζ*) from the upper half-*ζ*-plane onto the fluid domain in the *z*-plane, so that the boundary of the fluid domain in the *z*-plane is mapped onto the real axis of the *ζ*-plane, as shown in figure 2*d*. In particular, a bubble will be the image under *z*(*ζ*) of a corresponding interval on the real-*ζ*-axis (see below). For definiteness, we also choose to map the point to an arbitrary point on the upper edge between the upper left corner and the next bubble; see figure 2. In the most general configuration, one has one bubble on each corner and an arbitrary number of bubbles on each edge of the cell. To keep track of all the bubbles, let us denote by , *i*=1,…,4, the bubble on the *i*th corner, with the corners numbered anti-clockwise starting from the upper left corner. Under the map *z*(*ζ*) a corner bubble will be the image of the interval on the real-*ζ*-axis defined by2.13(The reason for distinguishing between even and odd corners is one of convenience, for with the above choice, one has that in the *W*-plane, the corners are the images of the points *ζ*=*α*_{i}, whereas in the -plane, the corners correspond to *ζ*=*β*_{i}; see figure 2.) Note that if we want a configuration with no bubble in the *i*th corner, for a given *i*, we must simply set *α*_{i}=*β*_{i}.

Similarly, we shall denote the bubbles on the edges of the unit cell by , where the superscript refers to the corresponding side, with *a*=l, r, u, d, indicating left, right, up and down, respectively, and *i*=1,…,*n*_{a}, where *n*_{a} is the number of bubbles on the respective edge. In the *ζ*-plane a given bubble will correspond to an interval on the real axis defined by2.14Owing to the three degrees of freedom allowed by Riemann’s mapping theorem, we could arbitrarily fix the values of three parameters, chosen among *α*_{i}, *β*_{i} or . At this stage, however, it is best to consider all such quantities as free parameters, deferring to the examples the moment when we shall need to fix the appropriate parameters.

With some abuse of notation, let us now writeso that the functions *W*(*ζ*) and can be viewed as the conformal mappings from the upper half-*ζ*-plane onto the flow domains in the *W*- and -planes, respectively. It then follows from equation (2.7) that the mapping *z*(*ζ*) can be written as2.15Once the mappings *W*(*ζ*) and are known, the bubble interface can then be readily computed from equation (2.15) by setting *ζ*=*s*, for .

Summarizing our procedure so far, we have reduced our original free-boundary problem to the much simpler problem of obtaining two conformal mappings, namely, *W*(*ζ*) and , from the upper half-*ζ*-plane onto respective rectangular domains in the *W*- and -planes. As we will see shortly, such mappings can be easily obtained from the Schwarz–Christoffel formula (Carrier *et al.* 1983). It should also be noted that in previous solutions which use conformal mapping (e.g. Tanveer 1987), once the mapping *W*(*ζ*) from the chosen domain in the *ζ*-plane to the flow domain in the *W*-plane is known, then the mapping *z*=*f*(*ζ*) is constructed explicitly so as to satisfy the appropriate boundary conditions. Our method has the advantage of making this step rather straighforward by identifying the function with the complex potential for the rotated problem, so that solving for the mapping then completes the solution. Before presenting our general solutions, however, it is worthwhile to point out some special features of the solutions for the case *U*=2*V* .

Solutions with *U*=2*V* are special in the sense that solutions for any *U*>*V* can be generated by a proper rescaling of the solutions for the former case (Vasconcelos 2001). To show this, we first note that from inspection of the flow domain in the *W*-plane (figure 2*b*), it follows that the mapping *W*_{U}(*ζ*), for given *a* and any *U*>*V* , can be obtained from the corresponding mapping *W*_{2V}(*ζ*), for *U*=2*V* , by the relation2.16In particular, the velocity (see figure 2*b*) for both cases are related by2.17Similarly, inspection of the flow domain in the -plane (see figure 2*c*) reveals that2.18with2.19Inserting equation (2.17) into equation (2.19) yields the half-period *L*_{U} in terms of the corresponding value *L*_{2V} for the case *U*=2*V* ,2.20

It then follows from equations (2.15), (2.16) and (2.18) that for any *U*>*V* , the mapping *z*_{U}(*ζ*) can be written in terms of the solutions for *U*=2*V* as2.21which can be more conveniently rewritten as2.22where2.23with *ρ*∈(−1,1). This result thus shows that solutions for any *U*>*V* can be obtained by an appropriate rescaling of the solution with *U*=2*V* (Vasconcelos 2001). In view of this property, we shall only consider the case *U*=2*V* in the remainder of the paper. Without loss of generality, we will also set *V* =1 and *a*=1 henceforth.

## 3. General solutions

As indicated in figure 2*b*, the flow domain in the *W*-plane is a rectangle with *n*_{l} (*n*_{r}) horizontal slits corresponding to the bubbles on the left (right) side of the unit cell. This domain can be viewed as a degenerate polygon and so the corresponding mapping *W*(*ζ*) can be obtained from the Schwarz–Christoffel formula (Carrier *et al.* 1983). One finds that *W*(*ζ*) is determined by3.1where *K* and are real-valued constants to be determined later. In the -plane, the flow domain is of the same form as that in the *W*-plane, namely, a rectangle with *n*_{d} (*n*_{u}) horizontal slits emanating from the right (left) side of the rectangle. By the same reasoning as before, we obtain that the map is determined by3.2where again and are real-valued constants to be determined. (Note that in equations (3.1) and (3.2), a given constant is the pre-image of the point on the bubble that lies furthest away from the respective edge; see figure 2.)

For calculation purposes, it is convenient to expand the numerator in equations (3.1) and (3.2), and so we write3.3and3.4where the *a*_{j} and *b*_{j} values are real-valued coefficients, with *a*_{nl+nr}=*b*_{nd+nu}=1. (The *a*_{j} and *b*_{j} values could of course be expressed in terms of the values, but it is more convenient to work directly with the mappings (3.3) and (3.4).) For a given set of parameters *α*_{i}, *β*_{i} and , the coefficients *a*_{j} and *b*_{j} are determined as follows.

First, let us define the quantities as3.5for *a*=l, r, with *i*=1,…,*n*_{a} and *j*=0,1,…,*n*_{l}+*n*_{r}. We then note from figure 2 that the mapping functions *W*(*ζ*) must satisfy the condition , for *a*=**l**,**r**, which in view of equations (3.3) and (3.5) yield3.6

Now recalling that *a*_{nl+nr}=1, we see that equation (3.6) gives a system of *n*_{l}+*n*_{r} linear equations for the *n*_{l}+*n*_{r} unknown coefficients *a*_{j}. Similarly, if we define the quantities3.7for *a*=d, u, with *i*=1,…,*n*_{a} and *j*=0,1,…,*n*_{d}+*n*_{u}, and use the fact that must satisfy the condition , for *a*=d, u, we obtain3.8which represent a system of *n*_{u}+*n*_{d} linear equations for the *n*_{u}+*n*_{d} coefficients *b*_{j}. Thus, once the set of parameters *α*_{i}, *β*_{i} and are given, we can readily obtain the constants *a*_{j} and *b*_{j} by solving the linear equations (3.6) and (3.8).

The remaining constants *K* and are determined by the boundary conditions (2.3) and (2.10), respectively. To see this, first note that *K* can be calculated from the requirement that as we go from point B to point F in figure 2*d*, the complex potential *W* in figure 2*b* must vary by i (recall that we have made *U*=2, *V* =1 and *a*=1), and so we write3.9which yields3.10The constant follows from a similar requirement,3.11which yields3.12

We have thus seen that the generic solution described in equations (3.3) and (3.4) is completely specified by prescribing the set of free parameters *α*_{i}, *β*_{i} and , from which all other constants appearing in the solution can be calculated. (Recall that on account of the three degrees of freedom allowed by the Riemann mapping theorem, we can fix any three of the above free parameters.) We note furthermore that all physical parameters of the solution can be calculated once the mathematical parameters are given. For example, the streamwise half-period *L* can be calculated from the requirement that as we go from point F to point M in figure 2*a*, the mapping *z*(*ζ*) must change by *L*,3.13which in view of equation (2.15) reads3.14so that from the knowledge of *W*(*ζ*) and , one readily obtains *L*. Similarly, the bubble centroids and areas can, in principle, be obtained by carrying out the appropriate calculations, but we shall not go into these details here. In the next section, we will give specific examples of the generic solutions described above.

## 4. Examples

In this section, we shall consider some particular cases of the generic solution given in the preceding section. We start discussing the situation where we can confine the solutions to a Hele-Shaw channel. After this, we will give an example of solutions that cannot be realized within the channel geometry and hence must be considered in an unbounded Hele-Shaw cell.

### (a) Periodic bubbles in a channel

Here, we consider the particular and experimentally relevant case of a periodic array of bubbles moving in a Hele-Shaw channel. We suppose the channel walls are at *y*=±1, so that *y*=0 corresponds to the channel centreline. For this case, we set in the generic solution given in equations (3.3) and (3.4) the following conditions:4.1since there can be no bubble attached to the channel walls. Thus, for the case of a periodic array of bubbles in the channel geometry, the solution is completely specified by prescribing the following set of free parameters: {*α*_{i}}_{i=2,3}, {*β*_{i}}_{i=2,3} and , with *a*=l, r, d. (Recall that from the third degree of freedom allowed by the Riemann mapping theorem, one of these parameters can be kept fixed.) We also note that the generic solutions for multiple bubbles in a Hele-Shaw channel presented by Vasconcelos (2001) can be obtained as a special case of our periodic solution above in the limit . This can be accomplished by simply setting *β*_{1}=*β*_{2}=−1 and *β*_{3}=*β*_{4}=1 in figure 2*d*. In what follows, we shall discuss some specific examples of periodic solutions in a Hele-Shaw channel.

We start with the case of a single bubble per unit cell. In the notation of figure 1, this is achieved by placing a bubble at the lower left corner of the reduced unit cell and no other bubble anywhere. In this case, we reproduce the exact solution for an infinite stream of bubbles in a Hele-Shaw cell originally obtained by Burgess & Tanveer (1991) using a rather different method. For completeness, we show in figure 3 an example of the Burgess–Tanveer solution. Note, in particular, that in this case, the bubble is placed at the centre of the unit cell and therefore has both centreline and fore-and-aft symmetry.

As the next example, we consider the situation in which there are bubbles only either at the lower edge or at the left edge of our reduced unit cell. Examples of these two cases are shown in figure 4. In such cases, we recover the class of periodic solutions with centreline symmetry or fore-and-aft symmetry, respectively, obtained earlier by Vasconcelos (1994). It should be noted, however, that in the work by Vasconcelos (1994), no specific example was given, only the general method of solution was presented. In this sense, this is the first time that explicit solutions for such symmetrical configurations have been calculated.

Next, we note that novel solutions can be generated if we place bubbles simultaneously at the left and at the right edges of the reduced unit cell (and nowhere else). For example, in figure 5, we show a solution with two bubbles per unit cell, where each bubble has fore-and-aft symmetry but they lie in different equipotentials of the flow. Seen in the laboratory frame, this solution corresponds to a staggered two-file array of bubbles moving steadily down the channel. It is interesting to note that similar zipper-like arrangements are observed in the flow of red cells in capillaries (Sugihara-Seki & Fu 2005).

Another type of new solution includes the case of ‘mixed symmetry’ where some bubbles have centreline symmetry and some others have fore-and-aft symmetry. An example of such a case is the multi-file array of bubbles shown in figure 6, where there is one file of bubbles symmetrical about the channel centreline, with the other files of bubbles having fore-and-aft symmetry. Another solution with mixed symmetry is shown in figure 7 where, in this case, there is a symmetrical bubble at the centre of the unit cell.

### (b) Solutions in an unbounded cell

Here, we briefly discuss the case in which the periodic array of bubbles cannot be physically realized within the channel geometry and hence must be considered in an unbounded Hele-Shaw cell. We begin by noting that if there are bubbles placed on the upper and lower edges of our reduced unit cell, then the solutions cannot be reduced to the channel geometry as before, and so one must consider an unbounded Hele-Shaw cell. An example of such a case is given in figure 8. It is worth mentioning that this type of solution represents the most general solution for a periodic array of (symmetrical) bubbles in a Hele-Shaw cell, in the sense that the other types of solution can be obtained by either ‘decimating’ or ‘adding’ bubbles along any of the edges of the reduced unit cell.

## 5. Conclusions and discussion

We have presented an exact solution for a doubly periodic array of steadily moving bubbles in a Hele-Shaw cell with an arbitrary number of bubbles per unit cell. Our solution represents the most general periodic solution for bubbles in a Hele-Shaw cell known to date, in the sense that the previously known solutions for a stream of bubbles in a Hele-Shaw channel are particular cases of the solutions reported here. In addition, our solution includes novel cases of periodic arrays of bubbles in a Hele-Shaw channel, such as multi-file flow of symmetrical bubbles and configurations with mixed symmetry, where some bubbles have centreline symmetry while others have fore-and-aft symmetry. Examples of a doubly periodic array of bubbles in an unbounded cell that cannot be restricted to the channel geometry were also discussed.

In constructing our solutions, we assumed that the bubbles are either symmetrical with respect to the cell centreline or have fore-and-aft symmetry (or both), so that the relevant flow domain can be reduced to a simply connected unit cell. In this context, the existence of periodic solutions in a Hele-Shaw cell without any symmetry requirement, where the unit cell is no longer simply connected, remains an open and interesting question. Recently, a new methodology to treat Hele-Shaw flows in multiply connected domains has been developed by Richardson (2001*a*) and Crowdy (2009*a*,*b*). In particular, solutions for a finite number of asymmetric steadily moving bubbles have been found both in an infinite Hele-Shaw cell (Crowdy 2009*a*) and in the channel geometry (Crowdy 2009*b*). It is therefore probable that our solutions can be generalized to include the case of a periodic array of asymmetric bubbles. This interesting (and more difficult) problem clearly deserves to be studied further.

Another interesting point to be investigated is the effect of surface tension on the shape and velocity of our periodic array of bubbles. From previous work on the related problem for a single bubble (Tanveer 1986, 1987; Combescot & Dombre 1988) as well for the Burgess–Tanveer periodic solution (Burgess & Tanveer 1991), one expects that the ‘degeneracy’ displayed by our solutions, meaning that fixing the geometrical parameters of the solution does not fix the bubble velocity *U*, will most certainly be removed by surface tension. However, it is not clear which particular bubble arrangements, among our large class of solutions, will have counterparts in the case of non-zero surface tension. For instance, in the case of a single bubble in a Hele-Shaw channel, Tanveer (1987) was unable to find solutions without symmetry with respect to the channel centreline, even though such non-symmetric solutions exist when surface tension is neglected. In light of these results, it appears that periodic solutions with centreline symmetry are more likely to ‘survive’ in the non-zero surface-tension scenario.

Unfortunately, however, the search for solutions with surface tension has to be performed numerically—a task that is likely to become increasingly more demanding as the number of bubble grows.

## Acknowledgements

This work was supported in part by the Brazilian agencies FINEP, CNPq and FACEPE, and by the special programs PRONEX and CT-PETRO.

## Footnotes

- Received April 30, 2010.
- Accepted June 16, 2010.

- © 2010 The Royal Society